Changeset 13621


Ignore:
Timestamp:
2020-10-16T14:16:38+02:00 (6 months ago)
Author:
techene
Message:

#2385 ln_dynvor=T should be ok

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90

    r13295 r13621  
    2121   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 
    2222   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation 
     23   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity 
    2324   !!---------------------------------------------------------------------- 
    2425 
    2526   !!---------------------------------------------------------------------- 
    2627   !!   dyn_vor       : Update the momentum trend with the vorticity trend 
     28   !!       vor_enT   : energy conserving scheme at T-pt  (ln_dynvor_enT=T) 
     29   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T) 
    2730   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T) 
    28    !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T) 
    2931   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T) 
     32   !!       vor_eeT   : energy conserving at T-pt         (ln_dynvor_eeT=T) 
    3033   !!   dyn_vor_init  : set and control of the different vorticity option 
    3134   !!---------------------------------------------------------------------- 
     
    8083   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation 
    8184   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -  
    82    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2v)/(2*e1e2f)  used in F-point metric term calculation 
    83    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       -  
     85   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation 
     86   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -  
    8487    
    8588   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
     
    218221      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
    219222      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace 
    220       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace 
     223      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwz      ! 3D workspace 
    221224      !!---------------------------------------------------------------------- 
    222225      ! 
     
    228231      ! 
    229232      ! 
    230       SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
    231       CASE ( np_RVO )                           !* relative vorticity 
    232          DO jk = 1, jpkm1                                 ! Horizontal slab 
     233      SELECT CASE( kvor )                 !== relative vorticity considered  ==! 
     234      ! 
     235      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used 
     236         ALLOCATE( zwz(jpi,jpj,jpk) ) 
     237         DO jk = 1, jpkm1                                ! Horizontal slab 
    233238            DO_2D( 1, 0, 1, 0 ) 
    234239               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
    235240                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    236241            END_2D 
    237             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     242            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity  
    238243               DO_2D( 1, 0, 1, 0 ) 
    239244                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     
    241246            ENDIF 
    242247         END DO 
    243  
    244          CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    245  
    246       CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    247          DO jk = 1, jpkm1                                 ! Horizontal slab 
    248             DO_2D( 1, 0, 1, 0 ) 
    249                zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
    250                   &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    251             END_2D 
    252             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
    253                DO_2D( 1, 0, 1, 0 ) 
    254                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    255                END_2D 
    256             ENDIF 
    257          END DO 
    258  
    259          CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    260  
     248         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. ) 
     249         ! 
    261250      END SELECT 
    262251 
    263252      !                                                ! =============== 
    264253      DO jk = 1, jpkm1                                 ! Horizontal slab 
    265       !                                                ! =============== 
    266  
     254         !                                             ! =============== 
     255         ! 
    267256         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
     257         ! 
    268258         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    269259            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 
    270260         CASE ( np_RVO )                           !* relative vorticity 
    271261            DO_2D( 0, 1, 0, 1 ) 
    272                zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    273                   &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 
    274                   &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     262               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)       & 
     263                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk)   )  & 
     264                  &              * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    275265            END_2D 
    276266         CASE ( np_MET )                           !* metric term 
    277267            DO_2D( 0, 1, 0, 1 ) 
    278                zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    279                   &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) & 
    280                   &             * e3t(ji,jj,jk,Kmm) 
     268               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)       & 
     269                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   )   & 
     270                  &       * e3t(ji,jj,jk,Kmm) 
    281271            END_2D 
    282272         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    283273            DO_2D( 0, 1, 0, 1 ) 
    284                zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    285                   &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) & 
    286                   &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     274               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)        & 
     275                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  )   & 
     276                  &       * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    287277            END_2D 
    288278         CASE ( np_CME )                           !* Coriolis + metric 
    289279            DO_2D( 0, 1, 0, 1 ) 
    290                zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    291                     &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    292                     &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) & 
    293                     &          * e3t(ji,jj,jk,Kmm) 
     280               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                               & 
     281                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)      & 
     282                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  )   & 
     283                    &     * e3t(ji,jj,jk,Kmm) 
    294284            END_2D 
    295285         CASE DEFAULT                                             ! error 
    296             CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
     286            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor') 
    297287         END SELECT 
    298288         ! 
     
    310300      END DO                                           !   End of slab 
    311301      !                                                ! =============== 
     302      ! 
     303      SELECT CASE( kvor )        ! deallocate zwz if necessary 
     304      CASE ( np_RVO , np_CRV )   ;   DEALLOCATE( zwz ) 
     305      END SELECT 
     306      ! 
    312307   END SUBROUTINE vor_enT 
    313308 
     
    362357                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    363358            END_2D 
     359            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     360               DO_2D( 1, 0, 1, 0 ) 
     361                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     362               END_2D 
     363            ENDIF 
    364364         CASE ( np_MET )                           !* metric term 
    365365            DO_2D( 1, 0, 1, 0 ) 
     
    372372                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    373373            END_2D 
     374            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term) 
     375               DO_2D( 1, 0, 1, 0 ) 
     376                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
     377               END_2D 
     378            ENDIF 
    374379         CASE ( np_CME )                           !* Coriolis + metric 
    375380            DO_2D( 1, 0, 1, 0 ) 
     
    381386         END SELECT 
    382387         ! 
    383          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    384             DO_2D( 1, 0, 1, 0 ) 
    385                zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    386             END_2D 
    387          ENDIF 
    388  
    389          IF( ln_sco ) THEN 
    390             zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    391             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    392             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    393          ELSE 
    394             zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
    395             zwy(:,:) = e1v(:,:) * pv(:,:,jk) 
    396          ENDIF 
     388         !                                   !==  horizontal fluxes and potential vorticity ==! 
     389         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     390         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     391         zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     392         ! 
    397393         !                                   !==  compute and add the vorticity term trend  =! 
    398394         DO_2D( 0, 0, 0, 0 ) 
     
    458454                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    459455            END_2D 
     456            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     457               DO_2D( 1, 0, 1, 0 ) 
     458                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     459               END_2D 
     460            ENDIF 
    460461         CASE ( np_MET )                           !* metric term 
    461462            DO_2D( 1, 0, 1, 0 ) 
     
    468469                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    469470            END_2D 
     471            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term) 
     472               DO_2D( 1, 0, 1, 0 ) 
     473                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
     474               END_2D 
     475            ENDIF 
    470476         CASE ( np_CME )                           !* Coriolis + metric 
    471477            DO_2D( 1, 0, 1, 0 ) 
     
    477483         END SELECT 
    478484         ! 
    479          IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==! 
    480             DO_2D( 1, 0, 1, 0 ) 
    481                zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    482             END_2D 
    483          ENDIF 
    484          ! 
    485          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    486             zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    487             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    488             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    489          ELSE 
    490             zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
    491             zwy(:,:) = e1v(:,:) * pv(:,:,jk) 
    492          ENDIF 
     485         ! 
     486         !                                   !==  horizontal fluxes and potential vorticity ==! 
     487         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     488         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     489         zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     490         ! 
    493491         !                                   !==  compute and add the vorticity term trend  =! 
    494492         DO_2D( 0, 0, 0, 0 ) 
     
    574572         ! 
    575573         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     574         ! 
    576575         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    577576            DO_2D( 1, 0, 1, 0 ) 
     
    583582                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 
    584583            END_2D 
     584            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     585               DO_2D( 1, 0, 1, 0 ) 
     586                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     587               END_2D 
     588            ENDIF 
    585589         CASE ( np_MET )                           !* metric term 
    586590            DO_2D( 1, 0, 1, 0 ) 
     
    594598                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    595599            END_2D 
     600            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     601               DO_2D( 1, 0, 1, 0 ) 
     602                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)  
     603               END_2D 
     604            ENDIF 
    596605         CASE ( np_CME )                           !* Coriolis + metric 
    597606            DO_2D( 1, 0, 1, 0 ) 
     
    602611            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    603612         END SELECT 
    604          ! 
    605          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    606             DO_2D( 1, 0, 1, 0 ) 
    607                zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    608             END_2D 
    609          ENDIF 
     613         !                                             ! =============== 
    610614      END DO                                           !   End of slab 
    611          ! 
     615      !                                                ! =============== 
     616      ! 
    612617      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    613  
     618      ! 
     619      !                                                ! =============== 
    614620      DO jk = 1, jpkm1                                 ! Horizontal slab 
     621         !                                             ! =============== 
    615622         ! 
    616623         !                                   !==  horizontal fluxes  ==! 
    617624         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    618625         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    619  
     626         ! 
    620627         !                                   !==  compute and add the vorticity term trend  =! 
    621          jj = 2 
    622          ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
    623          DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
    624                ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) 
    625                ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) 
    626                ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) 
    627                ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) 
    628          END DO 
    629          DO jj = 3, jpj 
    630             DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3 
    631                ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) 
    632                ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) 
    633                ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) 
    634                ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) 
    635             END DO 
    636          END DO 
     628         DO_2D( 0, 1, 0, 1 ) 
     629            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) 
     630            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) 
     631            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) 
     632            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) 
     633         END_2D 
     634         ! 
    637635         DO_2D( 0, 0, 0, 0 ) 
    638636            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     
    704702                  &          * r1_e1e2f(ji,jj) 
    705703            END_2D 
     704            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     705               DO_2D( 1, 0, 1, 0 ) 
     706                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     707               END_2D 
     708            ENDIF 
    706709         CASE ( np_MET )                           !* metric term 
    707710            DO_2D( 1, 0, 1, 0 ) 
     
    715718                  &                         * r1_e1e2f(ji,jj)    ) 
    716719            END_2D 
     720            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     721               DO_2D( 1, 0, 1, 0 ) 
     722                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)  
     723               END_2D 
     724            ENDIF 
    717725         CASE ( np_CME )                           !* Coriolis + metric 
    718726            DO_2D( 1, 0, 1, 0 ) 
     
    724732         END SELECT 
    725733         ! 
    726          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    727             DO_2D( 1, 0, 1, 0 ) 
    728                zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    729             END_2D 
    730          ENDIF 
    731       END DO 
     734         !                                             ! =============== 
     735      END DO                                           !   End of slab 
     736      !                                                ! =============== 
    732737      ! 
    733738      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    734739      ! 
     740      !                                                ! =============== 
    735741      DO jk = 1, jpkm1                                 ! Horizontal slab 
    736  
    737       !                                   !==  horizontal fluxes  ==! 
     742         !                                             ! =============== 
     743         ! 
     744         !                                   !==  horizontal fluxes  ==! 
    738745         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    739746         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    740747 
    741748         !                                   !==  compute and add the vorticity term trend  =! 
    742          jj = 2 
    743          ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
    744          DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
    745                z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    746                ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
    747                ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
    748                ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t 
    749                ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t 
    750          END DO 
    751          DO jj = 3, jpj 
    752             DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3 
    753                z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    754                ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
    755                ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
    756                ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t 
    757                ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t 
    758             END DO 
    759          END DO 
     749         DO_2D( 0, 1, 0, 1 ) 
     750            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
     751            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
     752            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
     753            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t 
     754            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t 
     755         END_2D 
     756         ! 
    760757         DO_2D( 0, 0, 0, 0 ) 
    761758            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     
    809806         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
    810807      ENDIF 
    811  
    812       IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available') 
    813808 
    814809!!gm  this should be removed when choosing a unique strategy for fmask at the coast 
Note: See TracChangeset for help on using the changeset viewer.