source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdftke.F90 @ 13006

Last change on this file since 13006 was 13006, checked in by cetlod, 5 months ago

branch clem_dan_fixcpl : add mixing length parameterization depending on sea-ice thickness, see ticket #2476

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