New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12377 for NEMO/trunk/src/OCE/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years 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/zdftke.F90

    r11536 r12377  
    8989 
    9090   !! * Substitutions 
    91 #  include "vectopt_loop_substitute.h90" 
     91#  include "do_loop_substitute.h90" 
    9292   !!---------------------------------------------------------------------- 
    9393   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    109109 
    110110 
    111    SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 
     111   SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    112112      !!---------------------------------------------------------------------- 
    113113      !!                   ***  ROUTINE zdf_tke  *** 
    … …  
    155155      !!---------------------------------------------------------------------- 
    156156      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
     157      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    157158      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    158159      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    159160      !!---------------------------------------------------------------------- 
    160161      ! 
    161       CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )   ! now tke (en) 
    162       ! 
    163       CALL tke_avn( gdepw_n, e3t_n, e3w_n,        p_avm, p_avt )   ! now avt, avm, dissl 
     162      CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )   ! now tke (en) 
     163      ! 
     164      CALL tke_avn( Kbb, Kmm,        p_avm, p_avt )   ! now avt, avm, dissl 
    164165      ! 
    165166  END SUBROUTINE zdf_tke 
    166167 
    167168 
    168    SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt ) 
     169   SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    169170      !!---------------------------------------------------------------------- 
    170171      !!                   ***  ROUTINE tke_tke  *** 
    … …  
    186187      USE zdf_oce , ONLY : en   ! ocean vertical physics 
    187188      !! 
    188       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! depth of w-points 
    189       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points) 
     189      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    190190      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term 
    191191      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    … …  
    215215      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    216216       
    217       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    218          DO ji = fs_2, fs_jpim1   ! vector opt. 
    219             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    220          END DO 
    221       END DO 
     217      DO_2D_00_00 
     218         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     219      END_2D 
    222220      IF ( ln_isfcav ) THEN 
    223          DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin 
    224             DO ji = fs_2, fs_jpim1   ! vector opt. 
    225                en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 
    226             END DO 
    227          END DO 
     221         DO_2D_00_00 
     222            en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 
     223         END_2D 
    228224      ENDIF 
    229225      ! 
    … …  
    238234      IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE 
    239235         ! 
    240          DO jj = 2, jpjm1           ! bottom friction 
    241             DO ji = fs_2, fs_jpim1     ! vector opt. 
    242                zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    243                zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
    244                !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
    245                zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  & 
    246                   &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
    247                en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
    248             END DO 
    249          END DO 
     236         DO_2D_00_00 
     237            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     238            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
     239            !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
     240            zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  & 
     241               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) 
     242            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
     243         END_2D 
    250244         IF( ln_isfcav ) THEN       ! top friction 
    251             DO jj = 2, jpjm1 
    252                DO ji = fs_2, fs_jpim1   ! vector opt. 
    253                   zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    254                   zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
    255                   !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
    256                   zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    257                      &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
    258                   en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface 
    259                END DO 
    260             END DO 
     245            DO_2D_00_00 
     246               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
     247               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
     248               !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
     249               zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  & 
     250                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  ) 
     251               en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface 
     252            END_2D 
    261253         ENDIF 
    262254         ! 
    … …  
    268260         ! 
    269261         !                        !* total energy produce by LC : cumulative sum over jk 
    270          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 
     262         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 
    271263         DO jk = 2, jpk 
    272             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 
     264            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 
    273265         END DO 
    274266         !                        !* finite Langmuir Circulation depth 
    275267         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    276268         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    277          DO jk = jpkm1, 2, -1 
    278             DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us  
    279                DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1) 
    280                   zus  = zcof * taum(ji,jj) 
    281                   IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
    282                END DO 
    283             END DO 
    284          END DO 
     269         DO_3DS_11_11( jpkm1, 2, -1 ) 
     270            zus  = zcof * taum(ji,jj) 
     271            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
     272         END_3D 
    285273         !                               ! finite LC depth 
    286          DO jj = 1, jpj  
    287             DO ji = 1, jpi 
    288                zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 
    289             END DO 
    290          END DO 
     274         DO_2D_11_11 
     275            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 
     276         END_2D 
    291277         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    292          DO jj = 2, jpjm1 
    293             DO ji = fs_2, fs_jpim1   ! vector opt. 
    294                zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    295                zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    296                IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
    297             END DO 
    298          END DO          
    299          DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    300             DO jj = 2, jpjm1 
    301                DO ji = fs_2, fs_jpim1   ! vector opt. 
    302                   IF ( zfr_i(ji,jj) /= 0. ) THEN                
    303                      ! vertical velocity due to LC    
    304                      IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    305                         !                                           ! vertical velocity due to LC 
    306                         zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
    307                         !                                           ! TKE Langmuir circulation source term 
    308                         en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    309                      ENDIF 
    310                   ENDIF 
    311                END DO 
    312             END DO 
    313          END DO 
     278         DO_2D_00_00 
     279            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
     280            zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
     281            IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     282         END_2D 
     283         DO_3D_00_00( 2, jpkm1 ) 
     284            IF ( zfr_i(ji,jj) /= 0. ) THEN                
     285               ! vertical velocity due to LC    
     286               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
     287                  !                                           ! vertical velocity due to LC 
     288                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     289                  !                                           ! TKE Langmuir circulation source term 
     290                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     291               ENDIF 
     292            ENDIF 
     293         END_3D 
    314294         ! 
    315295      ENDIF 
    … …  
    323303      ! 
    324304      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri ) 
    325          DO jk = 2, jpkm1 
    326             DO jj = 2, jpjm1 
    327                DO ji = 2, jpim1 
    328                   !                             ! local Richardson number 
    329                   zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 
    330                   !                             ! inverse of Prandtl number 
    331                   apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
    332                END DO 
    333             END DO 
    334          END DO 
     305         DO_3D_00_00( 2, jpkm1 ) 
     306            !                             ! local Richardson number 
     307            zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 
     308            !                             ! inverse of Prandtl number 
     309            apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
     310         END_3D 
    335311      ENDIF 
    336312      !          
    337       DO jk = 2, jpkm1           !* Matrix and right hand side in en 
    338          DO jj = 2, jpjm1 
    339             DO ji = fs_2, fs_jpim1   ! vector opt. 
    340                zcof   = zfact1 * tmask(ji,jj,jk) 
    341                !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 
    342                !                                   ! eddy coefficient (ensure numerical stability) 
    343                zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
    344                   &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  ) 
    345                zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
    346                   &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  ) 
    347                ! 
    348                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    349                zd_lw(ji,jj,jk) = zzd_lw 
    350                zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 
    351                ! 
    352                !                                   ! right hand side in en 
    353                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear 
    354                   &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
    355                   &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
    356                   &                                ) * wmask(ji,jj,jk) 
    357             END DO 
    358          END DO 
    359       END DO 
     313      DO_3D_00_00( 2, jpkm1 ) 
     314         zcof   = zfact1 * tmask(ji,jj,jk) 
     315         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 
     316         !                                   ! eddy coefficient (ensure numerical stability) 
     317         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
     318            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
     319         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
     320            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
     321         ! 
     322         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     323         zd_lw(ji,jj,jk) = zzd_lw 
     324         zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 
     325         ! 
     326         !                                   ! right hand side in en 
     327         en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear 
     328            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
     329            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
     330            &                                ) * wmask(ji,jj,jk) 
     331      END_3D 
    360332      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    361       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    362          DO jj = 2, jpjm1 
    363             DO ji = fs_2, fs_jpim1    ! vector opt. 
    364                zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    365             END DO 
    366          END DO 
    367       END DO 
    368       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    369          DO ji = fs_2, fs_jpim1   ! vector opt. 
    370             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    371          END DO 
    372       END DO 
    373       DO jk = 3, jpkm1 
    374          DO jj = 2, jpjm1 
    375             DO ji = fs_2, fs_jpim1    ! vector opt. 
    376                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) 
    377             END DO 
    378          END DO 
    379       END DO 
    380       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    381          DO ji = fs_2, fs_jpim1   ! vector opt. 
    382             en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    383          END DO 
    384       END DO 
    385       DO jk = jpk-2, 2, -1 
    386          DO jj = 2, jpjm1 
    387             DO ji = fs_2, fs_jpim1    ! vector opt. 
    388                en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    389             END DO 
    390          END DO 
    391       END DO 
    392       DO jk = 2, jpkm1                             ! set the minimum value of tke 
    393          DO jj = 2, jpjm1 
    394             DO ji = fs_2, fs_jpim1   ! vector opt. 
    395                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    396             END DO 
    397          END DO 
    398       END DO 
     333      DO_3D_00_00( 3, jpkm1 ) 
     334         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
     335      END_3D 
     336      DO_2D_00_00 
     337         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     338      END_2D 
     339      DO_3D_00_00( 3, jpkm1 ) 
     340         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) 
     341      END_3D 
     342      DO_2D_00_00 
     343         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
     344      END_2D 
     345      DO_3DS_00_00( jpk-2, 2, -1 ) 
     346         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
     347      END_3D 
     348      DO_3D_00_00( 2, jpkm1 ) 
     349         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
     350      END_3D 
    399351      ! 
    400352      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    402354      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    403355!!gm BUG : in the exp  remove the depth of ssh !!! 
    404 !!gm       i.e. use gde3w in argument (pdepw) 
     356!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 
    405357       
    406358       
    407359      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    408          DO jk = 2, jpkm1                       ! rn_eice =0 ON below sea-ice, =4 OFF when ice fraction > 0.25 
    409             DO jj = 2, jpjm1 
    410                DO ji = fs_2, fs_jpim1   ! vector opt. 
    411                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    412                      &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    413                END DO 
    414             END DO 
    415          END DO 
     360         DO_3D_00_00( 2, jpkm1 ) 
     361            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
     362               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     363         END_3D 
    416364      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction) 
    417          DO jj = 2, jpjm1 
    418             DO ji = fs_2, fs_jpim1   ! vector opt. 
    419                jk = nmln(ji,jj) 
    420                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    421                   &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    422             END DO 
    423          END DO 
     365         DO_2D_00_00 
     366            jk = nmln(ji,jj) 
     367            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
     368               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     369         END_2D 
    424370      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability) 
    425          DO jk = 2, jpkm1 
    426             DO jj = 2, jpjm1 
    427                DO ji = fs_2, fs_jpim1   ! vector opt. 
    428                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    429                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    430                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
    431                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    432                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    433                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    434                      &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    435                END DO 
    436             END DO 
    437          END DO 
     371         DO_3D_00_00( 2, jpkm1 ) 
     372            ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
     373            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
     374            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
     375            zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
     376            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
     377            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
     378               &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     379         END_3D 
    438380      ENDIF 
    439381      ! 
    … …  
    441383 
    442384 
    443    SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 
     385   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt ) 
    444386      !!---------------------------------------------------------------------- 
    445387      !!                   ***  ROUTINE tke_avn  *** 
    … …  
    477419      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics 
    478420      !! 
    479       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! depth (w-points) 
    480       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points) 
     421      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    481422      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    482423      ! 
    … …  
    500441      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    501442         zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    502          DO jj = 2, jpjm1 
    503             DO ji = fs_2, fs_jpim1 
    504                zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    505             END DO 
    506          END DO 
     443         DO_2D_00_00 
     444            zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
     445         END_2D 
    507446      ELSE  
    508447         zmxlm(:,:,1) = rn_mxl0 
    509448      ENDIF 
    510449      ! 
    511       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    512          DO jj = 2, jpjm1 
    513             DO ji = fs_2, fs_jpim1   ! vector opt. 
    514                zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    515                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    516             END DO 
    517          END DO 
    518       END DO 
     450      DO_3D_00_00( 2, jpkm1 ) 
     451         zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
     452         zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
     453      END_3D 
    519454      ! 
    520455      !                     !* Physical limits for the mixing length 
    … …  
    526461      ! 
    527462 !!gm Not sure of that coding for ISF.... 
    528       ! where wmask = 0 set zmxlm == p_e3w 
     463      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 
    529464      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    530          DO jk = 2, jpkm1 
    531             DO jj = 2, jpjm1 
    532                DO ji = fs_2, fs_jpim1   ! vector opt. 
    533                   zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    534                   &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 
    535                   ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
    536                   zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
    537                   zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
    538                END DO 
    539             END DO 
    540          END DO 
     465         DO_3D_00_00( 2, jpkm1 ) 
     466            zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   & 
     467            &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) 
     468            ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
     469            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 
     470            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 
     471         END_3D 
    541472         ! 
    542473      CASE ( 1 )           ! bounded by the vertical scale factor 
    543          DO jk = 2, jpkm1 
    544             DO jj = 2, jpjm1 
    545                DO ji = fs_2, fs_jpim1   ! vector opt. 
    546                   zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    547                   zmxlm(ji,jj,jk) = zemxl 
    548                   zmxld(ji,jj,jk) = zemxl 
    549                END DO 
    550             END DO 
    551          END DO 
     474         DO_3D_00_00( 2, jpkm1 ) 
     475            zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 
     476            zmxlm(ji,jj,jk) = zemxl 
     477            zmxld(ji,jj,jk) = zemxl 
     478         END_3D 
    552479         ! 
    553480      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    554          DO jk = 2, jpkm1         ! from the surface to the bottom : 
    555             DO jj = 2, jpjm1 
    556                DO ji = fs_2, fs_jpim1   ! vector opt. 
    557                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    558                END DO 
    559             END DO 
    560          END DO 
    561          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
    562             DO jj = 2, jpjm1 
    563                DO ji = fs_2, fs_jpim1   ! vector opt. 
    564                   zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    565                   zmxlm(ji,jj,jk) = zemxl 
    566                   zmxld(ji,jj,jk) = zemxl 
    567                END DO 
    568             END DO 
    569          END DO 
     481         DO_3D_00_00( 2, jpkm1 ) 
     482            zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
     483         END_3D 
     484         DO_3DS_00_00( jpkm1, 2, -1 ) 
     485            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
     486            zmxlm(ji,jj,jk) = zemxl 
     487            zmxld(ji,jj,jk) = zemxl 
     488         END_3D 
    570489         ! 
    571490      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    572          DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
    573             DO jj = 2, jpjm1 
    574                DO ji = fs_2, fs_jpim1   ! vector opt. 
    575                   zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    576                END DO 
    577             END DO 
    578          END DO 
    579          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
    580             DO jj = 2, jpjm1 
    581                DO ji = fs_2, fs_jpim1   ! vector opt. 
    582                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    583                END DO 
    584             END DO 
    585          END DO 
    586          DO jk = 2, jpkm1 
    587             DO jj = 2, jpjm1 
    588                DO ji = fs_2, fs_jpim1   ! vector opt. 
    589                   zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) ) 
    590                   zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 
    591                   zmxlm(ji,jj,jk) = zemlm 
    592                   zmxld(ji,jj,jk) = zemlp 
    593                END DO 
    594             END DO 
    595          END DO 
     491         DO_3D_00_00( 2, jpkm1 ) 
     492            zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
     493         END_3D 
     494         DO_3DS_00_00( jpkm1, 2, -1 ) 
     495            zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
     496         END_3D 
     497         DO_3D_00_00( 2, jpkm1 ) 
     498            zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) ) 
     499            zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 
     500            zmxlm(ji,jj,jk) = zemlm 
     501            zmxld(ji,jj,jk) = zemlp 
     502         END_3D 
    596503         ! 
    597504      END SELECT 
    … …  
    600507      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt) 
    601508      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    602       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points 
    603          DO jj = 2, jpjm1 
    604             DO ji = fs_2, fs_jpim1   ! vector opt. 
    605                zsqen = SQRT( en(ji,jj,jk) ) 
    606                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    607                p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
    608                p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    609                dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    610             END DO 
    611          END DO 
    612       END DO 
     509      DO_3D_00_00( 1, jpkm1 ) 
     510         zsqen = SQRT( en(ji,jj,jk) ) 
     511         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     512         p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     513         p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     514         dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
     515      END_3D 
    613516      ! 
    614517      ! 
    615518      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    616          DO jk = 2, jpkm1 
    617             DO jj = 2, jpjm1 
    618                DO ji = fs_2, fs_jpim1   ! vector opt. 
    619                   p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    620               END DO 
    621             END DO 
    622          END DO 
    623       ENDIF 
    624       ! 
    625       IF(ln_ctl) THEN 
     519         DO_3D_00_00( 2, jpkm1 ) 
     520            p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     521         END_3D 
     522      ENDIF 
     523      ! 
     524      IF(sn_cfctl%l_prtctl) THEN 
    626525         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk) 
    627526         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk ) 
    … …  
    631530 
    632531 
    633    SUBROUTINE zdf_tke_init 
     532   SUBROUTINE zdf_tke_init( Kmm ) 
    634533      !!---------------------------------------------------------------------- 
    635534      !!                  ***  ROUTINE zdf_tke_init  *** 
    … …  
    647546      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag 
    648547      !! 
    649       INTEGER ::   ji, jj, jk   ! dummy loop indices 
    650       INTEGER ::   ios 
     548      INTEGER, INTENT(in) ::   Kmm          ! time level index 
     549      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     550      INTEGER             ::   ios 
    651551      !! 
    652552      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,          & 
    … …  
    656556      !!---------------------------------------------------------------------- 
    657557      ! 
    658       REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy 
    659558      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) 
    660559901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' ) 
    661560 
    662       REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy 
    663561      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 ) 
    664562902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' ) 
    … …  
    725623      ENDIF 
    726624       
    727       IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     625      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln  
    728626 
    729627      !                               !* depth of penetration of surface tke 
Note: See TracChangeset for help on using the changeset viewer.