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

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

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