New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdftke.F90 in NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF – NEMO

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

Last change on this file since 13249 was 13249, checked in by clem, 4 years ago

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