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.
trczdf_iso_vopt.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trczdf_iso_vopt.F90 @ 340

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

nemo_v1_update_022 : CE + RB + CT : add print control possibility
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 24.6 KB
Line 
1MODULE trczdf_iso_vopt
2   !!==============================================================================
3   !!                 ***  MODULE  trczdf_iso_vopt  ***
4   !! Ocean passive tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6#if defined key_passivetrc && ( defined key_ldfslp   ||   defined key_esopa )
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
9   !!----------------------------------------------------------------------
10   !!   trc_zdf_iso_vopt : Update the tracer trend with the vertical part of
11   !!                  the isopycnal or geopotential s-coord. operator and
12   !!                  the vertical diffusion. vector optimization, use
13   !!                  k-j-i loops.
14   !!   trc_zdf_iso  :
15   !!   trc_zdf_zdf  :
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE oce_trc         ! ocean dynamics and tracers variables
19   USE trc             ! ocean passive tracers variables
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE trctrp_lec      ! passive tracers transport
22   USE prtctl_trc          ! Print control for debbuging
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC trc_zdf_iso_vopt   !  routine called by step.F90
29
30   !! * Module variables
31   REAL(wp), DIMENSION(jpk) ::  &
32      rdttrc                          ! vertical profile of 2 x time-step
33
34   !! * Substitutions
35#  include "passivetrc_substitute.h90"
36   !!----------------------------------------------------------------------
37
38CONTAINS
39   
40   SUBROUTINE trc_zdf_iso_vopt( kt )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE trc_zdf_iso_vopt  ***
43      !!
44      !! ** Purpose :
45      !! ** Method  :
46      !! ** Action  :
47      !!
48      !! History :
49      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
50      !!   9.0  !  04-03  (C. Ethe)   adapted for passive tracers
51      !!---------------------------------------------------------------------
52      !! * Arguments
53      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
54      CHARACTER (len=22) :: charout
55      !!---------------------------------------------------------------------
56      !!  OPA 8.5, LODYC-IPSL (2002)
57      !!---------------------------------------------------------------------
58
59      IF( kt == nittrc000 ) THEN
60         IF(lwp)WRITE(numout,*)
61         IF(lwp)WRITE(numout,*) 'trc_zdf_iso_vopt : vertical mixing computation'
62         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~~~~~   is  iso-neutral diffusion : implicit vertical time stepping'
63#if defined key_trcldf_eiv && defined key_diaeiv 
64         w_trc_eiv(:,:,:) = 0.e0
65#endif
66      ENDIF
67
68
69      ! I. vertical extra-diagonal part of the rotated tensor
70      ! -----------------------------------------------------
71
72      CALL trc_zdf_iso
73
74      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
75         WRITE(charout, FMT="('zdf - 1')")
76         CALL prt_ctl_trc_info(charout)
77         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
78      ENDIF
79
80      ! II. vertical diffusion (including the vertical diagonal part of the rotated tensor)
81      ! ----------------------
82
83      CALL trc_zdf_zdf( kt )
84
85      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
86         WRITE(charout, FMT="('zdf - 2')")
87         CALL prt_ctl_trc_info(charout)
88         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
89      ENDIF
90
91   END SUBROUTINE trc_zdf_iso_vopt
92
93
94   SUBROUTINE trc_zdf_zdf( kt )
95      !!----------------------------------------------------------------------
96      !!                 ***  ROUTINE trc_zdf_zdf  ***
97      !!                   
98      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
99      !!     including the vertical component of lateral mixing (only for 2nd
100      !!     order operator, for fourth order it is already computed and add
101      !!     to the general trend in traldf.F) and add it to the general trend
102      !!     of the tracer equations.
103      !!
104      !! ** Method  :   The vertical component of the lateral diffusive trends
105      !!      is provided by a 2nd order operator rotated along neural or geo-
106      !!      potential surfaces to which an eddy induced advection can be
107      !!      added. It is computed using before fields (forward in time) and
108      !!      isopycnal or geopotential slopes computed in routine ldfslp.
109      !!
110      !!      Second part: vertical trend associated with the vertical physics
111      !!      ===========  (including the vertical flux proportional to dk[t]
112      !!                  associated with the lateral mixing, through the
113      !!                  update of avt)
114      !!      The vertical diffusion of tracers  is given by:
115      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
116      !!      It is computed using a backward time scheme (t=tra).
117      !!      Surface and bottom boundary conditions: no diffusive flux on
118      !!      both tracers (bottom, applied through the masked field avt).
119      !!      Add this trend to the general trend tra :
120      !!         tra = tra + dz( avt dz(t) )
121      !!         (tra = tra + dz( avs dz(t) ) if lk_trc_zdfddm=T )
122      !!
123      !!      Third part: recover avt resulting from the vertical physics
124      !!      ==========  alone, for further diagnostics (for example to
125      !!                  compute the turbocline depth in diamld).
126      !!         avt = zavt
127      !!         (avs = zavs if lk_trc_zdfddm=T )
128      !!
129      !!      'key_trdtra' defined: trend saved for futher diagnostics.
130      !!
131      !!      macro-tasked on vertical slab (jj-loop)
132      !!
133      !! ** Action  : - Update tra with before vertical diffusion trend
134      !!              - Save the trend in trtrd  ('key_trc_diatrd')
135      !!
136      !! History :
137      !!   6.0  !  90-10  (B. Blanke)  Original code
138      !!   7.0  !  91-11 (G. Madec)
139      !!        !  92-06 (M. Imbard) correction on tracer trend loops
140      !!        !  96-01 (G. Madec) statement function for e3
141      !!        !  97-05 (G. Madec) vertical component of isopycnal
142      !!        !  97-07 (G. Madec) geopotential diffusion in s-coord
143      !!        !  98-03  (L. Bopp MA Foujols) passive tracer generalisation
144      !!        !  00-05  (MA Foujols) add lbc for tracer trends
145      !!        !  00-06  (O Aumont)  correct isopycnal scheme suppress
146      !!        !                     avt multiple correction
147      !!        !  00-08  (G. Madec)  double diffusive mixing
148      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
149      !!   9.0  !  04-03  (C. Ethe )  adapted for passive tracers
150      !!---------------------------------------------------------------------
151      !! * Modules used
152      USE oce_trc, ONLY :   zwd   => ua,  &  ! ua, va used as
153                            zws   => va      ! workspace
154      !! * Arguments
155      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
156
157      !! * Local declarations
158      INTEGER ::   ji, jj, jk,jn                ! dummy loop indices
159      REAL(wp) ::   &
160         zavi, zrhs                          ! temporary scalars
161      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
162         zwi, zwt, zavsi                     ! temporary workspace arrays
163      REAL(wp) ::    ztra              !temporary scalars
164      !!---------------------------------------------------------------------
165      !!  OPA 8.5, LODYC-IPSL (2002)
166      !!---------------------------------------------------------------------
167
168
169      ! I. Local constant initialization
170      ! --------------------------------
171      ! ... time step = 2 rdttra ex
172      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
173         ! time step = 2 rdttra with Arakawa or TVD advection scheme
174         IF( neuler == 0 .AND. kt == nittrc000 ) THEN
175            rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)             ! restarting with Euler time stepping
176         ELSEIF( kt <= nittrc000 + 1 ) THEN
177            rdttrc(:) = 2. * rdttra(:) * FLOAT(ndttrc)         ! leapfrog
178         ENDIF
179      ELSE
180         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)     
181      ENDIF
182
183      DO jn = 1, jptra
184         
185         zwd( 1 ,:,:)=0.e0     ;     zwd(jpi,:,:)=0.e0
186         zws( 1 ,:,:)=0.e0     ;     zws(jpi,:,:)=0.e0
187         zwi( 1 ,:,:)=0.e0     ;     zwi(jpi,:,:)=0.e0
188
189         zwt( 1 ,:,:)=0.e0     ;     zwt(jpi,:,:)=0.e0
190         zwt(  :,:,1)=0.e0     ;     zwt(:,:,jpk)= 0.e0
191         zavsi( 1 ,:,:)=0.e0   ;     zavsi(jpi,:,:)=0.e0 
192         zavsi(  :,:,1)=0.e0   ;     zavsi(:,:,jpk)=0.e0
193
194
195         ! II. Vertical trend associated with the vertical physics
196         ! =======================================================
197         !     (including the vertical flux proportional to dk[t] associated
198         !      with the lateral mixing, through the avt update)
199         !     dk[ avt dk[ (t,s) ] ] diffusive trends
200
201
202         ! II.0 Matrix construction
203         ! ------------------------       
204         ! update and save of avt (and avs if double diffusive mixing)
205         DO jk = 2, jpkm1
206            DO jj = 2, jpjm1
207               DO ji = fs_2, fs_jpim1   ! vector opt.
208                  zavi = fsahtw(ji,jj,jk) * (                 &   ! vertical mixing coef. due to lateral mixing
209                     &                           wslpi(ji,jj,jk) * wslpi(ji,jj,jk)      &
210                     &                         + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
211                  zavsi(ji,jj,jk) = fstravs(ji,jj,jk) + zavi        ! dd mixing: zavsi = total vertical mixing coef. on tracer
212
213               END DO
214            END DO
215         END DO
216
217
218         ! II.2 Vertical diffusion on tracer
219         ! ---------------------------========
220
221         ! Rebuild the Matrix as avt /= avs
222
223         ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
224         DO jk = 1, jpkm1
225            DO jj = 2, jpjm1
226               DO ji = fs_2, fs_jpim1   ! vector opt.
227                  zwi(ji,jj,jk) = - rdttrc(jk) * zavsi(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
228                  zws(ji,jj,jk) = - rdttrc(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
229                  zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
230               END DO
231            END DO
232         END DO
233
234         ! Surface boudary conditions
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1   ! vector opt.
237               zwi(ji,jj,1) = 0.e0
238               zwd(ji,jj,1) = 1. - zws(ji,jj,1)
239            END DO
240         END DO
241
242
243         !! Matrix inversion from the first level
244         !!----------------------------------------------------------------------
245         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
246         !
247         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
248         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
249         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
250         !        (        ...               )( ...  ) ( ...  )
251         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
252         !
253         !   m is decomposed in the product of an upper and lower triangular
254         !   matrix
255         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
256         !   The second member is in 2d array zwy
257         !   The solution is in 2d array zwx
258         !   The 3d arry zwt is a work space array
259         !   zwy is used and then used as a work space array : its value is modified!
260
261         ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
262         DO jj = 2, jpjm1
263            DO ji = fs_2, fs_jpim1
264               zwt(ji,jj,1) = zwd(ji,jj,1)
265            END DO
266         END DO
267         DO jk = 2, jpkm1
268            DO jj = 2, jpjm1
269               DO ji = fs_2, fs_jpim1
270                  zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
271               END DO
272            END DO
273         END DO
274
275         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
276         DO jj = 2, jpjm1
277            DO ji = fs_2, fs_jpim1
278               tra(ji,jj,1,jn) = trb(ji,jj,1,jn) + rdttrc(1) * tra(ji,jj,1,jn)
279            END DO
280         END DO
281         DO jk = 2, jpkm1
282            DO jj = 2, jpjm1
283               DO ji = fs_2, fs_jpim1
284                  zrhs = trb(ji,jj,jk,jn) + rdttrc(jk) * tra(ji,jj,jk,jn)   ! zrhs=right hand side
285                  tra(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * tra(ji,jj,jk-1,jn)
286               END DO
287            END DO
288         END DO
289
290         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
291         ! Save the masked passive tracer after in tra
292         ! (c a u t i o n: passive tracer not its trend, Leap-frog scheme done it will not be done in tranxt)
293         DO jj = 2, jpjm1
294            DO ji = fs_2, fs_jpim1
295               tra(ji,jj,jpkm1,jn) = tra(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
296            END DO
297         END DO
298         DO jk = jpk-2, 1, -1
299            DO jj = 2, jpjm1
300               DO ji = fs_2, fs_jpim1
301                  tra(ji,jj,jk,jn) = ( tra(ji,jj,jk,jn) - zws(ji,jj,jk) * tra(ji,jj,jk+1,jn) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
302               END DO
303            END DO
304         END DO
305
306#if defined key_trc_diatrd
307         ! Compute and save the vertical diffusive passive tracer trends
308#  if defined key_trc_ldfiso 
309         DO jk = 1, jpkm1
310            DO jj = 2, jpjm1
311               DO ji = fs_2, fs_jpim1   ! vector opt.
312                  ztra = ( tra(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
313                  trtrd(ji,jj,jk,jn,6) = ztra - tra(ji,jj,jk,jn) + trtrd(ji,jj,jk,jn,6)
314               END DO
315            END DO
316         END DO
317#  else
318         DO jk = 1, jpkm1
319            DO jj = 2, jpjm1
320               DO ji = fs_2, fs_jpim1   ! vector opt.
321                  ztra = ( tra(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
322                  trtrd(ji,jj,jk,jn,6) = ztra - tra(ji,jj,jk,jn)
323               END DO
324            END DO
325         END DO
326#  endif
327#endif
328
329      END DO
330
331   END SUBROUTINE trc_zdf_zdf
332
333
334   SUBROUTINE trc_zdf_iso
335      !!----------------------------------------------------------------------
336      !!                  ***  ROUTINE trc_zdf_iso  ***
337      !!
338      !! ** Purpose :
339      !!     Compute the trend due to the vertical tracer diffusion inclu-
340      !!     ding the vertical component of lateral mixing (only for second
341      !!     order operator, for fourth order it is already computed and
342      !!     add to the general trend in traldf.F) and add it to the general
343      !!     trend of the tracer equations.
344      !!
345      !! ** Method :
346      !!         The vertical component of the lateral diffusive trends is
347      !!      provided by a 2nd order operator rotated along neural or geopo-
348      !!      tential surfaces to which an eddy induced advection can be added
349      !!      It is computed using before fields (forward in time) and isopyc-
350      !!      nal or geopotential slopes computed in routine ldfslp.
351      !!
352      !!      First part: vertical trends associated with the lateral mixing
353      !!      ==========  (excluding the vertical flux proportional to dk[t] )
354      !!      vertical fluxes associated with the rotated lateral mixing:
355      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(trb)) ]
356      !!                     + e1t*wslpj dj[ mj(mk(trb)) ]  }
357      !!      save avt coef. resulting from vertical physics alone in zavt:
358      !!         zavt = avt
359      !!      update and save in zavt the vertical eddy viscosity coefficient:
360      !!         avt = avt + wslpi^2+wslj^2
361      !!      add vertical Eddy Induced advective fluxes (lk_traldf_eiv=T):
362      !!         zftw = zftw + { di[aht e2u mi(wslpi)]
363      !!                    +dj[aht e1v mj(wslpj)] } mk(trb)
364      !!      take the horizontal divergence of the fluxes:
365      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
366      !!      Add this trend to the general trend tra :
367      !!         tra = tra + difft
368      !!
369      !! ** Action :
370      !!         Update tra arrays with the before vertical diffusion trend
371      !!         Save in trtrd arrays the trends if 'key_trc_diatrd' defined
372      !!
373      !! History :
374      !!   6.0  !  90-10  (B. Blanke)  Original code
375      !!   7.0  !  91-11  (G. Madec)
376      !!        !  92-06  (M. Imbard) correction on tracer trend loops
377      !!        !  96-01  (G. Madec) statement function for e3
378      !!        !  97-05  (G. Madec) vertical component of isopycnal
379      !!        !  97-07  (G. Madec) geopotential diffusion in s-coord
380      !!        !  98-03  (L. Bopp MA Foujols) passive tracer generalisation
381      !!        !  00-05  (MA Foujols) add lbc for tracer trends
382      !!        !  00-06  (O Aumont)  correct isopycnal scheme suppress
383      !!        !                     avt multiple correction
384      !!        !  00-08  (G. Madec)  double diffusive mixing
385      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
386      !!   9.0  !  04-03  (C. Ethe )  adapted for passive tracers
387      !!---------------------------------------------------------------------
388      !! * Modules used
389      USE oce_trc, ONLY :   zwx => ua,  &  ! use ua, va as
390                            zwy => va      ! workspace arrays
391
392      !! * Local declarations
393      INTEGER ::   ji, jj, jk,jn       ! dummy loop indices
394#if defined key_partial_steps
395      INTEGER ::   iku, ikv
396#endif
397      REAL(wp) ::   &
398         ztavg,                  &  ! temporary scalars
399         zcoef0, zcoef3,         &  !    "         "
400         zcoef4,                 &  !    "         "
401         zbtr, zmku, zmkv,       &  !    "         "
402#if defined key_trcldf_eiv
403         zcoeg3,                 &  !    "         "
404         zuwki, zvwki,           &  !    "         "
405         zuwk, zvwk,             &  !    "         "
406#endif
407         ztav
408      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
409         zwz, zwt, ztfw             ! temporary workspace arrays
410
411      !!---------------------------------------------------------------------
412      !!  OPA 8.5, LODYC-IPSL (2002)
413      !!---------------------------------------------------------------------
414
415      DO jn = 1, jptra
416
417         ! 0. Local constant initialization
418         ! --------------------------------
419         ztavg = 0.e0
420
421         zwx( 1 ,:,:)=0.e0     ;     zwx(jpi,:,:)=0.e0
422         zwy( 1 ,:,:)=0.e0     ;     zwy(jpi,:,:)=0.e0
423         zwz( 1 ,:,:)=0.e0     ;     zwz(jpi,:,:)=0.e0
424         zwt( 1 ,:,:)=0.e0     ;     zwt(jpi,:,:)=0.e0
425         ztfw( 1 ,:,:)=0.e0    ;     ztfw(jpi,:,:)=0.e0
426
427         ! I. Vertical trends associated with lateral mixing
428         ! -------------------------------------------------
429         !    (excluding the vertical flux proportional to dk[t] )
430
431
432         ! I.1 horizontal tracer gradient
433         ! ------------------------------
434
435         DO jk = 1, jpkm1
436            DO jj = 1, jpjm1
437               DO ji = 1, fs_jpim1   ! vector opt.
438                  ! i-gradient of passive tracer at ji
439                  zwx (ji,jj,jk) = ( trb(ji+1,jj,jk,jn)-trb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
440                  ! j-gradient of passive tracer at jj
441                  zwy (ji,jj,jk) = ( trb(ji,jj+1,jk,jn)-trb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
442               END DO
443            END DO
444         END DO
445#  if defined key_partial_steps
446         ! partial steps correction at the bottom ocean level
447         DO jj = 1, jpjm1
448            DO ji = 1, fs_jpim1   ! vector opt.
449               ! last ocean level
450               iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
451               ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
452               ! i-gradient of passive tracer
453               zwx (ji,jj,iku) = gtru(ji,jj,jn)
454               ! j-gradient of passive tracer
455               zwy (ji,jj,ikv) = gtrv(ji,jj,jn) 
456            END DO
457         END DO
458#endif
459
460
461         ! I.2 Vertical fluxes
462         ! -------------------
463
464         ! Surface and bottom vertical fluxes set to zero
465         ztfw(:,:, 1 ) = 0.e0
466         ztfw(:,:,jpk) = 0.e0
467
468         ! interior (2=<jk=<jpk-1)
469         DO jk = 2, jpkm1
470            DO jj = 2, jpjm1
471               DO ji = fs_2, fs_jpim1   ! vector opt.
472                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
473
474                  zmku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
475                     &           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
476
477                  zmkv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
478                     &           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
479
480                  zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
481                  zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
482
483                  ztfw(ji,jj,jk) = zcoef3 * (   zwx(ji  ,jj  ,jk-1) + zwx(ji-1,jj  ,jk)      &
484                     &                        + zwx(ji-1,jj  ,jk-1) + zwx(ji  ,jj  ,jk)  )   &
485                     &           + zcoef4 * (   zwy(ji  ,jj  ,jk-1) + zwy(ji  ,jj-1,jk)      &
486                     &                        + zwy(ji  ,jj-1,jk-1) + zwy(ji  ,jj  ,jk)  )
487               END DO
488            END DO
489         END DO
490
491#if defined key_trcldf_eiv
492         !                              ! ---------------------------------------!
493         !                              ! Eddy induced vertical advective fluxes !
494         !                              ! ---------------------------------------!
495         zwx(:,:, 1 ) = 0.e0
496         zwx(:,:,jpk) = 0.e0
497
498         DO jk = 2, jpkm1
499            DO jj = 2, jpjm1
500               DO ji = fs_2, fs_jpim1   ! vector opt.
501#   if defined key_traldf_c2d || defined key_traldf_c3d
502                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
503                     &  * fsaeitru(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
504                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
505                     &  * fsaeitru(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
506                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
507                     &  * fsaeitrv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
508                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
509                     &  * fsaeitrv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
510
511                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
512#   else
513                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
514                     &  * e2u(ji-1,jj) * umask(ji-1,jj,jk)
515                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
516                     &  * e2u(ji  ,jj) * umask(ji  ,jj,jk)
517                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
518                     &  * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
519                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
520                     &  * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
521
522                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk)   &
523                     &            * ( zuwk - zuwki + zvwk - zvwki )
524#   endif
525                  zwx(ji,jj,jk) = + zcoeg3 * ( trb(ji,jj,jk,jn) + trb(ji,jj,jk-1,jn) )
526
527                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + zwx(ji,jj,jk)
528#   if defined key_diaeiv
529                  w_trc_eiv(ji,jj,jk) = -2. * zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
530#   endif
531               END DO
532            END DO
533         END DO
534#endif
535
536         ! I.5 Divergence of vertical fluxes added to the general tracer trend
537         ! -------------------------------------------------------------------
538
539         DO jk = 1, jpkm1
540            DO jj = 2, jpjm1
541               DO ji = fs_2, fs_jpim1   ! vector opt.
542                  zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
543                  ztav = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
544                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztav
545#if defined key_trc_diatrd
546#   if defined key_trcldf_eiv
547                  ztavg = ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * zbtr
548                  !  WARNING trtrd(ji,jj,jk,7) used for vertical gent velocity trend  not for damping !!!
549                  trtrd(ji,jj,jk,jn,7) = ztavg
550#   endif
551                  trtrd(ji,jj,jk,jn,6) = ztav - ztavg
552#endif
553               END DO
554            END DO
555         END DO
556
557      END DO
558
559   END SUBROUTINE trc_zdf_iso
560
561#else
562   !!----------------------------------------------------------------------
563   !!   Dummy module :             NO rotation of the lateral mixing tensor
564   !!----------------------------------------------------------------------
565CONTAINS
566   SUBROUTINE trc_zdf_iso_vopt( kt )              ! empty routine
567      WRITE(*,*) 'trc_zdf_iso_vopt: You should not have seen this print! error?', kt
568   END SUBROUTINE trc_zdf_iso_vopt
569#endif
570
571   !!==============================================================================
572END MODULE trczdf_iso_vopt
Note: See TracBrowser for help on using the repository browser.