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

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

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