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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 9546

Last change on this file since 9546 was 9546, checked in by gm, 6 years ago

dev_merge_2017: add a TKE parameter to control the activation of nn_etau below sea-ice

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