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 tags/nemo_v1_06/NEMO/TOP_SRC/TRP – NEMO

source: tags/nemo_v1_06/NEMO/TOP_SRC/TRP/trczdf_iso_vopt.F90 @ 2915

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

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