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 13472 for NEMO/trunk/src/OCE/ZDF – NEMO

Ignore:
Timestamp:
2020-09-16T15:05:19+02:00 (4 years ago)
Author:
smasson
Message:

trunk: commit changes from r4.0-HEAD from 13284 to 13449, see #2523

Location:
NEMO/trunk/src/OCE/ZDF
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/ZDF/zdfdrg.F90

    r13461 r13472  
    3232   USE lib_mpp        ! distributed memory computing 
    3333   USE prtctl         ! Print control 
     34   USE sbc_oce , ONLY : nn_ice  
    3435 
    3536   IMPLICIT NONE 
     
    4647   LOGICAL          ::   ln_loglayer  ! logarithmic drag: Cd = vkarmn/log(z/z0) 
    4748   LOGICAL , PUBLIC ::   ln_drgimp    ! implicit top/bottom friction flag 
    48  
     49   LOGICAL , PUBLIC ::   ln_drgice_imp ! implicit ice-ocean drag  
    4950   !                                 !!* Namelist namdrg_top & _bot: TOP or BOTTOM coefficient namelist * 
    5051   REAL(wp)         ::   rn_Cd0       !: drag coefficient                                           [ - ] 
     
    226227      INTEGER   ::   ios, ioptio   ! local integers 
    227228      !! 
    228       NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp 
     229      NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp, ln_drgice_imp 
    229230      !!---------------------------------------------------------------------- 
    230231      ! 
     
    247248         WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer 
    248249         WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp 
     250         WRITE(numout,*) '      implicit ice-ocean drag                   ln_drgice_imp  =', ln_drgice_imp 
    249251      ENDIF 
    250252      ! 
     
    257259      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' ) 
    258260      ! 
     261      IF ( ln_drgice_imp.AND.(.NOT.ln_drgimp) ) &  
     262         &                CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires ln_drgimp=T' ) 
     263      ! 
     264      IF ( ln_drgice_imp.AND.( nn_ice /=2 ) ) & 
     265         &  CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires si3' ) 
    259266      ! 
    260267      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor) 
     
    267274      !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities) 
    268275      ! 
    269       IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting 
    270          ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) ) 
     276      IF( ln_isfcav.OR.ln_drgice_imp ) THEN              ! Ocean cavities: top friction setting 
     277         ALLOCATE( rCdU_top(jpi,jpj) ) 
     278      ENDIF 
     279      ! 
     280      IF( ln_isfcav ) THEN 
     281         ALLOCATE( rCd0_top(jpi,jpj)) 
    271282         CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in 
    272283            &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out 
  • NEMO/trunk/src/OCE/ZDF/zdfgls.F90

    r13461 r13472  
    5454   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1) 
    5555   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation 
     56   INTEGER  ::   nn_z0_ice         ! Roughness accounting for sea ice 
    5657   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2) 
    5758   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen 
     
    6263   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
    6364   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
     65   REAL(wp) ::   rn_hsri           ! Ice ocean roughness 
    6466   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)  
    6567 
     
    153155      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves  
    154156      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves) 
     157      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice 
    155158      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before 
    156159      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before 
     
    168171      ustar2_bot (:,:) = 0._wp 
    169172 
     173      SELECT CASE ( nn_z0_ice ) 
     174      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     175      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     176      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     177      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     178      END SELECT 
     179       
    170180      ! Compute surface, top and bottom friction at T-points 
    171181      DO_2D( 0, 0, 0, 0 ) 
     
    207217      END SELECT 
    208218      ! 
     219      ! adapt roughness where there is sea ice 
     220      zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro 
     221      ! 
    209222      DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 
    210223         eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
     
    291304      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)  
    292305      ! First level 
    293       en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  ) 
     306      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  ) 
    294307      zd_lw(:,:,1) = en(:,:,1) 
    295308      zd_up(:,:,1) = 0._wp 
     
    297310      !  
    298311      ! One level below 
    299       en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm))  & 
    300          &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   ) 
     312      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 
     313         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   ) 
    301314      zd_lw(:,:,2) = 0._wp  
    302315      zd_up(:,:,2) = 0._wp 
     
    307320      ! 
    308321      ! Dirichlet conditions at k=1 
    309       en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  ) 
     322      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  ) 
    310323      zd_lw(:,:,1) = en(:,:,1) 
    311324      zd_up(:,:,1) = 0._wp 
     
    317330      zd_lw(:,:,2) = 0._wp 
    318331      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 
    319       zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
     332      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    320333          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
    321334!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 
     
    530543         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 
    531544         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 
    532          zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     545         zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) & 
     546            &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    533547         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 
    534548            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 
     
    753767      REAL(wp)::   zcr   ! local scalar 
    754768      !! 
    755       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    756          &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    757          &            rn_crban, rn_charn, rn_frac_hs,        & 
    758          &            nn_bc_surf, nn_bc_bot, nn_z0_met,     & 
     769      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       & 
     770         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   & 
     771         &            rn_crban, rn_charn, rn_frac_hs,              & 
     772         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, & 
    759773         &            nn_stab_func, nn_clos 
    760774      !!---------------------------------------------------------- 
     
    782796         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    783797         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     798         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice 
     799         SELECT CASE( nn_z0_ice ) 
     800         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking' 
     801         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     802         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)' 
     803         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     804         CASE DEFAULT 
     805            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3') 
     806         END SELECT 
    784807         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs 
    785808         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func 
    786809         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    787810         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
     811         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri 
    788812         WRITE(numout,*) 
    789813         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 
  • NEMO/trunk/src/OCE/ZDF/zdfphy.F90

    r13226 r13472  
    2828   USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test) 
    2929   USE sbcrnf         ! surface boundary condition: runoff variables 
     30   USE sbc_ice        ! sea ice drag 
    3031#if defined key_agrif 
    3132   USE agrif_oce_interp   ! interpavm 
     
    253254      ENDIF 
    254255      ! 
     256#if defined key_si3 
     257      IF ( ln_drgice_imp) THEN 
     258         IF ( ln_isfcav ) THEN 
     259            rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:) 
     260         ELSE 
     261            rCdU_top(:,:) = rCdU_ice(:,:) 
     262         ENDIF 
     263      ENDIF 
     264#endif 
     265      !  
    255266      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k) 
    256267      ! 
  • NEMO/trunk/src/OCE/ZDF/zdftke.F90

    r13461 r13472  
    5252#endif 
    5353   ! 
     54#if defined key_si3 
     55   USE ice, ONLY: hm_i, h_i 
     56#endif 
     57#if defined key_cice 
     58   USE sbc_ice, ONLY: h_i 
     59#endif 
    5460   USE in_out_manager ! I/O manager 
    5561   USE iom            ! I/O manager library 
     
    6874   !                      !!** Namelist  namzdf_tke  ** 
    6975   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not 
     76   INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 
     77   REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice 
    7078   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3) 
    7179   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) 
    7480   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1) 
    7581   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
     
    8288   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
    8389   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
    84    REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4    
    8590   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8691   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     92   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)    
    8793 
    8894   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
     
    199205      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    200206      ! 
    201       INTEGER ::   ji, jj, jk              ! dummy loop arguments 
     207      INTEGER ::   ji, jj, jk                  ! dummy loop arguments 
    202208      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars 
    203209      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3 
    204210      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient 
    205       REAL(wp) ::   zbbrau, 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             !   -         - 
     211      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars 
     212      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      - 
     213      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      - 
     214      REAL(wp) ::   ztau  , zdif               !   -      - 
     215      REAL(wp) ::   zus   , zwlc  , zind       !   -      - 
     216      REAL(wp) ::   zzd_up, zzd_lw             !   -      - 
    211217      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
    212       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i 
     218      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3 
    213219      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    214220      !!-------------------------------------------------------------------- 
    215221      ! 
    216       zbbrau = rn_ebb / rho0       ! Local constant initialisation 
    217       zfact1 = -.5_wp * rn_Dt  
    218       zfact2 = 1.5_wp * rn_Dt * rn_ediss 
    219       zfact3 = 0.5_wp       * rn_ediss 
     222      zbbrau  = rn_ebb / rho0       ! Local constant initialisation 
     223      zbbirau = 3.75_wp / rho0 
     224      zfact1  = -.5_wp * rn_Dt  
     225      zfact2  = 1.5_wp * rn_Dt * rn_ediss 
     226      zfact3  = 0.5_wp         * rn_ediss 
     227      ! 
     228      ! ice fraction considered for attenuation of langmuir & wave breaking 
     229      SELECT CASE ( nn_eice ) 
     230      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     231      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     232      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     233      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     234      END SELECT 
    220235      ! 
    221236      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    222237      !                     !  Surface/top/bottom boundary condition on tke 
    223238      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    224       !  
     239      ! 
    225240      DO_2D( 0, 0, 0, 0 ) 
     241!! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 
     242!!       one way around would be to increase zbbirau  
     243!!          en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 
     244!!             &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 
    226245         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    227246      END_2D 
     
    273292         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    274293         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    275          DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 
    276             zus  = zcof * taum(ji,jj) 
     294         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! Last w-level at which zpelc>=0.5*us*us  
     295            zus = zcof * taum(ji,jj)          !      with us=0.016*wind(starting from jpk-1) 
    277296            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
    278297         END_3D 
     
    284303         DO_2D( 0, 0, 0, 0 ) 
    285304            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. 
     305            zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    288306         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    
     307         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en 
     308            IF ( zus3(ji,jj) /= 0._wp ) THEN                
    292309               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    293310                  !                                           ! 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 
     311                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) 
    295312                  !                                           ! 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) 
     313                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    297314               ENDIF 
    298315            ENDIF 
     
    326343         !                                   ! eddy coefficient (ensure numerical stability) 
    327344         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)  ) 
     345            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    330346         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)  ) 
     347            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    333348         ! 
    334349         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     
    370385       
    371386      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    372          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     387         DO_3D( 0, 0, 0, 0, 2, jpkm1 )     ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF  
    373388            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) 
     389               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    375390         END_3D 
    376391      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction) 
     
    378393            jk = nmln(ji,jj) 
    379394            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) 
     395               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    381396         END_2D 
    382397      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability) 
     
    388403            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    389404            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) 
     405               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    391406         END_3D 
    392407      ENDIF 
     
    450465      zmxlm(:,:,:)  = rmxl_min     
    451466      zmxld(:,:,:)  = rmxl_min 
    452       ! 
     467      !  
    453468     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
    454469         ! 
     
    468483         CASE( 1 )                           ! scaling with constant sea-ice thickness 
    469484            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) 
     485               zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     486                  &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1) 
    471487            END_2D 
    472488            ! 
     
    474490            DO_2D( 0, 0, 0, 0 ) 
    475491#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) 
     492               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     493                  &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 
    477494#elif defined key_cice 
    478495               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) 
     496               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     497                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    480498#endif 
    481499            END_2D 
     
    484502            DO_2D( 0, 0, 0, 0 ) 
    485503               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) 
     504               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     505                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    487506            END_2D 
    488507            ! 
     
    610629         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             & 
    611630         &                 nn_pdl  , ln_lc    , rn_lc    ,             & 
    612          &                 nn_etau , nn_htau  , rn_efr   , rn_eice   
     631         &                 nn_etau , nn_htau  , rn_efr   , nn_eice   
    613632      !!---------------------------------------------------------------------- 
    614633      ! 
     
    642661         ENDIF          
    643662         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0 
     663         IF( ln_mxl0 ) THEN 
     664            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice 
     665            IF( nn_mxlice == 1 ) & 
     666            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice 
     667            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     668            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice' 
     669            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness' 
     670            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness' 
     671            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness' 
     672            CASE DEFAULT 
     673               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 
     674            END SELECT 
     675         ENDIF 
    644676         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc 
    645677         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc 
     
    647679         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau 
    648680         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr 
    649          WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice 
    650          WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   ' 
     681         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice 
     682         SELECT CASE( nn_eice )  
     683         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking' 
     684         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     685         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)' 
     686         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     687         CASE DEFAULT 
     688            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 
     689         END SELECT       
    651690         IF( .NOT.ln_drg_OFF ) THEN 
    652691            WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.