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

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

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