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

Last change on this file since 106 was 106, checked in by opalod, 20 years ago

CT : UPDATE067 : Add control indices nictl, njctl used in SUM function output to compare mono versus multi procs runs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.6 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      !!  OPA8.5, LODYC-IPSL (2002)
84      !!--------------------------------------------------------------------
85
86
87      ! 0. Initialization
88      !    --------------
89      IF( kt == nit000  )   CALL zdf_tke_init
90
91      ! Local constant
92      zmlmin = 1.e-8
93      zbbrau =  .5 * ebb / rau0
94      zfact1 = -.5 * rdt * efave
95      zfact2 = 1.5 * rdt * ediss
96      zfact3 = 0.5 * rdt * ediss
97
98
99      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
100      ! I.  Mixing length
101      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
102
103      !                                                ! ===============
104      DO jj = 2, jpjm1                                 !  Vertical slab
105         !                                             ! ===============
106         ! Buoyancy length scale: l=sqrt(2*e/n**2)
107         ! ---------------------
108         zmxlm(:,jj, 1 ) = zmlmin   ! surface set to the minimum value
109         zmxlm(:,jj,jpk) = zmlmin   ! bottom  set to the minimum value
110!CDIR NOVERRCHK
111         DO jk = 2, jpkm1
112!CDIR NOVERRCHK
113            DO ji = 2, jpim1
114               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
115               zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  )
116            END DO
117         END DO
118
119
120         ! Physical limits for the mixing length
121         ! -------------------------------------
122         zmxld(:,jj, 1 ) = zmlmin   ! surface set to the minimum value
123         zmxld(:,jj,jpk) = zmlmin   ! bottom  set to the minimum value
124
125         SELECT CASE ( nmxl )
126
127         CASE ( 0 )           ! bounded by the distance to surface and bottom
128
129            DO jk = 2, jpkm1
130               DO ji = 2, jpim1
131                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
132                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
133                  zmxlm(ji,jj,jk) = zemxl
134                  zmxld(ji,jj,jk) = zemxl
135               END DO
136            END DO
137
138         CASE ( 1 )           ! bounded by the vertical scale factor
139
140            DO jk = 2, jpkm1
141               DO ji = 2, jpim1
142                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
143                  zmxlm(ji,jj,jk) = zemxl
144                  zmxld(ji,jj,jk) = zemxl
145               END DO
146            END DO
147
148         CASE ( 2 )           ! |dk[xml]| bounded by e3t :
149
150            DO jk = 2, jpk           ! from the surface to the bottom :
151               DO ji = 2, jpim1
152                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
153               END DO
154            END DO
155            DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
156               DO ji = 2, jpim1
157                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
158                  zmxlm(ji,jj,jk) = zemxl
159                  zmxld(ji,jj,jk) = zemxl
160               END DO
161            END DO
162
163         CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
164
165            DO jk = 2, jpk           ! from the surface to the bottom : lup
166               DO ji = 2, jpim1
167                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
168               END DO
169            END DO
170            DO jk = jpkm1, 1, -1     ! from the bottom to the surface : ldown
171               DO ji = 2, jpim1
172                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
173               END DO
174            END DO
175!CDIR NOVERRCHK
176            DO jk = 1, jpk
177!CDIR NOVERRCHK
178               DO ji = 2, jpim1
179                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
180                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
181                  zmxlm(ji,jj,jk) = zemlm
182                  zmxld(ji,jj,jk) = zemlp
183               END DO
184            END DO
185
186         END SELECT
187
188
189         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
190         ! II  Tubulent kinetic energy time stepping
191         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
192
193
194         ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
195         ! ---------------------------------------------------------------------
196!CDIR NOVERRCHK
197         DO jk = 2, jpkm1
198!CDIR NOVERRCHK
199            DO ji = 2, jpim1
200               zsqen = SQRT( en(ji,jj,jk) )
201               zav   = ediff * zmxlm(ji,jj,jk) * zsqen
202               avt  (ji,jj,jk) = MAX( zav, avtb(jk) ) * tmask(ji,jj,jk)
203               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
204               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
205            END DO
206         END DO
207
208
209         ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
210         ! -------------------------------------------------
211         ! en(1)   = ebb sqrt(taux^2+tauy^2) / rau0  (min value emin0)
212         ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
213!CDIR NOVERRCHK
214         DO ji = 2, jpim1
215            ztx2 = taux(ji-1,jj  ) + taux(ji,jj)
216            zty2 = tauy(ji  ,jj-1) + tauy(ji,jj)
217            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
218            en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
219            zmxlm(ji,jj,1  ) = avmb(1) * tmask(ji,jj,1)
220            zmxlm(ji,jj,jpk) = 0.e0
221         END DO
222
223
224         ! 3. Now Turbulent kinetic energy (output in en)
225         ! -------------------------------
226         ! Resolution of a tridiagonal linear system by a "methode de chasse"
227         ! computation from level 2 to jpkm1  (e(1) already computed and
228         ! e(jpk)=0 ).
229
230         SELECT CASE ( npdl )
231
232         CASE ( 0 )           ! No Prandtl number
233            DO jk = 2, jpkm1
234               DO ji = 2, jpim1
235                  ! zesh2 = eboost * (du/dz)^2 - N^2
236                  zcoef = 0.5 / fse3w(ji,jj,jk)
237                  ! shear
238                  zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &
239                  &                - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
240                  zdkv = zcoef * (   vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
241                  &                - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
242                  ! coefficient (zesh2)
243                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)
244
245                  ! Matrix
246                  zcof = zfact1 * tmask(ji,jj,jk)
247                  ! lower diagonal
248                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
249                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
250                  ! upper diagonal
251                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
252                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
253                  ! diagonal
254                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
255                  ! right hand side in en
256                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
257                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
258               END DO
259            END DO
260
261         CASE ( 1 )           ! Prandtl number
262            DO jk = 2, jpkm1
263               DO ji = 2, jpim1
264                  ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2
265                  zcoef = 0.5 / fse3w(ji,jj,jk)
266                  ! shear
267                  zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) &
268                  &                - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )   )
269                  zdkv = zcoef * (   vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1) &
270                  &                - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )   )
271                  ! square of vertical shear
272                  zsh2 = zdku * zdku + zdkv * zdkv
273                  ! Prandtl number
274                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
275                  zpdl = 1.0
276                  IF( zri >= 0.2 ) zpdl = 0.2 / zri
277                  zpdl = MAX( 0.1, zpdl )
278                  ! coefficient (esh2)
279                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
280
281                  ! Matrix
282                  zcof = zfact1 * tmask(ji,jj,jk)
283                  ! lower diagonal
284                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
285                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
286                  ! upper diagonal
287                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
288                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
289                  ! diagonal
290                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
291                  ! right hand side in en
292                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
293                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
294                  ! save masked Prandlt number in zmxlm array
295                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
296               END DO
297            END DO
298
299         END SELECT
300
301
302         ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
303         !---------------------------------
304
305         ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
306         DO jk = 3, jpkm1
307            DO ji = 2, jpim1
308               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
309            END DO
310         END DO
311
312         ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
313         DO ji = 2, jpim1
314            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
315         END DO
316         DO jk = 3, jpkm1
317            DO ji = 2, jpim1
318               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
319            END DO
320         END DO
321
322         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
323         DO ji = 2, jpim1
324            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
325         END DO
326         DO jk = jpk-2, 2, -1
327            DO ji = 2, jpim1
328               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
329            END DO
330         END DO
331
332         ! Save the result in en and set minimum value of tke : emin
333         DO jk = 2, jpkm1
334            DO ji = 2, jpim1
335               en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
336            END DO
337         END DO
338         !                                             ! ===============
339      END DO                                           !   End of slab
340      !                                                ! ===============
341
342      ! Lateral boundary conditions on ( avt, en )  (sign unchanged)
343      ! --------------------------------=========
344      CALL lbc_lnk( avt, 'W', 1. )   ;   CALL lbc_lnk( en , 'W', 1. )
345
346
347      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
348      ! III.  Before vertical eddy vicosity and diffusivity coefficients
349      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
350
351      !                                                ! ===============
352      DO jk = 2, jpkm1                                 ! Horizontal slab
353         !                                             ! ===============
354         SELECT CASE ( nave )
355           
356         CASE ( 0 )                ! no horizontal average
357
358            ! Vertical eddy viscosity
359
360            DO jj = 2, jpjm1
361               DO ji = fs_2, fs_jpim1   ! vector opt.
362                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
363                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
364                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
365                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
366               END DO
367            END DO
368
369           
370         CASE ( 1 )                ! horizontal average
371
372            !                                                ( 1/2  1/2 )
373            ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1    1   )
374            !                      ( 1/2  1 1/2 )            ( 1/2  1/2 )
375            !           avmv = 1/4 ( 1/2  1 1/2 )     
376           
377!! caution vectopt_memory change the solution (last digit of the solver stat)
378#  if defined key_vectopt_memory
379            DO jj = 2, jpjm1
380               DO ji = fs_2, fs_jpim1   ! vector opt.
381                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
382                  &                 +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
383                  &                      +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
384
385                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
386                  &                 +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
387                  &                      +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
388               END DO
389            END DO
390#  else
391            DO jj = 2, jpjm1
392               DO ji = fs_2, fs_jpim1   ! vector opt.
393                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
394                  &              +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
395                  &                   +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)  &
396                  &       / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
397                  &              +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
398                  &                   +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
399
400                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
401                  &              +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
402                  &                   +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)  &
403                  &      /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
404                  &              +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
405                  &                   +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
406               END DO
407            END DO
408#  endif
409         END SELECT
410         !                                             ! ===============
411      END DO                                           !   End of slab
412      !                                                ! ===============
413
414      ! Lateral boundary conditions (avmu,avmv)  (sign unchanged)
415      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
416
417      !                                                ! ===============
418      DO jk = 2, jpkm1                                 ! Horizontal slab
419         !                                             ! ===============
420         SELECT CASE ( nave )
421           
422         CASE ( 1 )                ! horizontal average
423
424            ! Vertical eddy diffusivity
425            ! ------------------------------
426            !                                (1 2 1)
427            ! horizontal average  avt = 1/16 (2 4 2)
428            !                                (1 2 1)
429!! caution vectopt_memory change the solution (last digit of the solver stat)
430#  if defined key_vectopt_memory
431            DO jj = 2, jpjm1
432               DO ji = fs_2, fs_jpim1   ! vector opt.
433                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)    &
434                  &               + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
435               END DO
436            END DO
437#  else
438            DO jj = 2, jpjm1
439               DO ji = fs_2, fs_jpim1   ! vector opt.
440                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)   &
441                  &               + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
442                  &     / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
443                  &               + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
444               END DO
445            END DO
446#  endif
447         END SELECT
448
449
450         ! multiplied by the Prandtl number (npdl>1)
451         ! ----------------------------------------
452         IF( npdl == 1 ) THEN
453            DO jj = 2, jpjm1
454               DO ji = fs_2, fs_jpim1   ! vector opt.
455                  zpdl = zmxld(ji,jj,jk)
456                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
457               END DO
458            END DO
459         ENDIF
460
461         ! Minimum value on the eddy viscosity
462         ! ----------------------------------------
463         DO jj = 1, jpj
464            DO ji = 1, jpi
465               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
466               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
467            END DO
468         END DO
469         !                                             ! ===============
470      END DO                                           !   End of slab
471      !                                                ! ===============
472
473
474      ! Lateral boundary conditions on avt   (W-point (=T), sign unchanged)
475      ! ------------------------------=====
476      CALL lbc_lnk( avt, 'W', 1. )
477
478      IF(l_ctl) THEN
479         WRITE(numout,*) ' tke  e : ', SUM( en  (1:nictl+1,1:njctl+1,:) ), ' t : ', SUM( avt (1:nictl+1,1:njctl+1,:) )
480         WRITE(numout,*) '      u : ', SUM( avmu(1:nictl+1,1:njctl+1,:) ), ' v : ', SUM( avmv(1:nictl+1,1:njctl+1,:) )
481      ENDIF
482
483
484   END SUBROUTINE zdf_tke
Note: See TracBrowser for help on using the repository browser.