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/UKMO/NEMO_4.0.4_GO6_mixing/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_GO6_mixing/src/OCE/ZDF/zdftke.F90 @ 15649

Last change on this file since 15649 was 15649, checked in by dancopsey, 2 years ago

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