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

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

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