source: branches/TAM_V3_2_2/NEMO/OPA_SRC/ZDF/zdftke_old.F90 @ 2578

Last change on this file since 2578 was 2578, checked in by rblod, 11 years ago

first import of NEMOTAM 3.2.2

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