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

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

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