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

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

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