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/releases/r4.0/r4.0-HEAD/src/OCE/ZDF – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdftke.F90 @ 12697

Last change on this file since 12697 was 12697, checked in by mathiot, 4 years ago

ticket #2406: fix ticket #2406 for r4.0-HEAD

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