Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (14 months ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
  • NEMO/trunk/src/OCE/ZDF/zdfgls.F90

    r11536 r12377  
    104104 
    105105   !! * Substitutions 
    106 #  include "vectopt_loop_substitute.h90" 
     106#  include "do_loop_substitute.h90" 
    107107   !!---------------------------------------------------------------------- 
    108108   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    124124 
    125125 
    126    SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt ) 
     126   SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    127127      !!---------------------------------------------------------------------- 
    128128      !!                   ***  ROUTINE zdf_gls  *** 
     
    134134      !! 
    135135      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
     136      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    136137      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    137138      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
     
    166167 
    167168      ! Compute surface, top and bottom friction at T-points 
    168       DO jj = 2, jpjm1           
    169          DO ji = fs_2, fs_jpim1   ! vector opt.          
    170             ! 
    171             ! surface friction 
    172             ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
    173             !    
     169      DO_2D_00_00 
     170         ! 
     171         ! surface friction 
     172         ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
     173         !    
    174174!!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that... 
    175           ! bottom friction (explicit before friction) 
    176           zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    177           zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0) 
    178           ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  & 
    179              &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
    180          END DO 
    181       END DO 
     175       ! bottom friction (explicit before friction) 
     176       zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     177       zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0) 
     178       ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  & 
     179          &                                         + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) 
     180      END_2D 
    182181      IF( ln_isfcav ) THEN       !top friction 
    183          DO jj = 2, jpjm1 
    184             DO ji = fs_2, fs_jpim1   ! vector opt. 
    185                zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    186                zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0) 
    187                ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    188                   &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
    189             END DO 
    190          END DO 
     182         DO_2D_00_00 
     183            zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
     184            zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0) 
     185            ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  & 
     186               &                                         + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  ) 
     187         END_2D 
    191188      ENDIF 
    192189    
     
    206203      END SELECT 
    207204      ! 
    208       DO jk = 2, jpkm1              !==  Compute dissipation rate  ==! 
    209          DO jj = 1, jpjm1 
    210             DO ji = 1, jpim1 
    211                eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
    212             END DO 
    213          END DO 
    214       END DO 
     205      DO_3D_10_10( 2, jpkm1 ) 
     206         eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
     207      END_3D 
    215208 
    216209      ! Save tke at before time step 
     
    219212 
    220213      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
    221          DO jk = 2, jpkm1 
    222             DO jj = 2, jpjm1  
    223                DO ji = fs_2, fs_jpim1   ! vector opt. 
    224                   zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    225                   zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 
    226                   zcoef = ( zup / MAX( zdown, rsmall ) ) 
    227                   zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 
    228                END DO 
    229             END DO 
    230          END DO 
     214         DO_3D_00_00( 2, jpkm1 ) 
     215            zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
     216            zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) 
     217            zcoef = ( zup / MAX( zdown, rsmall ) ) 
     218            zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 
     219         END_3D 
    231220      ENDIF 
    232221 
     
    244233      ! Warning : after this step, en : right hand side of the matrix 
    245234 
    246       DO jk = 2, jpkm1 
    247          DO jj = 2, jpjm1 
    248             DO ji = 2, jpim1 
    249                ! 
    250                buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction 
    251                ! 
    252                diss = eps(ji,jj,jk)                         ! dissipation 
    253                ! 
    254                zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
    255                ! 
    256                zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk)          ! production term 
    257                zdiss = zdir*(diss/en(ji,jj,jk))   +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
     235      DO_3D_00_00( 2, jpkm1 ) 
     236         ! 
     237         buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction 
     238         ! 
     239         diss = eps(ji,jj,jk)                         ! dissipation 
     240         ! 
     241         zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     242         ! 
     243         zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk)          ! production term 
     244         zdiss = zdir*(diss/en(ji,jj,jk))   +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
    258245!!gm better coding, identical results 
    259246!               zesh2 =   p_sh2(ji,jj,jk) + zdir*buoy               ! production term 
    260247!               zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 
    261248!!gm 
    262                ! 
    263                ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 
    264                ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 
    265                ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 
    266                ! Otherwise, this should be rsc_psi/rsc_psi0 
    267                IF( ln_sigpsi ) THEN 
    268                   zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1. 
    269                   zwall_psi(ji,jj,jk) = rsc_psi /   &  
    270                      &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  ) 
    271                ELSE 
    272                   zwall_psi(ji,jj,jk) = 1._wp 
    273                ENDIF 
    274                ! 
    275                ! building the matrix 
    276                zcof = rfact_tke * tmask(ji,jj,jk) 
    277                !                                        ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 
    278                zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 
    279                !                                        ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 
    280                zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
    281                !                                        ! diagonal 
    282                zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk)  
    283                !                                        ! right hand side in en 
    284                en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
    285             END DO 
    286          END DO 
    287       END DO 
     249         ! 
     250         ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 
     251         ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 
     252         ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 
     253         ! Otherwise, this should be rsc_psi/rsc_psi0 
     254         IF( ln_sigpsi ) THEN 
     255            zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1. 
     256            zwall_psi(ji,jj,jk) = rsc_psi /   &  
     257               &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  ) 
     258         ELSE 
     259            zwall_psi(ji,jj,jk) = 1._wp 
     260         ENDIF 
     261         ! 
     262         ! building the matrix 
     263         zcof = rfact_tke * tmask(ji,jj,jk) 
     264         !                                        ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 
     265         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 
     266         !                                        ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 
     267         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) ) 
     268         !                                        ! diagonal 
     269         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk)  
     270         !                                        ! right hand side in en 
     271         en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
     272      END_3D 
    288273      ! 
    289274      zdiag(:,:,jpk) = 1._wp 
     
    306291      !  
    307292      ! One level below 
    308       en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   & 
     293      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm))   & 
    309294         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   ) 
    310295      zd_lw(:,:,2) = 0._wp  
     
    325310      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag 
    326311      zd_lw(:,:,2) = 0._wp 
    327       zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
     312      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 
    328313      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    329           &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
    330 !!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf) 
    331       en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
     314          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
     315!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 
     316      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 
    332317      ! 
    333318      ! 
     
    342327         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 
    343328         !                      ! Balance between the production and the dissipation terms 
    344          DO jj = 2, jpjm1 
    345             DO ji = fs_2, fs_jpim1   ! vector opt. 
     329         DO_2D_00_00 
    346330!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ???? 
    347331!!   With thick deep ocean level thickness, this may be quite large, no ??? 
    348332!!   in particular in ocean cavities where top stratification can be large... 
    349                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    350                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     333            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     334            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     335            ! 
     336            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 
     337            ! 
     338            ! Dirichlet condition applied at:  
     339            !     Bottom level (ibot)      &      Just above it (ibotm1)    
     340            zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp 
     341            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp 
     342            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp 
     343            en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en 
     344         END_2D 
     345         ! 
     346         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     347            DO_2D_00_00 
     348               itop   = mikt(ji,jj)       ! k   top w-point 
     349               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     350               !                                                ! mask at the ocean surface points 
     351               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
    351352               ! 
    352                z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 
    353                ! 
     353 !!gm TO BE VERIFIED !!! 
    354354               ! Dirichlet condition applied at:  
    355                !     Bottom level (ibot)      &      Just above it (ibotm1)    
    356                zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp 
    357                zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp 
    358                zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp 
    359                en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en 
    360             END DO 
    361          END DO 
    362          ! 
    363          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
    364             DO jj = 2, jpjm1 
    365                DO ji = fs_2, fs_jpim1   ! vector opt. 
    366                   itop   = mikt(ji,jj)       ! k   top w-point 
    367                   itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
    368                   !                                                ! mask at the ocean surface points 
    369                   z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
    370                   ! 
    371  !!gm TO BE VERIFIED !!! 
    372                   ! Dirichlet condition applied at:  
    373                   !     top level (itop)         &      Just below it (itopp1)    
    374                   zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp 
    375                   zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp 
    376                   zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp 
    377                   en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en 
    378                END DO 
    379             END DO 
     355               !     top level (itop)         &      Just below it (itopp1)    
     356               zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp 
     357               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp 
     358               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp 
     359               en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en 
     360            END_2D 
    380361         ENDIF 
    381362         ! 
    382363      CASE ( 1 )             ! Neumman boundary condition 
    383364         !                       
    384          DO jj = 2, jpjm1 
    385             DO ji = fs_2, fs_jpim1   ! vector opt. 
    386                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    387                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    388                ! 
    389                z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 
     365         DO_2D_00_00 
     366            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     367            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     368            ! 
     369            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 
     370            ! 
     371            ! Bottom level Dirichlet condition: 
     372            !     Bottom level (ibot)      &      Just above it (ibotm1)    
     373            !         Dirichlet            !         Neumann 
     374            zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag 
     375            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 
     376            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp 
     377         END_2D 
     378         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     379            DO_2D_00_00 
     380               itop   = mikt(ji,jj)       ! k   top w-point 
     381               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     382               !                                                ! mask at the ocean surface points 
     383               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
    390384               ! 
    391385               ! Bottom level Dirichlet condition: 
    392386               !     Bottom level (ibot)      &      Just above it (ibotm1)    
    393387               !         Dirichlet            !         Neumann 
    394                zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag 
    395                zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 
    396                zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp 
    397             END DO 
    398          END DO 
    399          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
    400             DO jj = 2, jpjm1 
    401                DO ji = fs_2, fs_jpim1   ! vector opt. 
    402                   itop   = mikt(ji,jj)       ! k   top w-point 
    403                   itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
    404                   !                                                ! mask at the ocean surface points 
    405                   z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
    406                   ! 
    407                   ! Bottom level Dirichlet condition: 
    408                   !     Bottom level (ibot)      &      Just above it (ibotm1)    
    409                   !         Dirichlet            !         Neumann 
    410                   zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag 
    411                   zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 
    412                   zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp 
    413                END DO 
    414             END DO 
     388               zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag 
     389               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 
     390               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp 
     391            END_2D 
    415392         ENDIF 
    416393         ! 
     
    420397      ! ---------------------------------------------------------- 
    421398      ! 
    422       DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    423          DO jj = 2, jpjm1 
    424             DO ji = fs_2, fs_jpim1    ! vector opt. 
    425                zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    426             END DO 
    427          END DO 
    428       END DO 
    429       DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    430          DO jj = 2, jpjm1 
    431             DO ji = fs_2, fs_jpim1    ! vector opt. 
    432                zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
    433             END DO 
    434          END DO 
    435       END DO 
    436       DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    437          DO jj = 2, jpjm1 
    438             DO ji = fs_2, fs_jpim1    ! vector opt. 
    439                en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    440             END DO 
    441          END DO 
    442       END DO 
     399      DO_3D_00_00( 2, jpkm1 ) 
     400         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
     401      END_3D 
     402      DO_3D_00_00( 2, jpk ) 
     403         zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
     404      END_3D 
     405      DO_3DS_00_00( jpk-1, 2, -1 ) 
     406         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
     407      END_3D 
    443408      !                                            ! set the minimum value of tke  
    444409      en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
     
    453418      ! 
    454419      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    455          DO jk = 2, jpkm1 
    456             DO jj = 2, jpjm1 
    457                DO ji = fs_2, fs_jpim1   ! vector opt. 
    458                   psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 
    459                END DO 
    460             END DO 
    461          END DO 
     420         DO_3D_00_00( 2, jpkm1 ) 
     421            psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 
     422         END_3D 
    462423         ! 
    463424      CASE( 1 )               ! k-eps 
    464          DO jk = 2, jpkm1 
    465             DO jj = 2, jpjm1 
    466                DO ji = fs_2, fs_jpim1   ! vector opt. 
    467                   psi(ji,jj,jk)  = eps(ji,jj,jk) 
    468                END DO 
    469             END DO 
    470          END DO 
     425         DO_3D_00_00( 2, jpkm1 ) 
     426            psi(ji,jj,jk)  = eps(ji,jj,jk) 
     427         END_3D 
    471428         ! 
    472429      CASE( 2 )               ! k-w 
    473          DO jk = 2, jpkm1 
    474             DO jj = 2, jpjm1 
    475                DO ji = fs_2, fs_jpim1   ! vector opt. 
    476                   psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 
    477                END DO 
    478             END DO 
    479          END DO 
     430         DO_3D_00_00( 2, jpkm1 ) 
     431            psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 
     432         END_3D 
    480433         ! 
    481434      CASE( 3 )               ! generic 
    482          DO jk = 2, jpkm1 
    483             DO jj = 2, jpjm1 
    484                DO ji = fs_2, fs_jpim1   ! vector opt. 
    485                   psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn  
    486                END DO 
    487             END DO 
    488          END DO 
     435         DO_3D_00_00( 2, jpkm1 ) 
     436            psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn  
     437         END_3D 
    489438         ! 
    490439      END SELECT 
     
    497446      ! Warning : after this step, en : right hand side of the matrix 
    498447 
    499       DO jk = 2, jpkm1 
    500          DO jj = 2, jpjm1 
    501             DO ji = fs_2, fs_jpim1   ! vector opt. 
    502                ! 
    503                ! psi / k 
    504                zratio = psi(ji,jj,jk) / eb(ji,jj,jk)  
    505                ! 
    506                ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 
    507                zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 
    508                ! 
    509                rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 
    510                ! 
    511                ! shear prod. - stratif. destruction 
    512                prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 
    513                ! 
    514                ! stratif. destruction 
    515                buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
    516                ! 
    517                ! shear prod. - stratif. destruction 
    518                diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 
    519                ! 
    520                zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
    521                ! 
    522                zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term 
    523                zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
    524                !                                                         
    525                ! building the matrix 
    526                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
    527                !                                               ! lower diagonal 
    528                zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 
    529                !                                               ! upper diagonal 
    530                zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
    531                !                                               ! diagonal 
    532                zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 
    533                !                                               ! right hand side in psi 
    534                psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
    535             END DO 
    536          END DO 
    537       END DO 
     448      DO_3D_00_00( 2, jpkm1 ) 
     449         ! 
     450         ! psi / k 
     451         zratio = psi(ji,jj,jk) / eb(ji,jj,jk)  
     452         ! 
     453         ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 
     454         zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 
     455         ! 
     456         rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 
     457         ! 
     458         ! shear prod. - stratif. destruction 
     459         prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 
     460         ! 
     461         ! stratif. destruction 
     462         buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
     463         ! 
     464         ! shear prod. - stratif. destruction 
     465         diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 
     466         ! 
     467         zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     468         ! 
     469         zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term 
     470         zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
     471         !                                                         
     472         ! building the matrix 
     473         zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
     474         !                                               ! lower diagonal 
     475         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 
     476         !                                               ! upper diagonal 
     477         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) ) 
     478         !                                               ! diagonal 
     479         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 
     480         !                                               ! right hand side in psi 
     481         psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
     482      END_3D 
    538483      ! 
    539484      zdiag(:,:,jpk) = 1._wp 
     
    554499         ! 
    555500         ! One level below 
    556          zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 
    557          zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 
     501         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 
     502         zdep    (:,:)   = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 
    558503         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    559504         zd_lw(:,:,2) = 0._wp 
     
    575520         ! 
    576521         ! Set psi vertical flux at the surface: 
    577          zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    578          zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
     522         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 
     523         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 
    579524         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    580525         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 
    581             &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 
     526            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 
    582527         zflxs(:,:)   = zdep(:,:) * zflxs(:,:) 
    583          psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
     528         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 
    584529         ! 
    585530      END SELECT 
     
    596541         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 
    597542         !                      ! Balance between the production and the dissipation terms 
    598          DO jj = 2, jpjm1 
    599             DO ji = fs_2, fs_jpim1   ! vector opt. 
    600                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    601                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    602                zdep(ji,jj) = vkarmn * r_z0_bot 
    603                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    604                zd_lw(ji,jj,ibot) = 0._wp 
    605                zd_up(ji,jj,ibot) = 0._wp 
    606                zdiag(ji,jj,ibot) = 1._wp 
    607                ! 
    608                ! Just above last level, Dirichlet condition again (GOTM like) 
    609                zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 
    610                psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
    611                zd_lw(ji,jj,ibotm1) = 0._wp 
    612                zd_up(ji,jj,ibotm1) = 0._wp 
    613                zdiag(ji,jj,ibotm1) = 1._wp 
    614             END DO 
    615          END DO 
     543         DO_2D_00_00 
     544            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     545            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     546            zdep(ji,jj) = vkarmn * r_z0_bot 
     547            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
     548            zd_lw(ji,jj,ibot) = 0._wp 
     549            zd_up(ji,jj,ibot) = 0._wp 
     550            zdiag(ji,jj,ibot) = 1._wp 
     551            ! 
     552            ! Just above last level, Dirichlet condition again (GOTM like) 
     553            zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) ) 
     554            psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
     555            zd_lw(ji,jj,ibotm1) = 0._wp 
     556            zd_up(ji,jj,ibotm1) = 0._wp 
     557            zdiag(ji,jj,ibotm1) = 1._wp 
     558         END_2D 
    616559         ! 
    617560      CASE ( 1 )             ! Neumman boundary condition 
    618561         !                       
    619          DO jj = 2, jpjm1 
    620             DO ji = fs_2, fs_jpim1   ! vector opt. 
    621                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    622                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    623                ! 
    624                ! Bottom level Dirichlet condition: 
    625                zdep(ji,jj) = vkarmn * r_z0_bot 
    626                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    627                ! 
    628                zd_lw(ji,jj,ibot) = 0._wp 
    629                zd_up(ji,jj,ibot) = 0._wp 
    630                zdiag(ji,jj,ibot) = 1._wp 
    631                ! 
    632                ! Just above last level: Neumann condition with flux injection 
    633                zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 
    634                zd_up(ji,jj,ibotm1) = 0. 
    635                ! 
    636                ! Set psi vertical flux at the bottom: 
    637                zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 
    638                zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   & 
    639                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
    640                psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 
    641             END DO 
    642          END DO 
     562         DO_2D_00_00 
     563            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     564            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     565            ! 
     566            ! Bottom level Dirichlet condition: 
     567            zdep(ji,jj) = vkarmn * r_z0_bot 
     568            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
     569            ! 
     570            zd_lw(ji,jj,ibot) = 0._wp 
     571            zd_up(ji,jj,ibot) = 0._wp 
     572            zdiag(ji,jj,ibot) = 1._wp 
     573            ! 
     574            ! Just above last level: Neumann condition with flux injection 
     575            zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 
     576            zd_up(ji,jj,ibotm1) = 0. 
     577            ! 
     578            ! Set psi vertical flux at the bottom: 
     579            zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm) 
     580            zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   & 
     581               &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
     582            psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm) 
     583         END_2D 
    643584         ! 
    644585      END SELECT 
     
    647588      ! ---------------- 
    648589      ! 
    649       DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    650          DO jj = 2, jpjm1 
    651             DO ji = fs_2, fs_jpim1    ! vector opt. 
    652                zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    653             END DO 
    654          END DO 
    655       END DO 
    656       DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    657          DO jj = 2, jpjm1 
    658             DO ji = fs_2, fs_jpim1    ! vector opt. 
    659                zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
    660             END DO 
    661          END DO 
    662       END DO 
    663       DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    664          DO jj = 2, jpjm1 
    665             DO ji = fs_2, fs_jpim1    ! vector opt. 
    666                psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    667             END DO 
    668          END DO 
    669       END DO 
     590      DO_3D_00_00( 2, jpkm1 ) 
     591         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
     592      END_3D 
     593      DO_3D_00_00( 2, jpk ) 
     594         zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
     595      END_3D 
     596      DO_3DS_00_00( jpk-1, 2, -1 ) 
     597         psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
     598      END_3D 
    670599 
    671600      ! Set dissipation 
     
    675604      ! 
    676605      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    677          DO jk = 1, jpkm1 
    678             DO jj = 2, jpjm1 
    679                DO ji = fs_2, fs_jpim1   ! vector opt. 
    680                   eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 
    681                END DO 
    682             END DO 
    683          END DO 
     606         DO_3D_00_00( 1, jpkm1 ) 
     607            eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 
     608         END_3D 
    684609         ! 
    685610      CASE( 1 )               ! k-eps 
    686          DO jk = 1, jpkm1 
    687             DO jj = 2, jpjm1 
    688                DO ji = fs_2, fs_jpim1   ! vector opt. 
    689                   eps(ji,jj,jk) = psi(ji,jj,jk) 
    690                END DO 
    691             END DO 
    692          END DO 
     611         DO_3D_00_00( 1, jpkm1 ) 
     612            eps(ji,jj,jk) = psi(ji,jj,jk) 
     613         END_3D 
    693614         ! 
    694615      CASE( 2 )               ! k-w 
    695          DO jk = 1, jpkm1 
    696             DO jj = 2, jpjm1 
    697                DO ji = fs_2, fs_jpim1   ! vector opt. 
    698                   eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)  
    699                END DO 
    700             END DO 
    701          END DO 
     616         DO_3D_00_00( 1, jpkm1 ) 
     617            eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)  
     618         END_3D 
    702619         ! 
    703620      CASE( 3 )               ! generic 
     
    705622         zex1  =      ( 1.5_wp + rmm/rnn ) 
    706623         zex2  = -1._wp / rnn 
    707          DO jk = 1, jpkm1 
    708             DO jj = 2, jpjm1 
    709                DO ji = fs_2, fs_jpim1   ! vector opt. 
    710                   eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 
    711                END DO 
    712             END DO 
    713          END DO 
     624         DO_3D_00_00( 1, jpkm1 ) 
     625            eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 
     626         END_3D 
    714627         ! 
    715628      END SELECT 
     
    717630      ! Limit dissipation rate under stable stratification 
    718631      ! -------------------------------------------------- 
    719       DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 
    720          DO jj = 2, jpjm1 
    721             DO ji = fs_2, fs_jpim1    ! vector opt. 
    722                ! limitation 
    723                eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
    724                hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
    725                ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
    726                zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    727                IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 
    728             END DO 
    729          END DO 
    730       END DO  
     632      DO_3D_00_00( 1, jpkm1 ) 
     633         ! limitation 
     634         eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
     635         hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
     636         ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
     637         zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
     638         IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 
     639      END_3D 
    731640 
    732641      ! 
     
    737646      ! 
    738647      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions 
    739          DO jk = 2, jpkm1 
    740             DO jj = 2, jpjm1 
    741                DO ji = fs_2, fs_jpim1   ! vector opt. 
    742                   ! zcof =  l²/q² 
    743                   zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
    744                   ! Gh = -N²l²/q² 
    745                   gh = - rn2(ji,jj,jk) * zcof 
    746                   gh = MIN( gh, rgh0   ) 
    747                   gh = MAX( gh, rghmin ) 
    748                   ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 
    749                   sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 
    750                   sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 
    751                   ! 
    752                   ! Store stability function in zstt and zstm 
    753                   zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
    754                   zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
    755                END DO 
    756             END DO 
    757          END DO 
     648         DO_3D_00_00( 2, jpkm1 ) 
     649            ! zcof =  l²/q² 
     650            zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
     651            ! Gh = -N²l²/q² 
     652            gh = - rn2(ji,jj,jk) * zcof 
     653            gh = MIN( gh, rgh0   ) 
     654            gh = MAX( gh, rghmin ) 
     655            ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 
     656            sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 
     657            sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 
     658            ! 
     659            ! Store stability function in zstt and zstm 
     660            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
     661            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
     662         END_3D 
    758663         ! 
    759664      CASE ( 2, 3 )               ! Canuto stability functions 
    760          DO jk = 2, jpkm1 
    761             DO jj = 2, jpjm1 
    762                DO ji = fs_2, fs_jpim1   ! vector opt. 
    763                   ! zcof =  l²/q² 
    764                   zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
    765                   ! Gh = -N²l²/q² 
    766                   gh = - rn2(ji,jj,jk) * zcof 
    767                   gh = MIN( gh, rgh0   ) 
    768                   gh = MAX( gh, rghmin ) 
    769                   gh = gh * rf6 
    770                   ! Gm =  M²l²/q² Shear number 
    771                   shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 
    772                   gm = MAX( shr * zcof , 1.e-10 ) 
    773                   gm = gm * rf6 
    774                   gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 
    775                   ! Stability functions from Canuto 
    776                   rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 
    777                   sm = (rs0 - rs1*gh + rs2*gm) / rcff 
    778                   sh = (rs4 - rs5*gh + rs6*gm) / rcff 
    779                   ! 
    780                   ! Store stability function in zstt and zstm 
    781                   zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
    782                   zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
    783                END DO 
    784             END DO 
    785          END DO 
     665         DO_3D_00_00( 2, jpkm1 ) 
     666            ! zcof =  l²/q² 
     667            zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
     668            ! Gh = -N²l²/q² 
     669            gh = - rn2(ji,jj,jk) * zcof 
     670            gh = MIN( gh, rgh0   ) 
     671            gh = MAX( gh, rghmin ) 
     672            gh = gh * rf6 
     673            ! Gm =  M²l²/q² Shear number 
     674            shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 
     675            gm = MAX( shr * zcof , 1.e-10 ) 
     676            gm = gm * rf6 
     677            gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 
     678            ! Stability functions from Canuto 
     679            rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 
     680            sm = (rs0 - rs1*gh + rs2*gm) / rcff 
     681            sh = (rs4 - rs5*gh + rs6*gm) / rcff 
     682            ! 
     683            ! Store stability function in zstt and zstm 
     684            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
     685            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
     686         END_3D 
    786687         ! 
    787688      END SELECT 
     
    794695      ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 
    795696      zstm(:,:,jpk) = 0.   
    796       DO jj = 2, jpjm1                ! update bottom with good values 
    797          DO ji = fs_2, fs_jpim1   ! vector opt. 
    798             zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 
    799          END DO 
    800       END DO 
     697      DO_2D_00_00 
     698         zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 
     699      END_2D 
    801700 
    802701      zstt(:,:,  1) = wmask(:,:,  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 
     
    811710      !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 
    812711      !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 
    813       DO jk = 1, jpk 
    814          DO jj = 2, jpjm1 
    815             DO ji = fs_2, fs_jpim1   ! vector opt. 
    816                zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
    817                zavt  = zsqen * zstt(ji,jj,jk) 
    818                zavm  = zsqen * zstm(ji,jj,jk) 
    819                p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
    820                p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
    821             END DO 
    822          END DO 
    823       END DO 
     712      DO_3D_00_00( 1, jpk ) 
     713         zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
     714         zavt  = zsqen * zstt(ji,jj,jk) 
     715         zavm  = zsqen * zstm(ji,jj,jk) 
     716         p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
     717         p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
     718      END_3D 
    824719      p_avt(:,:,1) = 0._wp 
    825720      ! 
    826       IF(ln_ctl) THEN 
     721      IF(sn_cfctl%l_prtctl) THEN 
    827722         CALL prt_ctl( tab3d_1=en   , clinfo1=' gls  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk) 
    828723         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls  - m: ', kdim=jpk ) 
     
    857752      !!---------------------------------------------------------- 
    858753      ! 
    859       REWIND( numnam_ref )              ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme 
    860754      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901) 
    861755901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_gls in reference namelist' ) 
    862756 
    863       REWIND( numnam_cfg )              ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme 
    864757      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 ) 
    865758902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist' ) 
Note: See TracChangeset for help on using the changeset viewer.