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

source: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 1146

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago

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