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.
zdftke.F90 in NEMO/trunk/src/OCE/ZDF – NEMO

source: NEMO/trunk/src/OCE/ZDF/zdftke.F90 @ 13530

Last change on this file since 13530 was 13530, checked in by clem, 3 years ago

remove duplicated lines

  • Property svn:keywords set to Id
File size: 44.5 KB
RevLine 
[1531]1MODULE zdftke
[1239]2   !!======================================================================
[1531]3   !!                       ***  MODULE  zdftke  ***
[1239]4   !! Ocean physics:  vertical mixing coefficient computed from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
[1492]7   !! History :  OPA  !  1991-03  (b. blanke)  Original code
8   !!            7.0  !  1991-11  (G. Madec)   bug fix
9   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
10   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
11   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
12   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
13   !!            8.0  !  1997-07  (G. Madec)   lbc
14   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
15   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
16   !!             -   !  2004-10  (C. Ethe )  1D configuration
17   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
18   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
19   !!                 !           - tke penetration (wind steering)
20   !!                 !           - suface condition for tke & mixing length
21   !!                 !           - Langmuir cells
22   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
23   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
26   !!                 !                                + cleaning of the parameters + bugs correction
[2528]27   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
[5120]28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
[9019]29   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
[13461]30   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition
[1239]31   !!----------------------------------------------------------------------
[9019]32
[1239]33   !!----------------------------------------------------------------------
[3625]34   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
35   !!   tke_tke       : tke time stepping: update tke at now time step (en)
36   !!   tke_avn       : compute mixing length scale and deduce avm and avt
37   !!   zdf_tke_init  : initialization, namelist read, and parameters control
38   !!   tke_rst       : read/write tke restart in ocean restart file
[1239]39   !!----------------------------------------------------------------------
[2528]40   USE oce            ! ocean: dynamics and active tracers variables
41   USE phycst         ! physical constants
42   USE dom_oce        ! domain: ocean
43   USE domvvl         ! domain: variable volume layer
[1492]44   USE sbc_oce        ! surface boundary condition: ocean
[9019]45   USE zdfdrg         ! vertical physics: top/bottom drag coef.
[2528]46   USE zdfmxl         ! vertical physics: mixed layer
[13007]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
[9019]53   !
[1492]54   USE in_out_manager ! I/O manager
55   USE iom            ! I/O manager library
[2715]56   USE lib_mpp        ! MPP library
[9019]57   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
58   USE prtctl         ! Print control
[3625]59   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[1239]60
61   IMPLICIT NONE
62   PRIVATE
63
[2528]64   PUBLIC   zdf_tke        ! routine called in step module
65   PUBLIC   zdf_tke_init   ! routine called in opa module
66   PUBLIC   tke_rst        ! routine called in step module
[1239]67
[4147]68   !                      !!** Namelist  namzdf_tke  **
69   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
[13472]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
[4147]72   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
73   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
74   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
75   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
76   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
77   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
78   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
79   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
80   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
81   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
[9019]82   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
83   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
[4147]84   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
[9019]85   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells
[13472]86   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)   
[1239]87
[4147]88   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
89   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
[2528]90   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
91   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
[1239]92
[9019]93   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau)
94   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation
95   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation
[1492]96
[1239]97   !! * Substitutions
[12377]98#  include "do_loop_substitute.h90"
[13237]99#  include "domzgr_substitute.h90"
[1239]100   !!----------------------------------------------------------------------
[9598]101   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]102   !! $Id$
[10068]103   !! Software governed by the CeCILL license (see ./LICENSE)
[1239]104   !!----------------------------------------------------------------------
105CONTAINS
106
[2715]107   INTEGER FUNCTION zdf_tke_alloc()
108      !!----------------------------------------------------------------------
109      !!                ***  FUNCTION zdf_tke_alloc  ***
110      !!----------------------------------------------------------------------
[9019]111      ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) ,   STAT= zdf_tke_alloc )
112      !
[10425]113      CALL mpp_sum ( 'zdftke', zdf_tke_alloc )
114      IF( zdf_tke_alloc /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_alloc: failed to allocate arrays' )
[2715]115      !
116   END FUNCTION zdf_tke_alloc
117
118
[12377]119   SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
[1239]120      !!----------------------------------------------------------------------
[1531]121      !!                   ***  ROUTINE zdf_tke  ***
[1239]122      !!
123      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
[1492]124      !!              coefficients using a turbulent closure scheme (TKE).
[1239]125      !!
[1492]126      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
127      !!              is computed from a prognostic equation :
128      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
129      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
130      !!                  + avt N^2                      ! stratif. destruc.
131      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
[1239]132      !!      with the boundary conditions:
[1695]133      !!         surface: en = max( rn_emin0, rn_ebb * taum )
[1239]134      !!         bottom : en = rn_emin
[1492]135      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
136      !!
137      !!        The now Turbulent kinetic energy is computed using the following
138      !!      time stepping: implicit for vertical diffusion term, linearized semi
139      !!      implicit for kolmogoroff dissipation term, and explicit forward for
140      !!      both buoyancy and shear production terms. Therefore a tridiagonal
141      !!      linear system is solved. Note that buoyancy and shear terms are
142      !!      discretized in a energy conserving form (Bruchard 2002).
143      !!
144      !!        The dissipative and mixing length scale are computed from en and
145      !!      the stratification (see tke_avn)
146      !!
147      !!        The now vertical eddy vicosity and diffusivity coefficients are
148      !!      given by:
149      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
150      !!              avt = max( avmb, pdl * avm                 ) 
[1239]151      !!              eav = max( avmb, avm )
[1492]152      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
153      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
[1239]154      !!
155      !! ** Action  :   compute en (now turbulent kinetic energy)
[9019]156      !!                update avt, avm (before vertical eddy coef.)
[1239]157      !!
158      !! References : Gaspar et al., JGR, 1990,
159      !!              Blanke and Delecluse, JPO, 1991
160      !!              Mellor and Blumberg, JPO 2004
161      !!              Axell, JGR, 2002
[1492]162      !!              Bruchard OM 2002
[1239]163      !!----------------------------------------------------------------------
[9019]164      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
[12377]165      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
[9019]166      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
167      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
[1492]168      !!----------------------------------------------------------------------
[1481]169      !
[12377]170      CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )   ! now tke (en)
[5656]171      !
[12377]172      CALL tke_avn( Kbb, Kmm,        p_avm, p_avt )   ! now avt, avm, dissl
[3632]173      !
[5656]174  END SUBROUTINE zdf_tke
[1239]175
[1492]176
[12377]177   SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )
[1239]178      !!----------------------------------------------------------------------
[1492]179      !!                   ***  ROUTINE tke_tke  ***
180      !!
181      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
182      !!
183      !! ** Method  : - TKE surface boundary condition
[2528]184      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
[9019]185      !!              - source term due to shear (= Kz dz[Ub] * dz[Un] )
[1492]186      !!              - Now TKE : resolution of the TKE equation by inverting
187      !!                a tridiagonal linear system by a "methode de chasse"
188      !!              - increase TKE due to surface and internal wave breaking
[9019]189      !!             NB: when sea-ice is present, both LC parameterization
190      !!                 and TKE penetration are turned off when the ice fraction
191      !!                 is smaller than 0.25
[1492]192      !!
193      !! ** Action  : - en : now turbulent kinetic energy)
[1239]194      !! ---------------------------------------------------------------------
[9019]195      USE zdf_oce , ONLY : en   ! ocean vertical physics
196      !!
[12377]197      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
[9019]198      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term
199      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
200      !
[13472]201      INTEGER ::   ji, jj, jk                  ! dummy loop arguments
[9019]202      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars
203      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3
204      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient
[13472]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             !   -      -
[9019]211      INTEGER , DIMENSION(jpi,jpj)     ::   imlc
[13472]212      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3
[9019]213      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
[1239]214      !!--------------------------------------------------------------------
[1492]215      !
[13472]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
[1492]221      !
[13472]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
229      !
[1492]230      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[9019]231      !                     !  Surface/top/bottom boundary condition on tke
[1492]232      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[13472]233      !
[13497]234      DO_2D( 0, 0, 0, 0 )         ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
[13472]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)
[12377]239         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
240      END_2D
[9019]241      !
[1492]242      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
243      !                     !  Bottom boundary condition on tke
244      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1719]245      !
[12489]246      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
[9019]247      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
248      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
[1492]249      !
[13461]250      IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE
[9019]251         !
[13461]252         DO_2D( 0, 0, 0, 0 )        ! bottom friction
[12377]253            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
254            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
[12489]255            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
[12377]256            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  &
257               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
258            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
259         END_2D
[13461]260         IF( ln_isfcav ) THEN
261            DO_2D( 0, 0, 0, 0 )     ! top friction
[12377]262               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
263               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
[12489]264               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
[12377]265               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  &
266                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
[12702]267               ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
268               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &
269                  &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
[12377]270            END_2D
[9019]271         ENDIF
272         !
273      ENDIF
274      !
[1492]275      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[9019]276      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
[1492]277         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[1239]278         !
[1492]279         !                        !* total energy produce by LC : cumulative sum over jk
[12377]280         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm)
[1239]281         DO jk = 2, jpk
[13237]282            zpelc(:,:,jk)  = zpelc(:,:,jk-1) +   &
283               &        MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
[1239]284         END DO
[1492]285         !                        !* finite Langmuir Circulation depth
[1705]286         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
[7753]287         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
[13472]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)
[12377]290            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
291         END_3D
[1492]292         !                               ! finite LC depth
[13295]293         DO_2D( 1, 1, 1, 1 )
[12377]294            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
295         END_2D
[1705]296         zcof = 0.016 / SQRT( zrhoa * zcdrag )
[13295]297         DO_2D( 0, 0, 0, 0 )
[12377]298            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
[13472]299            zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
[12377]300         END_2D
[13472]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               
[12377]303               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
304                  !                                           ! vertical velocity due to LC
[13472]305                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )
[12377]306                  !                                           ! TKE Langmuir circulation source term
[13472]307                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
[12377]308               ENDIF
309            ENDIF
310         END_3D
[1239]311         !
312      ENDIF
[1492]313      !
314      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
315      !                     !  Now Turbulent kinetic energy (output in en)
316      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
317      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
318      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
319      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
320      !
[13497]321      IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
[13295]322         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]323            !                             ! local Richardson number
[13226]324            IF (rn2b(ji,jj,jk) <= 0.0_wp) then
325                zri = 0.0_wp
326            ELSE
327                zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
328            ENDIF
[12377]329            !                             ! inverse of Prandtl number
330            apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
331         END_3D
[5656]332      ENDIF
[5836]333      !         
[13497]334      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en
[12377]335         zcof   = zfact1 * tmask(ji,jj,jk)
336         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
337         !                                   ! eddy coefficient (ensure numerical stability)
338         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
[13472]339            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
[12377]340         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
[13472]341            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
[12377]342         !
343         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
344         zd_lw(ji,jj,jk) = zzd_lw
345         zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
346         !
347         !                                   ! right hand side in en
[12489]348         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                        &   ! shear
[12377]349            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
350            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
351            &                                ) * wmask(ji,jj,jk)
352      END_3D
[5120]353      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
[13497]354      DO_3D( 0, 0, 0, 0, 3, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
[12377]355         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
356      END_3D
[13497]357      DO_2D( 0, 0, 0, 0 )                          ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
[12377]358         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
359      END_2D
[13295]360      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
[12377]361         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)
362      END_3D
[13497]363      DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
[12377]364         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
365      END_2D
[13295]366      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
[12377]367         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
368      END_3D
[13497]369      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke
[12377]370         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
371      END_3D
[9019]372      !
[1492]373      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
374      !                            !  TKE due to surface and internal wave breaking
375      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[6140]376!!gm BUG : in the exp  remove the depth of ssh !!!
[12377]377!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
[13530]378      !
379      ! penetration is partly switched off below sea-ice if nn_eice/=0
380      !
[2528]381      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
[13530]382         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
[12377]383            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
[13472]384               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[12377]385         END_3D
[2528]386      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
[13295]387         DO_2D( 0, 0, 0, 0 )
[12377]388            jk = nmln(ji,jj)
389            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
[13472]390               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[12377]391         END_2D
[2528]392      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
[13530]393         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]394            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
395            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
396            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
397            zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
398            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
399            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
[13472]400               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[12377]401         END_3D
[1239]402      ENDIF
[1492]403      !
[1239]404   END SUBROUTINE tke_tke
405
[1492]406
[12377]407   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
[1239]408      !!----------------------------------------------------------------------
[1492]409      !!                   ***  ROUTINE tke_avn  ***
[1239]410      !!
[1492]411      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
412      !!
413      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
414      !!              the tke_tke routine). First, the now mixing lenth is
415      !!      computed from en and the strafification (N^2), then the mixings
416      !!      coefficients are computed.
417      !!              - Mixing length : a first evaluation of the mixing lengh
418      !!      scales is:
419      !!                      mxl = sqrt(2*en) / N 
420      !!      where N is the brunt-vaisala frequency, with a minimum value set
[2528]421      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
[1492]422      !!        The mixing and dissipative length scale are bound as follow :
423      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
424      !!                        zmxld = zmxlm = mxl
425      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
426      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
427      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
428      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
429      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
430      !!                    the surface to obtain ldown. the resulting length
431      !!                    scales are:
432      !!                        zmxld = sqrt( lup * ldown )
433      !!                        zmxlm = min ( lup , ldown )
434      !!              - Vertical eddy viscosity and diffusivity:
435      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
436      !!                      avt = max( avmb, pdlr * avm ) 
437      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
438      !!
[9019]439      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
[1239]440      !!----------------------------------------------------------------------
[9019]441      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics
442      !!
[12377]443      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
[9019]444      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
445      !
[2715]446      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[9019]447      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
448      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
[13007]449      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
[9019]450      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
[1239]451      !!--------------------------------------------------------------------
[3294]452      !
[1492]453      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
454      !                     !  Mixing length
455      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
456      !
457      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
458      !
[5120]459      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
[7753]460      zmxlm(:,:,:)  = rmxl_min   
461      zmxld(:,:,:)  = rmxl_min
[13472]462      !
[13012]463     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
[13007]464         !
[13012]465         zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
[13007]466#if ! defined key_si3 && ! defined key_cice
[13497]467         DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
[13007]468            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
[12377]469         END_2D
[13007]470#else
471         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
472         !
473         CASE( 0 )                      ! No scaling under sea-ice
[13295]474            DO_2D( 0, 0, 0, 0 )
[13007]475               zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
476            END_2D
477            !
[13497]478         CASE( 1 )                      ! scaling with constant sea-ice thickness
[13295]479            DO_2D( 0, 0, 0, 0 )
[13472]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)
[13007]482            END_2D
483            !
[13497]484         CASE( 2 )                      ! scaling with mean sea-ice thickness
[13295]485            DO_2D( 0, 0, 0, 0 )
[13007]486#if defined key_si3
[13472]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)
[13007]489#elif defined key_cice
490               zmaxice = MAXVAL( h_i(ji,jj,:) )
[13472]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)
[13007]493#endif
494            END_2D
495            !
[13497]496         CASE( 3 )                      ! scaling with max sea-ice thickness
[13295]497            DO_2D( 0, 0, 0, 0 )
[13007]498               zmaxice = MAXVAL( h_i(ji,jj,:) )
[13472]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)
[13007]501            END_2D
502            !
503         END SELECT
504#endif
505         !
[13295]506         DO_2D( 0, 0, 0, 0 )
[13007]507            zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
508         END_2D
509         !
510      ELSE
[7753]511         zmxlm(:,:,1) = rn_mxl0
[1239]512      ENDIF
[13007]513
[1239]514      !
[13295]515      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]516         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
517         zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
518      END_3D
[1492]519      !
520      !                     !* Physical limits for the mixing length
521      !
[7753]522      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
523      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
[1492]524      !
[1239]525      SELECT CASE ( nn_mxl )
526      !
[5836]527 !!gm Not sure of that coding for ISF....
[12377]528      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
[1239]529      CASE ( 0 )           ! bounded by the distance to surface and bottom
[13295]530         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]531            zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
532            &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
533            ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
[13237]534            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
535               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
536            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
537               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
[12377]538         END_3D
[1239]539         !
540      CASE ( 1 )           ! bounded by the vertical scale factor
[13295]541         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]542            zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
543            zmxlm(ji,jj,jk) = zemxl
544            zmxld(ji,jj,jk) = zemxl
545         END_3D
[1239]546         !
547      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
[13497]548         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom :
[13237]549            zmxlm(ji,jj,jk) =   &
550               &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
[12377]551         END_3D
[13497]552         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface :
[12377]553            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
554            zmxlm(ji,jj,jk) = zemxl
555            zmxld(ji,jj,jk) = zemxl
556         END_3D
[1239]557         !
558      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
[13497]559         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup
[13237]560            zmxld(ji,jj,jk) =    &
561               &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
[12377]562         END_3D
[13497]563         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
[13237]564            zmxlm(ji,jj,jk) =   &
565               &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
[12377]566         END_3D
[13295]567         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]568            zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
569            zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
570            zmxlm(ji,jj,jk) = zemlm
571            zmxld(ji,jj,jk) = zemlp
572         END_3D
[1239]573         !
574      END SELECT
[1492]575      !
576      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[9019]577      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
[1492]578      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[13497]579      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
[12377]580         zsqen = SQRT( en(ji,jj,jk) )
581         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
582         p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
583         p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
584         dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
585      END_3D
[1492]586      !
587      !
[13497]588      IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
[13295]589         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12698]590            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)
[12377]591         END_3D
[1239]592      ENDIF
[9019]593      !
[12377]594      IF(sn_cfctl%l_prtctl) THEN
[9440]595         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
596         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk )
[1239]597      ENDIF
598      !
[1492]599   END SUBROUTINE tke_avn
[1239]600
[1492]601
[12377]602   SUBROUTINE zdf_tke_init( Kmm )
[1239]603      !!----------------------------------------------------------------------
[2528]604      !!                  ***  ROUTINE zdf_tke_init  ***
[1239]605      !!                     
606      !! ** Purpose :   Initialization of the vertical eddy diffivity and
[1492]607      !!              viscosity when using a tke turbulent closure scheme
[1239]608      !!
[1601]609      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
[1492]610      !!              called at the first timestep (nit000)
[1239]611      !!
[1601]612      !! ** input   :   Namlist namzdf_tke
[1239]613      !!
614      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
615      !!----------------------------------------------------------------------
[9019]616      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
617      !!
[12377]618      INTEGER, INTENT(in) ::   Kmm          ! time level index
619      INTEGER             ::   ji, jj, jk   ! dummy loop indices
620      INTEGER             ::   ios
[1239]621      !!
[13007]622      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  &
623         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
624         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
[13461]625         &                 nn_pdl  , ln_lc    , rn_lc    ,             &
[13472]626         &                 nn_etau , nn_htau  , rn_efr   , nn_eice 
[1239]627      !!----------------------------------------------------------------------
[2715]628      !
[4147]629      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
[11536]630901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' )
[4147]631
632      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
[11536]633902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' )
[4624]634      IF(lwm) WRITE ( numond, namzdf_tke )
[2715]635      !
[2528]636      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
[2715]637      !
[1492]638      IF(lwp) THEN                    !* Control print
[1239]639         WRITE(numout,*)
[2528]640         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
641         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]642         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
[1705]643         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
644         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
645         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
646         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
647         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
[9019]648         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
[1705]649         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
650         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
[9019]651         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
[13472]653         IF( ln_mxl0 ) THEN
654            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice
655            IF( nn_mxlice == 1 ) &
656            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice
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
[9019]666         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
667         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
[1705]668         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
[9019]669         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
670         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
[13472]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     
[13461]680         IF( .NOT.ln_drg_OFF ) THEN
[9169]681            WRITE(numout,*)
[9019]682            WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
683            WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top
684            WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot
685         ENDIF
686         WRITE(numout,*)
[9190]687         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
[9019]688         WRITE(numout,*)
[1239]689      ENDIF
[2715]690      !
[9019]691      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
692         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
693         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
[9190]694         IF(lwp) WRITE(numout,*) '   ==>>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
[9019]695      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
696         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
[9190]697         IF(lwp) WRITE(numout,*) '   ==>>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
[9019]698      ENDIF
699      !
[2715]700      !                              ! allocate tke arrays
701      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
702      !
[1492]703      !                               !* Check of some namelist values
[4990]704      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
705      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
706      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
[5407]707      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
[9019]708      !
[2528]709      IF( ln_mxl0 ) THEN
[9169]710         IF(lwp) WRITE(numout,*)
[9190]711         IF(lwp) WRITE(numout,*) '   ==>>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
[2528]712         rn_mxl0 = rmxl_min
713      ENDIF
714     
[12377]715      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln
[1239]716
[1492]717      !                               !* depth of penetration of surface tke
718      IF( nn_etau /= 0 ) THEN     
[1601]719         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
[2528]720         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
[7753]721            htau(:,:) = 10._wp
[2528]722         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
[7753]723            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
[1492]724         END SELECT
725      ENDIF
[9019]726      !                                !* read or initialize all required files
727      CALL tke_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, dissl)
[1239]728      !
[9367]729      IF( lwxios ) THEN
730         CALL iom_set_rstw_var_active('en')
731         CALL iom_set_rstw_var_active('avt_k')
732         CALL iom_set_rstw_var_active('avm_k')
733         CALL iom_set_rstw_var_active('dissl')
734      ENDIF
[2528]735   END SUBROUTINE zdf_tke_init
[1239]736
737
[1531]738   SUBROUTINE tke_rst( kt, cdrw )
[9019]739      !!---------------------------------------------------------------------
740      !!                   ***  ROUTINE tke_rst  ***
741      !!                     
742      !! ** Purpose :   Read or write TKE file (en) in restart file
743      !!
744      !! ** Method  :   use of IOM library
745      !!                if the restart does not contain TKE, en is either
746      !!                set to rn_emin or recomputed
747      !!----------------------------------------------------------------------
748      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
749      !!
750      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
751      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
752      !
753      INTEGER ::   jit, jk              ! dummy loop indices
754      INTEGER ::   id1, id2, id3, id4   ! local integers
755      !!----------------------------------------------------------------------
756      !
757      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
758         !                                   ! ---------------
759         IF( ln_rstart ) THEN                   !* Read the restart file
760            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
761            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
762            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
763            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
764            !
765            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
[13286]766               CALL iom_get( numror, jpdom_auto, 'en'   , en   , ldxios = lrxios )
767               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k, ldxios = lrxios )
768               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k, ldxios = lrxios )
769               CALL iom_get( numror, jpdom_auto, 'dissl', dissl, ldxios = lrxios )
[9019]770            ELSE                                          ! start TKE from rest
[9169]771               IF(lwp) WRITE(numout,*)
[9190]772               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without TKE scheme, set en to background values'
[9019]773               en   (:,:,:) = rn_emin * wmask(:,:,:)
774               dissl(:,:,:) = 1.e-12_wp
775               ! avt_k, avm_k already set to the background value in zdf_phy_init
776            ENDIF
777         ELSE                                   !* Start from rest
[9169]778            IF(lwp) WRITE(numout,*)
[9190]779            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set en to the background value'
[9019]780            en   (:,:,:) = rn_emin * wmask(:,:,:)
781            dissl(:,:,:) = 1.e-12_wp
782            ! avt_k, avm_k already set to the background value in zdf_phy_init
783         ENDIF
784         !
785      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
786         !                                   ! -------------------
[9169]787         IF(lwp) WRITE(numout,*) '---- tke_rst ----'
[9367]788         IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
789         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en   , ldxios = lwxios )
790         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios )
791         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios )
792         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl, ldxios = lwxios )
793         IF( lwxios ) CALL iom_swap(      cxios_context          )
[9019]794         !
795      ENDIF
796      !
[1531]797   END SUBROUTINE tke_rst
[1239]798
799   !!======================================================================
[1531]800END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.