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.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trczdf_iso.F90 @ 724

Last change on this file since 724 was 724, checked in by cetlod, 17 years ago

Update modules for passive tracers transport trends computation, see ticket:13
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.4 KB
Line 
1MODULE trczdf_iso
2   !!==============================================================================
3   !!                    ***  MODULE  trczdf_iso  ***
4   !! Ocean passive tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6#if defined key_passivetrc && ( defined key_ldfslp   ||   defined key_esopa )
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
9   !!----------------------------------------------------------------------
10   !!   trc_zdf_iso  : update the tracer trend with the vertical part of
11   !!                  the isopycnal or geopotential s-coord. operator and
12   !!                  the vertical diffusion
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce_trc          ! ocean dynamics and tracers variables
16   USE trc              ! ocean passive tracers variables
17   USE lbclnk           ! ocean lateral boundary conditions (or mpp link)
18   USE trctrp_lec       ! passive tracers transport
19   USE prtctl_trc          ! Print control for debbuging
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! * Accessibility
25   PUBLIC trc_zdf_iso    ! routine called by step.F90
26
27   !! * Module variable
28   REAL(wp), DIMENSION(jpk) ::   &
29      rdttrc                     ! vertical profile of 2 x tracer time-step
30
31   !! * Substitutions
32#  include "passivetrc_substitute.h90"
33   !!----------------------------------------------------------------------
34   !!   TOP 1.0 , LOCEAN-IPSL (2005)
35   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trczdf_iso.F90,v 1.13 2007/10/12 09:26:30 opalod Exp $
36   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE trc_zdf_iso( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE trc_zdf_iso  ***
44      !!
45      !! ** Purpose :
46      !!     Compute the trend due to the vertical tracer diffusion inclu-
47      !!     ding the vertical component of lateral mixing (only for second
48      !!     order operator, for fourth order it is already computed and
49      !!     add to the general trend in trcldf.F) and add it to the general
50      !!     trend of the tracer equations.
51      !!
52      !! ** Method :
53      !!         The vertical component of the lateral diffusive trends is
54      !!      provided by a 2nd order operator rotated along neural or geopo-
55      !!      tential surfaces to which an eddy induced advection can be added
56      !!      It is computed using before fields (forward in time) and isopyc-
57      !!      nal or geopotential slopes computed in routine ldfslp.
58      !!
59      !!      First part: vertical trends associated with the lateral mixing
60      !!      ==========  (excluding the vertical flux proportional to dk[t] )
61      !!      vertical fluxes associated with the rotated lateral mixing:
62      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(trb)) ]
63      !!                     + e1t*wslpj dj[ mj(mk(trb)) ]  }
64      !!      save avt coef. resulting from vertical physics alone in zavt:
65      !!         zavt = avt
66      !!      update and save in zavt the vertical eddy viscosity coefficient:
67      !!         avt = avt + wslpi^2+wslj^2
68      !!      add vertical Eddy Induced advective fluxes ('lk_trcldf_eiv=T):
69      !!         zftw = zftw + { di[aht e2u mi(wslpi)]
70      !!                    +dj[aht e1v mj(wslpj)] } mk(trb)
71      !!      take the horizontal divergence of the fluxes:
72      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
73      !!      Add this trend to the general trend tra :
74      !!         tra = tra + difft
75      !!
76      !!      Second part: vertical trend associated with the vertical physics
77      !!      ===========  (including the vertical flux proportional to dk[t]
78      !!                  associated with the lateral mixing, through the
79      !!                  update of avt)
80      !!      The vertical diffusion of tracers tra  is given by:
81      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
82      !!      It is computed using a backward time scheme, t=ta.
83      !!      Surface and bottom boundary conditions: no diffusive flux on
84      !!      both tracers (bottom, applied through the masked field avt).
85      !!      Add this trend to the general trend tra  :
86      !!         tra = tra + dz( avt dz(t) )
87      !!         (tra = tra + dz( avs dz(t) ) if lk_trc_zdfddm=T )
88      !!
89      !!      Third part: recover avt resulting from the vertical physics
90      !!      ==========  alone, for further diagnostics (for example to
91      !!                  compute the turbocline depth in diamld).
92      !!         avt = zavt
93      !!         (avs = zavs if lk_trc_zdfddm=T )
94      !!
95      !!      'key_trc_diatrd' defined: trend saved for futher diagnostics.
96      !!
97      !!      macro-tasked on vertical slab (jj-loop)
98      !!
99      !! ** Action :
100      !!         Update tra arrays with the before vertical diffusion trend
101      !!         Save in trtrd arrays the trends if 'key_trc_diatrd' defined
102      !!
103      !! History :
104      !!   7.0  !  91-11  (G. Madec)  Original code
105      !!        !  92-06  (M. Imbard)  correction on tracer trend loops
106      !!        !  96-01  (G. Madec)  statement function for e3
107      !!        !  97-05  (G. Madec)  vertical component of isopycnal
108      !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord
109      !!        !  98-03  (L. Bopp MA Foujols) passive tracer generalisation
110      !!        !  00-05  (MA Foujols) add lbc for tracer trends
111      !!        !  00-06  (O Aumont)  correct isopycnal scheme suppress
112      !!        !                     avt multiple correction
113      !!        !  00-08  (G. Madec)  double diffusive mixing
114      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
115      !!   9.0  !  04-03  (C. Ethe )  adapted for passive tracers
116      !!---------------------------------------------------------------------
117      !! * Modules used
118      USE oce_trc               ,   &
119         zavs => va
120
121      !! * Arguments
122      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
123
124      !! * Local declarations
125      INTEGER ::   ji, jj, jk, jn                 ! dummy loop indices
126      INTEGER ::   ikst, ikenm2, ikstp1       ! temporary integers
127      INTEGER ::   iku, ikv, ikv1             ! temporary integers
128
129      REAL(wp) ::   ztra
130      REAL(wp) ::   &
131         ztavg,                 &  ! ???
132         zcoef0, zcoef3,        &  ! ???
133         zcoef4, zavi,          &  ! ???
134         zbtr, zmku, zmkv,      &  !
135         ztav
136      REAL(wp), DIMENSION(jpi,jpk) ::   &
137         zwd, zws, zwi,         &  ! ???
138         zwx, zwy, zwz, zwt        ! ???
139      REAL(wp), DIMENSION(jpi,jpk) ::   &
140         ztfw, zdit, zdjt, zdj1t
141#if defined key_trcldf_eiv   ||   defined key_esopa
142      REAL(wp), DIMENSION(jpi,jpk) ::   &
143         ztfwg
144
145      REAL(wp) ::         &
146         zcoeg3,          &
147         zuwk, zvwk,      &
148         zuwki, zvwki
149#endif
150      CHARACTER (len=22) :: charout
151      !!---------------------------------------------------------------------
152
153      IF( kt == nittrc000 ) THEN
154         IF(lwp) WRITE(numout,*)
155         IF(lwp) WRITE(numout,*) 'trc_zdf_iso : vertical mixing (including isopycnal component)'
156         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
157#if defined key_trcldf_eiv && defined key_diaeiv
158         w_trc_eiv(:,:,:) = 0.e0
159#endif
160      ENDIF
161
162      ! 0.0  Local constant initialization
163      ! --------------------------------
164      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
165         ! time step = 2 rdttra with Arakawa or TVD advection scheme
166         IF( neuler == 0 .AND. kt == nittrc000 ) THEN
167            rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)             ! restarting with Euler time stepping
168         ELSEIF( kt <= nittrc000 + ndttrc ) THEN
169            rdttrc(:) = 2. * rdttra(:) * FLOAT(ndttrc)         ! leapfrog
170         ENDIF
171      ELSE
172         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)     
173      ENDIF
174
175
176
177      ! 0.1 Save avs in zavs to recover avs in output files
178      !---------------------------------------------------
179      zavs(:,:,:) = fstravs(:,:,:)
180
181
182
183      DO jn = 1, jptra
184
185         ztavg = 0.e0
186
187         !                                                ! ===============
188         DO jj = 2, jpjm1                                 !  Vertical slab
189            !                                             ! ===============
190
191            ! I. vertical trends associated with the lateral mixing
192            ! =====================================================
193            !  (excluding the vertical flux proportional to dk[t]
194
195
196            ! I.1 horizontal tracer gradient
197            ! ------------------------------
198
199            DO jk = 1, jpkm1
200               DO ji = 1, jpim1
201                  ! i-gradient of passive tracer at jj
202                  zdit (ji,jk) = ( trb(ji+1,jj,jk,jn)-trb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
203                  ! j-gradient of passive tracer at jj
204                  zdjt (ji,jk) = ( trb(ji,jj+1,jk,jn)-trb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
205                  ! j-gradient of passive tracer at jj+1
206                  zdj1t(ji,jk) = ( trb(ji,jj,jk,jn)-trb(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
207               END DO
208            END DO
209
210            IF( ln_zps ) THEN
211               ! partial steps correction at the bottom ocean level
212               DO ji = 1, jpim1
213                  ! last ocean level
214                  iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
215                  ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
216                  ikv1 = MIN( mbathy(ji,jj), mbathy(ji  ,jj-1) ) - 1
217                  ! i-gradient of of passive tracer at jj
218                  zdit (ji,iku) = gtru(ji,jj,jn)
219                  ! j-gradient of of passive tracer at jj
220                  zdjt (ji,ikv) = gtrv(ji,jj,jn) 
221                  ! j-gradient of of passive tracer at jj+1
222                  zdj1t(ji,ikv1)= gtrv(ji,jj-1,jn)
223               END DO
224            ENDIF
225
226            ! I.2 Vertical fluxes
227            ! -------------------
228
229            ! Surface and bottom vertical fluxes set to zero
230            ztfw(:, 1 ) = 0.e0
231            ztfw(:,jpk) = 0.e0
232
233#if defined key_trcldf_eiv
234            ztfwg(:, 1 ) = 0.e0
235            ztfwg(:,jpk) = 0.e0
236#endif
237
238            ! interior (2=<jk=<jpk-1)
239            DO jk = 2, jpkm1
240               DO ji = 2, jpim1
241                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
242
243                  zmku = 1./MAX( umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)   &
244                     &          +umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1. )
245
246                  zmkv = 1./MAX( vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)   &
247                     &          +vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1. )
248
249                  zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
250                  zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
251
252                  ztfw(ji,jk) = zcoef3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
253                     &                    +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
254                     &        + zcoef4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
255                     &                    +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )
256
257               END DO
258            END DO
259
260
261            ! I.3  update and save of avt (and avs if double diffusive mixing)
262            ! ---------------------------
263
264            DO jk = 2, jpkm1
265               DO ji = 2, jpim1
266
267                  zavi = fsahtw(ji,jj,jk)*( wslpi(ji,jj,jk)*wslpi(ji,jj,jk)   &
268                     &                     +wslpj(ji,jj,jk)*wslpj(ji,jj,jk) )
269
270                  ! add isopycnal vertical coeff. to avs
271                  fstravs(ji,jj,jk) = fstravs(ji,jj,jk) + zavi
272
273               END DO
274            END DO
275
276#if defined key_trcldf_eiv
277            !                              ! ---------------------------------------!
278            !                              ! Eddy induced vertical advective fluxes !
279            !                              ! ---------------------------------------!
280#if defined key_traldf_c2d || defined key_traldf_c3d || defined key_off_degrad
281            DO jk = 2, jpkm1
282               DO ji = 2, jpim1
283                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
284                     &  * fsaeitru(ji-1,jj,jk) * e2u(ji-1,jj)*umask(ji-1,jj,jk)
285                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
286                     &  * fsaeitru(ji  ,jj,jk) * e2u(ji  ,jj)*umask(ji  ,jj,jk)
287                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
288                     &  * fsaeitrv(ji,jj-1,jk) * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
289                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
290                     &  * fsaeitrv(ji,jj  ,jk) * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
291
292                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
293
294                  ztfwg(ji,jk) = + zcoeg3 * ( trb(ji,jj,jk,jn) + trb(ji,jj,jk-1,jn) ) 
295                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jk)
296
297# if defined key_diaeiv
298                  w_trc_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
299# endif
300               END DO
301            END DO
302
303#else
304            DO jk = 2, jpkm1
305               DO ji = 2, jpim1
306                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
307                     &  * e2u(ji-1,jj)*umask(ji-1,jj,jk)
308                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
309                     &  * e2u(ji  ,jj)*umask(ji  ,jj,jk)
310                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
311                     &  * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
312                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
313                     &  * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
314
315                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeitrw(ji,jj,jk)   &
316                     &            * ( zuwk - zuwki + zvwk - zvwki )
317
318                  ztfwg(ji,jk) = + zcoeg3 * ( trb(ji,jj,jk,jn) + trb(ji,jj,jk-1,jn) )
319                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jk)
320
321# if defined key_diaeiv
322                  w_trc_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
323# endif
324               END DO
325            END DO
326#endif
327
328#endif
329
330
331            ! I.5 Divergence of vertical fluxes added to the general tracer trend
332            ! -------------------------------------------------------------------
333
334            DO jk = 1, jpkm1
335               DO ji = 2, jpim1
336                  zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
337                  ztav = (  ztfw(ji,jk) - ztfw(ji,jk+1)  ) * zbtr
338                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztav
339
340#if defined key_trc_diatrd
341#   if defined key_trcldf_eiv
342                  ztavg = ( ztfwg(ji,jk) - ztfwg(ji,jk+1) ) * zbtr
343                  !  WARNING trtrd(ji,jj,jk,6) used for vertical gent velocity trend
344                  !                           not for damping !!!
345                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),9) = ztavg
346#   endif
347                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztav - ztavg
348#endif
349               END DO
350            END DO
351            !                                             ! ===============
352         END DO                                           !   End of slab
353         !                                                ! ===============
354
355      END DO
356
357      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
358         WRITE(charout, FMT="('zdf - 1')")
359         CALL prt_ctl_trc_info(charout)
360         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
361      ENDIF
362
363      DO jn = 1, jptra
364         !                                                ! ===============
365         DO jj = 2, jpjm1                                 !  Vertical slab
366            !                                             ! ===============
367
368            ! II. Vertical trend associated with the vertical physics
369            ! =======================================================
370            !     (including the vertical flux proportional to dk[t] associated
371            !      with the lateral mixing, through the avt update)
372            !     dk[ avt dk[ (t,s) ] ] diffusive trends
373
374
375            ! Diagonal, inferior, superior
376            ! (including the bottom boundary condition via avs masked)
377            DO jk = 1, jpkm1
378               DO ji = 2, jpim1
379                  zwi(ji,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  )   &
380                     /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
381                  zws(ji,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1)   &
382                     /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
383                  zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
384               END DO
385            END DO
386
387            ! Surface boudary conditions
388            DO ji = 2, jpim1
389               zwi(ji,1) = 0.e0
390               zwd(ji,1) = 1. - zws(ji,1)
391            END DO
392
393            ! Second member construction
394            DO jk = 1, jpkm1
395               DO ji = 2, jpim1
396                  zwy(ji,jk) = trb(ji,jj,jk,jn) + rdttrc(jk) * tra(ji,jj,jk,jn)
397               END DO
398            END DO
399
400            ! Matrix inversion from the first level
401            ikst = 1
402#   include "zdf.matrixsolver.h90"
403#if defined key_trc_diatrd
404            ! Compute and save the vertical diffusive of tracers trends
405#  if defined key_trcldf_iso
406            DO jk = 1, jpkm1
407               DO ji = 2, jpim1
408                  ztra = ( zwx(ji,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
409                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztra - tra(ji,jj,jk,jn) + trtrd(ji,jj,jk,ikeep(jn),6)
410               END DO
411            END DO
412#  else
413            DO jk = 1, jpkm1
414               DO ji = 2, jpim1
415                  ztra = ( zwx(ji,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
416                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztra - tra(ji,jj,jk,jn)
417               END DO
418            END DO
419#  endif
420#endif
421            ! Save the masked passive tracer after in tra
422            ! (c a u t i o n:  tracer not its trend, Leap-frog scheme done
423            !                  it will not be done in tranxt)
424            DO jk = 1, jpkm1
425               DO ji = 2, jpim1
426                  tra(ji,jj,jk,jn) = zwx(ji,jk)  * tmask(ji,jj,jk)
427               END DO
428            END DO
429            !                                             ! ===============
430         END DO                                           !   End of slab
431         !                                                ! ===============
432
433      END DO
434
435
436
437      ! III. recover the avt (avs) resulting from vertical physics only
438      !---------------------------------------------------------------
439      fstravs(:,:,:) = zavs(:,:,:)
440
441     
442      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
443         WRITE(charout, FMT="('zdf - 2')")
444         CALL prt_ctl_trc_info(charout)
445         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
446      ENDIF
447
448   END SUBROUTINE trc_zdf_iso
449
450#else
451   !!----------------------------------------------------------------------
452   !!   Dummy module               NO rotation of the lateral mixing tensor
453   !!----------------------------------------------------------------------
454CONTAINS
455   SUBROUTINE trc_zdf_iso( kt )              ! empty routine
456      WRITE(*,*) 'trc_zdf_iso: You should not have seen this print! error?', kt
457   END SUBROUTINE trc_zdf_iso
458#endif
459
460   !!==============================================================================
461END MODULE trczdf_iso
Note: See TracBrowser for help on using the repository browser.