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 branches/dev_001_SBC/NEMO/OPA_SRC/ZDF – NEMO

source: branches/dev_001_SBC/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 762

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

phasing DYN directory, solisl.F90, restart.F90 and zdftke.F90 modules with the trunk, to get changeset 746
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
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 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   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   !! $Id$
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(utau^2 + vtau^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(utau^2+vtau^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 = utau(ji-1,jj  ) + utau(ji,jj)
310            zty2 = vtau(ji  ,jj-1) + vtau(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', ldstop = .FALSE. ) > 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', ldstop = .FALSE. ) > 0 )   &
805                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
806              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
807              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en set by iterative loop'
808              IF( lwp                )   WRITE(numout,*) ' =======             ========='
809              en (:,:,:) = emin * tmask(:,:,:)
810              DO jit = 2, nitke+1
811                 CALL zdf_tke( jit )
812              END DO
813           ENDIF
814        ELSE
815           en(:,:,:) = emin * tmask(:,:,:)      ! no restart: en set to emin
816        ENDIF
817     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
818        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
819     ENDIF
820     !
821   END SUBROUTINE tke_rst
822
823#else
824   !!----------------------------------------------------------------------
825   !!   Dummy module :                                        NO TKE scheme
826   !!----------------------------------------------------------------------
827   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
828CONTAINS
829   SUBROUTINE zdf_tke( kt )          ! Empty routine
830      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
831   END SUBROUTINE zdf_tke
832#endif
833
834   !!======================================================================
835END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.