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

Last change on this file since 258 was 258, checked in by opalod, 19 years ago

nemo_v1_update_004 : CT : Integration of the control print option for debugging work
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.2 KB
Line 
1MODULE zdftke
2   !!======================================================================
3   !!                       ***  MODULE  zdftke  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7#if defined key_zdftke   ||   defined key_esopa
8   !!----------------------------------------------------------------------
9   !!   'key_zdftke'                                             TKE scheme
10   !!----------------------------------------------------------------------
11   !!   zdf_tke      : update momentum and tracer Kz from a tke scheme
12   !!   zdf_tke_init : initialization, namelist read, and parameters control
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and active tracers
16   USE dom_oce         ! ocean space and time domain
17   USE zdf_oce         ! ocean vertical physics
18   USE in_out_manager  ! I/O manager
19   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
20   USE phycst          ! physical constants
21   USE taumod          ! surface stress
22   USE prtctl          ! Print control
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC zdf_tke   ! routine called by step.F90
29
30   !! * Share Module variables
31   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.    !: TKE vertical mixing flag
32   LOGICAL, PUBLIC ::         & !!: ** tke namelist (namtke) **
33     ln_rstke = .FALSE.          !: =T restart with tke from a run without tke with
34     !                           !  a none zero initial value for en
35   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
36      en                         !: now turbulent kinetic energy
37
38   !! * Module variables
39   INTEGER ::                 & !!! ** tke namelist (namtke) **
40      nitke = 50 ,            &  ! number of restart iterative loops
41      nmxl  =  2 ,            &  ! = 0/1/2/3 flag for the type of mixing length used
42      npdl  =  1 ,            &  ! = 0/1/2 flag on prandtl number on vert. eddy coeff.
43      nave  =  1 ,            &  ! = 0/1 flag for horizontal average on avt, avmu, avmv
44      navb  =  0                 ! = 0/1 flag for constant or profile background avt
45   REAL(wp) ::                & !!! ** tke namlist (namtke) **
46      ediff = 0.1_wp       ,  &  ! coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e)
47      ediss = 0.7_wp       ,  &  ! coef. of the Kolmogoroff dissipation
48      ebb   = 3.75_wp      ,  &  ! coef. of the surface input of tke
49      efave = 1._wp        ,  &  ! coef. for the tke vert. diff. coeff.; avtke=efave*avm
50      emin  = 0.7071e-6_wp ,  &  ! minimum value of tke (m2/s2)
51      emin0 = 1.e-4_wp     ,  &  ! surface minimum value of tke (m2/s2)
52      ri_c  = 2._wp / 9._wp      ! critic Richardson number
53   REAL(wp) ::   & 
54      eboost        ! multiplicative coeff of the shear product.
55
56   !! caution vectopt_memory change the solution (last digit of the solver stat)
57# if defined key_vectopt_memory
58   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
59      etmean,    &  ! coefficient used for horizontal smoothing
60      eumean,    &  ! at t-, u- and v-points
61      evmean        !
62# endif
63
64# if defined key_cfg_1d
65   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   
66      e_dis,    &   ! dissipation turbulent lengh scale
67      e_mix,    &   ! mixing turbulent lengh scale
68      e_pdl,    &   ! prandl number
69      e_ric         ! local Richardson number
70#endif
71
72   !! * Substitutions
73#  include "domzgr_substitute.h90"
74#  include "vectopt_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !!   OPA 9.0 , LOCEAN-IPSL (2005)
77   !! $Header$
78   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
79   !!----------------------------------------------------------------------
80
81CONTAINS
82
83# if defined key_autotasking
84   !!----------------------------------------------------------------------
85   !!   'key_autotasking' :                             j-k-i loop (j-slab)
86   !!----------------------------------------------------------------------
87#  include "zdftke_atsk.h90"
88
89# else
90   !!----------------------------------------------------------------------
91   !!   default option :                                         k-j-i loop
92   !!----------------------------------------------------------------------
93
94   SUBROUTINE zdf_tke ( kt )
95      !!----------------------------------------------------------------------
96      !!                   ***  ROUTINE zdf_tke  ***
97      !!
98      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
99      !!      coefficients using a 1.5 turbulent closure scheme.
100      !!
101      !! ** Method  :   The time evolution of the turbulent kinetic energy
102      !!      (tke) is computed from a prognostic equation :
103      !!         d(en)/dt = eboost eav (d(u)/dz)**2       ! shear production
104      !!                  + d( efave eav d(en)/dz )/dz    ! diffusion of tke
105      !!                  + grav/rau0 pdl eav d(rau)/dz   ! stratif. destruc.
106      !!                  - ediss / emxl en**(2/3)        ! dissipation
107      !!      with the boundary conditions:
108      !!         surface: en = max( emin0,ebb sqrt(taux^2 + tauy^2) )
109      !!         bottom : en = emin
110      !!      -1- The dissipation and mixing turbulent lengh scales are computed
111      !!      from the usual diagnostic buoyancy length scale: 
112      !!         mxl= 1/(sqrt(en)/N)  WHERE N is the brunt-vaisala frequency
113      !!      Four cases :
114      !!         nmxl=0 : mxl bounded by the distance to surface and bottom.
115      !!                  zmxld = zmxlm = mxl
116      !!         nmxl=1 : mxl bounded by the vertical scale factor.
117      !!                  zmxld = zmxlm = mxl
118      !!         nmxl=2 : mxl bounded such that the vertical derivative of mxl
119      !!                  is less than 1 (|d/dz(xml)|<1).
120      !!                  zmxld = zmxlm = mxl
121      !!         nmxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface
122      !!                        to the bottom
123      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom
124      !!                        to the surface
125      !!                  zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
126      !!      -2- Compute the now Turbulent kinetic energy. The time differencing
127      !!      is implicit for vertical diffusion term, linearized for kolmo-
128      !!      goroff dissipation term, and explicit forward for both buoyancy
129      !!      and dynamic production terms. Thus a tridiagonal linear system is
130      !!      solved.
131      !!         Note that - the shear production is multiplied by eboost in order
132      !!      to set the critic richardson number to ri_c (namelist parameter)
133      !!                   - the destruction by stratification term is multiplied
134      !!      by the Prandtl number (defined by an empirical funtion of the local
135      !!      Richardson number) if npdl=1 (namelist parameter)
136      !!      coefficient (zesh2):
137      !!      -3- Compute the now vertical eddy vicosity and diffusivity
138      !!      coefficients from en (before the time stepping) and zmxlm:
139      !!              avm = max( avtb, ediff*zmxlm*en^1/2 )
140      !!              avt = max( avmb, pdl*avm )  (pdl=1 if npdl=0)
141      !!              eav = max( avmb, avm )
142      !!      avt and avm are horizontally averaged to avoid numerical insta-
143      !!      bilities.
144      !!        N.B. The computation is done from jk=2 to jpkm1 except for
145      !!      en. Surface value of avt avmu avmv are set once a time to
146      !!      their background value in routine zdf_tke_init.
147      !!
148      !! ** Action  :   compute en (now turbulent kinetic energy)
149      !!                update avt, avmu, avmv (before vertical eddy coef.)
150      !!
151      !! References :
152      !!      Gaspar et al., jgr, 95, 1990,
153      !!      Blanke and Delecluse, jpo, 1991
154      !! History :
155      !!   6.0  !  91-03  (b. blanke)  Original code
156      !!   7.0  !  91-11  (G. Madec)   bug fix
157      !!   7.1  !  92-10  (G. Madec)   new mixing length and eav
158      !!   7.2  !  93-03  (M. Guyon)   symetrical conditions
159      !!   7.3  !  94-08  (G. Madec, M. Imbard)   npdl flag
160      !!   7.5  !  96-01  (G. Madec)   s-coordinates
161      !!   8.0  !  97-07  (G. Madec)   lbc
162      !!   8.1  !  99-01  (E. Stretta) new option for the mixing length
163      !!   8.5  !  02-08  (G. Madec)  ri_c and Free form, F90
164      !!   9.0  !  04-10  (C. Ethe )  1D configuration
165      !!----------------------------------------------------------------------
166      !! * Modules used
167      USE oce     , zwd   => ua,  &  ! use ua as workspace
168         &          zmxlm => ta,  &  ! use ta as workspace
169         &          zmxld => sa      ! use sa as workspace
170
171      !! * arguments
172      INTEGER, INTENT( in  ) ::   kt ! ocean time step
173
174      !! * local declarations
175      INTEGER ::   ji, jj, jk        ! dummy loop arguments
176      REAL(wp) ::   &
177         zmlmin, zbbrau,          &  ! temporary scalars
178         zfact1, zfact2, zfact3,  &  !
179         zrn2, zesurf,            &  !
180         ztx2, zty2, zav,         &  !
181         zcoef, zcof, zsh2,       &  !
182         zdku, zdkv, zpdl, zri,   &  !
183         zsqen, zesh2,            &  !
184         zemxl, zemlm, zemlp
185      !!--------------------------------------------------------------------
186
187      ! Initialization (first time-step only)
188      ! --------------
189      IF( kt == nit000  )   CALL zdf_tke_init
190
191      ! Local constant initialization
192      zmlmin = 1.e-8
193      zbbrau =  .5 * ebb / rau0
194      zfact1 = -.5 * rdt * efave
195      zfact2 = 1.5 * rdt * ediss
196      zfact3 = 0.5 * rdt * ediss
197
198
199      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
200      ! I.  Mixing length
201      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
202
203
204      ! Buoyancy length scale: l=sqrt(2*e/n**2)
205      ! ---------------------
206      zmxlm(:,:, 1 ) = zmlmin   ! surface set to the minimum value
207      zmxlm(:,:,jpk) = zmlmin   ! bottom  set to the minimum value
208!CDIR NOVERRCHK
209      DO jk = 2, jpkm1
210!CDIR NOVERRCHK
211         DO jj = 2, jpjm1
212!CDIR NOVERRCHK
213            DO ji = fs_2, fs_jpim1   ! vector opt.
214               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
215               zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  )
216            END DO
217         END DO
218      END DO
219
220
221      ! Physical limits for the mixing length
222      ! -------------------------------------
223      zmxld(:,:, 1 ) = zmlmin   ! surface set to the minimum value
224      zmxld(:,:,jpk) = zmlmin   ! bottom  set to the minimum value
225
226      SELECT CASE ( nmxl )
227
228      CASE ( 0 )           ! bounded by the distance to surface and bottom
229         DO jk = 2, jpkm1
230            DO jj = 2, jpjm1
231               DO ji = fs_2, fs_jpim1   ! vector opt.
232                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
233                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
234                  zmxlm(ji,jj,jk) = zemxl
235                  zmxld(ji,jj,jk) = zemxl
236               END DO
237            END DO
238         END DO
239
240      CASE ( 1 )           ! bounded by the vertical scale factor
241         DO jk = 2, jpkm1
242            DO jj = 2, jpjm1
243               DO ji = fs_2, fs_jpim1   ! vector opt.
244                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(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 ( 2 )           ! |dk[xml]| bounded by e3t :
252         DO jk = 2, jpkm1         ! from the surface to the bottom :
253            DO jj = 2, jpjm1
254               DO ji = fs_2, fs_jpim1   ! vector opt.
255                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
256               END DO
257            END DO
258         END DO
259         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
260            DO jj = 2, jpjm1
261               DO ji = fs_2, fs_jpim1   ! vector opt.
262                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
263                  zmxlm(ji,jj,jk) = zemxl
264                  zmxld(ji,jj,jk) = zemxl
265               END DO
266            END DO
267         END DO
268
269      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
270         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1   ! vector opt.
273                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
274               END DO
275            END DO
276         END DO
277         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
278            DO jj = 2, jpjm1
279               DO ji = fs_2, fs_jpim1   ! vector opt.
280                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
281               END DO
282            END DO
283         END DO
284!CDIR NOVERRCHK
285         DO jk = 2, jpkm1
286!CDIR NOVERRCHK
287            DO jj = 2, jpjm1
288!CDIR NOVERRCHK
289               DO ji = fs_2, fs_jpim1   ! vector opt.
290                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
291                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
292                  zmxlm(ji,jj,jk) = zemlm
293                  zmxld(ji,jj,jk) = zemlp
294               END DO
295            END DO
296         END DO
297
298      END SELECT
299
300# if defined key_cfg_1d
301      ! save mixing and dissipation turbulent length scales
302      e_dis(:,:,:) = zmxld(:,:,:)
303      e_mix(:,:,:) = zmxlm(:,:,:)
304# endif
305
306
307      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
308      ! II  Tubulent kinetic energy time stepping
309      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
310
311
312      ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
313      ! ---------------------------------------------------------------------
314!CDIR NOVERRCHK
315      DO jk = 2, jpkm1
316!CDIR NOVERRCHK
317         DO jj = 2, jpjm1
318!CDIR NOVERRCHK
319            DO ji = fs_2, fs_jpim1   ! vector opt.
320               zsqen = SQRT( en(ji,jj,jk) )
321               zav   = ediff * zmxlm(ji,jj,jk) * zsqen
322               avt  (ji,jj,jk) = MAX( zav, avtb(jk) ) * tmask(ji,jj,jk)
323               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
324               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
325            END DO
326         END DO
327      END DO
328
329      ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
330      ! -------------------------------------------------
331      ! en(1)   = ebb sqrt(taux^2+tauy^2) / rau0  (min value emin0)
332      ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
333!CDIR NOVERRCHK
334      DO jj = 2, jpjm1
335!CDIR NOVERRCHK
336         DO ji = fs_2, fs_jpim1   ! vector opt.
337            ztx2 = taux(ji-1,jj  ) + taux(ji,jj)
338            zty2 = tauy(ji  ,jj-1) + tauy(ji,jj)
339            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
340            en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
341            zmxlm(ji,jj,1  ) = avmb(1) * tmask(ji,jj,1)
342            zmxlm(ji,jj,jpk) = 0.e0
343         END DO
344      END DO
345
346      ! 3. Now Turbulent kinetic energy (output in en)
347      ! -------------------------------
348      ! Resolution of a tridiagonal linear system by a "methode de chasse"
349      ! computation from level 2 to jpkm1  (e(1) already computed and
350      ! e(jpk)=0 ).
351
352      SELECT CASE ( npdl )
353
354      CASE ( 0 )           ! No Prandtl number
355         DO jk = 2, jpkm1
356            DO jj = 2, jpjm1
357               DO ji = fs_2, fs_jpim1   ! vector opt.
358                  ! zesh2 = eboost * (du/dz)^2 - N^2
359                  zcoef = 0.5 / fse3w(ji,jj,jk)
360                  ! shear
361                  zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &
362                  &                - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
363                  zdkv = zcoef * (   vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
364                  &                - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
365                  ! coefficient (zesh2)
366                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)
367
368                  ! Matrix
369                  zcof = zfact1 * tmask(ji,jj,jk)
370                  ! lower diagonal
371                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
372                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
373                  ! upper diagonal
374                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
375                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
376                  ! diagonal
377                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
378                  ! right hand side in en
379                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
380                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
381               END DO
382            END DO
383         END DO
384
385      CASE ( 1 )           ! Prandtl number
386         DO jk = 2, jpkm1
387            DO jj = 2, jpjm1
388               DO ji = fs_2, fs_jpim1   ! vector opt.
389                  ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2
390                  zcoef = 0.5 / fse3w(ji,jj,jk)
391                  ! shear
392                  zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) &
393                  &                - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )   )
394                  zdkv = zcoef * (   vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1) &
395                  &                - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )   )
396                  ! square of vertical shear
397                  zsh2 = zdku * zdku + zdkv * zdkv
398                  ! local Richardson number
399                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
400# if defined key_cfg_1d
401                  ! save masked local Richardson number in zmxlm array
402                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)
403# endif
404                  ! Prandtl number
405                  zpdl = 1.0
406                  IF( zri >= 0.2 ) zpdl = 0.2 / zri
407                  zpdl = MAX( 0.1, zpdl )
408                  ! coefficient (esh2)
409                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
410
411                  ! Matrix
412                  zcof = zfact1 * tmask(ji,jj,jk)
413                  ! lower diagonal
414                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
415                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
416                  ! upper diagonal
417                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
418                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
419                  ! diagonal
420                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
421                  ! right hand side in en
422                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
423                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
424                  ! save masked Prandlt number in zmxlm array
425                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
426               END DO
427            END DO
428         END DO
429
430      END SELECT
431
432# if defined key_cfg_1d
433      !  save masked Prandlt number
434      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)
435      e_pdl(:,:,      1) = e_pdl(:,:,      2)
436      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)     
437# endif
438
439      ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
440      !!--------------------------------
441      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
442      DO jk = 3, jpkm1
443         DO jj = 2, jpjm1
444            DO ji = fs_2, fs_jpim1    ! vector opt.
445               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
446            END DO
447         END DO
448      END DO
449
450      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
451      DO jj = 2, jpjm1
452         DO ji = fs_2, fs_jpim1    ! vector opt.
453            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
454         END DO
455      END DO
456      DO jk = 3, jpkm1
457         DO jj = 2, jpjm1
458            DO ji = fs_2, fs_jpim1    ! vector opt.
459               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
460            END DO
461         END DO
462      END DO
463
464      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
465      DO jj = 2, jpjm1
466         DO ji = fs_2, fs_jpim1    ! vector opt.
467            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
468         END DO
469      END DO
470      DO jk = jpk-2, 2, -1
471         DO jj = 2, jpjm1
472            DO ji = fs_2, fs_jpim1    ! vector opt.
473               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
474            END DO
475         END DO
476      END DO
477
478      ! Save the result in en and set minimum value of tke : emin
479      DO jk = 2, jpkm1
480         DO jj = 2, jpjm1
481            DO ji = fs_2, fs_jpim1   ! vector opt.
482               en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
483            END DO
484         END DO
485      END DO
486
487      ! Lateral boundary conditions on ( avt, en )   (sign unchanged)
488      CALL lbc_lnk( en , 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
489
490
491      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
492      ! III.  Before vertical eddy vicosity and diffusivity coefficients
493      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
494
495      SELECT CASE ( nave )
496         
497      CASE ( 0 )                ! no horizontal average
498
499         ! Vertical eddy viscosity
500
501         DO jk = 2, jpkm1                                 ! Horizontal slab
502            DO jj = 2, jpjm1
503               DO ji = fs_2, fs_jpim1   ! vector opt.
504                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
505                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
506                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
507                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
508               END DO
509            END DO
510         END DO
511
512         ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
513         CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
514         
515      CASE ( 1 )                ! horizontal average
516
517         !                                                ( 1/2  1/2 )
518         ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1    1   )
519         !                      ( 1/2  1 1/2 )            ( 1/2  1/2 )
520         !           avmv = 1/4 ( 1/2  1 1/2 )     
521         
522!! caution vectopt_memory change the solution (last digit of the solver stat)
523#  if defined key_vectopt_memory
524         DO jk = 2, jpkm1                                 ! Horizontal slab
525            DO jj = 2, jpjm1
526               DO ji = fs_2, fs_jpim1   ! vector opt.
527                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
528                  &                 +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
529                  &                      +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
530
531                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
532                  &                 +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
533                  &                      +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
534               END DO
535            END DO
536         END DO
537#  else
538         DO jk = 2, jpkm1                                 ! Horizontal slab
539            DO jj = 2, jpjm1
540               DO ji = fs_2, fs_jpim1   ! vector opt.
541                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
542                  &              +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
543                  &                   +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)  &
544                  &       / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
545                  &              +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
546                  &                   +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
547
548                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
549                  &              +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
550                  &                   +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)  &
551                  &      /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
552                  &              +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
553                  &                   +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
554               END DO
555            END DO
556         END DO
557#  endif
558
559         ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
560         CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
561
562         ! Vertical eddy diffusivity
563         ! ------------------------------
564         !                                (1 2 1)
565         ! horizontal average  avt = 1/16 (2 4 2)
566         !                                (1 2 1)
567         DO jk = 2, jpkm1                                 ! Horizontal slab
568#  if defined key_vectopt_memory
569            DO jj = 2, jpjm1
570               DO ji = fs_2, fs_jpim1   ! vector opt.
571                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)    &
572                  &               + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
573               END DO
574            END DO
575#  else
576            DO jj = 2, jpjm1
577               DO ji = fs_2, fs_jpim1   ! vector opt.
578                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)   &
579                  &               + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
580                  &     / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
581                  &               + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
582               END DO
583            END DO
584#  endif
585         END DO
586
587      END SELECT
588
589      ! multiplied by the Prandtl number (npdl>1)
590      ! ----------------------------------------
591      IF( npdl == 1 ) THEN
592         DO jk = 2, jpkm1                                 ! Horizontal slab
593            DO jj = 2, jpjm1
594               DO ji = fs_2, fs_jpim1   ! vector opt.
595                  zpdl = zmxld(ji,jj,jk)
596                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
597               END DO
598            END DO
599         END DO
600      ENDIF
601
602      ! Minimum value on the eddy viscosity
603      ! ----------------------------------------
604      DO jk = 2, jpkm1                                 ! Horizontal slab
605         DO jj = 1, jpj
606            DO ji = 1, jpi
607               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
608               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
609            END DO
610         END DO
611      END DO
612
613      ! Lateral boundary conditions on avt  (sign unchanged)
614      ! ------------------------------=====
615      CALL lbc_lnk( avt, 'W', 1. )
616
617      IF(ln_ctl) THEN
618         CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
619         CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
620      ENDIF
621
622   END SUBROUTINE zdf_tke
623
624# endif
625
626   SUBROUTINE zdf_tke_init
627      !!----------------------------------------------------------------------
628      !!                  ***  ROUTINE zdf_tke_init  ***
629      !!                     
630      !! ** Purpose :   Initialization of the vertical eddy diffivity and
631      !!      viscosity when using a tke turbulent closure scheme
632      !!
633      !! ** Method  :   Read the namtke namelist and check the parameters
634      !!      called at the first timestep (nit000)
635      !!
636      !! ** input   :   Namlist namtke
637      !!
638      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
639      !!
640      !! history :
641      !!  8.5  ! 02-06 (G. Madec) original code
642      !!----------------------------------------------------------------------
643      !! * Module used
644      USE dynzdf_exp
645      USE trazdf_exp
646
647      !! * local declarations
648      !! caution vectopt_memory change the solution (last digit of the solver stat)
649# if defined key_vectopt_memory
650      INTEGER ::   ji, jj, jk, jit   ! dummy loop indices
651# else
652      INTEGER ::           jk, jit   ! dummy loop indices
653# endif
654
655      NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
656         ri_c, nitke, nmxl, npdl, nave, navb
657      !!----------------------------------------------------------------------
658
659      ! Read Namelist namtke : Turbulente Kinetic Energy
660      ! --------------------
661      REWIND ( numnam )
662      READ   ( numnam, namtke )
663
664      ! Compute boost associated with the Richardson critic
665      !     (control values: ri_c = 0.3   ==> eboost=1.25 for npdl=1 or 2)
666      !     (                ri_c = 0.222 ==> eboost=1.                  )
667      eboost = ri_c * ( 2. + ediss / ediff ) / 2.
668
669
670      ! Parameter control and print
671      ! ---------------------------
672      ! Control print
673      IF(lwp) THEN
674         WRITE(numout,*)
675         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme'
676         WRITE(numout,*) '~~~~~~~~~~~~'
677         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
678         WRITE(numout,*) '             restart with tke from no tke ln_rstke = ', ln_rstke
679         WRITE(numout,*) '             coef. to compute avt           ediff  = ', ediff
680         WRITE(numout,*) '             Kolmogoroff dissipation coef.  ediss  = ', ediss
681         WRITE(numout,*) '             tke surface input coef.        ebb    = ', ebb
682         WRITE(numout,*) '             tke diffusion coef.            efave  = ', efave
683         WRITE(numout,*) '             minimum value of tke           emin   = ', emin
684         WRITE(numout,*) '             surface minimum value of tke   emin0  = ', emin0
685         WRITE(numout,*) '             number of restart iter loops   nitke  = ', nitke
686         WRITE(numout,*) '             mixing length type             nmxl   = ', nmxl
687         WRITE(numout,*) '             prandl number flag             npdl   = ', npdl
688         WRITE(numout,*) '             horizontal average flag        nave   = ', nave
689         WRITE(numout,*) '             critic Richardson nb           ri_c   = ', ri_c
690         WRITE(numout,*) '                and its associated coeff.   eboost = ', eboost
691         WRITE(numout,*) '             constant background or profile navb   = ', navb
692         WRITE(numout,*)
693      ENDIF
694
695      ! Check nmxl and npdl values
696      IF( nmxl < 0 .OR. nmxl > 3 ) THEN
697         IF(lwp) WRITE(numout,cform_err)
698         IF(lwp) WRITE(numout,*) '          bad flag: nmxl is < 0 or > 3 '
699         nstop = nstop + 1
700      ENDIF
701      IF ( npdl < 0 .OR. npdl > 1 ) THEN
702         IF(lwp) WRITE(numout,cform_err)
703         IF(lwp) WRITE(numout,*) '          bad flag: npdl is < 0 or > 1 '
704         nstop = nstop + 1
705      ENDIF
706
707
708      ! Horizontal average : initialization of weighting arrays
709      ! -------------------
710     
711      SELECT CASE ( nave )
712
713      CASE ( 0 )                ! no horizontal average
714         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
715         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
716!! caution vectopt_memory change the solution (last digit of the solver stat)
717# if defined key_vectopt_memory
718         ! weighting mean arrays etmean, eumean and evmean
719         !           ( 1  1 )                                          ( 1 )
720         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
721         !                         
722         etmean(:,:,:) = 0.e0
723         eumean(:,:,:) = 0.e0
724         evmean(:,:,:) = 0.e0
725         
726         DO jk = 1, jpkm1
727            DO jj = 2, jpjm1
728               DO ji = fs_2, fs_jpim1   ! vector opt.
729                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
730                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
731                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
732                 
733                  eumean(ji,jj,jk) = umask(ji,jj,jk)                     &
734                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk)  )
735
736                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                     &
737                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  )
738               END DO
739            END DO
740         END DO
741# endif
742
743      CASE ( 1 )                ! horizontal average
744         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
745!! caution vectopt_memory change the solution (last digit of the solver stat)
746# if defined key_vectopt_memory
747         ! weighting mean arrays etmean, eumean and evmean
748         !           ( 1  1 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
749         ! avt = 1/4 ( 1  1 )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
750         !                                 ( 1/2  1/2 )
751         etmean(:,:,:) = 0.e0
752         eumean(:,:,:) = 0.e0
753         evmean(:,:,:) = 0.e0
754         
755         DO jk = 1, jpkm1
756            DO jj = 2, jpjm1
757               DO ji = fs_2, fs_jpim1   ! vector opt.
758                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
759                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
760                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
761                 
762                  eumean(ji,jj,jk) = umask(ji,jj,jk)                        &
763                  &  / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
764                  &       +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
765                  &              +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
766
767                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                        &
768                  &  / MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
769                  &       +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
770                  &              +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
771               END DO
772            END DO
773         END DO
774# endif
775
776      CASE DEFAULT
777         IF(lwp) WRITE(numout,cform_err)
778         IF(lwp) WRITE(numout,*) '          bad flag value for nave = ', nave
779         nstop = nstop + 1
780
781      END SELECT
782
783
784      ! Background eddy viscosity and diffusivity profil
785      ! ------------------------------------------------
786      IF( navb == 0 ) THEN
787         !   Define avmb, avtb from namelist parameter
788         avmb(:) = avm0
789         avtb(:) = avt0
790      ELSE
791         !   Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
792         avmb(:) = avm0
793         avtb(:) = 1.e-5 + 2.8e-8 * gdepw(:)   ! m2/s
794      ENDIF
795
796      !   Increase the background in the surface layers
797      avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
798      avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
799      avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
800      avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
801
802
803      ! Initialization of vertical eddy coef. to the background value
804      ! -------------------------------------------------------------
805      DO jk = 1, jpk
806         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
807         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
808         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
809      END DO
810
811
812      ! Initialization of turbulent kinetic energy ( en )
813      ! -------------------------------------------------
814      IF( ln_rstart ) THEN
815         ! no en field in the restart file, en set by iterative loop
816         IF( ln_rstke ) THEN
817            en (:,:,:) = emin * tmask(:,:,:)
818            DO jit = 2, nitke+1
819               CALL zdf_tke( jit )
820            END DO
821         ENDIF
822         ! otherwise en is already read in dtrlec called by inidtr
823      ELSE
824         ! no restart: en set to emin
825         en(:,:,:) = emin * tmask(:,:,:)
826      ENDIF
827
828   END SUBROUTINE zdf_tke_init
829
830#else
831   !!----------------------------------------------------------------------
832   !!   Dummy module :                                        NO TKE scheme
833   !!----------------------------------------------------------------------
834   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
835CONTAINS
836   SUBROUTINE zdf_tke( kt )          ! Empty routine
837      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
838   END SUBROUTINE zdf_tke
839#endif
840
841   !!======================================================================
842END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.