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

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

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