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

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

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