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_jki.F90 in trunk/NEMO/OPA_SRC/ZDF – NEMO

source: trunk/NEMO/OPA_SRC/ZDF/zdftke_jki.F90 @ 467

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

nemo_v1_update_054:RB: take into account tracers and dynamics reorganization, change atsk to jki
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.9 KB
Line 
1MODULE zdftke_jki
2   !!======================================================================
3   !!                       ***  MODULE  zdftke_jki  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7#if defined key_zdftke   ||   defined key_esopa
8   !!----------------------------------------------------------------------
9   !!   'key_zdftke'                                             TKE scheme
10   !!----------------------------------------------------------------------
11   !!   zdf_tke_jki  : update momentum and tracer Kz from a tke scheme
12   !!                  j-k-i loops for NEC autotasking or OpenMP
13   !!                  (used if 'key_mpp_omp' activated)
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce             ! ocean dynamics and active tracers
17   USE dom_oce         ! ocean space and time domain
18   USE zdf_oce         ! ocean vertical physics
19   USE zdftke          ! ???
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE phycst          ! physical constants
23   USE taumod          ! surface stress
24   USE prtctl          ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC zdf_tke_jki   ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !!   OPA 9.0 , LOCEAN-IPSL (2005)
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE zdf_tke_jki( kt )
42      !!----------------------------------------------------------------------
43      !!                   ***  ROUTINE zdf_tke  ***
44      !!                     
45      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
46      !!      coefficients using a 1.5 turbulent closure scheme.
47      !!
48      !! ** Method  :   The time evolution of the turbulent kinetic energy
49      !!      (tke) is computed from a prognostic equation :
50      !!         d(en)/dt = eboost eav (d(u)/dz)**2       ! shear production
51      !!                  + d( efave eav d(en)/dz )/dz    ! diffusion of tke
52      !!                  + g/rau0 pdl eav d(rau)/dz      ! stratif. destruc.
53      !!                  - ediss / emxl en**(2/3)        ! dissipation
54      !!      with the boundary conditions:
55      !!         surface: en = max( emin0,ebb sqrt(taux^2 + tauy^2) )
56      !!         bottom : en = emin
57      !!      -1- The dissipation and mixing turbulent lengh scales are computed
58      !!      from the usual diagnostic buoyancy length scale: 
59      !!         mxl= 1/(sqrt(en)/N)  WHERE N is the brunt-vaisala frequency
60      !!      Four cases :
61      !!         nmxl=0 : mxl bounded by the distance to surface and bottom.
62      !!                  zmxld = zmxlm = mxl
63      !!         nmxl=1 : mxl bounded by the vertical scale factor.
64      !!                  zmxld = zmxlm = mxl
65      !!         nmxl=2 : mxl bounded such that the vertical derivative of mxl
66      !!                  is less than 1 (|d/dz(xml)|<1).
67      !!                  zmxld = zmxlm = mxl
68      !!         nmxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface
69      !!                        to the bottom
70      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom
71      !!                        to the surface
72      !!                  zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
73      !!      -2- Compute the now Turbulent kinetic energy. The time differencing
74      !!      is implicit for vertical diffusion term, linearized for kolmo-
75      !!      goroff dissipation term, and explicit forward for both buoyancy
76      !!      and dynamic production terms. Thus a tridiagonal linear system is
77      !!      solved.
78      !!         Note that - the shear production is multiplied by eboost in order
79      !!      to set the critic richardson number to ri_c (namelist parameter)
80      !!                   - the destruction by stratification term is multiplied
81      !!      by the Prandtl number (defined by an empirical funtion of the local
82      !!      Richardson number) if npdl=1 (namelist parameter)
83      !!      coefficient (zesh2):
84      !!      -3- Compute the now vertical eddy vicosity and diffusivity
85      !!      coefficients from en (before the time stepping) and zmxlm:
86      !!              avm = max( avtb, ediff*zmxlm*en^1/2 )
87      !!              avt = max( avmb, pdl*avm )  (pdl=1 if npdl=0)
88      !!              eav = max( avmb, avm )
89      !!      avt and avm are horizontally averaged to avoid numerical insta-
90      !!      bilities.
91      !!        N.B. The computation is done from jk=2 to jpkm1 except for
92      !!      en. Surface value of avt avmu avmv are set once a time to
93      !!      their background value in routine zdftke_init.
94      !!
95      !! ** Action :   compute en (now turbulent kinetic energy)
96      !!               update avt, avmu, avmv (before vertical eddy coeff.)
97      !!
98      !! References :
99      !!      Gaspar et al., jgr, 95, 1990,
100      !!      Blanke and Delecluse, jpo, 1991
101      !! History :
102      !!   9.0  !  02-08  (G. Madec)  autotasking optimization
103      !!----------------------------------------------------------------------
104      !! * Modules used
105      USE oce       , zwd   => ua,  &  ! use ua as workspace
106                      zmxlm => ta,  &  ! use ta as workspace
107                      zmxld => sa      ! use sa as workspace
108      !! * arguments
109      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
110
111      !! * local declarations
112      INTEGER ::   ji, jj, jk          ! dummy loop arguments
113      REAL(wp) ::                   &
114         zmlmin, zbbrau,            &  ! temporary scalars
115         zfact1, zfact2, zfact3,    &  !
116         zrn2, zesurf,              &  !
117         ztx2, zty2, zav,           &  !
118         zcoef, zcof, zsh2,         &  !
119         zdku, zdkv, zpdl, zri,     &  !
120         zsqen, zesh2,              &  !
121         zemxl, zemlm, zemlp
122      !!--------------------------------------------------------------------
123      !!  OPA 9.0 , LOCEAN-IPSL (2005)
124      !! $Header$
125      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
126      !!--------------------------------------------------------------------
127
128
129      ! 0. Initialization
130      !    --------------
131      IF( kt == nit000  )   CALL zdf_tke_init
132
133      ! Local constant
134      zmlmin = 1.e-8
135      zbbrau =  .5 * ebb / rau0
136      zfact1 = -.5 * rdt * efave
137      zfact2 = 1.5 * rdt * ediss
138      zfact3 = 0.5 * rdt * ediss
139
140
141      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
142      ! I.  Mixing length
143      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
144
145      !                                                ! ===============
146      DO jj = 2, jpjm1                                 !  Vertical slab
147         !                                             ! ===============
148         ! Buoyancy length scale: l=sqrt(2*e/n**2)
149         ! ---------------------
150         zmxlm(:,jj, 1 ) = zmlmin   ! surface set to the minimum value
151         zmxlm(:,jj,jpk) = zmlmin   ! bottom  set to the minimum value
152!CDIR NOVERRCHK
153         DO jk = 2, jpkm1
154!CDIR NOVERRCHK
155            DO ji = 2, jpim1
156               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
157               zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  )
158            END DO
159         END DO
160
161
162         ! Physical limits for the mixing length
163         ! -------------------------------------
164         zmxld(:,jj, 1 ) = zmlmin   ! surface set to the minimum value
165         zmxld(:,jj,jpk) = zmlmin   ! bottom  set to the minimum value
166
167         SELECT CASE ( nmxl )
168
169         CASE ( 0 )           ! bounded by the distance to surface and bottom
170
171            DO jk = 2, jpkm1
172               DO ji = 2, jpim1
173                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
174                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
175                  zmxlm(ji,jj,jk) = zemxl
176                  zmxld(ji,jj,jk) = zemxl
177               END DO
178            END DO
179
180         CASE ( 1 )           ! bounded by the vertical scale factor
181
182            DO jk = 2, jpkm1
183               DO ji = 2, jpim1
184                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
185                  zmxlm(ji,jj,jk) = zemxl
186                  zmxld(ji,jj,jk) = zemxl
187               END DO
188            END DO
189
190         CASE ( 2 )           ! |dk[xml]| bounded by e3t :
191
192            DO jk = 2, jpk           ! from the surface to the bottom :
193               DO ji = 2, jpim1
194                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
195               END DO
196            END DO
197            DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
198               DO ji = 2, jpim1
199                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
200                  zmxlm(ji,jj,jk) = zemxl
201                  zmxld(ji,jj,jk) = zemxl
202               END DO
203            END DO
204
205         CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
206
207            DO jk = 2, jpk           ! from the surface to the bottom : lup
208               DO ji = 2, jpim1
209                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
210               END DO
211            END DO
212            DO jk = jpkm1, 1, -1     ! from the bottom to the surface : ldown
213               DO ji = 2, jpim1
214                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
215               END DO
216            END DO
217!CDIR NOVERRCHK
218            DO jk = 1, jpk
219!CDIR NOVERRCHK
220               DO ji = 2, jpim1
221                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
222                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
223                  zmxlm(ji,jj,jk) = zemlm
224                  zmxld(ji,jj,jk) = zemlp
225               END DO
226            END DO
227
228         END SELECT
229
230
231         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
232         ! II  Tubulent kinetic energy time stepping
233         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
234
235
236         ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
237         ! ---------------------------------------------------------------------
238!CDIR NOVERRCHK
239         DO jk = 2, jpkm1
240!CDIR NOVERRCHK
241            DO ji = 2, jpim1
242               zsqen = SQRT( en(ji,jj,jk) )
243               zav   = ediff * zmxlm(ji,jj,jk) * zsqen
244               avt  (ji,jj,jk) = MAX( zav, avtb(jk) ) * tmask(ji,jj,jk)
245               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
246               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
247            END DO
248         END DO
249
250
251         ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
252         ! -------------------------------------------------
253         ! en(1)   = ebb sqrt(taux^2+tauy^2) / rau0  (min value emin0)
254         ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
255!CDIR NOVERRCHK
256         DO ji = 2, jpim1
257            ztx2 = taux(ji-1,jj  ) + taux(ji,jj)
258            zty2 = tauy(ji  ,jj-1) + tauy(ji,jj)
259            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
260            en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
261            zmxlm(ji,jj,1  ) = avmb(1) * tmask(ji,jj,1)
262            zmxlm(ji,jj,jpk) = 0.e0
263         END DO
264
265
266         ! 3. Now Turbulent kinetic energy (output in en)
267         ! -------------------------------
268         ! Resolution of a tridiagonal linear system by a "methode de chasse"
269         ! computation from level 2 to jpkm1  (e(1) already computed and
270         ! e(jpk)=0 ).
271
272         SELECT CASE ( npdl )
273
274         CASE ( 0 )           ! No Prandtl number
275            DO jk = 2, jpkm1
276               DO ji = 2, jpim1
277                  ! zesh2 = eboost * (du/dz)^2 - N^2
278                  zcoef = 0.5 / fse3w(ji,jj,jk)
279                  ! shear
280                  zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &
281                  &                - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
282                  zdkv = zcoef * (   vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
283                  &                - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
284                  ! coefficient (zesh2)
285                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)
286
287                  ! Matrix
288                  zcof = zfact1 * tmask(ji,jj,jk)
289                  ! lower diagonal
290                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
291                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
292                  ! upper diagonal
293                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
294                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
295                  ! diagonal
296                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
297                  ! right hand side in en
298                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
299                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
300               END DO
301            END DO
302
303         CASE ( 1 )           ! Prandtl number
304            DO jk = 2, jpkm1
305               DO ji = 2, jpim1
306                  ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2
307                  zcoef = 0.5 / fse3w(ji,jj,jk)
308                  ! shear
309                  zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) &
310                  &                - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )   )
311                  zdkv = zcoef * (   vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1) &
312                  &                - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )   )
313                  ! square of vertical shear
314                  zsh2 = zdku * zdku + zdkv * zdkv
315                  ! Prandtl number
316                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
317                  zpdl = 1.0
318                  IF( zri >= 0.2 ) zpdl = 0.2 / zri
319                  zpdl = MAX( 0.1, zpdl )
320                  ! coefficient (esh2)
321                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
322
323                  ! Matrix
324                  zcof = zfact1 * tmask(ji,jj,jk)
325                  ! lower diagonal
326                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
327                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
328                  ! upper diagonal
329                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
330                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
331                  ! diagonal
332                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
333                  ! right hand side in en
334                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
335                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
336                  ! save masked Prandlt number in zmxlm array
337                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
338               END DO
339            END DO
340
341         END SELECT
342
343
344         ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
345         !---------------------------------
346
347         ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
348         DO jk = 3, jpkm1
349            DO ji = 2, jpim1
350               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
351            END DO
352         END DO
353
354         ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
355         DO ji = 2, jpim1
356            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
357         END DO
358         DO jk = 3, jpkm1
359            DO ji = 2, jpim1
360               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
361            END DO
362         END DO
363
364         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
365         DO ji = 2, jpim1
366            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
367         END DO
368         DO jk = jpk-2, 2, -1
369            DO ji = 2, jpim1
370               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
371            END DO
372         END DO
373
374         ! Save the result in en and set minimum value of tke : emin
375         DO jk = 2, jpkm1
376            DO ji = 2, jpim1
377               en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
378            END DO
379         END DO
380         !                                             ! ===============
381      END DO                                           !   End of slab
382      !                                                ! ===============
383
384      ! Lateral boundary conditions on ( avt, en )  (sign unchanged)
385      ! --------------------------------=========
386      CALL lbc_lnk( avt, 'W', 1. )   ;   CALL lbc_lnk( en , 'W', 1. )
387
388
389      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
390      ! III.  Before vertical eddy vicosity and diffusivity coefficients
391      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
392
393      !                                                ! ===============
394      DO jk = 2, jpkm1                                 ! Horizontal slab
395         !                                             ! ===============
396         SELECT CASE ( nave )
397           
398         CASE ( 0 )                ! no horizontal average
399
400            ! Vertical eddy viscosity
401
402            DO jj = 2, jpjm1
403               DO ji = fs_2, fs_jpim1   ! vector opt.
404                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
405                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
406                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
407                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
408               END DO
409            END DO
410
411           
412         CASE ( 1 )                ! horizontal average
413
414            !                                                ( 1/2  1/2 )
415            ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1    1   )
416            !                      ( 1/2  1 1/2 )            ( 1/2  1/2 )
417            !           avmv = 1/4 ( 1/2  1 1/2 )     
418           
419!! caution vectopt_memory change the solution (last digit of the solver stat)
420#  if defined key_vectopt_memory
421            DO jj = 2, jpjm1
422               DO ji = fs_2, fs_jpim1   ! vector opt.
423                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
424                  &                 +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
425                  &                      +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
426
427                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
428                  &                 +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
429                  &                      +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
430               END DO
431            END DO
432#  else
433            DO jj = 2, jpjm1
434               DO ji = fs_2, fs_jpim1   ! vector opt.
435                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
436                  &              +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
437                  &                   +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)  &
438                  &       / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
439                  &              +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
440                  &                   +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
441
442                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
443                  &              +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
444                  &                   +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)  &
445                  &      /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
446                  &              +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
447                  &                   +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
448               END DO
449            END DO
450#  endif
451         END SELECT
452         !                                             ! ===============
453      END DO                                           !   End of slab
454      !                                                ! ===============
455
456      ! Lateral boundary conditions (avmu,avmv)  (sign unchanged)
457      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
458
459      !                                                ! ===============
460      DO jk = 2, jpkm1                                 ! Horizontal slab
461         !                                             ! ===============
462         SELECT CASE ( nave )
463           
464         CASE ( 1 )                ! horizontal average
465
466            ! Vertical eddy diffusivity
467            ! ------------------------------
468            !                                (1 2 1)
469            ! horizontal average  avt = 1/16 (2 4 2)
470            !                                (1 2 1)
471!! caution vectopt_memory change the solution (last digit of the solver stat)
472#  if defined key_vectopt_memory
473            DO jj = 2, jpjm1
474               DO ji = fs_2, fs_jpim1   ! vector opt.
475                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)    &
476                  &               + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
477               END DO
478            END DO
479#  else
480            DO jj = 2, jpjm1
481               DO ji = fs_2, fs_jpim1   ! vector opt.
482                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)   &
483                  &               + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
484                  &     / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
485                  &               + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
486               END DO
487            END DO
488#  endif
489         END SELECT
490
491
492         ! multiplied by the Prandtl number (npdl>1)
493         ! ----------------------------------------
494         IF( npdl == 1 ) THEN
495            DO jj = 2, jpjm1
496               DO ji = fs_2, fs_jpim1   ! vector opt.
497                  zpdl = zmxld(ji,jj,jk)
498                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
499               END DO
500            END DO
501         ENDIF
502
503         ! Minimum value on the eddy viscosity
504         ! ----------------------------------------
505         DO jj = 1, jpj
506            DO ji = 1, jpi
507               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
508               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
509            END DO
510         END DO
511         !                                             ! ===============
512      END DO                                           !   End of slab
513      !                                                ! ===============
514
515
516      ! Lateral boundary conditions on avt   (W-point (=T), sign unchanged)
517      ! ------------------------------=====
518      CALL lbc_lnk( avt, 'W', 1. )
519
520      IF(ln_ctl) THEN
521         CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
522         CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
523      ENDIF
524
525   END SUBROUTINE zdf_tke_jki
526
527#else
528   !!----------------------------------------------------------------------
529   !!   Dummy module :                                        NO TKE scheme
530   !!----------------------------------------------------------------------
531CONTAINS
532   SUBROUTINE zdf_tke_jki( kt )          ! Empty routine
533      WRITE(*,*) 'zdf_tke_jki: You should not have seen this print! error?', kt
534   END SUBROUTINE zdf_tke_jki
535#endif
536
537   !!======================================================================
538END MODULE zdftke_jki
Note: See TracBrowser for help on using the repository browser.