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

Last change on this file since 1175 was 1175, checked in by cetlod, 16 years ago

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