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/branches/2020/dev_r12953_ENHANCE-10_acc_fix_traqsr/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r12953_ENHANCE-10_acc_fix_traqsr/src/OCE/ZDF/zdftke.F90 @ 13121

Last change on this file since 13121 was 13121, checked in by acc, 4 years ago

2020/dev_r12953_ENHANCE-10_acc_fix_traqsr. Merge in trunk changes from 12953 to current HEAD (13115). Fully SETTE tested

  • Property svn:keywords set to Id
File size: 41.4 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
30   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition (ln_drg)
[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
[13121]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
70   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
71   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
[13121]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)
[4147]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)
[9019]81   LOGICAL  ::   ln_drg    ! top/bottom friction forcing flag
[4147]82   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
[9019]83   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
84   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
[9546]85   REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4   
[4147]86   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
[9019]87   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells
[1239]88
[4147]89   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
90   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
[2528]91   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
92   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
[1239]93
[9019]94   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau)
95   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation
96   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation
[1492]97
[1239]98   !! * Substitutions
[12377]99#  include "do_loop_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      !
201      INTEGER ::   ji, jj, jk              ! dummy loop arguments
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
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      INTEGER , DIMENSION(jpi,jpj)     ::   imlc
[10425]212      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i
[9019]213      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
[1239]214      !!--------------------------------------------------------------------
[1492]215      !
[12489]216      zbbrau = rn_ebb / rho0       ! Local constant initialisation
217      zfact1 = -.5_wp * rn_Dt 
218      zfact2 = 1.5_wp * rn_Dt * rn_ediss
[2528]219      zfact3 = 0.5_wp       * rn_ediss
[1492]220      !
221      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[9019]222      !                     !  Surface/top/bottom boundary condition on tke
[1492]223      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[12698]224      !
[12377]225      DO_2D_00_00
226         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
227      END_2D
[9019]228      !
[1492]229      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
230      !                     !  Bottom boundary condition on tke
231      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1719]232      !
[12489]233      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
[9019]234      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
235      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
[1492]236      !
[9019]237      IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE
238         !
[12377]239         DO_2D_00_00
240            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
241            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
[12489]242            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
[12377]243            zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
244               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
245            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
246         END_2D
[9019]247         IF( ln_isfcav ) THEN       ! top friction
[12377]248            DO_2D_00_00
249               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
250               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
[12489]251               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
[12377]252               zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
253                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
[12702]254               ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
255               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &
256                  &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
[12377]257            END_2D
[9019]258         ENDIF
259         !
260      ENDIF
261      !
[1492]262      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[9019]263      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
[1492]264         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[1239]265         !
[1492]266         !                        !* total energy produce by LC : cumulative sum over jk
[12377]267         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm)
[1239]268         DO jk = 2, jpk
[12377]269            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
[1239]270         END DO
[1492]271         !                        !* finite Langmuir Circulation depth
[1705]272         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
[7753]273         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
[12377]274         DO_3DS_11_11( jpkm1, 2, -1 )
275            zus  = zcof * taum(ji,jj)
276            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
277         END_3D
[1492]278         !                               ! finite LC depth
[12377]279         DO_2D_11_11
280            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
281         END_2D
[1705]282         zcof = 0.016 / SQRT( zrhoa * zcdrag )
[12377]283         DO_2D_00_00
284            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
285            zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
286            IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0.
287         END_2D
288         DO_3D_00_00( 2, jpkm1 )
289            IF ( zfr_i(ji,jj) /= 0. ) THEN               
290               ! vertical velocity due to LC   
291               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
292                  !                                           ! vertical velocity due to LC
293                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i
294                  !                                           ! TKE Langmuir circulation source term
[12489]295                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
[12377]296               ENDIF
297            ENDIF
298         END_3D
[1239]299         !
300      ENDIF
[1492]301      !
302      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
303      !                     !  Now Turbulent kinetic energy (output in en)
304      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
305      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
306      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
307      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
308      !
[9019]309      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
[12377]310         DO_3D_00_00( 2, jpkm1 )
311            !                             ! local Richardson number
312            zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
313            !                             ! inverse of Prandtl number
314            apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
315         END_3D
[5656]316      ENDIF
[5836]317      !         
[12377]318      DO_3D_00_00( 2, jpkm1 )
319         zcof   = zfact1 * tmask(ji,jj,jk)
320         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
321         !                                   ! eddy coefficient (ensure numerical stability)
322         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
323            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
324         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
325            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
326         !
327         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
328         zd_lw(ji,jj,jk) = zzd_lw
329         zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
330         !
331         !                                   ! right hand side in en
[12489]332         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                        &   ! shear
[12377]333            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
334            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
335            &                                ) * wmask(ji,jj,jk)
336      END_3D
[5120]337      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
[12377]338      DO_3D_00_00( 3, jpkm1 )
339         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
340      END_3D
341      DO_2D_00_00
342         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
343      END_2D
344      DO_3D_00_00( 3, jpkm1 )
345         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)
346      END_3D
347      DO_2D_00_00
348         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
349      END_2D
350      DO_3DS_00_00( jpk-2, 2, -1 )
351         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
352      END_3D
353      DO_3D_00_00( 2, jpkm1 )
354         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
355      END_3D
[9019]356      !
[1492]357      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
358      !                            !  TKE due to surface and internal wave breaking
359      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[6140]360!!gm BUG : in the exp  remove the depth of ssh !!!
[12377]361!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
[6140]362     
363     
[2528]364      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
[12377]365         DO_3D_00_00( 2, jpkm1 )
366            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
367               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
368         END_3D
[2528]369      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
[12377]370         DO_2D_00_00
371            jk = nmln(ji,jj)
372            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
373               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
374         END_2D
[2528]375      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
[12377]376         DO_3D_00_00( 2, jpkm1 )
377            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
378            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
379            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
380            zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
381            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
382            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
383               &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
384         END_3D
[1239]385      ENDIF
[1492]386      !
[1239]387   END SUBROUTINE tke_tke
388
[1492]389
[12377]390   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
[1239]391      !!----------------------------------------------------------------------
[1492]392      !!                   ***  ROUTINE tke_avn  ***
[1239]393      !!
[1492]394      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
395      !!
396      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
397      !!              the tke_tke routine). First, the now mixing lenth is
398      !!      computed from en and the strafification (N^2), then the mixings
399      !!      coefficients are computed.
400      !!              - Mixing length : a first evaluation of the mixing lengh
401      !!      scales is:
402      !!                      mxl = sqrt(2*en) / N 
403      !!      where N is the brunt-vaisala frequency, with a minimum value set
[2528]404      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
[1492]405      !!        The mixing and dissipative length scale are bound as follow :
406      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
407      !!                        zmxld = zmxlm = mxl
408      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
409      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
410      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
411      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
412      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
413      !!                    the surface to obtain ldown. the resulting length
414      !!                    scales are:
415      !!                        zmxld = sqrt( lup * ldown )
416      !!                        zmxlm = min ( lup , ldown )
417      !!              - Vertical eddy viscosity and diffusivity:
418      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
419      !!                      avt = max( avmb, pdlr * avm ) 
420      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
421      !!
[9019]422      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
[1239]423      !!----------------------------------------------------------------------
[9019]424      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics
425      !!
[12377]426      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
[9019]427      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
428      !
[2715]429      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[9019]430      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
431      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
[13121]432      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
[9019]433      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
[1239]434      !!--------------------------------------------------------------------
[3294]435      !
[1492]436      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
437      !                     !  Mixing length
438      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
439      !
440      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
441      !
[5120]442      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
[7753]443      zmxlm(:,:,:)  = rmxl_min   
444      zmxld(:,:,:)  = rmxl_min
[5120]445      !
[13121]446     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
447         !
[12489]448         zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
[13121]449#if ! defined key_si3 && ! defined key_cice
[12377]450         DO_2D_00_00
[13121]451            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
[12377]452         END_2D
[13121]453#else
454         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
455         !
456         CASE( 0 )                      ! No scaling under sea-ice
457            DO_2D_00_00
458               zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
459            END_2D
460            !
461         CASE( 1 )                           ! scaling with constant sea-ice thickness
462            DO_2D_00_00
463               zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1)
464            END_2D
465            !
466         CASE( 2 )                                 ! scaling with mean sea-ice thickness
467            DO_2D_00_00
468#if defined key_si3
469               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)
470#elif defined key_cice
471               zmaxice = MAXVAL( h_i(ji,jj,:) )
472               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
473#endif
474            END_2D
475            !
476         CASE( 3 )                                 ! scaling with max sea-ice thickness
477            DO_2D_00_00
478               zmaxice = MAXVAL( h_i(ji,jj,:) )
479               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
480            END_2D
481            !
482         END SELECT
483#endif
484         !
485         DO_2D_00_00
486            zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
487         END_2D
488         !
489      ELSE
[7753]490         zmxlm(:,:,1) = rn_mxl0
[1239]491      ENDIF
[13121]492
[1239]493      !
[12377]494      DO_3D_00_00( 2, jpkm1 )
495         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
496         zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
497      END_3D
[1492]498      !
499      !                     !* Physical limits for the mixing length
500      !
[7753]501      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
502      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
[1492]503      !
[1239]504      SELECT CASE ( nn_mxl )
505      !
[5836]506 !!gm Not sure of that coding for ISF....
[12377]507      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
[1239]508      CASE ( 0 )           ! bounded by the distance to surface and bottom
[12377]509         DO_3D_00_00( 2, jpkm1 )
510            zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
511            &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
512            ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
513            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
514            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
515         END_3D
[1239]516         !
517      CASE ( 1 )           ! bounded by the vertical scale factor
[12377]518         DO_3D_00_00( 2, jpkm1 )
519            zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
520            zmxlm(ji,jj,jk) = zemxl
521            zmxld(ji,jj,jk) = zemxl
522         END_3D
[1239]523         !
524      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
[12377]525         DO_3D_00_00( 2, jpkm1 )
526            zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
527         END_3D
528         DO_3DS_00_00( jpkm1, 2, -1 )
529            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
530            zmxlm(ji,jj,jk) = zemxl
531            zmxld(ji,jj,jk) = zemxl
532         END_3D
[1239]533         !
534      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
[12377]535         DO_3D_00_00( 2, jpkm1 )
536            zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
537         END_3D
538         DO_3DS_00_00( jpkm1, 2, -1 )
539            zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
540         END_3D
541         DO_3D_00_00( 2, jpkm1 )
542            zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
543            zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
544            zmxlm(ji,jj,jk) = zemlm
545            zmxld(ji,jj,jk) = zemlp
546         END_3D
[1239]547         !
548      END SELECT
[1492]549      !
550      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[9019]551      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
[1492]552      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[12377]553      DO_3D_00_00( 1, jpkm1 )
554         zsqen = SQRT( en(ji,jj,jk) )
555         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
556         p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
557         p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
558         dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
559      END_3D
[1492]560      !
561      !
562      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
[12377]563         DO_3D_00_00( 2, jpkm1 )
[12698]564            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]565         END_3D
[1239]566      ENDIF
[9019]567      !
[12377]568      IF(sn_cfctl%l_prtctl) THEN
[9440]569         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
570         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk )
[1239]571      ENDIF
572      !
[1492]573   END SUBROUTINE tke_avn
[1239]574
[1492]575
[12377]576   SUBROUTINE zdf_tke_init( Kmm )
[1239]577      !!----------------------------------------------------------------------
[2528]578      !!                  ***  ROUTINE zdf_tke_init  ***
[1239]579      !!                     
580      !! ** Purpose :   Initialization of the vertical eddy diffivity and
[1492]581      !!              viscosity when using a tke turbulent closure scheme
[1239]582      !!
[1601]583      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
[1492]584      !!              called at the first timestep (nit000)
[1239]585      !!
[1601]586      !! ** input   :   Namlist namzdf_tke
[1239]587      !!
588      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
589      !!----------------------------------------------------------------------
[9019]590      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
591      !!
[12377]592      INTEGER, INTENT(in) ::   Kmm          ! time level index
593      INTEGER             ::   ji, jj, jk   ! dummy loop indices
594      INTEGER             ::   ios
[1239]595      !!
[13121]596      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  &
597         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
598         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
599         &                 nn_pdl  , ln_drg   , ln_lc    , rn_lc,      &
600         &                 nn_etau , nn_htau  , rn_efr   , rn_eice 
[1239]601      !!----------------------------------------------------------------------
[2715]602      !
[4147]603      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
[11536]604901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' )
[4147]605
606      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
[11536]607902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' )
[4624]608      IF(lwm) WRITE ( numond, namzdf_tke )
[2715]609      !
[2528]610      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
[2715]611      !
[1492]612      IF(lwp) THEN                    !* Control print
[1239]613         WRITE(numout,*)
[2528]614         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
615         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]616         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
[1705]617         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
618         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
619         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
620         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
621         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
[9019]622         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
[1705]623         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
624         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
[9019]625         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
[13121]626         IF( ln_mxl0 ) THEN
627            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice
628            IF( nn_mxlice == 1 ) &
629            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice
630         ENDIF         
[9019]631         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
632         WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg
633         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
634         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
[1705]635         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
[9019]636         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
637         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
[9546]638         WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice
639         WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   '
[9019]640         IF( ln_drg ) THEN
[9169]641            WRITE(numout,*)
[9019]642            WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
643            WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top
644            WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot
645         ENDIF
646         WRITE(numout,*)
[9190]647         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
[9019]648         WRITE(numout,*)
[1239]649      ENDIF
[2715]650      !
[9019]651      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
652         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
653         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
[9190]654         IF(lwp) WRITE(numout,*) '   ==>>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
[9019]655      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
656         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
[9190]657         IF(lwp) WRITE(numout,*) '   ==>>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
[9019]658      ENDIF
659      !
[2715]660      !                              ! allocate tke arrays
661      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
662      !
[1492]663      !                               !* Check of some namelist values
[4990]664      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
665      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
666      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
[5407]667      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
[9019]668      !
[2528]669      IF( ln_mxl0 ) THEN
[9169]670         IF(lwp) WRITE(numout,*)
[9190]671         IF(lwp) WRITE(numout,*) '   ==>>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
[2528]672         rn_mxl0 = rmxl_min
673      ENDIF
674     
[12377]675      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln
[1239]676
[1492]677      !                               !* depth of penetration of surface tke
678      IF( nn_etau /= 0 ) THEN     
[1601]679         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
[2528]680         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
[7753]681            htau(:,:) = 10._wp
[2528]682         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
[7753]683            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
[1492]684         END SELECT
685      ENDIF
[9019]686      !                                !* read or initialize all required files
687      CALL tke_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, dissl)
[1239]688      !
[9367]689      IF( lwxios ) THEN
690         CALL iom_set_rstw_var_active('en')
691         CALL iom_set_rstw_var_active('avt_k')
692         CALL iom_set_rstw_var_active('avm_k')
693         CALL iom_set_rstw_var_active('dissl')
694      ENDIF
[2528]695   END SUBROUTINE zdf_tke_init
[1239]696
697
[1531]698   SUBROUTINE tke_rst( kt, cdrw )
[9019]699      !!---------------------------------------------------------------------
700      !!                   ***  ROUTINE tke_rst  ***
701      !!                     
702      !! ** Purpose :   Read or write TKE file (en) in restart file
703      !!
704      !! ** Method  :   use of IOM library
705      !!                if the restart does not contain TKE, en is either
706      !!                set to rn_emin or recomputed
707      !!----------------------------------------------------------------------
708      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
709      !!
710      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
711      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
712      !
713      INTEGER ::   jit, jk              ! dummy loop indices
714      INTEGER ::   id1, id2, id3, id4   ! local integers
715      !!----------------------------------------------------------------------
716      !
717      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
718         !                                   ! ---------------
719         IF( ln_rstart ) THEN                   !* Read the restart file
720            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
721            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
722            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
723            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
724            !
725            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
[9367]726               CALL iom_get( numror, jpdom_autoglo, 'en'   , en   , ldxios = lrxios )
727               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios )
728               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios )
729               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl, ldxios = lrxios )
[9019]730            ELSE                                          ! start TKE from rest
[9169]731               IF(lwp) WRITE(numout,*)
[9190]732               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without TKE scheme, set en to background values'
[9019]733               en   (:,:,:) = rn_emin * wmask(:,:,:)
734               dissl(:,:,:) = 1.e-12_wp
735               ! avt_k, avm_k already set to the background value in zdf_phy_init
736            ENDIF
737         ELSE                                   !* Start from rest
[9169]738            IF(lwp) WRITE(numout,*)
[9190]739            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set en to the background value'
[9019]740            en   (:,:,:) = rn_emin * wmask(:,:,:)
741            dissl(:,:,:) = 1.e-12_wp
742            ! avt_k, avm_k already set to the background value in zdf_phy_init
743         ENDIF
744         !
745      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
746         !                                   ! -------------------
[9169]747         IF(lwp) WRITE(numout,*) '---- tke_rst ----'
[9367]748         IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
749         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en   , ldxios = lwxios )
750         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios )
751         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios )
752         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl, ldxios = lwxios )
753         IF( lwxios ) CALL iom_swap(      cxios_context          )
[9019]754         !
755      ENDIF
756      !
[1531]757   END SUBROUTINE tke_rst
[1239]758
759   !!======================================================================
[1531]760END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.