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

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

nemo_v1_update_071:RB: add iom for restart and reorganization of restart

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