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

Last change on this file since 13497 was 13497, checked in by techene, 7 months ago

re-introduce comments that have been erased by loops transformation see #2525

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