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 branches/dev_004_VVL/NEMO/OPA_SRC/ZDF – NEMO

source: branches/dev_004_VVL/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 1361

Last change on this file since 1361 was 1268, checked in by ctlod, 15 years ago

include the missing Von Karman constant in surface mixing length computation and set to 0.1 the default value for the rn_lmin parameter, see ticket: #291
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
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 zdf_tke_init routine
16   !!              -   !  2002-08  (G. Madec)  rn_cri and Free form, F90
17   !!              -   !  2004-10  (C. Ethe )  1D configuration
18   !!             2.0  !  2006-07  (S. Masson)  distributed restart using iom
19   !!             3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
20   !!                              - tke penetration (wind steering)
21   !!                              - suface condition for tke & mixing length
22   !!                              - Langmuir cells
23   !!              -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
24   !!              -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
25   !!----------------------------------------------------------------------
26#if defined key_zdftke   ||   defined key_esopa
27   !!----------------------------------------------------------------------
28   !!   'key_zdftke'                                   TKE vertical physics
29   !!----------------------------------------------------------------------
30   !!----------------------------------------------------------------------
31   !!   zdf_tke      : update momentum and tracer Kz from a tke scheme
32   !!   zdf_tke_init : initialization, namelist read, and parameters control
33   !!   tke_rst      : read/write tke restart in ocean restart file
34   !!----------------------------------------------------------------------
35   USE oce             ! ocean dynamics and active tracers
36   USE dom_oce         ! ocean space and time domain
37   USE zdf_oce         ! ocean vertical physics
38   USE sbc_oce         ! surface boundary condition: ocean
39   USE phycst          ! physical constants
40   USE zdfmxl          ! mixed layer
41   USE restart         ! only for lrst_oce
42   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
43   USE prtctl          ! Print control
44   USE in_out_manager  ! I/O manager
45   USE iom             ! I/O manager library
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   zdf_tke    ! routine called in step module
51
52   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
53   REAL(wp), PUBLIC                         ::   eboost              !: multiplicative coeff of the shear product.
54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy
55# if defined key_vectopt_memory
56   !                                                                !!! key_vectopt_memory
57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etmean              !: coefficient used for horizontal smoothing
58   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   eumean, evmean      !: at t-, u- and v-points
59# endif
60#if defined key_c1d
61   !                                                                !!! 1D cfg only
62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   hlc                 !: save finite Langmuir Circulation depth
65#endif
66
67   !                                       !!! ** Namelist  namtke  **
68   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke
69   LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not
70   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
71   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops
72   INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3)
73   INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
74   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1)
75   INTEGER  ::   nn_avb   =  0              ! constant or profile background on avt (=0/1)
76   REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
77   REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation
78   REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke
79   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm
80   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke (m2/s2)
81   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke (m2/s2)
82   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number
83   INTEGER  ::   nn_havtb = 1               ! horizontal shape or not for avtb (=0/1)
84   INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2)
85   INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1/2)
86   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length
87   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length
88   REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
89   REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells
90
91   REAL(wp), DIMENSION (jpi,jpj) ::   avtb_2d   ! set in tke_init, for other modif than ice
92
93   !! * Substitutions
94#  include "domzgr_substitute.h90"
95#  include "vectopt_loop_substitute.h90"
96   !!----------------------------------------------------------------------
97   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
98   !! $Id$
99   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
100   !!----------------------------------------------------------------------
101
102CONTAINS
103
104   SUBROUTINE zdf_tke( kt )
105      !!----------------------------------------------------------------------
106      !!                   ***  ROUTINE zdf_tke  ***
107      !!
108      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
109      !!      coefficients using a 1.5 turbulent closure scheme.
110      !!
111      !! ** Method  :   The time evolution of the turbulent kinetic energy
112      !!      (tke) is computed from a prognostic equation :
113      !!         d(en)/dt = eboost eav (d(u)/dz)**2         ! shear production
114      !!                  + d( rn_efave eav d(en)/dz )/dz   ! diffusion of tke
115      !!                  + grav/rau0 pdl eav d(rau)/dz     ! stratif. destruc.
116      !!                  - rn_ediss / emxl en**(2/3)       ! dissipation
117      !!      with the boundary conditions:
118      !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) )
119      !!         bottom : en = rn_emin
120      !!      -1- The dissipation and mixing turbulent lengh scales are computed
121      !!         from the usual diagnostic buoyancy length scale: 
122      !!         mxl= sqrt(2*en)/N  where N is the brunt-vaisala frequency
123      !!         with mxl = rn_lmin at the bottom minimum value of 0.4
124      !!      Four cases :
125      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
126      !!                  zmxld = zmxlm = mxl
127      !!         nn_mxl=1 : mxl bounded by the vertical scale factor.
128      !!                  zmxld = zmxlm = mxl
129      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl
130      !!                  is less than 1 (|d/dz(xml)|<1).
131      !!                  zmxld = zmxlm = mxl
132      !!         nn_mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface
133      !!                        to the bottom
134      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom
135      !!                        to the surface
136      !!                  zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
137      !!      -2- Compute the now Turbulent kinetic energy. The time differencing
138      !!      is implicit for vertical diffusion term, linearized for kolmo-
139      !!      goroff dissipation term, and explicit forward for both buoyancy
140      !!      and dynamic production terms. Thus a tridiagonal linear system is
141      !!      solved.
142      !!         Note that - the shear production is multiplied by eboost in order
143      !!      to set the critic richardson number to rn_cri (namelist parameter)
144      !!                   - the destruction by stratification term is multiplied
145      !!      by the Prandtl number (defined by an empirical funtion of the local
146      !!      Richardson number) if nn_pdl=1 (namelist parameter)
147      !!      coefficient (zesh2):
148      !!      -3- Compute the now vertical eddy vicosity and diffusivity
149      !!      coefficients from en (before the time stepping) and zmxlm:
150      !!              avm = max( avtb, rn_ediff*zmxlm*en^1/2 )
151      !!              avt = max( avmb, pdl*avm )  (pdl=1 if nn_pdl=0)
152      !!              eav = max( avmb, avm )
153      !!      avt and avm are horizontally averaged to avoid numerical insta-
154      !!      bilities.
155      !!        N.B. The computation is done from jk=2 to jpkm1 except for
156      !!      en. Surface value of avt avmu avmv are set once a time to
157      !!      their background value in routine zdf_tke_init.
158      !!
159      !! ** Action  :   compute en (now turbulent kinetic energy)
160      !!                update avt, avmu, avmv (before vertical eddy coef.)
161      !!
162      !! References : Gaspar et al., JGR, 1990,
163      !!              Blanke and Delecluse, JPO, 1991
164      !!              Mellor and Blumberg, JPO 2004
165      !!              Axell, JGR, 2002
166      !!----------------------------------------------------------------------
167      USE oce,     zwd    =>   ua   ! use ua as workspace
168      USE oce,     zmxlm  =>   va   ! use va as workspace
169      USE oce,     zmxld  =>   ta   ! use ta as workspace
170      USE oce,     ztkelc =>   sa   ! use sa as workspace
171      !
172      INTEGER, INTENT(in) ::   kt   ! ocean time step
173      !
174      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
175      REAL(wp) ::   zbbrau, zrn2, zesurf            ! temporary scalars
176      REAL(wp) ::   zfact1, ztx2, zdku              !    -         -
177      REAL(wp) ::   zfact2, zty2, zdkv              !    -         -
178      REAL(wp) ::   zfact3, zcoef, zcof, zav        !    -         -
179      REAL(wp) ::   zsh2, zpdl, zri, zsqen, zesh2   !    -         -
180      REAL(wp) ::   zemxl, zemlm, zemlp             !    -         -
181      REAL(wp) ::   zraug, zus, zwlc, zind          !    -         -
182      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
183      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   !  -      -
184      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
185      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
186      !!--------------------------------------------------------------------
187
188      IF( kt == nit000  )   CALL zdf_tke_init      ! Initialization (first time-step only)
189
190      !                                            ! Local constant initialization
191      zbbrau =  .5 * rn_ebb / rau0
192      zfact1 = -.5 * rdt * rn_efave
193      zfact2 = 1.5 * rdt * rn_ediss
194      zfact3 = 0.5 * rdt * rn_ediss
195
196      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
197      ! I.  Mixing length
198      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
199
200      ! Buoyancy length scale: l=sqrt(2*e/n**2)
201      ! ---------------------
202      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
203!!gm  this should be useless
204         zmxlm(:,:,1) = 0.e0
205!!gm end
206         zraug = 0.5 * vkarmn * 2.e5 / ( rau0 * grav )
207         DO jj = 2, jpjm1
208!CDIR NOVERRCHK
209            DO ji = fs_2, fs_jpim1   ! vector opt.
210               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
211               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
212               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )  )
213            END DO
214         END DO
215      ELSE                       ! surface set to the minimum value
216         zmxlm(:,:,1) = rn_lmin0
217      ENDIF
218      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value
219      !
220!CDIR NOVERRCHK
221      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2)
222!CDIR NOVERRCHK
223         DO jj = 2, jpjm1
224!CDIR NOVERRCHK
225            DO ji = fs_2, fs_jpim1   ! vector opt.
226               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
227               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
228            END DO
229         END DO
230      END DO
231
232      ! Physical limits for the mixing length
233      ! -------------------------------------
234      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
235      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
236
237      SELECT CASE ( nn_mxl )
238      !
239      CASE ( 0 )           ! bounded by the distance to surface and bottom
240         DO jk = 2, jpkm1
241            DO jj = 2, jpjm1
242               DO ji = fs_2, fs_jpim1   ! vector opt.
243                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
244                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
245                  zmxlm(ji,jj,jk) = zemxl
246                  zmxld(ji,jj,jk) = zemxl
247               END DO
248            END DO
249         END DO
250         !
251      CASE ( 1 )           ! bounded by the vertical scale factor
252         DO jk = 2, jpkm1
253            DO jj = 2, jpjm1
254               DO ji = fs_2, fs_jpim1   ! vector opt.
255                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
256                  zmxlm(ji,jj,jk) = zemxl
257                  zmxld(ji,jj,jk) = zemxl
258               END DO
259            END DO
260         END DO
261         !
262      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
263         DO jk = 2, jpkm1         ! from the surface to the bottom :
264            DO jj = 2, jpjm1
265               DO ji = fs_2, fs_jpim1   ! vector opt.
266                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
267               END DO
268            END DO
269         END DO
270         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1   ! vector opt.
273                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
274                  zmxlm(ji,jj,jk) = zemxl
275                  zmxld(ji,jj,jk) = zemxl
276               END DO
277            END DO
278         END DO
279         !
280      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
281         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
282            DO jj = 2, jpjm1
283               DO ji = fs_2, fs_jpim1   ! vector opt.
284                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
285               END DO
286            END DO
287         END DO
288         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
289            DO jj = 2, jpjm1
290               DO ji = fs_2, fs_jpim1   ! vector opt.
291                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
292               END DO
293            END DO
294         END DO
295!CDIR NOVERRCHK
296         DO jk = 2, jpkm1
297!CDIR NOVERRCHK
298            DO jj = 2, jpjm1
299!CDIR NOVERRCHK
300               DO ji = fs_2, fs_jpim1   ! vector opt.
301                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
302                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
303                  zmxlm(ji,jj,jk) = zemlm
304                  zmxld(ji,jj,jk) = zemlp
305               END DO
306            END DO
307         END DO
308         !
309      END SELECT
310
311# if defined key_c1d
312      ! c1d configuration : save mixing and dissipation turbulent length scales
313      e_dis(:,:,:) = zmxld(:,:,:)
314      e_mix(:,:,:) = zmxlm(:,:,:)
315# endif
316
317      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
318      ! II  TKE Langmuir circulation source term
319      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
320      IF( ln_lc ) THEN
321         !
322         ! Computation of total energy produce by LC : cumulative sum over jk
323         zpelc(:,:,1) =  MAX( rn2(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
324         DO jk = 2, jpk
325            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
326         END DO
327         !
328         ! Computation of finite Langmuir Circulation depth
329         ! Initialization to the number of w ocean point mbathy
330         imlc(:,:) = mbathy(:,:)
331         DO jk = jpkm1, 2, -1
332            DO jj = 1, jpj
333               DO ji = 1, jpi
334                  ! Last w-level at which zpelc>=0.5*us*us with us=0.016*wind(starting from jpk-1)
335                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
336                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
337               END DO
338            END DO
339         END DO
340         !
341         ! finite LC depth
342         DO jj = 1, jpj
343            DO ji = 1, jpi
344               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
345            END DO
346         END DO
347         !
348# if defined key_c1d
349         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
350# endif
351         !
352         ! TKE Langmuir circulation source term
353!CDIR NOVERRCHK
354         DO jk = 2, jpkm1
355!CDIR NOVERRCHK
356            DO jj = 2, jpjm1
357!CDIR NOVERRCHK
358               DO ji = fs_2, fs_jpim1   ! vector opt.
359                  ! Stokes drift
360                  zus  = 0.016 * wndm(ji,jj)
361                  ! computation of vertical velocity due to LC
362                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
363                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
364                  ! TKE Langmuir circulation source term
365                  ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
366               END DO
367            END DO
368         END DO
369         !
370      ENDIF
371
372      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
373      ! III  Tubulent kinetic energy time stepping
374      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
375
376      ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
377      ! ---------------------------------------------------------------------
378!CDIR NOVERRCHK
379      DO jk = 2, jpkm1
380!CDIR NOVERRCHK
381         DO jj = 2, jpjm1
382!CDIR NOVERRCHK
383            DO ji = fs_2, fs_jpim1   ! vector opt.
384               zsqen = SQRT( en(ji,jj,jk) )
385               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
386               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
387               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
388               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
389            END DO
390         END DO
391      END DO
392
393      ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
394      ! -------------------------------------------------
395      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
396      ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
397!CDIR NOVERRCHK
398      DO jj = 2, jpjm1
399!CDIR NOVERRCHK
400         DO ji = fs_2, fs_jpim1   ! vector opt.
401            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
402            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
403            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
404            en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
405            zav  =  rn_ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )
406            zmxlm(ji,jj,1  ) = MAX( zav, avmb(1)                  ) * tmask(ji,jj,1)
407            avt  (ji,jj,1  ) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1)
408            zmxlm(ji,jj,jpk) = 0.e0
409         END DO
410      END DO
411
412      ! 3. Now Turbulent kinetic energy (output in en)
413      ! -------------------------------
414      ! Resolution of a tridiagonal linear system by a "methode de chasse"
415      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
416
417      SELECT CASE ( nn_pdl )
418      !
419      CASE ( 0 )           ! No Prandtl number
420         DO jk = 2, jpkm1
421            DO jj = 2, jpjm1
422               DO ji = fs_2, fs_jpim1   ! vector opt.
423                  !                                                             ! shear prod. - stratif. destruction
424                  zcoef = 0.5 / fse3w(ji,jj,jk)
425                  zdku = zcoef * (  ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &          ! shear
426                     &            - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
427                  zdkv = zcoef * (  vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
428                     &            - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
429                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)        ! coefficient (zesh2)
430                  !
431                  !                                                             ! Matrix
432                  zcof = zfact1 * tmask(ji,jj,jk)
433                  !                                                                  ! lower diagonal
434                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
435                     &                  / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
436                  !                                                                  ! upper diagonal
437                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
438                     &                  / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
439                  !                                                                  ! diagonal
440                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
441                  !
442                  !                                                             ! right hand side in en
443                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
444                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
445                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
446                  ELSE                                                               ! Langmuir cell source term
447                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
448                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2          &
449                        &                        +   rdt  * ztkelc(ji,jj,jk)
450                  ENDIF
451               END DO
452            END DO
453         END DO
454         !
455      CASE ( 1 )           ! Prandtl number
456         DO jk = 2, jpkm1
457            DO jj = 2, jpjm1
458               DO ji = fs_2, fs_jpim1   ! vector opt.
459                  !                                                             ! shear prod. - stratif. destruction
460                  zcoef = 0.5 / fse3w(ji,jj,jk)                                     
461                  zdku = zcoef * (  ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1)   &          ! shear
462                  &               - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )  )
463                  zdkv = zcoef * (  vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
464                  &               - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )  )
465                  zsh2 = zdku * zdku + zdkv * zdkv                                   ! square of shear
466                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number
467# if defined key_c1d
468                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
469# endif
470                  zpdl = 1.0                                                         ! Prandtl number
471                  IF( zri >= 0.2 )   zpdl = 0.2 / zri
472                  zpdl = MAX( 0.1, zpdl )
473                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)                       ! coefficient (esh2)
474                  !
475                  !                                                             ! Matrix
476                  zcof = zfact1 * tmask(ji,jj,jk)
477                  !                                                                 ! lower diagonal
478                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
479                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
480                  !                                                                 ! upper diagonal
481                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
482                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
483                  !                                                                 ! diagonal
484                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
485                  !
486                  !                                                             ! right hand side in en
487                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
488                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
489                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
490                  ELSE                                                               ! Langmuir cell source term
491                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
492                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2             &
493                        &                        +   rdt  * ztkelc(ji,jj,jk)
494                  ENDIF
495                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)                      ! store masked Prandlt number in zmxld array
496               END DO
497            END DO
498         END DO
499         !
500      END SELECT
501
502# if defined key_c1d
503      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)      ! c1d configuration : save masked Prandlt number
504      e_pdl(:,:,      1) = e_pdl(:,:,      2)
505      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)     
506# endif
507
508      ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
509      !!--------------------------------
510      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
511         DO jj = 2, jpjm1
512            DO ji = fs_2, fs_jpim1    ! vector opt.
513               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
514            END DO
515         END DO
516      END DO
517      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
518         DO ji = fs_2, fs_jpim1    ! vector opt.
519            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
520         END DO
521      END DO
522      DO jk = 3, jpkm1
523         DO jj = 2, jpjm1
524            DO ji = fs_2, fs_jpim1    ! vector opt.
525               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
526            END DO
527         END DO
528      END DO
529      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
530         DO ji = fs_2, fs_jpim1    ! vector opt.
531            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
532         END DO
533      END DO
534      DO jk = jpk-2, 2, -1
535         DO jj = 2, jpjm1
536            DO ji = fs_2, fs_jpim1    ! vector opt.
537               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
538            END DO
539         END DO
540      END DO
541      DO jk = 2, jpkm1                             ! set the minimum value of tke
542         DO jj = 2, jpjm1
543            DO ji = fs_2, fs_jpim1   ! vector opt.
544               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
545            END DO
546         END DO
547      END DO
548     
549      ! 5. Add extra TKE due to surface and internal wave breaking (nn_etau /= 0)
550      !!----------------------------------------------------------
551      IF( nn_etau /= 0 ) THEN        ! extra tke : en = en + rn_efr * en(1) * exp( -z/zhtau )
552         !
553         SELECT CASE( nn_htau )           ! Choice of the depth of penetration
554         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
555            DO jj = 2, jpjm1
556               DO ji = fs_2, fs_jpim1   ! vector opt.
557                  zhtau(ji,jj) = 10.
558               END DO
559            END DO
560         CASE( 1 )                                    ! meridional profile  1
561            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 40 m at high lat.)
562               DO ji = fs_2, fs_jpim1   ! vector opt.
563                  zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
564               END DO
565            END DO
566         CASE( 2 )                                    ! meridional profile 2
567            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 60 m at high lat.)
568               DO ji = fs_2, fs_jpim1   ! vector opt.
569                  zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
570               END DO
571            END DO
572         END SELECT
573         !
574         IF( nn_etau == 1 ) THEN          ! extra term throughout the water column
575            DO jk = 2, jpkm1
576               DO jj = 2, jpjm1
577                  DO ji = fs_2, fs_jpim1   ! vector opt.
578                     en(ji,jj,jk) = en(ji,jj,jk)   &
579                        &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
580                        &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
581                  END DO
582               END DO
583            END DO
584         ELSEIF( nn_etau == 2 ) THEN      ! extra term only at the base of the mixed layer
585            DO jj = 2, jpjm1
586               DO ji = fs_2, fs_jpim1   ! vector opt.
587                  jk = nmln(ji,jj)
588                  en(ji,jj,jk) = en(ji,jj,jk)    &
589                     &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
590                     &                   * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
591               END DO
592            END DO
593         ENDIF
594         !
595      ENDIF
596
597
598      ! Lateral boundary conditions (sign unchanged)
599      CALL lbc_lnk( en, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
600
601
602      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
603      ! IV.  Before vertical eddy vicosity and diffusivity coefficients
604      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
605      !
606      SELECT CASE ( nn_ave )
607      CASE ( 0 )                                     ! no horizontal average
608         DO jk = 2, jpkm1                                 ! only vertical eddy viscosity
609            DO jj = 2, jpjm1
610               DO ji = fs_2, fs_jpim1   ! vector opt.
611                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
612                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
613                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
614                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
615               END DO
616            END DO
617         END DO
618         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
619         !
620      CASE ( 1 )                                     ! horizontal average                         ( 1/2  1/2 )
621         !                                                ! Vertical eddy viscosity    avmu = 1/4 ( 1    1   )
622         !                                                !                                       ( 1/2  1/2 )
623         !                                                !
624         !                                                !                                       ( 1/2  1   1/2 )           
625         !                                                !                            avmv = 1/4 ( 1/2  1   1/2 )     
626         DO jk = 2, jpkm1                               
627            DO jj = 2, jpjm1
628               DO ji = fs_2, fs_jpim1   ! vector opt.
629# if defined key_vectopt_memory
630                  !                                       ! caution : vectopt_memory change the solution
631                  !                                       !           (last digit of the solver stat)
632                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
633                     &              +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
634                     &                   +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
635
636                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
637                     &              +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
638                     &                   +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
639# else
640                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
641                     &           +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
642                     &                +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)   &
643                     &    / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
644                     &           +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
645                     &                +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
646
647                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
648                     &           +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
649                     &                +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)   &
650                     &   /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
651                     &           +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
652                     &                +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
653# endif
654               END DO
655            END DO
656         END DO
657         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
658         !
659         !                                                ! Vertical eddy diffusivity             (1 2 1)
660         !                                                !                            avt = 1/16 (2 4 2)
661         !                                                !                                       (1 2 1)
662         DO jk = 2, jpkm1         
663            DO jj = 2, jpjm1
664               DO ji = fs_2, fs_jpim1   ! vector opt.
665#  if defined key_vectopt_memory
666                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)                          &
667                     &            + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
668#  else
669                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)                        &
670                     &            + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
671                     &  / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)                        &
672                     &            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
673#  endif
674               END DO
675            END DO
676         END DO
677         !
678      END SELECT
679      !
680      IF( nn_pdl == 1 ) THEN                            ! Ponderation by the Prandtl number (nn_pdl=1)
681         DO jk = 2, jpkm1 
682            DO jj = 2, jpjm1
683               DO ji = fs_2, fs_jpim1   ! vector opt.
684                  zpdl = zmxld(ji,jj,jk)
685                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
686               END DO
687            END DO
688         END DO
689      ENDIF
690      !
691      DO jk = 2, jpkm1                                  ! Minimum value on the eddy viscosity
692         avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk)
693         avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk)
694      END DO
695
696      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
697
698      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )      ! write en in restart file
699
700      IF(ln_ctl) THEN
701         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
702         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
703            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
704      ENDIF
705      !
706   END SUBROUTINE zdf_tke
707
708
709   SUBROUTINE zdf_tke_init
710      !!----------------------------------------------------------------------
711      !!                  ***  ROUTINE zdf_tke_init  ***
712      !!                     
713      !! ** Purpose :   Initialization of the vertical eddy diffivity and
714      !!      viscosity when using a tke turbulent closure scheme
715      !!
716      !! ** Method  :   Read the namtke namelist and check the parameters
717      !!      called at the first timestep (nit000)
718      !!
719      !! ** input   :   Namlist namtke
720      !!
721      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
722      !!
723      !!----------------------------------------------------------------------
724      USE dynzdf_exp
725      USE trazdf_exp
726      !
727# if defined key_vectopt_memory
728      INTEGER ::   ji, jj, jk   ! dummy loop indices
729# else
730      INTEGER ::           jk   ! dummy loop indices
731# endif
732      !!
733      NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb  , rn_efave, rn_emin,   &
734         &             rn_emin0, rn_cri  , nn_itke , nn_mxl  , nn_pdl  , nn_ave ,   &
735         &             nn_avb  , ln_mxl0 , rn_lmin , rn_lmin0, nn_havtb, nn_etau,   &
736         &             nn_htau , rn_efr  , ln_lc   , rn_lc 
737      !!----------------------------------------------------------------------
738
739      ! Read Namelist namtke : Turbulente Kinetic Energy
740      ! --------------------
741      REWIND ( numnam )
742      READ   ( numnam, namtke )
743
744      ! Compute boost associated with the Richardson critic
745      !     (control values: rn_cri = 0.3   ==> eboost=1.25 for nn_pdl=1)
746      !     (                rn_cri = 0.222 ==> eboost=1.               )
747      eboost = rn_cri * ( 2. + rn_ediss / rn_ediff ) / 2.
748
749
750
751      ! Parameter control and print
752      ! ---------------------------
753      ! Control print
754      IF(lwp) THEN
755         WRITE(numout,*)
756         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme'
757         WRITE(numout,*) '~~~~~~~~~~~~'
758         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
759         WRITE(numout,*) '             restart with tke from no tke              ln_rstke = ', ln_rstke
760         WRITE(numout,*) '             coef. to compute avt                      rn_ediff = ', rn_ediff
761         WRITE(numout,*) '             Kolmogoroff dissipation coef.             rn_ediss = ', rn_ediss
762         WRITE(numout,*) '             tke surface input coef.                   rn_ebb   = ', rn_ebb
763         WRITE(numout,*) '             tke diffusion coef.                       rn_efave = ', rn_efave
764         WRITE(numout,*) '             minimum value of tke                      rn_emin  = ', rn_emin
765         WRITE(numout,*) '             surface minimum value of tke              rn_emin0 = ', rn_emin0
766         WRITE(numout,*) '             number of restart iter loops              nn_itke  = ', nn_itke
767         WRITE(numout,*) '             mixing length type                        nn_mxl   = ', nn_mxl
768         WRITE(numout,*) '             prandl number flag                        nn_pdl   = ', nn_pdl
769         WRITE(numout,*) '             horizontal average flag                   nn_ave   = ', nn_ave
770         WRITE(numout,*) '             critic Richardson nb                      rn_cri   = ', rn_cri
771         WRITE(numout,*) '                and its associated coeff.              eboost   = ', eboost
772         WRITE(numout,*) '             constant background or profile            nn_avb   = ', nn_avb
773         WRITE(numout,*) '             surface mixing length = F(stress) or not  ln_mxl0  = ', ln_mxl0
774         WRITE(numout,*) '             surface  mixing length minimum value      rn_lmin0 = ', rn_lmin0
775         WRITE(numout,*) '             interior mixing length minimum value      rn_lmin0 = ', rn_lmin0
776         WRITE(numout,*) '             horizontal variation for avtb             nn_havtb = ', nn_havtb
777         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau  = ', nn_etau
778         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau  = ', nn_htau
779         WRITE(numout,*) '             % of rn_emin0 which pene. the thermocline rn_efr   = ', rn_efr
780         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc    = ', ln_lc
781         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc    = ', rn_lc
782         WRITE(numout,*)
783      ENDIF
784
785      ! Check of some namelist values
786      IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
787      IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
788      IF( nn_ave  < 0    .OR. nn_ave  > 1   )   CALL ctl_stop( 'bad flag: nn_ave is  0 or 1    ' )
789      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
790      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
791
792      IF( nn_etau == 2 )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
793
794      ! Horizontal average : initialization of weighting arrays
795      ! -------------------
796      !
797      SELECT CASE ( nn_ave )
798      ! Notice: mean arrays have not to by defined at local domain boundaries due to the vertical nature of TKE
799      !
800      CASE ( 0 )                ! no horizontal average
801         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
802         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
803# if defined key_vectopt_memory
804         ! caution vectopt_memory change the solution (last digit of the solver stat)
805         ! weighting mean arrays etmean, eumean and evmean
806         !           ( 1  1 )                                          ( 1 )
807         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
808         !                         
809         etmean(:,:,:) = 0.e0
810         eumean(:,:,:) = 0.e0
811         evmean(:,:,:) = 0.e0
812         !
813         DO jk = 1, jpkm1
814            DO jj = 2, jpjm1
815               DO ji = fs_2, fs_jpim1   ! vector opt.
816                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
817                  &                                            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk) )
818                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
819                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
820               END DO
821            END DO
822         END DO
823# endif
824         !
825      CASE ( 1 )                ! horizontal average
826         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
827# if defined key_vectopt_memory
828         ! caution vectopt_memory change the solution (last digit of the solver stat)
829         ! weighting mean arrays etmean, eumean and evmean
830         !           ( 1  1 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
831         ! avt = 1/4 ( 1  1 )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
832         !                                 ( 1/2  1/2 )
833         etmean(:,:,:) = 0.e0
834         eumean(:,:,:) = 0.e0
835         evmean(:,:,:) = 0.e0
836         !
837         DO jk = 1, jpkm1
838            DO jj = 2, jpjm1
839               DO ji = fs_2, fs_jpim1   ! vector opt.
840                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,   umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk)   &
841                  &                                       +       vmask(ji  ,jj-1,jk) + vmask(ji  ,jj  ,jk) )
842                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
843                  &                                       +.5 * ( tmask(ji  ,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
844                  &                                             + tmask(ji  ,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
845                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
846                  &                                       +.5 * ( tmask(ji-1,jj  ,jk) + tmask(ji-1,jj+1,jk)   &
847                  &                                             + tmask(ji+1,jj  ,jk) + tmask(ji+1,jj+1,jk) )  )
848               END DO
849            END DO
850         END DO
851# endif
852         !
853      END SELECT
854
855
856      ! Background eddy viscosity and diffusivity profil
857      ! ------------------------------------------------
858      IF( nn_avb == 0 ) THEN      ! Define avmb, avtb from namelist parameter
859         avmb(:) = avm0
860         avtb(:) = avt0
861      ELSE                      ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
862         avmb(:) = avm0
863         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
864         IF(ln_sco .AND. lwp)   CALL ctl_warn( '          avtb profile not valid in sco' )
865      ENDIF
866      !
867      !                         ! 2D shape of the avtb
868      avtb_2d(:,:) = 1.e0             ! uniform
869      !
870      IF( nn_havtb == 1 ) THEN      ! decrease avtb in the equatorial band
871           !  -15S -5S : linear decrease from avt0 to avt0/10.
872           !  -5S  +5N : cst value avt0/10.
873           !   5N  15N : linear increase from avt0/10, to avt0
874           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
875           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
876           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
877      ENDIF
878
879      ! Initialization of vertical eddy coef. to the background value
880      ! -------------------------------------------------------------
881      DO jk = 1, jpk
882         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
883         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
884         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
885      END DO
886
887
888      ! read or initialize turbulent kinetic energy ( en )
889      ! -------------------------------------------------
890      CALL tke_rst( nit000, 'READ' )
891      !
892   END SUBROUTINE zdf_tke_init
893
894
895   SUBROUTINE tke_rst( kt, cdrw )
896     !!---------------------------------------------------------------------
897     !!                   ***  ROUTINE ts_rst  ***
898     !!                     
899     !! ** Purpose :   Read or write TKE file (en) in restart file
900     !!
901     !! ** Method  :   use of IOM library
902     !!                if the restart does not contain TKE, en is either
903     !!                set to rn_emin or recomputed (nn_itke/=0)
904     !!----------------------------------------------------------------------
905     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
906     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
907     !
908     INTEGER ::   jit   ! dummy loop indices
909     !!----------------------------------------------------------------------
910     !
911     IF( TRIM(cdrw) == 'READ' ) THEN
912        IF( ln_rstart ) THEN
913           IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN
914              CALL iom_get( numror, jpdom_autoglo, 'en', en )
915           ELSE
916              IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 )   &
917                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
918              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
919              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop'
920              en (:,:,:) = rn_emin * tmask(:,:,:)
921              DO jit = 2, nn_itke + 1
922                 CALL zdf_tke( jit )
923              END DO
924           ENDIF
925        ELSE
926           en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value
927        ENDIF
928     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
929        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
930     ENDIF
931     !
932   END SUBROUTINE tke_rst
933
934#else
935   !!----------------------------------------------------------------------
936   !!   Dummy module :                                        NO TKE scheme
937   !!----------------------------------------------------------------------
938   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
939CONTAINS
940   SUBROUTINE zdf_tke( kt )          ! Empty routine
941      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
942   END SUBROUTINE zdf_tke
943#endif
944
945   !!======================================================================
946END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.