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

source: trunk/NEMO/OPA_SRC/ZDF/zdftke_atsk.h90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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