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 tags/nemo_v1_08/NEMO/TOP_SRC/TRP – NEMO

source: tags/nemo_v1_08/NEMO/TOP_SRC/TRP/trczdf_iso.F90 @ 9353

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

nemo_v1_update_031 : CT : change header names
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.2 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$
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#if defined key_partial_steps
128      INTEGER ::   iku, ikv, ikv1             ! temporary integers
129#endif
130      REAL(wp) ::   ztra
131      REAL(wp) ::   &
132         ztavg,                 &  ! ???
133         zcoef0, zcoef3,        &  ! ???
134         zcoef4, zavi,          &  ! ???
135         zbtr, zmku, zmkv,      &  !
136         ztav
137      REAL(wp), DIMENSION(jpi,jpk) ::   &
138         zwd, zws, zwi,         &  ! ???
139         zwx, zwy, zwz, zwt        ! ???
140      REAL(wp), DIMENSION(jpi,jpk) ::   &
141         ztfw, zdit, zdjt, zdj1t
142#if defined key_trcldf_eiv   ||   defined key_esopa
143      REAL(wp), DIMENSION(jpi,jpk) ::   &
144         ztfwg
145
146      REAL(wp) ::         &
147         zcoeg3,          &
148         zuwk, zvwk,      &
149         zuwki, zvwki
150#endif
151      CHARACTER (len=22) :: charout
152      !!---------------------------------------------------------------------
153
154      IF( kt == nittrc000 ) THEN
155         IF(lwp) WRITE(numout,*)
156         IF(lwp) WRITE(numout,*) 'trc_zdf_iso : vertical mixing (including isopycnal component)'
157         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
158#if defined key_trcldf_eiv && defined key_diaeiv
159         w_trc_eiv(:,:,:) = 0.e0
160#endif
161      ENDIF
162
163      ! 0. Local constant initialization
164      ! --------------------------------
165      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
166         ! time step = 2 rdttra with Arakawa or TVD advection scheme
167         IF( neuler == 0 .AND. kt == nittrc000 ) THEN
168            rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)             ! restarting with Euler time stepping
169         ELSEIF( kt <= nittrc000 + 1 ) THEN
170            rdttrc(:) = 2. * rdttra(:) * FLOAT(ndttrc)         ! leapfrog
171         ENDIF
172      ELSE
173         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)     
174      ENDIF
175
176
177      DO jn = 1, jptra
178
179         ztavg = 0.e0
180
181         !                                                ! ===============
182         DO jj = 2, jpjm1                                 !  Vertical slab
183            !                                             ! ===============
184
185            ! I. vertical trends associated with the lateral mixing
186            ! =====================================================
187            !  (excluding the vertical flux proportional to dk[t]
188
189
190            ! I.1 horizontal tracer gradient
191            ! ------------------------------
192
193            DO jk = 1, jpkm1
194               DO ji = 1, jpim1
195                  ! i-gradient of passive tracer at jj
196                  zdit (ji,jk) = ( trb(ji+1,jj,jk,jn)-trb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
197                  ! j-gradient of passive tracer at jj
198                  zdjt (ji,jk) = ( trb(ji,jj+1,jk,jn)-trb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
199                  ! j-gradient of passive tracer at jj+1
200                  zdj1t(ji,jk) = ( trb(ji,jj,jk,jn)-trb(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
201               END DO
202            END DO
203#  if defined key_partial_steps
204            ! partial steps correction at the bottom ocean level
205            DO ji = 1, jpim1
206               ! last ocean level
207               iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
208               ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
209               ikv1 = MIN( mbathy(ji,jj), mbathy(ji  ,jj-1) ) - 1
210               ! i-gradient of of passive tracer at jj
211               zdit (ji,iku) = gtru(ji,jj,jn)
212               ! j-gradient of of passive tracer at jj
213               zdjt (ji,ikv) = gtrv(ji,jj,jn) 
214               ! j-gradient of of passive tracer at jj+1
215               zdj1t(ji,ikv1)= gtrv(ji,jj-1,jn)
216            END DO
217#endif
218
219
220            ! I.2 Vertical fluxes
221            ! -------------------
222
223            ! Surface and bottom vertical fluxes set to zero
224            ztfw(:, 1 ) = 0.e0
225            ztfw(:,jpk) = 0.e0
226
227#if defined key_trcldf_eiv
228            ztfwg(:, 1 ) = 0.e0
229            ztfwg(:,jpk) = 0.e0
230#endif
231
232            ! interior (2=<jk=<jpk-1)
233            DO jk = 2, jpkm1
234               DO ji = 2, jpim1
235                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
236
237                  zmku = 1./MAX( umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)   &
238                     &          +umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1. )
239
240                  zmkv = 1./MAX( vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)   &
241                     &          +vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1. )
242
243                  zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
244                  zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
245
246                  ztfw(ji,jk) = zcoef3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
247                     &                    +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
248                     &        + zcoef4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
249                     &                    +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )
250
251               END DO
252            END DO
253
254
255            ! I.3  update and save of avt (and avs if double diffusive mixing)
256            ! ---------------------------
257
258            DO jk = 2, jpkm1
259               DO ji = 2, jpim1
260
261                  zavi = fsahtw(ji,jj,jk)*( wslpi(ji,jj,jk)*wslpi(ji,jj,jk)   &
262                     &                     +wslpj(ji,jj,jk)*wslpj(ji,jj,jk) )
263
264                  ! save avs in zavs to recover avs in output files
265                  zavs(ji,jj,jk) = fstravs(ji,jj,jk)
266                  ! add isopycnal vertical coeff. to avs
267                  fstravs(ji,jj,jk) = fstravs(ji,jj,jk) + zavi
268
269               END DO
270            END DO
271
272#if defined key_trcldf_eiv
273            !                              ! ---------------------------------------!
274            !                              ! Eddy induced vertical advective fluxes !
275            !                              ! ---------------------------------------!
276#if defined key_traldf_c2d || defined key_traldf_c3d
277            DO jk = 2, jpkm1
278               DO ji = 2, jpim1
279                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
280                     &  * fsaeitru(ji-1,jj,jk) * e2u(ji-1,jj)*umask(ji-1,jj,jk)
281                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
282                     &  * fsaeitru(ji  ,jj,jk) * e2u(ji  ,jj)*umask(ji  ,jj,jk)
283                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
284                     &  * fsaeitrv(ji,jj-1,jk) * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
285                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
286                     &  * fsaeitrv(ji,jj  ,jk) * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
287
288                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
289
290                  ztfwg(ji,jk) = + zcoeg3 * ( trb(ji,jj,jk,jn) + trb(ji,jj,jk-1,jn) ) 
291                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jk)
292
293# if defined key_diaeiv
294                  w_trc_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
295# endif
296               END DO
297            END DO
298
299#else
300            DO jk = 2, jpkm1
301               DO ji = 2, jpim1
302                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
303                     &  * e2u(ji-1,jj)*umask(ji-1,jj,jk)
304                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
305                     &  * e2u(ji  ,jj)*umask(ji  ,jj,jk)
306                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
307                     &  * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
308                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
309                     &  * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
310
311                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeitrw(ji,jj,jk)   &
312                     &            * ( zuwk - zuwki + zvwk - zvwki )
313
314                  ztfwg(ji,jk) = + zcoeg3 * ( trb(ji,jj,jk,jn) + trb(ji,jj,jk-1,jn) )
315                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jk)
316
317# if defined key_diaeiv
318                  w_trc_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
319# endif
320               END DO
321            END DO
322#endif
323
324#endif
325
326
327            ! I.5 Divergence of vertical fluxes added to the general tracer trend
328            ! -------------------------------------------------------------------
329
330            DO jk = 1, jpkm1
331               DO ji = 2, jpim1
332                  zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
333                  ztav = (  ztfw(ji,jk) - ztfw(ji,jk+1)  ) * zbtr
334                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztav
335
336#if defined key_trc_diatrd
337#   if defined key_trcldf_eiv
338                  ztavg = ( ztfwg(ji,jk) - ztfwg(ji,jk+1) ) * zbtr
339                  !  WARNING trtrd(ji,jj,jk,6) used for vertical gent velocity trend
340                  !                           not for damping !!!
341                  trtrd(ji,jj,jk,jn,6) = ztavg
342#   endif
343                  trtrd(ji,jj,jk,jn,6) = ztav - ztavg
344#endif
345               END DO
346            END DO
347            !                                             ! ===============
348         END DO                                           !   End of slab
349         !                                                ! ===============
350
351      END DO
352
353      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
354         WRITE(charout, FMT="('zdf - 1')")
355         CALL prt_ctl_trc_info(charout)
356         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
357      ENDIF
358
359      DO jn = 1, jptra
360         !                                                ! ===============
361         DO jj = 2, jpjm1                                 !  Vertical slab
362            !                                             ! ===============
363
364            ! II. Vertical trend associated with the vertical physics
365            ! =======================================================
366            !     (including the vertical flux proportional to dk[t] associated
367            !      with the lateral mixing, through the avt update)
368            !     dk[ avt dk[ (t,s) ] ] diffusive trends
369
370
371            ! Diagonal, inferior, superior
372            ! (including the bottom boundary condition via avs masked)
373            DO jk = 1, jpkm1
374               DO ji = 2, jpim1
375                  zwi(ji,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  )   &
376                     /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
377                  zws(ji,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1)   &
378                     /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
379                  zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
380               END DO
381            END DO
382
383            ! Surface boudary conditions
384            DO ji = 2, jpim1
385               zwi(ji,1) = 0.e0
386               zwd(ji,1) = 1. - zws(ji,1)
387            END DO
388
389            ! Second member construction
390            DO jk = 1, jpkm1
391               DO ji = 2, jpim1
392                  zwy(ji,jk) = trb(ji,jj,jk,jn) + rdttrc(jk) * tra(ji,jj,jk,jn)
393               END DO
394            END DO
395
396            ! Matrix inversion from the first level
397            ikst = 1
398#   include "zdf.matrixsolver.h90"
399#if defined key_trc_diatrd
400            ! Compute and save the vertical diffusive of tracers trends
401#  if defined key_trc_ldfiso
402            DO jk = 1, jpkm1
403               DO ji = 2, jpim1
404                  ztra = ( zwx(ji,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
405                  trtrd(ji,jj,jk,jn,6) = ztra - tra(ji,jj,jk,jn) + trtrd(ji,jj,jk,jn,6)
406               END DO
407            END DO
408#  else
409            DO jk = 1, jpkm1
410               DO ji = 2, jpim1
411                  ztra = ( zwx(ji,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
412                  trtrd(ji,jj,jk,jn,6) = ztra - tra(ji,jj,jk,jn)
413               END DO
414            END DO
415#  endif
416#endif
417            ! Save the masked passive tracer after in tra
418            ! (c a u t i o n:  tracer not its trend, Leap-frog scheme done
419            !                  it will not be done in tranxt)
420            DO jk = 1, jpkm1
421               DO ji = 2, jpim1
422                  tra(ji,jj,jk,jn) = zwx(ji,jk)  * tmask(ji,jj,jk)
423               END DO
424            END DO
425
426
427            ! III. recover the avt (avs) resulting from vertical physics only
428            ! ===============================================================
429
430            DO jk = 2, jpkm1
431               DO ji = 2, jpim1
432                  fstravs(ji,jj,jk) = zavs(ji,jj,jk)
433               END DO
434            END DO
435
436            !                                             ! ===============
437         END DO                                           !   End of slab
438         !                                                ! ===============
439
440      END DO
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.