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

Last change on this file since 1264 was 1264, checked in by cetlod, 15 years ago

clean TOP model routines to avoid warning when compiling, see ticket:303
  • 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
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_trc, 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
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_trc, ONLY :   zwx => ua,  &  ! use ua, va as
411                            zwy => va      ! workspace arrays
412
413      INTEGER ::   ji, jj, jk, jn                ! dummy loop indices
414      INTEGER ::   iku, ikv                     
415      REAL(wp) ::   ztavg                        ! temporary scalars
416      REAL(wp) ::   zcoef0, zcoef3               !    "         "
417      REAL(wp) ::   zcoef4                       !    "         "
418      REAL(wp) ::   zbtr, zmku, zmkv             !    "         "
419#if defined key_trcldf_eiv                       
420      REAL(wp) ::   zcoeg3, z_hdivn_z            !    "         "
421      REAL(wp) ::   zuwki, zvwki                 !    "         "
422      REAL(wp) ::   zuwk, zvwk                   !    "         "
423#endif
424      REAL(wp) ::   ztav
425      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz  ! temporary workspace arrays
426      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwt
427      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztfw
428      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
429      !!---------------------------------------------------------------------
430
431
432      IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
433
434      !                                                          ! ===========
435      DO jn = 1, jptra                                           ! tracer loop
436         !                                                       ! ===========
437
438         ! 0. Local constant initialization
439         ! --------------------------------
440         zwx (1,:,:) = 0.e0    ;     zwx (jpi,:,:) = 0.e0
441         zwy (1,:,:) = 0.e0    ;     zwy (jpi,:,:) = 0.e0
442         zwz (1,:,:) = 0.e0    ;     zwz (jpi,:,:) = 0.e0
443         zwt (1,:,:) = 0.e0    ;     zwt (jpi,:,:) = 0.e0
444         ztfw(1,:,:) = 0.e0    ;     ztfw(jpi,:,:) = 0.e0
445
446         IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)          ! save trends
447
448         ztavg = 0.e0
449
450         ! I. Vertical trends associated with lateral mixing
451         ! -------------------------------------------------
452         !    (excluding the vertical flux proportional to dk[t] )
453
454         ! I.1 horizontal tracer gradient
455         ! ------------------------------
456
457         DO jk = 1, jpkm1
458            DO jj = 1, jpjm1
459               DO ji = 1, fs_jpim1   ! vector opt.
460                  ! i-gradient of passive tracer at ji
461                  zwx (ji,jj,jk) = ( trb(ji+1,jj,jk,jn)-trb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
462                  ! j-gradient of passive tracer at jj
463                  zwy (ji,jj,jk) = ( trb(ji,jj+1,jk,jn)-trb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
464               END DO
465            END DO
466         END DO
467         IF( ln_zps ) THEN
468            ! partial steps correction at the bottom ocean level
469            DO jj = 1, jpjm1
470               DO ji = 1, fs_jpim1   ! vector opt.
471                  ! last ocean level
472                  iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
473                  ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
474                  ! i-gradient of passive tracer
475                  zwx (ji,jj,iku) = gtru(ji,jj,jn)
476                  ! j-gradient of passive tracer
477                  zwy (ji,jj,ikv) = gtrv(ji,jj,jn) 
478               END DO
479            END DO
480         ENDIF
481
482         ! I.2 Vertical fluxes
483         ! -------------------
484
485         ! Surface and bottom vertical fluxes set to zero
486         ztfw(:,:, 1 ) = 0.e0
487         ztfw(:,:,jpk) = 0.e0
488
489         ! interior (2=<jk=<jpk-1)
490         DO jk = 2, jpkm1
491            DO jj = 2, jpjm1
492               DO ji = fs_2, fs_jpim1   ! vector opt.
493                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
494
495                  zmku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
496                     &           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
497
498                  zmkv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
499                     &           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
500
501                  zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
502                  zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
503
504                  ztfw(ji,jj,jk) = zcoef3 * (   zwx(ji  ,jj  ,jk-1) + zwx(ji-1,jj  ,jk)      &
505                     &                        + zwx(ji-1,jj  ,jk-1) + zwx(ji  ,jj  ,jk)  )   &
506                     &           + zcoef4 * (   zwy(ji  ,jj  ,jk-1) + zwy(ji  ,jj-1,jk)      &
507                     &                        + zwy(ji  ,jj-1,jk-1) + zwy(ji  ,jj  ,jk)  )
508               END DO
509            END DO
510         END DO
511
512#if defined key_trcldf_eiv
513         !                              ! ---------------------------------------!
514         !                              ! Eddy induced vertical advective fluxes !
515         !                              ! ---------------------------------------!
516         zwx(:,:, 1 ) = 0.e0
517         zwx(:,:,jpk) = 0.e0
518
519         DO jk = 2, jpkm1
520            DO jj = 2, jpjm1
521               DO ji = fs_2, fs_jpim1   ! vector opt.
522#   if defined key_traldf_c2d || defined key_traldf_c3d  || defined key_off_degrad
523                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
524                     &  * fsaeitru(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
525                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
526                     &  * fsaeitru(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
527                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
528                     &  * fsaeitrv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
529                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
530                     &  * fsaeitrv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
531
532                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
533#   else
534                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
535                     &  * e2u(ji-1,jj) * umask(ji-1,jj,jk)
536                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
537                     &  * e2u(ji  ,jj) * umask(ji  ,jj,jk)
538                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
539                     &  * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
540                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
541                     &  * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
542
543                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk)   &
544                     &            * ( zuwk - zuwki + zvwk - zvwki )
545#   endif
546                  zwx(ji,jj,jk) = + zcoeg3 * ( trb(ji,jj,jk,jn) + trb(ji,jj,jk-1,jn) )
547
548                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + zwx(ji,jj,jk)
549#   if defined key_diaeiv
550                  w_trc_eiv(ji,jj,jk) = -2. * zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
551#   endif
552               END DO
553            END DO
554         END DO
555#endif
556
557         ! I.3 Divergence of vertical fluxes added to the general tracer trend
558         ! -------------------------------------------------------------------
559
560         DO jk = 1, jpkm1
561            DO jj = 2, jpjm1
562               DO ji = fs_2, fs_jpim1   ! vector opt.
563                  zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
564                  ztav = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
565                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztav
566#if defined key_trc_diatrd
567#   if defined key_trcldf_eiv
568                  ztavg = ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * zbtr
569                  !  WARNING trtrd(ji,jj,jk,7) used for vertical gent velocity trend  not for damping !!!
570                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),7) = ztavg
571#   endif
572                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztav - ztavg
573#endif
574
575               END DO
576            END DO
577         END DO
578
579         ! II. Save the trends for diagnostics
580         ! -----------------------------------
581         IF( l_trdtrc )   THEN
582#if defined key_trcldf_eiv
583
584            ! II.1) Compute the eiv VERTICAL trend
585            DO jk = 1, jpkm1
586               DO jj = 2, jpjm1
587                  DO ji = fs_2, fs_jpim1   ! vector opt.
588                     
589                     !-- Compute the eiv vertical divergence : 1/e3t ( dk[w_eiv] )
590                     !   N.B. This is only possible if key_diaeiv is switched on.
591                     !     Else, the vertical eiv is not diagnosed,
592                     !     so we can only store the flux form trend d_z ( T * w_eiv )
593                     !     instead of w_eiv * d_z( T ). Then, ONLY THE SUM of zonal,
594                     !     meridional, and vertical trends are valid.
595#   if defined key_diaeiv
596                     z_hdivn_z = ( 1. / fse3t(ji,jj,jk) ) * ( w_trc_eiv(ji,jj,jk) - w_trc_eiv(ji,jj,jk+1) )
597#   else
598                     z_hdivn_z = 0.e0
599#   endif
600                     zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
601                     ztrcavg(ji,jj,jk,jn) = ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * zbtr &
602                           &                - trn(ji,jj,jk,jn) * z_hdivn_z
603                  END DO
604               END DO
605            END DO
606
607            ! II.2) save the trends for diagnostic
608            !       N.B. The other part of the computed trend is stored below for later
609            !         output (see trc_zdf_zdf) 
610            IF (luttrd(jn)) CALL trd_mod_trc( ztrcavg(:,:,:,jn), jn, jptrc_trd_zei, kt )
611
612#endif
613            !-- Retain only the vertical diff. trends due to the extra diagonal
614            !   part of the rotated tensor (i.e. remove vert. eiv from the trend)
615            !   N.B. ztrcavg is recycled for this purpose
616            ztrcavg(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:) - ztrcavg(:,:,:,jn)
617
618          END IF
619
620         !                                                       ! ===========
621      END DO                                                     ! tracer loop
622      !                                                          ! ===========
623
624      IF( l_trdtrc ) DEALLOCATE( ztrtrd )
625
626   END SUBROUTINE trc_zdf_iso
627
628#else
629   !!----------------------------------------------------------------------
630   !!   Dummy module :             NO rotation of the lateral mixing tensor
631   !!----------------------------------------------------------------------
632CONTAINS
633   SUBROUTINE trc_zdf_iso_vopt( kt )              ! empty routine
634      WRITE(*,*) 'trc_zdf_iso_vopt: You should not have seen this print! error?', kt
635   END SUBROUTINE trc_zdf_iso_vopt
636#endif
637
638   !!==============================================================================
639END MODULE trczdf_iso_vopt
Note: See TracBrowser for help on using the repository browser.