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

Last change on this file since 13530 was 13530, checked in by clem, 7 months ago

remove duplicated lines

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