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 @ 1050

Last change on this file since 1050 was 1050, checked in by ctlod, 16 years ago

trunk: update TKE physics, see ticket: #183
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.6 KB
Line 
1MODULE zdftke
2   !!======================================================================
3   !!                       ***  MODULE  zdftke  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !! History :   6.0  !  91-03  (b. blanke)  Original code
8   !!             7.0  !  91-11  (G. Madec)   bug fix
9   !!             7.1  !  92-10  (G. Madec)   new mixing length and eav
10   !!             7.2  !  93-03  (M. Guyon)   symetrical conditions
11   !!             7.3  !  94-08  (G. Madec, M. Imbard)   npdl flag
12   !!             7.5  !  96-01  (G. Madec)   s-coordinates
13   !!             8.0  !  97-07  (G. Madec)   lbc
14   !!             8.1  !  99-01  (E. Stretta) new option for the mixing length
15   !!             8.5  !  02-06  (G. Madec) add zdf_tke_init routine
16   !!             8.5  !  02-08  (G. Madec)  ri_c and Free form, F90
17   !!             9.0  !  04-10  (C. Ethe )  1D configuration
18   !!             9.0  !  02-08  (G. Madec)  autotasking optimization
19   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
20   !!              -   !  2005-07  (C. Ethe,  G.Madec) : update TKE physics:
21   !!                              - tke penetration (wind steering)
22   !!                              - suface condition for tke & mixing length
23   !!                              - Langmuir cells
24   !!             3.0  !  2007-11  (G. Madec)  avtb_2d, nn_avtb_2d
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 lbclnk          ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl          ! Print control
43   USE in_out_manager  ! I/O manager
44   USE iom
45   USE restart         ! only for lrst_oce
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   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etmean              !: coefficient used for horizontal smoothing
57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   eumean, evmean      !: at t-, u- and v-points
58# endif
59
60   !! * Namelist (namtke)
61   LOGICAL , PUBLIC ::   ln_rstke = .FALSE.         !: =T restart with tke from a run without tke with
62     !                                              !  a none zero initial value for en
63   INTEGER , PUBLIC ::   nitke = 50 ,            &  !: number of restart iterative loops
64      &                  nmxl  =  2 ,            &  !: = 0/1/2/3 flag for the type of mixing length used
65      &                  npdl  =  1 ,            &  !: = 0/1/2 flag on prandtl number on vert. eddy coeff.
66      &                  nave  =  1 ,            &  !: = 0/1 flag for horizontal average on avt, avmu, avmv
67      &                  navb  =  0                 !: = 0/1 flag for constant or profile background avt
68   REAL(wp), PUBLIC ::   ediff = 0.1_wp       ,  &  !: coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e)
69      &                  ediss = 0.7_wp       ,  &  !: coef. of the Kolmogoroff dissipation
70      &                  ebb   = 3.75_wp      ,  &  !: coef. of the surface input of tke
71      &                  efave = 1._wp        ,  &  !: coef. for the tke vert. diff. coeff.; avtke=efave*avm
72      &                  emin  = 0.7071e-6_wp ,  &  !: minimum value of tke (m2/s2)
73      &                  emin0 = 1.e-4_wp     ,  &  !: surface minimum value of tke (m2/s2)
74      &                  ri_c  = 2._wp / 9._wp      !: critic Richardson number
75   !                                      !!! ** Namelist (namtke) **
76   INTEGER  ::   nn_avtb_2d = 1            !: = 0/1 horizontal shape for avtb
77   INTEGER  ::   nn_etau    = 0            !: = 0/1/2 tke depth penetration
78   INTEGER  ::   nn_htau    = 0            !: = 0/1/2 type of tke profile of penetration
79   REAL(wp) ::   rn_lmin    = 0.4_wp       !: surface min value of mixing turbulent length scale
80   REAL(wp) ::   rn_efr     = 1.0_wp       !: fraction of TKE surface value which penetrates in the ocean
81   LOGICAL  ::   ln_lsfc    = .FALSE.      !: mixing length scale surface value as function of wind stress or not
82   LOGICAL  ::   ln_lc      = .FALSE.      !: Langmuir cells (LC) as a source term of TKE or not
83   REAL(wp) ::   rn_lc      = 0.15_wp      !: coef to compute vertical velocity of LC
84
85   REAL(wp), DIMENSION (jpi,jpj)  :: avtb_2d ! set in tke_init, for other modif than ice
86
87# if defined key_c1d
88   !                                                                   ! 1D cfg only
89   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix,      &  ! dissipation and mixing turbulent lengh scales
90      &                                          e_pdl, e_ric          ! prandl and local Richardson numbers
91#endif
92
93   !! * Substitutions
94#  include "domzgr_substitute.h90"
95#  include "vectopt_loop_substitute.h90"
96   !!----------------------------------------------------------------------
97   !!   OPA 9.0 , LOCEAN-IPSL (2006)
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( efave eav d(en)/dz )/dz    ! diffusion of tke
115      !!                  + grav/rau0 pdl eav d(rau)/dz   ! stratif. destruc.
116      !!                  - ediss / emxl en**(2/3)        ! dissipation
117      !!      with the boundary conditions:
118      !!         surface: en = max( emin0,ebb sqrt(utau^2 + vtau^2) )
119      !!         bottom : en = emin
120      !!      -1- The dissipation and mixing turbulent lengh scales are computed
121      !!      from the usual diagnostic buoyancy length scale: 
122      !!         mxl= 1/(sqrt(en)/N)  WHERE N is the brunt-vaisala frequency
123      !!      Four cases :
124      !!         nmxl=0 : mxl bounded by the distance to surface and bottom.
125      !!                  zmxld = zmxlm = mxl
126      !!         nmxl=1 : mxl bounded by the vertical scale factor.
127      !!                  zmxld = zmxlm = mxl
128      !!         nmxl=2 : mxl bounded such that the vertical derivative of mxl
129      !!                  is less than 1 (|d/dz(xml)|<1).
130      !!                  zmxld = zmxlm = mxl
131      !!         nmxl=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 ri_c (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 npdl=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, ediff*zmxlm*en^1/2 )
150      !!              avt = max( avmb, pdl*avm )  (pdl=1 if npdl=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, 95, 1990,
162      !!              Blanke and Delecluse, jpo, 1991
163      !!----------------------------------------------------------------------
164      USE oce     , zwd   => ua,  &  ! use ua as workspace
165         &          zmxlm => ta,  &  ! use ta as workspace
166         &          zmxld => sa      ! use sa as workspace
167      USE oce     , ztkelc => va     ! use va as workspace
168      !
169      INTEGER, INTENT(in) ::   kt ! ocean time step
170      !
171      INTEGER  ::   ji, jj, jk                  ! dummy loop arguments
172      REAL(wp) ::   zmlmin, zbbrau,          &  ! temporary scalars
173         &          zfact1, zfact2, zfact3,  &  !
174         &          zrn2, zesurf,            &  !
175         &          ztx2, zty2, zav,         &  !
176         &          zcoef, zcof, zsh2,       &  !
177         &          zdku, zdkv, zpdl, zri,   &  !
178         &          zsqen, zesh2,            &  !
179         &          zemxl, zemlm, zemlp
180      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
181      REAL(wp) ::   zraug, zus, zwlc, zind     ! temporary scalar
182      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   ! 2D workspace
183      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    ! 2D workspace
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      zmlmin = 0.4
191      zbbrau =  .5 * ebb / rau0
192      zfact1 = -.5 * rdt * efave
193      zfact2 = 1.5 * rdt * ediss
194      zfact3 = 0.5 * rdt * ediss
195
196      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
197      ! I.  Mixing length
198      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
199
200      ! Buoyancy length scale: l=sqrt(2*e/n**2)
201      ! ---------------------
202      IF( ln_lsfc ) THEN      ! lsurf is function of wind stress : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
203         zmxlm(:,:,1) = 0.e0
204         zraug = 0.5 * 2.e5 / ( rau0 * grav )
205         DO jj = 2, jpjm1
206!CDIR NOVERRCHK
207            DO ji = fs_2, fs_jpim1   ! vector opt.
208               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
209               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
210               zmxlm(ji,jj,1  ) = zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )
211               ! set the surface minimum value to rn_lmin
212               zmxlm(ji,jj,1  ) = MAX( zmxlm(ji,jj,1) , rn_lmin )
213            END DO
214         END DO
215      ELSE                    ! surface set to the minimum value
216         zmxlm(:,:,1) = rn_lmin 
217      ENDIF
218      zmxlm(:,:,jpk) = rn_lmin   ! bottom  set to the minimum value
219!CDIR NOVERRCHK
220      DO jk = 2, jpkm1
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( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  )
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 ( nmxl )
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      ! 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)      ! 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   = 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)   = ebb sqrt(utau^2+vtau^2) / rau0  (min value 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, emin0 ) * tmask(ji,jj,1)
404            zav  =  ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )* tmask(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
415      ! e(jpk)=0 ).
416
417      SELECT CASE ( npdl )
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                  ! zesh2 = eboost * (du/dz)^2 - N^2
424                  zcoef = 0.5 / fse3w(ji,jj,jk)
425                  ! shear
426                  zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &
427                  &                - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
428                  zdkv = zcoef * (   vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
429                  &                - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
430                  ! coefficient (zesh2)
431                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)
432
433                  ! Matrix
434                  zcof = zfact1 * tmask(ji,jj,jk)
435                  ! lower diagonal
436                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
437                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
438                  ! upper diagonal
439                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
440                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
441                  ! diagonal
442                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
443                  ! right hand side in en
444                  IF( .NOT. ln_lc ) THEN
445                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
446                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
447                  ELSE
448                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
449                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2          &
450                        &                        +   rdt  * ztkelc(ji,jj,jk)
451                  ENDIF
452               END DO
453            END DO
454         END DO
455
456      CASE ( 1 )           ! Prandtl number
457         DO jk = 2, jpkm1
458            DO jj = 2, jpjm1
459               DO ji = fs_2, fs_jpim1   ! vector opt.
460                  ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2
461                  zcoef = 0.5 / fse3w(ji,jj,jk)
462                  ! shear
463                  zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) &
464                  &                - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )   )
465                  zdkv = zcoef * (   vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1) &
466                  &                - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )   )
467                  ! square of vertical shear
468                  zsh2 = zdku * zdku + zdkv * zdkv
469                  ! local Richardson number
470                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
471# if defined key_c1d
472                  ! save masked local Richardson number in zmxlm array
473                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)
474# endif
475                  ! Prandtl number
476                  zpdl = 1.0
477                  IF( zri >= 0.2 ) zpdl = 0.2 / zri
478                  zpdl = MAX( 0.1, zpdl )
479                  ! coefficient (esh2)
480                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
481
482                  ! Matrix
483                  zcof = zfact1 * tmask(ji,jj,jk)
484                  ! lower diagonal
485                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
486                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
487                  ! upper diagonal
488                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
489                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
490                  ! diagonal
491                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
492                  ! right hand side in en
493                  IF( .NOT. ln_lc ) THEN
494                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
495                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
496                  ELSE
497                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
498                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2             &
499                        &                        +   rdt  * ztkelc(ji,jj,jk)
500                  ENDIF
501                  ! save masked Prandlt number in zmxld array
502                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
503               END DO
504            END DO
505         END DO
506
507      END SELECT
508
509# if defined key_c1d
510      !  save masked Prandlt number
511      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)
512      e_pdl(:,:,      1) = e_pdl(:,:,      2)
513      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)     
514# endif
515
516      ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
517      !!--------------------------------
518      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
519      DO jk = 3, jpkm1
520         DO jj = 2, jpjm1
521            DO ji = fs_2, fs_jpim1    ! vector opt.
522               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
523            END DO
524         END DO
525      END DO
526
527      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
528      DO jj = 2, jpjm1
529         DO ji = fs_2, fs_jpim1    ! vector opt.
530            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
531         END DO
532      END DO
533      DO jk = 3, jpkm1
534         DO jj = 2, jpjm1
535            DO ji = fs_2, fs_jpim1    ! vector opt.
536               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
537            END DO
538         END DO
539      END DO
540
541      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
542      DO jj = 2, jpjm1
543         DO ji = fs_2, fs_jpim1    ! vector opt.
544            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
545         END DO
546      END DO
547      DO jk = jpk-2, 2, -1
548         DO jj = 2, jpjm1
549            DO ji = fs_2, fs_jpim1    ! vector opt.
550               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
551            END DO
552         END DO
553      END DO
554
555      ! Save the result in en and set minimum value of tke : emin
556      DO jk = 2, jpkm1
557         DO jj = 2, jpjm1
558            DO ji = fs_2, fs_jpim1   ! vector opt.
559               en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
560            END DO
561         END DO
562      END DO
563     
564      ! Modify the value of TKE by an exponential assumption
565      ! en(ji,jj,jk)=en(ji,jj,jk)+en(ji,jj,1)*EXP(-fsdepw(ji,jj,jk)/ zhtau(ji,jj) )
566
567      SELECT CASE( nn_htau )      ! Choice of H value
568      !
569      CASE( 0 )
570         DO jj = 2, jpjm1
571            DO ji = fs_2, fs_jpim1   ! vector opt.
572               zhtau(ji,jj) = 10.
573            END DO
574         END DO
575         !
576      CASE( 1 )
577         DO jj = 2, jpjm1
578            DO ji = fs_2, fs_jpim1   ! vector opt.
579               zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
580            END DO
581         END DO
582         !
583      CASE( 2 )
584         DO jj = 2, jpjm1
585            DO ji = fs_2, fs_jpim1   ! vector opt.
586               zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
587            END DO
588         END DO
589         !
590      END SELECT
591
592      IF( nn_etau == 1 ) THEN
593         DO jk = 2, jpkm1
594            DO jj = 2, jpjm1
595               DO ji = fs_2, fs_jpim1   ! vector opt.
596                  en(ji,jj,jk) = en(ji,jj,jk)   &
597                     &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
598                     &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
599               END DO
600            END DO
601         END DO
602      ELSEIF( nn_etau == 2 ) THEN     !  only at the base of the mixed layer
603         DO jj = 2, jpjm1
604            DO ji = fs_2, fs_jpim1   ! vector opt.
605               jk = nmln(ji,jj)
606               en(ji,jj,jk) = en(ji,jj,jk)    &
607                  &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
608                  &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
609            END DO
610         END DO
611      ENDIF
612
613
614      ! Lateral boundary conditions on ( avt, en )   (sign unchanged)
615      CALL lbc_lnk( en , 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
616
617      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
618      ! IV.  Before vertical eddy vicosity and diffusivity coefficients
619      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
620
621      SELECT CASE ( nave )
622         
623      CASE ( 0 )                ! no horizontal average
624
625         ! Vertical eddy viscosity
626
627         DO jk = 2, jpkm1                                 ! Horizontal slab
628            DO jj = 2, jpjm1
629               DO ji = fs_2, fs_jpim1   ! vector opt.
630                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
631                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
632                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
633                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
634               END DO
635            END DO
636         END DO
637
638         ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
639         CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
640         
641      CASE ( 1 )                ! horizontal average
642
643         !                                                ( 1/2  1/2 )
644         ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1    1   )
645         !                      ( 1/2  1 1/2 )            ( 1/2  1/2 )
646         !           avmv = 1/4 ( 1/2  1 1/2 )     
647         
648!! caution vectopt_memory change the solution (last digit of the solver stat)
649#  if defined key_vectopt_memory
650         DO jk = 2, jpkm1                                 ! Horizontal slab
651            DO jj = 2, jpjm1
652               DO ji = fs_2, fs_jpim1   ! vector opt.
653                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
654                  &                 +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
655                  &                      +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
656
657                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
658                  &                 +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
659                  &                      +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
660               END DO
661            END DO
662         END DO
663#  else
664         DO jk = 2, jpkm1                                 ! Horizontal slab
665            DO jj = 2, jpjm1
666               DO ji = fs_2, fs_jpim1   ! vector opt.
667                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
668                  &              +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
669                  &                   +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)  &
670                  &       / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
671                  &              +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
672                  &                   +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
673
674                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
675                  &              +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
676                  &                   +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)  &
677                  &      /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
678                  &              +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
679                  &                   +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
680               END DO
681            END DO
682         END DO
683#  endif
684
685         ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
686         CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
687
688         ! Vertical eddy diffusivity
689         ! ------------------------------
690         !                                (1 2 1)
691         ! horizontal average  avt = 1/16 (2 4 2)
692         !                                (1 2 1)
693         DO jk = 2, jpkm1                                 ! Horizontal slab
694#  if defined key_vectopt_memory
695            DO jj = 2, jpjm1
696               DO ji = fs_2, fs_jpim1   ! vector opt.
697                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)    &
698                  &               + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
699               END DO
700            END DO
701#  else
702            DO jj = 2, jpjm1
703               DO ji = fs_2, fs_jpim1   ! vector opt.
704                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)   &
705                  &               + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
706                  &     / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
707                  &               + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
708               END DO
709            END DO
710#  endif
711         END DO
712
713      END SELECT
714
715      ! multiplied by the Prandtl number (npdl>1)
716      ! ----------------------------------------
717      IF( npdl == 1 ) THEN
718         DO jk = 2, jpkm1                                 ! Horizontal slab
719            DO jj = 2, jpjm1
720               DO ji = fs_2, fs_jpim1   ! vector opt.
721                  zpdl = zmxld(ji,jj,jk)
722                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
723               END DO
724            END DO
725         END DO
726      ENDIF
727
728      ! Minimum value on the eddy viscosity
729      ! ----------------------------------------
730      DO jk = 2, jpkm1                                 ! Horizontal slab
731         DO jj = 1, jpj
732            DO ji = 1, jpi
733               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
734               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
735            END DO
736         END DO
737      END DO
738
739      ! Lateral boundary conditions on avt  (sign unchanged)
740      ! ------------------------------=====
741      CALL lbc_lnk( avt, 'W', 1. )
742
743      ! write en in restart file
744      ! ------------------------
745      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
746
747      IF(ln_ctl) THEN
748         CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
749         CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask, &
750            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
751      ENDIF
752
753   END SUBROUTINE zdf_tke
754
755
756   SUBROUTINE zdf_tke_init
757      !!----------------------------------------------------------------------
758      !!                  ***  ROUTINE zdf_tke_init  ***
759      !!                     
760      !! ** Purpose :   Initialization of the vertical eddy diffivity and
761      !!      viscosity when using a tke turbulent closure scheme
762      !!
763      !! ** Method  :   Read the namtke namelist and check the parameters
764      !!      called at the first timestep (nit000)
765      !!
766      !! ** input   :   Namlist namtke
767      !!
768      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
769      !!
770      !!----------------------------------------------------------------------
771      USE dynzdf_exp
772      USE trazdf_exp
773      !
774# if defined key_vectopt_memory
775      ! caution vectopt_memory change the solution (last digit of the solver stat)
776      INTEGER ::   ji, jj, jk   ! dummy loop indices
777# else
778      INTEGER ::           jk   ! dummy loop indices
779# endif
780      !!
781      NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
782                       ri_c, nitke, nmxl, npdl, nave, navb,    &
783                       ln_lsfc, rn_lmin, nn_avtb_2d, nn_etau, nn_htau, rn_efr, &
784                       ln_lc , rn_lc 
785      !!----------------------------------------------------------------------
786
787      ! Read Namelist namtke : Turbulente Kinetic Energy
788      ! --------------------
789      REWIND ( numnam )
790      READ   ( numnam, namtke )
791
792      ! Compute boost associated with the Richardson critic
793      !     (control values: ri_c = 0.3   ==> eboost=1.25 for npdl=1 or 2)
794      !     (                ri_c = 0.222 ==> eboost=1.                  )
795      eboost = ri_c * ( 2. + ediss / ediff ) / 2.
796
797
798      ! Parameter control and print
799      ! ---------------------------
800      ! Control print
801      IF(lwp) THEN
802         WRITE(numout,*)
803         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme'
804         WRITE(numout,*) '~~~~~~~~~~~~'
805         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
806         WRITE(numout,*) '             restart with tke from no tke ln_rstke = ', ln_rstke
807         WRITE(numout,*) '             coef. to compute avt           ediff  = ', ediff
808         WRITE(numout,*) '             Kolmogoroff dissipation coef.  ediss  = ', ediss
809         WRITE(numout,*) '             tke surface input coef.        ebb    = ', ebb
810         WRITE(numout,*) '             tke diffusion coef.            efave  = ', efave
811         WRITE(numout,*) '             minimum value of tke           emin   = ', emin
812         WRITE(numout,*) '             surface minimum value of tke   emin0  = ', emin0
813         WRITE(numout,*) '             number of restart iter loops   nitke  = ', nitke
814         WRITE(numout,*) '             mixing length type             nmxl   = ', nmxl
815         WRITE(numout,*) '             prandl number flag             npdl   = ', npdl
816         WRITE(numout,*) '             horizontal average flag        nave   = ', nave
817         WRITE(numout,*) '             critic Richardson nb           ri_c   = ', ri_c
818         WRITE(numout,*) '                and its associated coeff.   eboost = ', eboost
819         WRITE(numout,*) '             constant background or profile navb   = ', navb
820         WRITE(numout,*) '             flag for compu.of bls as fun. of wind     ln_lsfc    = ', ln_lsfc
821         WRITE(numout,*) '             buoyancy lenght scale surface min value   rn_lmin    = ', rn_lmin
822         WRITE(numout,*) '             horizontal variation for avtb             nn_avtb_2d = ', nn_avtb_2d
823         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau    = ', nn_etau
824         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau    = ', nn_htau
825         WRITE(numout,*) '             % of emin which pene. the thermocline     rn_efr     = ', rn_efr
826         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc      = ', ln_lc
827         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc      = ', rn_lc
828         WRITE(numout,*)
829      ENDIF
830
831      ! Check nmxl and npdl values
832      IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( '          bad flag: nmxl is < 0 or > 3 ' )
833      IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( '          bad flag: npdl is < 0 or > 1 ' )
834      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is < 0 or > 3 ' )
835      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be > 0.15 and < 0.2 ' )
836
837      IF( nn_etau == 2 )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
838
839      ! Horizontal average : initialization of weighting arrays
840      ! -------------------
841     
842      SELECT CASE ( nave )
843
844      CASE ( 0 )                ! no horizontal average
845         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
846         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
847# if defined key_vectopt_memory
848         ! caution vectopt_memory change the solution (last digit of the solver stat)
849         ! weighting mean arrays etmean, eumean and evmean
850         !           ( 1  1 )                                          ( 1 )
851         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
852         !                         
853         etmean(:,:,:) = 0.e0
854         eumean(:,:,:) = 0.e0
855         evmean(:,:,:) = 0.e0
856         
857         DO jk = 1, jpkm1
858            DO jj = 2, jpjm1
859               DO ji = fs_2, fs_jpim1   ! vector opt.
860                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
861                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
862                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
863                 
864                  eumean(ji,jj,jk) = umask(ji,jj,jk)                     &
865                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk)  )
866
867                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                     &
868                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  )
869               END DO
870            END DO
871         END DO
872# endif
873
874      CASE ( 1 )                ! horizontal average
875         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
876# if defined key_vectopt_memory
877         ! caution vectopt_memory change the solution (last digit of the solver stat)
878         ! weighting mean arrays etmean, eumean and evmean
879         !           ( 1  1 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
880         ! avt = 1/4 ( 1  1 )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
881         !                                 ( 1/2  1/2 )
882         etmean(:,:,:) = 0.e0
883         eumean(:,:,:) = 0.e0
884         evmean(:,:,:) = 0.e0
885         
886         DO jk = 1, jpkm1
887            DO jj = 2, jpjm1
888               DO ji = fs_2, fs_jpim1   ! vector opt.
889                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
890                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
891                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
892                 
893                  eumean(ji,jj,jk) = umask(ji,jj,jk)                        &
894                  &  / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
895                  &       +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
896                  &              +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
897
898                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                        &
899                  &  / MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
900                  &       +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
901                  &              +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
902               END DO
903            END DO
904         END DO
905# endif
906
907      CASE DEFAULT
908         WRITE(ctmp1,*) '          bad flag value for nave = ', nave
909         CALL ctl_stop( ctmp1 )
910
911      END SELECT
912
913
914      ! Background eddy viscosity and diffusivity profil
915      ! ------------------------------------------------
916      IF( navb == 0 ) THEN      ! Define avmb, avtb from namelist parameter
917         avmb(:) = avm0
918         avtb(:) = avt0
919      ELSE                      ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
920         avmb(:) = avm0
921!!bug  this is not valide neither in scoord
922         IF(ln_sco .AND. lwp) WRITE(numout,cform_war)
923         IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile not valid in sco'
924
925         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
926      ENDIF
927
928      !                         ! 2D shape of the avtb
929      avtb_2d(:,:) = 1.e0             ! uniform
930      !
931      IF( nn_avtb_2d == 1 ) THEN      ! decrease avtb in the equatorial band
932           !  -15S -5S : linear decrease from avt0 to avt0/10.
933           !  -5S  +5N : cst value avt0/10.
934           !   5N  15N : linear increase from avt0/10, to avt0
935           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
936           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
937           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
938      ENDIF
939
940      !   Increase the background in the surface layers
941!!gm  the increase is only required when using cen2 advection scheme
942!!gm  for the equatorial upwelling.
943      avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
944      avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
945      avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
946      avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
947
948
949      ! Initialization of vertical eddy coef. to the background value
950      ! -------------------------------------------------------------
951      DO jk = 1, jpk
952         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
953         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
954         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
955      END DO
956
957
958      ! read or initialize turbulent kinetic energy ( en )
959      ! -------------------------------------------------
960      CALL tke_rst( nit000, 'READ' )
961      !
962   END SUBROUTINE zdf_tke_init
963
964
965   SUBROUTINE tke_rst( kt, cdrw )
966     !!---------------------------------------------------------------------
967     !!                   ***  ROUTINE ts_rst  ***
968     !!                     
969     !! ** Purpose : Read or write filtered free surface arrays in restart file
970     !!
971     !! ** Method  :
972     !!
973     !!----------------------------------------------------------------------
974     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
975     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
976     !
977     INTEGER ::   jit   ! dummy loop indices
978     !!----------------------------------------------------------------------
979     !
980     IF( TRIM(cdrw) == 'READ' ) THEN
981        IF( ln_rstart ) THEN
982           IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN
983              CALL iom_get( numror, jpdom_autoglo, 'en', en )
984           ELSE
985              IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 )   &
986                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
987              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
988              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en set by iterative loop'
989              IF( lwp                )   WRITE(numout,*) ' =======             ========='
990              en (:,:,:) = emin * tmask(:,:,:)
991              DO jit = 2, nitke+1
992                 CALL zdf_tke( jit )
993              END DO
994           ENDIF
995        ELSE
996           en(:,:,:) = emin * tmask(:,:,:)      ! no restart: en set to emin
997        ENDIF
998     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
999        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
1000     ENDIF
1001     !
1002   END SUBROUTINE tke_rst
1003
1004#else
1005   !!----------------------------------------------------------------------
1006   !!   Dummy module :                                        NO TKE scheme
1007   !!----------------------------------------------------------------------
1008   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
1009CONTAINS
1010   SUBROUTINE zdf_tke( kt )          ! Empty routine
1011      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
1012   END SUBROUTINE zdf_tke
1013#endif
1014
1015   !!======================================================================
1016END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.