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/releases/r4.0/r4.0-HEAD/src/OCE/ZDF – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdftke.F90 @ 13284

Last change on this file since 13284 was 13284, checked in by smasson, 4 years ago

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