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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

Location:
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • 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_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/ZDF/zdftke.F90

    r11405 r13463  
    2828   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    2929   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only  
    30    !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition (ln_drg) 
     30   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition 
    3131   !!---------------------------------------------------------------------- 
    3232 
     
    4545   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
    4646   USE zdfmxl         ! vertical physics: mixed layer 
     47#if defined key_si3 
     48   USE ice, ONLY: hm_i, h_i 
     49#endif 
     50#if defined key_cice 
     51   USE sbc_ice, ONLY: h_i 
     52#endif 
    4753   ! 
    4854   USE in_out_manager ! I/O manager 
     
    6470   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3) 
    6571   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m] 
     72   INTEGER  ::      nn_mxlice ! type of scaling under sea-ice 
     73   REAL(wp) ::      rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1) 
    6674   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1) 
    6775   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
     
    7179   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2] 
    7280   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it) 
    73    LOGICAL  ::   ln_drg    ! top/bottom friction forcing flag  
    7481   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3) 
    7582   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
     
    8996 
    9097   !! * Substitutions 
    91 #  include "vectopt_loop_substitute.h90" 
     98#  include "do_loop_substitute.h90" 
     99#  include "domzgr_substitute.h90" 
    92100   !!---------------------------------------------------------------------- 
    93101   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    109117 
    110118 
    111    SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 
     119   SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    112120      !!---------------------------------------------------------------------- 
    113121      !!                   ***  ROUTINE zdf_tke  *** 
     
    155163      !!---------------------------------------------------------------------- 
    156164      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
     165      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    157166      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    158167      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    159168      !!---------------------------------------------------------------------- 
    160169      ! 
    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 
     170      CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )   ! now tke (en) 
     171      ! 
     172      CALL tke_avn( Kbb, Kmm,        p_avm, p_avt )   ! now avt, avm, dissl 
    164173      ! 
    165174  END SUBROUTINE zdf_tke 
    166175 
    167176 
    168    SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt ) 
     177   SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    169178      !!---------------------------------------------------------------------- 
    170179      !!                   ***  ROUTINE tke_tke  *** 
     
    186195      USE zdf_oce , ONLY : en   ! ocean vertical physics 
    187196      !! 
    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) 
     197      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    190198      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term 
    191199      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
     
    206214      !!-------------------------------------------------------------------- 
    207215      ! 
    208       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
    209       zfact1 = -.5_wp * rdt  
    210       zfact2 = 1.5_wp * rdt * rn_ediss 
     216      zbbrau = rn_ebb / rho0       ! Local constant initialisation 
     217      zfact1 = -.5_wp * rn_Dt  
     218      zfact2 = 1.5_wp * rn_Dt * rn_ediss 
    211219      zfact3 = 0.5_wp       * rn_ediss 
    212220      ! 
     
    214222      !                     !  Surface/top/bottom boundary condition on tke 
    215223      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    216        
    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 
    222       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 
    228       ENDIF 
     224      !  
     225      DO_2D( 0, 0, 0, 0 ) 
     226         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     227      END_2D 
    229228      ! 
    230229      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    232231      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    233232      ! 
    234       !   en(bot)   = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
     233      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
    235234      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 
    236235      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 
    237236      ! 
    238       IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE 
    239          ! 
    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 
    250          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 
     237      IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE 
     238         ! 
     239         DO_2D( 0, 0, 0, 0 )        ! bottom friction 
     240            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     241            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
     242            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
     243            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  & 
     244               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) 
     245            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
     246         END_2D 
     247         IF( ln_isfcav ) THEN 
     248            DO_2D( 0, 0, 0, 0 )     ! top friction 
     249               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
     250               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
     251               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
     252               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  & 
     253                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  ) 
     254               ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present 
     255               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) & 
     256                  &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 
     257            END_2D 
    261258         ENDIF 
    262259         ! 
     
    268265         ! 
    269266         !                        !* total energy produce by LC : cumulative sum over jk 
    270          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 
     267         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 
    271268         DO jk = 2, jpk 
    272             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 
     269            zpelc(:,:,jk)  = zpelc(:,:,jk-1) +   & 
     270               &        MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 
    273271         END DO 
    274272         !                        !* finite Langmuir Circulation depth 
    275273         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    276274         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 
     275         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 
     276            zus  = zcof * taum(ji,jj) 
     277            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
     278         END_3D 
    285279         !                               ! 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 
     280         DO_2D( 1, 1, 1, 1 ) 
     281            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 
     282         END_2D 
    291283         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 
     284         DO_2D( 0, 0, 0, 0 ) 
     285            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
     286            zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
     287            IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     288         END_2D 
     289         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     290            IF ( zfr_i(ji,jj) /= 0. ) THEN                
     291               ! vertical velocity due to LC    
     292               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
     293                  !                                           ! vertical velocity due to LC 
     294                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     295                  !                                           ! TKE Langmuir circulation source term 
     296                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     297               ENDIF 
     298            ENDIF 
     299         END_3D 
    314300         ! 
    315301      ENDIF 
     
    323309      ! 
    324310      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 
     311         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     312            !                             ! local Richardson number 
     313            IF (rn2b(ji,jj,jk) <= 0.0_wp) then 
     314                zri = 0.0_wp 
     315            ELSE 
     316                zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 
     317            ENDIF 
     318            !                             ! inverse of Prandtl number 
     319            apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
     320         END_3D 
    335321      ENDIF 
    336322      !          
    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 
     323      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     324         zcof   = zfact1 * tmask(ji,jj,jk) 
     325         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 
     326         !                                   ! eddy coefficient (ensure numerical stability) 
     327         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
     328            &          /    (  e3t(ji,jj,jk  ,Kmm)   & 
     329            &                * e3w(ji,jj,jk  ,Kmm)  ) 
     330         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
     331            &          /    (  e3t(ji,jj,jk-1,Kmm)   & 
     332            &                * e3w(ji,jj,jk  ,Kmm)  ) 
     333         ! 
     334         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     335         zd_lw(ji,jj,jk) = zzd_lw 
     336         zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 
     337         ! 
     338         !                                   ! right hand side in en 
     339         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                        &   ! shear 
     340            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
     341            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
     342            &                                ) * wmask(ji,jj,jk) 
     343      END_3D 
    360344      !                          !* 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 
     345      DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 
     346         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
     347      END_3D 
     348      DO_2D( 0, 0, 0, 0 ) 
     349         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     350      END_2D 
     351      DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 
     352         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) 
     353      END_3D 
     354      DO_2D( 0, 0, 0, 0 ) 
     355         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
     356      END_2D 
     357      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 ) 
     358         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
     359      END_3D 
     360      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     361         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
     362      END_3D 
    399363      ! 
    400364      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    402366      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    403367!!gm BUG : in the exp  remove the depth of ssh !!! 
    404 !!gm       i.e. use gde3w in argument (pdepw) 
     368!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 
    405369       
    406370       
    407371      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 
     372         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     373            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
     374               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     375         END_3D 
    416376      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 
     377         DO_2D( 0, 0, 0, 0 ) 
     378            jk = nmln(ji,jj) 
     379            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
     380               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     381         END_2D 
    424382      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 
     383         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     384            ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
     385            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
     386            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
     387            zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
     388            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
     389            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
     390               &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     391         END_3D 
    438392      ENDIF 
    439393      ! 
     
    441395 
    442396 
    443    SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 
     397   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt ) 
    444398      !!---------------------------------------------------------------------- 
    445399      !!                   ***  ROUTINE tke_avn  *** 
     
    477431      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics 
    478432      !! 
    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) 
     433      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    481434      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    482435      ! 
     
    484437      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars 
    485438      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      - 
    486       REAL(wp) ::   zemxl, zemlm, zemlp       !   -      - 
     439      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      - 
    487440      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace 
    488441      !!-------------------------------------------------------------------- 
     
    498451      zmxld(:,:,:)  = rmxl_min 
    499452      ! 
    500       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    501          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 
    507       ELSE  
     453     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
     454         ! 
     455         zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 
     456#if ! defined key_si3 && ! defined key_cice 
     457         DO_2D( 0, 0, 0, 0 ) 
     458            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1) 
     459         END_2D 
     460#else 
     461         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     462         ! 
     463         CASE( 0 )                      ! No scaling under sea-ice 
     464            DO_2D( 0, 0, 0, 0 ) 
     465               zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 
     466            END_2D 
     467            ! 
     468         CASE( 1 )                           ! scaling with constant sea-ice thickness 
     469            DO_2D( 0, 0, 0, 0 ) 
     470               zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 
     471            END_2D 
     472            ! 
     473         CASE( 2 )                                 ! scaling with mean sea-ice thickness 
     474            DO_2D( 0, 0, 0, 0 ) 
     475#if defined key_si3 
     476               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1) 
     477#elif defined key_cice 
     478               zmaxice = MAXVAL( h_i(ji,jj,:) ) 
     479               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     480#endif 
     481            END_2D 
     482            ! 
     483         CASE( 3 )                                 ! scaling with max sea-ice thickness 
     484            DO_2D( 0, 0, 0, 0 ) 
     485               zmaxice = MAXVAL( h_i(ji,jj,:) ) 
     486               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     487            END_2D 
     488            ! 
     489         END SELECT 
     490#endif 
     491         ! 
     492         DO_2D( 0, 0, 0, 0 ) 
     493            zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 
     494         END_2D 
     495         ! 
     496      ELSE 
    508497         zmxlm(:,:,1) = rn_mxl0 
    509498      ENDIF 
    510       ! 
    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 
     499 
     500      ! 
     501      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     502         zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
     503         zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
     504      END_3D 
    519505      ! 
    520506      !                     !* Physical limits for the mixing length 
     
    526512      ! 
    527513 !!gm Not sure of that coding for ISF.... 
    528       ! where wmask = 0 set zmxlm == p_e3w 
     514      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 
    529515      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 
     516         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     517            zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   & 
     518            &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) 
     519            ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
     520            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   & 
     521               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 
     522            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   & 
     523               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 
     524         END_3D 
    541525         ! 
    542526      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 
     527         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     528            zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 
     529            zmxlm(ji,jj,jk) = zemxl 
     530            zmxld(ji,jj,jk) = zemxl 
     531         END_3D 
    552532         ! 
    553533      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 
     534         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     535            zmxlm(ji,jj,jk) =   & 
     536               &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
     537         END_3D 
     538         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) 
     539            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
     540            zmxlm(ji,jj,jk) = zemxl 
     541            zmxld(ji,jj,jk) = zemxl 
     542         END_3D 
    570543         ! 
    571544      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 
     545         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     546            zmxld(ji,jj,jk) =    & 
     547               &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
     548         END_3D 
     549         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) 
     550            zmxlm(ji,jj,jk) =   & 
     551               &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
     552         END_3D 
     553         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     554            zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) ) 
     555            zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 
     556            zmxlm(ji,jj,jk) = zemlm 
     557            zmxld(ji,jj,jk) = zemlp 
     558         END_3D 
    596559         ! 
    597560      END SELECT 
     
    600563      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt) 
    601564      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    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 
     565      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     566         zsqen = SQRT( en(ji,jj,jk) ) 
     567         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     568         p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     569         p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     570         dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
     571      END_3D 
    613572      ! 
    614573      ! 
    615574      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 
     575         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     576            p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     577         END_3D 
     578      ENDIF 
     579      ! 
     580      IF(sn_cfctl%l_prtctl) THEN 
    626581         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk) 
    627582         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk ) 
     
    631586 
    632587 
    633    SUBROUTINE zdf_tke_init 
     588   SUBROUTINE zdf_tke_init( Kmm ) 
    634589      !!---------------------------------------------------------------------- 
    635590      !!                  ***  ROUTINE zdf_tke_init  *** 
     
    647602      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag 
    648603      !! 
    649       INTEGER ::   ji, jj, jk   ! dummy loop indices 
    650       INTEGER ::   ios 
    651       !! 
    652       NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,          & 
    653          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,          & 
    654          &                 rn_mxl0 , nn_pdl   , ln_drg , ln_lc    , rn_lc,   & 
    655          &                 nn_etau , nn_htau  , rn_efr , rn_eice   
    656       !!---------------------------------------------------------------------- 
    657       ! 
    658       REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy 
     604      INTEGER, INTENT(in) ::   Kmm          ! time level index 
     605      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     606      INTEGER             ::   ios 
     607      !! 
     608      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  & 
     609         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  & 
     610         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             & 
     611         &                 nn_pdl  , ln_lc    , rn_lc    ,             & 
     612         &                 nn_etau , nn_htau  , rn_efr   , rn_eice   
     613      !!---------------------------------------------------------------------- 
     614      ! 
    659615      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) 
    660 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp ) 
    661  
    662       REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy 
     616901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' ) 
     617 
    663618      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 ) 
    664 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp ) 
     619902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' ) 
    665620      IF(lwm) WRITE ( numond, namzdf_tke ) 
    666621      ! 
     
    681636         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl 
    682637         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0 
     638         IF( ln_mxl0 ) THEN 
     639            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice 
     640            IF( nn_mxlice == 1 ) & 
     641            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice 
     642         ENDIF          
    683643         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0 
    684          WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg 
    685644         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc 
    686645         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc 
     
    690649         WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice 
    691650         WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   ' 
    692          IF( ln_drg ) THEN 
     651         IF( .NOT.ln_drg_OFF ) THEN 
    693652            WRITE(numout,*) 
    694653            WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 
     
    725684      ENDIF 
    726685       
    727       IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     686      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln  
    728687 
    729688      !                               !* depth of penetration of surface tke 
     
    777736            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist 
    778737               IF(lrxios) CALL iom_swap( TRIM(crxios_context) )  
    779                CALL iom_get( numror, jpdom_autoglo, 'en'   , en   , ldxios = lrxios ) 
    780                CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios ) 
    781                CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios ) 
    782                CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl, ldxios = lrxios ) 
     738               CALL iom_get( numror, jpdom_auto, 'en'   , en   , ldxios = lrxios ) 
     739               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k, ldxios = lrxios ) 
     740               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k, ldxios = lrxios ) 
     741               CALL iom_get( numror, jpdom_auto, 'dissl', dissl, ldxios = lrxios ) 
    783742               IF(lrxios) CALL iom_swap( TRIM(cxios_context) ) 
    784743            ELSE                                          ! start TKE from rest 
Note: See TracChangeset for help on using the changeset viewer.