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

source: tags/nemo_v1_05/NEMO/TOP_SRC/TRP/trczdf_iso.F90 @ 2007

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

nemo_v1_update_005:RB: update headers for the TOP component.
  • 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
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   !!  TOP 1.0,  LOCEAN-IPSL (2005)
34   !! $Header$
35   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE trc_zdf_iso( kt )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE trc_zdf_iso  ***
43      !!
44      !! ** Purpose :
45      !!     Compute the trend due to the vertical tracer diffusion inclu-
46      !!     ding the vertical component of lateral mixing (only for second
47      !!     order operator, for fourth order it is already computed and
48      !!     add to the general trend in trcldf.F) and add it to the general
49      !!     trend of the tracer equations.
50      !!
51      !! ** Method :
52      !!         The vertical component of the lateral diffusive trends is
53      !!      provided by a 2nd order operator rotated along neural or geopo-
54      !!      tential surfaces to which an eddy induced advection can be added
55      !!      It is computed using before fields (forward in time) and isopyc-
56      !!      nal or geopotential slopes computed in routine ldfslp.
57      !!
58      !!      First part: vertical trends associated with the lateral mixing
59      !!      ==========  (excluding the vertical flux proportional to dk[t] )
60      !!      vertical fluxes associated with the rotated lateral mixing:
61      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(trb)) ]
62      !!                     + e1t*wslpj dj[ mj(mk(trb)) ]  }
63      !!      save avt coef. resulting from vertical physics alone in zavt:
64      !!         zavt = avt
65      !!      update and save in zavt the vertical eddy viscosity coefficient:
66      !!         avt = avt + wslpi^2+wslj^2
67      !!      add vertical Eddy Induced advective fluxes ('lk_trcldf_eiv=T):
68      !!         zftw = zftw + { di[aht e2u mi(wslpi)]
69      !!                    +dj[aht e1v mj(wslpj)] } mk(trb)
70      !!      take the horizontal divergence of the fluxes:
71      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
72      !!      Add this trend to the general trend tra :
73      !!         tra = tra + difft
74      !!
75      !!      Second part: vertical trend associated with the vertical physics
76      !!      ===========  (including the vertical flux proportional to dk[t]
77      !!                  associated with the lateral mixing, through the
78      !!                  update of avt)
79      !!      The vertical diffusion of tracers tra  is given by:
80      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
81      !!      It is computed using a backward time scheme, t=ta.
82      !!      Surface and bottom boundary conditions: no diffusive flux on
83      !!      both tracers (bottom, applied through the masked field avt).
84      !!      Add this trend to the general trend tra  :
85      !!         tra = tra + dz( avt dz(t) )
86      !!         (tra = tra + dz( avs dz(t) ) if lk_trc_zdfddm=T )
87      !!
88      !!      Third part: recover avt resulting from the vertical physics
89      !!      ==========  alone, for further diagnostics (for example to
90      !!                  compute the turbocline depth in diamld).
91      !!         avt = zavt
92      !!         (avs = zavs if lk_trc_zdfddm=T )
93      !!
94      !!      'key_trc_diatrd' defined: trend saved for futher diagnostics.
95      !!
96      !!      macro-tasked on vertical slab (jj-loop)
97      !!
98      !! ** Action :
99      !!         Update tra arrays with the before vertical diffusion trend
100      !!         Save in trtrd arrays the trends if 'key_trc_diatrd' defined
101      !!
102      !! History :
103      !!   7.0  !  91-11  (G. Madec)  Original code
104      !!        !  92-06  (M. Imbard)  correction on tracer trend loops
105      !!        !  96-01  (G. Madec)  statement function for e3
106      !!        !  97-05  (G. Madec)  vertical component of isopycnal
107      !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord
108      !!        !  98-03  (L. Bopp MA Foujols) passive tracer generalisation
109      !!        !  00-05  (MA Foujols) add lbc for tracer trends
110      !!        !  00-06  (O Aumont)  correct isopycnal scheme suppress
111      !!        !                     avt multiple correction
112      !!        !  00-08  (G. Madec)  double diffusive mixing
113      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
114      !!   9.0  !  04-03  (C. Ethe )  adapted for passive tracers
115      !!---------------------------------------------------------------------
116      !! * Modules used
117      USE oce_trc               ,   &
118         zavs => va
119
120      !! * Arguments
121      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
122
123      !! * Local declarations
124      INTEGER ::   ji, jj, jk, jn                 ! dummy loop indices
125      INTEGER ::   ikst, ikenm2, ikstp1       ! temporary integers
126#if defined key_partial_steps
127      INTEGER ::   iku, ikv, ikv1             ! temporary integers
128#endif
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      !!---------------------------------------------------------------------
151      !!  OPA 8.5, LODYC-IPSL (2002)
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         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
352            ztra = SUM( tra(2:nictle,2:njctle,1:jpkm1,jn) * tmask(2:nictle,2:njctle,1:jpkm1) )
353            WRITE(numout,*) ' trc/zdf 1  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn)
354            tra_ctl(jn) = ztra 
355         ENDIF
356
357      END DO
358
359
360      DO jn = 1, jptra
361         !                                                ! ===============
362         DO jj = 2, jpjm1                                 !  Vertical slab
363            !                                             ! ===============
364
365            ! II. Vertical trend associated with the vertical physics
366            ! =======================================================
367            !     (including the vertical flux proportional to dk[t] associated
368            !      with the lateral mixing, through the avt update)
369            !     dk[ avt dk[ (t,s) ] ] diffusive trends
370
371
372            ! Diagonal, inferior, superior
373            ! (including the bottom boundary condition via avs masked)
374            DO jk = 1, jpkm1
375               DO ji = 2, jpim1
376                  zwi(ji,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  )   &
377                     /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
378                  zws(ji,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1)   &
379                     /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
380                  zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
381               END DO
382            END DO
383
384            ! Surface boudary conditions
385            DO ji = 2, jpim1
386               zwi(ji,1) = 0.e0
387               zwd(ji,1) = 1. - zws(ji,1)
388            END DO
389
390            ! Second member construction
391            DO jk = 1, jpkm1
392               DO ji = 2, jpim1
393                  zwy(ji,jk) = trb(ji,jj,jk,jn) + rdttrc(jk) * tra(ji,jj,jk,jn)
394               END DO
395            END DO
396
397            ! Matrix inversion from the first level
398            ikst = 1
399#   include "zdf.matrixsolver.h90"
400#if defined key_trc_diatrd
401            ! Compute and save the vertical diffusive of tracers trends
402#  if defined key_trc_ldfiso
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) + trtrd(ji,jj,jk,jn,6)
407               END DO
408            END DO
409#  else
410            DO jk = 1, jpkm1
411               DO ji = 2, jpim1
412                  ztra = ( zwx(ji,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
413                  trtrd(ji,jj,jk,jn,6) = ztra - tra(ji,jj,jk,jn)
414               END DO
415            END DO
416#  endif
417#endif
418            ! Save the masked passive tracer after in tra
419            ! (c a u t i o n:  tracer not its trend, Leap-frog scheme done
420            !                  it will not be done in tranxt)
421            DO jk = 1, jpkm1
422               DO ji = 2, jpim1
423                  tra(ji,jj,jk,jn) = zwx(ji,jk)  * tmask(ji,jj,jk)
424               END DO
425            END DO
426
427
428            ! III. recover the avt (avs) resulting from vertical physics only
429            ! ===============================================================
430
431            DO jk = 2, jpkm1
432               DO ji = 2, jpim1
433                  fstravs(ji,jj,jk) = zavs(ji,jj,jk)
434               END DO
435            END DO
436
437            !                                             ! ===============
438         END DO                                           !   End of slab
439         !                                                ! ===============
440
441
442         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
443            ztra = SUM( tra(2:nictle,2:njctle,1:jpkm1,jn) * tmask(2:nictle,2:njctle,1:jpkm1) )
444            WRITE(numout,*) ' trc/zdf 2  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn)
445            tra_ctl(jn) = ztra 
446         ENDIF
447
448      END DO
449
450   END SUBROUTINE trc_zdf_iso
451
452#else
453   !!----------------------------------------------------------------------
454   !!   Dummy module               NO rotation of the lateral mixing tensor
455   !!----------------------------------------------------------------------
456CONTAINS
457   SUBROUTINE trc_zdf_iso( kt )              ! empty routine
458      WRITE(*,*) 'trc_zdf_iso: You should not have seen this print! error?', kt
459   END SUBROUTINE trc_zdf_iso
460#endif
461
462   !!==============================================================================
463END MODULE trczdf_iso
Note: See TracBrowser for help on using the repository browser.