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

Last change on this file since 1271 was 1271, checked in by rblod, 15 years ago

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