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 13998 for NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2020-12-02T14:55:21+01:00 (3 years ago)
Author:
techene
Message:

branch updated with trunk 13787

Location:
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13292        sette 
         10^/utils/CI/sette@13559        sette 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdftke.F90

    r13295 r13998  
    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 
     
    6868   !                      !!** Namelist  namzdf_tke  ** 
    6969   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not 
     70   INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 
     71   REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice 
    7072   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3) 
    7173   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) 
    7474   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1) 
    7575   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
     
    7979   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2] 
    8080   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it) 
    81    LOGICAL  ::   ln_drg    ! top/bottom friction forcing flag  
    8281   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3) 
    8382   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
    8483   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
    85    REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4    
    8684   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8785   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     86   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)    
    8887 
    8988   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
     
    200199      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    201200      ! 
    202       INTEGER ::   ji, jj, jk              ! dummy loop arguments 
     201      INTEGER ::   ji, jj, jk                  ! dummy loop arguments 
    203202      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars 
    204203      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3 
    205204      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient 
    206       REAL(wp) ::   zbbrau, zri                ! local scalars 
    207       REAL(wp) ::   zfact1, zfact2, zfact3     !   -         - 
    208       REAL(wp) ::   ztx2  , zty2  , zcof       !   -         - 
    209       REAL(wp) ::   ztau  , zdif               !   -         - 
    210       REAL(wp) ::   zus   , zwlc  , zind       !   -         - 
    211       REAL(wp) ::   zzd_up, zzd_lw             !   -         - 
     205      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars 
     206      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      - 
     207      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      - 
     208      REAL(wp) ::   ztau  , zdif               !   -      - 
     209      REAL(wp) ::   zus   , zwlc  , zind       !   -      - 
     210      REAL(wp) ::   zzd_up, zzd_lw             !   -      - 
    212211      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
    213       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i 
     212      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3 
    214213      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    215214      !!-------------------------------------------------------------------- 
    216215      ! 
    217       zbbrau = rn_ebb / rho0       ! Local constant initialisation 
    218       zfact1 = -.5_wp * rn_Dt  
    219       zfact2 = 1.5_wp * rn_Dt * rn_ediss 
    220       zfact3 = 0.5_wp       * rn_ediss 
     216      zbbrau  = rn_ebb / rho0       ! Local constant initialisation 
     217      zbbirau = 3.75_wp / rho0 
     218      zfact1  = -.5_wp * rn_Dt  
     219      zfact2  = 1.5_wp * rn_Dt * rn_ediss 
     220      zfact3  = 0.5_wp         * rn_ediss 
     221      ! 
     222      ! ice fraction considered for attenuation of langmuir & wave breaking 
     223      SELECT CASE ( nn_eice ) 
     224      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     225      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     226      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     227      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     228      END SELECT 
    221229      ! 
    222230      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    223231      !                     !  Surface/top/bottom boundary condition on tke 
    224232      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    225       !  
    226       DO_2D( 0, 0, 0, 0 ) 
     233      ! 
     234      DO_2D( 0, 0, 0, 0 )         ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
     235!! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 
     236!!       one way around would be to increase zbbirau  
     237!!          en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 
     238!!             &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 
    227239         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    228240      END_2D 
     
    236248      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 
    237249      ! 
    238       IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE 
    239          ! 
    240          DO_2D( 0, 0, 0, 0 ) 
     250      IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE 
     251         ! 
     252         DO_2D( 0, 0, 0, 0 )        ! bottom friction 
    241253            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    242254            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
     
    246258            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
    247259         END_2D 
    248          IF( ln_isfcav ) THEN       ! top friction 
    249             DO_2D( 0, 0, 0, 0 ) 
     260         IF( ln_isfcav ) THEN 
     261            DO_2D( 0, 0, 0, 0 )     ! top friction 
    250262               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    251263               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
     
    274286         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    275287         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    276          DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 
    277             zus  = zcof * taum(ji,jj) 
     288         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! Last w-level at which zpelc>=0.5*us*us  
     289            zus = zcof * taum(ji,jj)          !      with us=0.016*wind(starting from jpk-1) 
    278290            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
    279291         END_3D 
     
    285297         DO_2D( 0, 0, 0, 0 ) 
    286298            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    287             zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    288             IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     299            zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    289300         END_2D 
    290          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    291             IF ( zfr_i(ji,jj) /= 0. ) THEN                
    292                ! vertical velocity due to LC    
     301         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en 
     302            IF ( zus3(ji,jj) /= 0._wp ) THEN                
    293303               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    294304                  !                                           ! vertical velocity due to LC 
    295                   zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     305                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) 
    296306                  !                                           ! TKE Langmuir circulation source term 
    297                   en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     307                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    298308               ENDIF 
    299309            ENDIF 
     
    309319      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 
    310320      ! 
    311       IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri ) 
     321      IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri ) 
    312322         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    313323            !                             ! local Richardson number 
     
    322332      ENDIF 
    323333      !          
    324       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     334      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en 
    325335         zcof   = zfact1 * tmask(ji,jj,jk) 
    326336         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 
    327337         !                                   ! eddy coefficient (ensure numerical stability) 
    328338         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
    329             &          /    (  e3t(ji,jj,jk  ,Kmm)   & 
    330             &                * e3w(ji,jj,jk  ,Kmm)  ) 
     339            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    331340         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
    332             &          /    (  e3t(ji,jj,jk-1,Kmm)   & 
    333             &                * e3w(ji,jj,jk  ,Kmm)  ) 
     341            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    334342         ! 
    335343         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     
    344352      END_3D 
    345353      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    346       DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 
     354      DO_3D( 0, 0, 0, 0, 3, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    347355         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    348356      END_3D 
    349       DO_2D( 0, 0, 0, 0 ) 
     357      DO_2D( 0, 0, 0, 0 )                          ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    350358         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    351359      END_2D 
     
    353361         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) 
    354362      END_3D 
    355       DO_2D( 0, 0, 0, 0 ) 
     363      DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    356364         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    357365      END_2D 
     
    359367         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    360368      END_3D 
    361       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     369      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke 
    362370         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    363371      END_3D 
     
    368376!!gm BUG : in the exp  remove the depth of ssh !!! 
    369377!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 
    370        
    371        
     378      ! 
     379      ! penetration is partly switched off below sea-ice if nn_eice/=0 
     380      ! 
    372381      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    373          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     382         DO_3D( 0, 0, 0, 0, 2, jpkm1 )  
    374383            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    375                &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     384               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    376385         END_3D 
    377386      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction) 
     
    379388            jk = nmln(ji,jj) 
    380389            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    381                &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     390               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    382391         END_2D 
    383392      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability) 
     
    389398            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    390399            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    391                &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     400               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    392401         END_3D 
    393402      ENDIF 
     
    451460      zmxlm(:,:,:)  = rmxl_min     
    452461      zmxld(:,:,:)  = rmxl_min 
    453       ! 
     462      !  
    454463     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
    455464         ! 
    456465         zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 
    457466#if ! defined key_si3 && ! defined key_cice 
    458          DO_2D( 0, 0, 0, 0 ) 
     467         DO_2D( 0, 0, 0, 0 )                  ! No sea-ice 
    459468            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1) 
    460469         END_2D 
     
    467476            END_2D 
    468477            ! 
    469          CASE( 1 )                           ! scaling with constant sea-ice thickness 
     478         CASE( 1 )                      ! scaling with constant sea-ice thickness 
    470479            DO_2D( 0, 0, 0, 0 ) 
    471                zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 
     480               zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     481                  &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1) 
    472482            END_2D 
    473483            ! 
    474          CASE( 2 )                                 ! scaling with mean sea-ice thickness 
     484         CASE( 2 )                      ! scaling with mean sea-ice thickness 
    475485            DO_2D( 0, 0, 0, 0 ) 
    476486#if defined key_si3 
    477                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) 
     487               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     488                  &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 
    478489#elif defined key_cice 
    479490               zmaxice = MAXVAL( h_i(ji,jj,:) ) 
    480                zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     491               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     492                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    481493#endif 
    482494            END_2D 
    483495            ! 
    484          CASE( 3 )                                 ! scaling with max sea-ice thickness 
     496         CASE( 3 )                      ! scaling with max sea-ice thickness 
    485497            DO_2D( 0, 0, 0, 0 ) 
    486498               zmaxice = MAXVAL( h_i(ji,jj,:) ) 
    487                zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     499               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     500                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    488501            END_2D 
    489502            ! 
     
    533546         ! 
    534547      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    535          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     548         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : 
    536549            zmxlm(ji,jj,jk) =   & 
    537550               &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
    538551         END_3D 
    539          DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) 
     552         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : 
    540553            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
    541554            zmxlm(ji,jj,jk) = zemxl 
     
    544557         ! 
    545558      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    546          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     559         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup 
    547560            zmxld(ji,jj,jk) =    & 
    548561               &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
    549562         END_3D 
    550          DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) 
     563         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown 
    551564            zmxlm(ji,jj,jk) =   & 
    552565               &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
     
    564577      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt) 
    565578      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    566       DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     579      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points 
    567580         zsqen = SQRT( en(ji,jj,jk) ) 
    568581         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     
    573586      ! 
    574587      ! 
    575       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
     588      IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt 
    576589         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    577590            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) 
     
    610623         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  & 
    611624         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             & 
    612          &                 nn_pdl  , ln_drg   , ln_lc    , rn_lc,      & 
    613          &                 nn_etau , nn_htau  , rn_efr   , rn_eice   
     625         &                 nn_pdl  , ln_lc    , rn_lc    ,             & 
     626         &                 nn_etau , nn_htau  , rn_efr   , nn_eice   
    614627      !!---------------------------------------------------------------------- 
    615628      ! 
     
    637650         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl 
    638651         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0 
     652         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0 
    639653         IF( ln_mxl0 ) THEN 
    640654            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice 
    641655            IF( nn_mxlice == 1 ) & 
    642656            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice 
    643          ENDIF          
    644          WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0 
    645          WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg 
     657            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     658            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice' 
     659            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness' 
     660            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness' 
     661            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness' 
     662            CASE DEFAULT 
     663               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 
     664            END SELECT 
     665         ENDIF 
    646666         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc 
    647667         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc 
     
    649669         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau 
    650670         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr 
    651          WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice 
    652          WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   ' 
    653          IF( ln_drg ) THEN 
    654             WRITE(numout,*) 
    655             WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 
    656             WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top 
    657             WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot 
    658          ENDIF 
     671         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice 
     672         SELECT CASE( nn_eice )  
     673         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking' 
     674         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     675         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)' 
     676         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     677         CASE DEFAULT 
     678            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 
     679         END SELECT       
    659680         WRITE(numout,*) 
    660681         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri 
Note: See TracChangeset for help on using the changeset viewer.