source: trunk/NEMO/TOP_SRC/TRP/trczdf_iso.F90 @ 941

Last change on this file since 941 was 941, checked in by cetlod, 13 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143

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