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

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

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