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

Last change on this file since 508 was 508, checked in by opalod, 18 years ago

nemo_v1_update_071:RB: add iom for restart and reorganization of restart
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.6 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   NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
72      &             ri_c, nitke, nmxl, npdl, nave, navb
73
74# if defined key_cfg_1d
75   !                                                                   ! 1D cfg only
76   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix,      &  ! dissipation and mixing turbulent lengh scales
77      &                                          e_pdl, e_ric          ! prandl and local Richardson numbers
78#endif
79
80   !! * Substitutions
81#  include "domzgr_substitute.h90"
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !!   OPA 9.0 , LOCEAN-IPSL (2006)
85   !! $Header$
86   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
87   !!----------------------------------------------------------------------
88
89CONTAINS
90
91   SUBROUTINE zdf_tke( kt )
92      !!----------------------------------------------------------------------
93      !!                   ***  ROUTINE zdf_tke  ***
94      !!
95      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
96      !!      coefficients using a 1.5 turbulent closure scheme.
97      !!
98      !! ** Method  :   The time evolution of the turbulent kinetic energy
99      !!      (tke) is computed from a prognostic equation :
100      !!         d(en)/dt = eboost eav (d(u)/dz)**2       ! shear production
101      !!                  + d( efave eav d(en)/dz )/dz    ! diffusion of tke
102      !!                  + grav/rau0 pdl eav d(rau)/dz   ! stratif. destruc.
103      !!                  - ediss / emxl en**(2/3)        ! dissipation
104      !!      with the boundary conditions:
105      !!         surface: en = max( emin0,ebb sqrt(taux^2 + tauy^2) )
106      !!         bottom : en = emin
107      !!      -1- The dissipation and mixing turbulent lengh scales are computed
108      !!      from the usual diagnostic buoyancy length scale: 
109      !!         mxl= 1/(sqrt(en)/N)  WHERE N is the brunt-vaisala frequency
110      !!      Four cases :
111      !!         nmxl=0 : mxl bounded by the distance to surface and bottom.
112      !!                  zmxld = zmxlm = mxl
113      !!         nmxl=1 : mxl bounded by the vertical scale factor.
114      !!                  zmxld = zmxlm = mxl
115      !!         nmxl=2 : mxl bounded such that the vertical derivative of mxl
116      !!                  is less than 1 (|d/dz(xml)|<1).
117      !!                  zmxld = zmxlm = mxl
118      !!         nmxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface
119      !!                        to the bottom
120      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom
121      !!                        to the surface
122      !!                  zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
123      !!      -2- Compute the now Turbulent kinetic energy. The time differencing
124      !!      is implicit for vertical diffusion term, linearized for kolmo-
125      !!      goroff dissipation term, and explicit forward for both buoyancy
126      !!      and dynamic production terms. Thus a tridiagonal linear system is
127      !!      solved.
128      !!         Note that - the shear production is multiplied by eboost in order
129      !!      to set the critic richardson number to ri_c (namelist parameter)
130      !!                   - the destruction by stratification term is multiplied
131      !!      by the Prandtl number (defined by an empirical funtion of the local
132      !!      Richardson number) if npdl=1 (namelist parameter)
133      !!      coefficient (zesh2):
134      !!      -3- Compute the now vertical eddy vicosity and diffusivity
135      !!      coefficients from en (before the time stepping) and zmxlm:
136      !!              avm = max( avtb, ediff*zmxlm*en^1/2 )
137      !!              avt = max( avmb, pdl*avm )  (pdl=1 if npdl=0)
138      !!              eav = max( avmb, avm )
139      !!      avt and avm are horizontally averaged to avoid numerical insta-
140      !!      bilities.
141      !!        N.B. The computation is done from jk=2 to jpkm1 except for
142      !!      en. Surface value of avt avmu avmv are set once a time to
143      !!      their background value in routine zdf_tke_init.
144      !!
145      !! ** Action  :   compute en (now turbulent kinetic energy)
146      !!                update avt, avmu, avmv (before vertical eddy coef.)
147      !!
148      !! References : Gaspar et al., jgr, 95, 1990,
149      !!              Blanke and Delecluse, jpo, 1991
150      !!----------------------------------------------------------------------
151      USE oce     , zwd   => ua,  &  ! use ua as workspace
152         &          zmxlm => ta,  &  ! use ta as workspace
153         &          zmxld => sa      ! use sa as workspace
154      !
155      INTEGER, INTENT(in) ::   kt ! ocean time step
156      !
157      INTEGER  ::   ji, jj, jk                  ! dummy loop arguments
158      REAL(wp) ::   zmlmin, zbbrau,          &  ! temporary scalars
159         &          zfact1, zfact2, zfact3,  &  !
160         &          zrn2, zesurf,            &  !
161         &          ztx2, zty2, zav,         &  !
162         &          zcoef, zcof, zsh2,       &  !
163         &          zdku, zdkv, zpdl, zri,   &  !
164         &          zsqen, zesh2,            &  !
165         &          zemxl, zemlm, zemlp
166      !!--------------------------------------------------------------------
167
168      IF( kt == nit000  )   CALL zdf_tke_init      ! Initialization (first time-step only)
169
170      !                                            ! Local constant initialization
171      zmlmin = 1.e-8
172      zbbrau =  .5 * ebb / rau0
173      zfact1 = -.5 * rdt * efave
174      zfact2 = 1.5 * rdt * ediss
175      zfact3 = 0.5 * rdt * ediss
176
177      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
178      ! I.  Mixing length
179      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
180
181      ! Buoyancy length scale: l=sqrt(2*e/n**2)
182      ! ---------------------
183      zmxlm(:,:, 1 ) = zmlmin   ! surface set to the minimum value
184      zmxlm(:,:,jpk) = zmlmin   ! bottom  set to the minimum value
185!CDIR NOVERRCHK
186      DO jk = 2, jpkm1
187!CDIR NOVERRCHK
188         DO jj = 2, jpjm1
189!CDIR NOVERRCHK
190            DO ji = fs_2, fs_jpim1   ! vector opt.
191               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
192               zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  )
193            END DO
194         END DO
195      END DO
196
197      ! Physical limits for the mixing length
198      ! -------------------------------------
199      zmxld(:,:, 1 ) = zmlmin   ! surface set to the minimum value
200      zmxld(:,:,jpk) = zmlmin   ! bottom  set to the minimum value
201
202      SELECT CASE ( nmxl )
203
204      CASE ( 0 )           ! bounded by the distance to surface and bottom
205         DO jk = 2, jpkm1
206            DO jj = 2, jpjm1
207               DO ji = fs_2, fs_jpim1   ! vector opt.
208                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
209                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
210                  zmxlm(ji,jj,jk) = zemxl
211                  zmxld(ji,jj,jk) = zemxl
212               END DO
213            END DO
214         END DO
215
216      CASE ( 1 )           ! bounded by the vertical scale factor
217         DO jk = 2, jpkm1
218            DO jj = 2, jpjm1
219               DO ji = fs_2, fs_jpim1   ! vector opt.
220                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
221                  zmxlm(ji,jj,jk) = zemxl
222                  zmxld(ji,jj,jk) = zemxl
223               END DO
224            END DO
225         END DO
226
227      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
228         DO jk = 2, jpkm1         ! from the surface to the bottom :
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
232               END DO
233            END DO
234         END DO
235         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
236            DO jj = 2, jpjm1
237               DO ji = fs_2, fs_jpim1   ! vector opt.
238                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
239                  zmxlm(ji,jj,jk) = zemxl
240                  zmxld(ji,jj,jk) = zemxl
241               END DO
242            END DO
243         END DO
244
245      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
246         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
247            DO jj = 2, jpjm1
248               DO ji = fs_2, fs_jpim1   ! vector opt.
249                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
250               END DO
251            END DO
252         END DO
253         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
254            DO jj = 2, jpjm1
255               DO ji = fs_2, fs_jpim1   ! vector opt.
256                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
257               END DO
258            END DO
259         END DO
260!CDIR NOVERRCHK
261         DO jk = 2, jpkm1
262!CDIR NOVERRCHK
263            DO jj = 2, jpjm1
264!CDIR NOVERRCHK
265               DO ji = fs_2, fs_jpim1   ! vector opt.
266                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
267                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
268                  zmxlm(ji,jj,jk) = zemlm
269                  zmxld(ji,jj,jk) = zemlp
270               END DO
271            END DO
272         END DO
273
274      END SELECT
275
276# if defined key_cfg_1d
277      ! save mixing and dissipation turbulent length scales
278      e_dis(:,:,:) = zmxld(:,:,:)
279      e_mix(:,:,:) = zmxlm(:,:,:)
280# endif
281
282      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
283      ! II  Tubulent kinetic energy time stepping
284      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
285
286      ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
287      ! ---------------------------------------------------------------------
288!CDIR NOVERRCHK
289      DO jk = 2, jpkm1
290!CDIR NOVERRCHK
291         DO jj = 2, jpjm1
292!CDIR NOVERRCHK
293            DO ji = fs_2, fs_jpim1   ! vector opt.
294               zsqen = SQRT( en(ji,jj,jk) )
295               zav   = ediff * zmxlm(ji,jj,jk) * zsqen
296               avt  (ji,jj,jk) = MAX( zav, avtb(jk) ) * tmask(ji,jj,jk)
297               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
298               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
299            END DO
300         END DO
301      END DO
302
303      ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
304      ! -------------------------------------------------
305      ! en(1)   = ebb sqrt(taux^2+tauy^2) / rau0  (min value emin0)
306      ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
307!CDIR NOVERRCHK
308      DO jj = 2, jpjm1
309!CDIR NOVERRCHK
310         DO ji = fs_2, fs_jpim1   ! vector opt.
311            ztx2 = taux(ji-1,jj  ) + taux(ji,jj)
312            zty2 = tauy(ji  ,jj-1) + tauy(ji,jj)
313            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
314            en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
315            zmxlm(ji,jj,1  ) = avmb(1) * tmask(ji,jj,1)
316            zmxlm(ji,jj,jpk) = 0.e0
317         END DO
318      END DO
319
320      ! 3. Now Turbulent kinetic energy (output in en)
321      ! -------------------------------
322      ! Resolution of a tridiagonal linear system by a "methode de chasse"
323      ! computation from level 2 to jpkm1  (e(1) already computed and
324      ! e(jpk)=0 ).
325
326      SELECT CASE ( npdl )
327
328      CASE ( 0 )           ! No Prandtl number
329         DO jk = 2, jpkm1
330            DO jj = 2, jpjm1
331               DO ji = fs_2, fs_jpim1   ! vector opt.
332                  ! zesh2 = eboost * (du/dz)^2 - N^2
333                  zcoef = 0.5 / fse3w(ji,jj,jk)
334                  ! shear
335                  zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &
336                  &                - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
337                  zdkv = zcoef * (   vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
338                  &                - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
339                  ! coefficient (zesh2)
340                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)
341
342                  ! Matrix
343                  zcof = zfact1 * tmask(ji,jj,jk)
344                  ! lower diagonal
345                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
346                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
347                  ! upper diagonal
348                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
349                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
350                  ! diagonal
351                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
352                  ! right hand side in en
353                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
354                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
355               END DO
356            END DO
357         END DO
358
359      CASE ( 1 )           ! Prandtl number
360         DO jk = 2, jpkm1
361            DO jj = 2, jpjm1
362               DO ji = fs_2, fs_jpim1   ! vector opt.
363                  ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2
364                  zcoef = 0.5 / fse3w(ji,jj,jk)
365                  ! shear
366                  zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) &
367                  &                - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )   )
368                  zdkv = zcoef * (   vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1) &
369                  &                - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )   )
370                  ! square of vertical shear
371                  zsh2 = zdku * zdku + zdkv * zdkv
372                  ! local Richardson number
373                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
374# if defined key_cfg_1d
375                  ! save masked local Richardson number in zmxlm array
376                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)
377# endif
378                  ! Prandtl number
379                  zpdl = 1.0
380                  IF( zri >= 0.2 ) zpdl = 0.2 / zri
381                  zpdl = MAX( 0.1, zpdl )
382                  ! coefficient (esh2)
383                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
384
385                  ! Matrix
386                  zcof = zfact1 * tmask(ji,jj,jk)
387                  ! lower diagonal
388                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
389                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
390                  ! upper diagonal
391                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
392                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
393                  ! diagonal
394                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
395                  ! right hand side in en
396                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
397                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
398                  ! save masked Prandlt number in zmxlm array
399                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
400               END DO
401            END DO
402         END DO
403
404      END SELECT
405
406# if defined key_cfg_1d
407      !  save masked Prandlt number
408      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)
409      e_pdl(:,:,      1) = e_pdl(:,:,      2)
410      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)     
411# endif
412
413      ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
414      !!--------------------------------
415      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
416      DO jk = 3, jpkm1
417         DO jj = 2, jpjm1
418            DO ji = fs_2, fs_jpim1    ! vector opt.
419               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
420            END DO
421         END DO
422      END DO
423
424      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
425      DO jj = 2, jpjm1
426         DO ji = fs_2, fs_jpim1    ! vector opt.
427            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
428         END DO
429      END DO
430      DO jk = 3, jpkm1
431         DO jj = 2, jpjm1
432            DO ji = fs_2, fs_jpim1    ! vector opt.
433               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
434            END DO
435         END DO
436      END DO
437
438      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
439      DO jj = 2, jpjm1
440         DO ji = fs_2, fs_jpim1    ! vector opt.
441            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
442         END DO
443      END DO
444      DO jk = jpk-2, 2, -1
445         DO jj = 2, jpjm1
446            DO ji = fs_2, fs_jpim1    ! vector opt.
447               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
448            END DO
449         END DO
450      END DO
451
452      ! Save the result in en and set minimum value of tke : emin
453      DO jk = 2, jpkm1
454         DO jj = 2, jpjm1
455            DO ji = fs_2, fs_jpim1   ! vector opt.
456               en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
457            END DO
458         END DO
459      END DO
460
461      ! Lateral boundary conditions on ( avt, en )   (sign unchanged)
462      CALL lbc_lnk( en , 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
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      ! write en in restart file
591      ! ------------------------
592      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
593
594      IF(ln_ctl) THEN
595         CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
596         CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
597      ENDIF
598
599   END SUBROUTINE zdf_tke
600
601
602   SUBROUTINE zdf_tke_init
603      !!----------------------------------------------------------------------
604      !!                  ***  ROUTINE zdf_tke_init  ***
605      !!                     
606      !! ** Purpose :   Initialization of the vertical eddy diffivity and
607      !!      viscosity when using a tke turbulent closure scheme
608      !!
609      !! ** Method  :   Read the namtke namelist and check the parameters
610      !!      called at the first timestep (nit000)
611      !!
612      !! ** input   :   Namlist namtke
613      !!
614      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
615      !!
616      !!----------------------------------------------------------------------
617      USE dynzdf_exp
618      USE trazdf_exp
619      !
620# if defined key_vectopt_memory
621      ! caution vectopt_memory change the solution (last digit of the solver stat)
622      INTEGER ::   ji, jj, jk   ! dummy loop indices
623# else
624      INTEGER ::           jk   ! dummy loop indices
625# endif
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' ) > 0 .AND. .NOT.(ln_rstke) ) THEN
800              CALL iom_get( numror, jpdom_local, 'en', en )
801           ELSE
802              IF(lwp .AND. iom_varid(numror,'en') > 0 ) WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
803              IF(lwp .AND. ln_rstke ) WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
804              IF(lwp) WRITE(numout,*) ' ===>>>> : en set by iterative loop'
805              IF(lwp) WRITE(numout,*) ' =======             ========='
806              en (:,:,:) = emin * tmask(:,:,:)
807              DO jit = 2, nitke+1
808                 CALL zdf_tke( jit )
809              END DO
810           ENDIF
811        ELSE
812           en(:,:,:) = emin * tmask(:,:,:)      ! no restart: en set to emin
813        ENDIF
814     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
815        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
816     ENDIF
817     !
818   END SUBROUTINE tke_rst
819
820#else
821   !!----------------------------------------------------------------------
822   !!   Dummy module :                                        NO TKE scheme
823   !!----------------------------------------------------------------------
824   PUBLIC, LOGICAL, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
825CONTAINS
826   SUBROUTINE zdf_tke( kt )          ! Empty routine
827      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
828   END SUBROUTINE zdf_tke
829#endif
830
831   !!======================================================================
832END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.