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.
trcadv_smolar.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trcadv_smolar.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: 22.1 KB
Line 
1MODULE trcadv_smolar
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_smolar  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_adv_smolar : update the passive tracer trend with the horizontal
9   !!                  and vertical advection trends using a Smolarkiewicz
10   !!                  FCT scheme
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce_trc             ! ocean dynamics and active tracers variables
14   USE trc                 ! ocean passive tracers variables
15   USE lbclnk              ! ocean lateral boundary conditions (or mpp link)
16   USE trcbbl              ! advective passive tracers in the BBL
17
18   IMPLICIT NONE
19   PRIVATE
20
21   !! * Accessibility
22   PUBLIC trc_adv_smolar    ! routine called by trcstp.F90
23
24   !! * Module variable
25   REAL(wp), DIMENSION(jpk) ::   &
26      rdttrc                     ! vertical profile of tracer time-step
27 
28   !! * Substitutions
29#  include "passivetrc_substitute.h90"
30   !!----------------------------------------------------------------------
31   !!   TOP 1.0 , LOCEAN-IPSL (2005)
32   !! $Header$
33   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE trc_adv_smolar( kt )
38      !!----------------------------------------------------------------------
39      !!                   ***  ROUTINE trc_adv_smolar  ***
40      !!
41      !! ** Purpose :   Compute the now trend due to total advection of passi-
42      !!      ve tracer using a Smolarkiewicz FCT (Flux Corrected Transport )
43      !!      scheme and add it to the general tracer trend.
44      !!
45      !! ** Method : Computation of not exactly the advection but the
46      !!             transport term, i.e.  div(u*tra).
47      !!             Computes the now horizontal and vertical advection with
48      !!             the complete 3d method.
49      !!
50      !!       note: - sc is an empirical factor to be used with care
51      !!             - this advection scheme needs an euler-forward time scheme
52      !!
53      !! ** Action : - update tra with the now advective tracer trends
54      !!             - save trends in trtrd ('key_trc_diatrd')
55      !!
56      !! References :               
57      !!     Piotr K. Smolarkiewicz, 1983,
58      !!       "A simple positive definit advection
59      !!        scheme with small IMPLICIT diffusion"
60      !!        Monthly Weather Review, pp 479-486
61      !!
62      !! History :
63      !!         !  87-06 (pa-dl) Original
64      !!         !  91-11 (G. Madec)
65      !!         !  94-08 (A. Czaja)
66      !!         !  95-09 (M. Levy) passive tracers
67      !!         !  98-03 (M.A. Foujols) lateral boundary conditions
68      !!         !  99-02 (M.A. Foujols) lbc in conjonction with ORCA
69      !!         !  00-05 (MA Foujols) add lbc for tracer trends
70      !!         !  00-10 (MA Foujols and E.Kestenare) INCLUDE instead of routine
71      !!         !  01-05 (E.Kestenare) fix bug in trtrd indexes
72      !!         !  02-05 (M-A Filiberti, and M.Levy) correction in trtrd computation
73      !!   9.0   !  03-04  (C. Ethe)  F90: Free form and module
74      !!----------------------------------------------------------------------
75      !! * modules used
76#if defined key_trcbbl_adv
77      USE oce_trc            , zun => ua,  &  ! use ua as workspace
78         &                     zvn => va      ! use va as workspace
79      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
80#else
81      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
82                               zvn => vn,  &  !              zvn == vn
83                               zwn => wn      !              zwn == wn
84#endif
85      !! * Arguments
86      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
87
88      !! * Local declarations
89      INTEGER :: ji, jj, jk,jt, jn            ! dummy loop indices
90
91      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &     
92         zti, ztj,              &
93         zaa, zbb, zcc,         &
94         zx , zy , zz ,         &
95         zkx, zky, zkz,         &
96         zbuf
97
98      REAL(wp) :: zbtr, ztra
99      REAL(wp) :: zfp_ui, zfp_vj, zfm_ui, zfm_vj, zfp_w, zfm_w
100
101      !!----------------------------------------------------------------------
102
103
104      IF( kt == nittrc000  .AND. lwp ) THEN
105         WRITE(numout,*)
106         WRITE(numout,*) 'trc_adv_smolar : SMOLARKIEWICZ advection scheme'
107         WRITE(numout,*) '~~~~~~~~~~~~~~~'
108         rdttrc(:) = rdttra(:) * FLOAT(ndttrc)
109      ENDIF
110
111
112#if defined key_trcbbl_adv       
113      ! Advective bottom boundary layer
114      ! -------------------------------
115      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
116      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
117      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl( :,:,:)
118#endif
119
120      ! tracer loop parallelized (macrotasking)
121      ! =======================================
122     
123      DO jn = 1, jptra
124         
125         ! 1. tracer flux in the 3 directions
126         ! ----------------------------------
127         
128         ! 1.1 mass flux at u v and t-points and initialization
129
130        DO jk = 1,jpk
131
132           DO jj = 1,jpj
133              DO ji = 1,jpi
134                 zaa(ji,jj,jk) = e2u(ji,jj)*fse3u(ji,jj,jk) * zun(ji,jj,jk)
135                 zbb(ji,jj,jk) = e1v(ji,jj)*fse3v(ji,jj,jk) * zvn(ji,jj,jk)
136                 zcc(ji,jj,jk) = e1t(ji,jj)*e2t(ji,jj)      * zwn(ji,jj,jk)
137                 zbuf(ji,jj,jk) = 0.
138                 ztj(ji,jj,jk) = 0.
139                 zx(ji,jj,jk) = 0.
140                 zy(ji,jj,jk) = 0.
141                 zz(ji,jj,jk) = 0.
142                 zti(ji,jj,jk) = trn(ji,jj,jk,jn)
143#if defined key_trc_diatrd
144                 trtrd(ji,jj,jk,jn,1) = 0.
145                 trtrd(ji,jj,jk,jn,2) = 0.
146                 trtrd(ji,jj,jk,jn,3) = 0.
147#endif
148              END DO
149           END DO
150           
151           ! 1.2 calcul of intermediate field with an upstream advection scheme
152           !     and mass fluxes calculated above
153           
154           ! calcul of tracer flux in the i and j direction
155           
156           DO jj=1,jpj
157              zkx(  1,jj,jk)=0.
158              zkx(jpi,jj,jk)=0.
159           END DO
160           
161           DO ji=1,jpi
162              zky(ji,  1,jk)=0.
163              zky(ji,jpj,jk)=0.
164           END DO
165           
166           DO jj = 2,jpjm1
167              DO ji = 2,jpim1
168                 zfp_ui = 0.5 * ( zaa(ji,jj,jk) + ABS( zaa(ji,jj,jk) ) )
169                 zfp_vj = 0.5 * ( zbb(ji,jj,jk) + ABS( zbb(ji,jj,jk) ) )
170                 zfm_ui = 0.5 * ( zaa(ji,jj,jk) - ABS( zaa(ji,jj,jk) ) )
171                 zfm_vj = 0.5 * ( zbb(ji,jj,jk) - ABS( zbb(ji,jj,jk) ) )           
172                 zkx(ji,jj,jk) = zfp_ui * zti(ji,jj,jk) + zfm_ui * zti(ji+1,jj  ,jk) 
173                 zky(ji,jj,jk) = zfp_vj * zti(ji,jj,jk) + zfm_vj * zti(ji    ,jj+1,jk)             
174              END DO
175           END DO
176
177        END DO
178
179         ! II. Vertical advection
180         ! ----------------------
181
182         ! Surface value
183         IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN        ! free surface-constant volume
184            zkz(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn) * tmask(ji,jj,1)
185         ELSE                                                   ! rigid lid : flux set to zero
186            zkz(:,:, 1 ) = 0.e0 
187         ENDIF
188
189        DO jk = 2,jpk
190          DO jj = 1,jpj
191            DO ji = 1,jpi
192               zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
193               zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
194               zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
195            END DO
196          END DO
197        END DO
198
199! ... Lateral boundary conditions on zk[xy]
200      CALL lbc_lnk( zkx, 'U', -1. )
201      CALL lbc_lnk( zky, 'V', -1. )
202
203
204! 2. calcul of after field using an upstream advection scheme
205! -----------------------------------------------------------
206
207        DO jk = 1,jpkm1
208          DO jj = 2,jpjm1
209            DO ji = 2,jpim1
210              zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
211              ztj(ji,jj,jk) = -zbtr*    &
212     &            ( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk)  &
213     &            + zky(ji,jj,jk) - zky(ji,jj - 1,jk)  &
214     &            + zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
215#if defined key_trc_diatrd
216              trtrd(ji,jj,jk,jn,1) = trtrd(ji,jj,jk,jn,1) -  &
217     &                       zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )
218
219              trtrd(ji,jj,jk,jn,2) = trtrd(ji,jj,jk,jn,2) -  &
220     &            zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
221
222              trtrd(ji,jj,jk,jn,3) = trtrd(ji,jj,jk,jn,3) -  &
223     &            zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
224#endif
225            END DO
226          END DO
227        END DO
228
229! 2.1 start of antidiffusive correction loop
230
231        DO jt = 1,ncortrc
232
233! 2.2 calcul of intermediary field zti
234
235          DO jk = 1,jpkm1
236            DO jj = 2,jpjm1
237              DO ji = 2,jpim1
238                zti(ji,jj,jk) = zti(ji,jj,jk)+rdttrc(jk)*ztj(ji,jj,jk)
239                zbuf(ji,jj,jk) = zbuf(ji,jj,jk) + ztj(ji,jj,jk)
240              END DO
241            END DO
242          END DO
243
244! ... Lateral boundary conditions on zti
245      CALL lbc_lnk( zti, 'T', 1. )
246
247
248! 2.3 calcul of the antidiffusive flux
249
250          DO jk = 1,jpkm1
251            DO jj = 2,jpjm1
252              DO ji = 2,jpim1
253                zx(ji,jj,jk) = ( abs(zaa(ji,jj,jk)) - rdttrc(jk)       &
254     &              *zaa(ji,jj,jk)**2/                          &
255     &              (e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) )  &
256     &              *(zti(ji + 1,jj,jk) - zti( ji ,jj,jk))      &
257     &              /(zti( ji ,jj,jk) + zti(ji + 1,jj,jk) + rtrn)    &
258     &              * rsc
259
260                zy(ji,jj,jk) = ( abs(zbb(ji,jj,jk)) - rdttrc(jk)       &
261     &              *zbb(ji,jj,jk)**2/                          &
262     &              (e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) )  &
263     &              *(zti(ji,jj + 1,jk) - zti(ji, jj ,jk))      &
264     &              /(zti(ji, jj ,jk) + zti(ji,jj + 1,jk) + rtrn)    &
265     &              * rsc
266              END DO
267            END DO
268          END DO
269
270          DO jk = 2,jpkm1
271            DO jj = 2,jpjm1
272              DO ji = 2,jpim1
273                zz(ji,jj,jk) = ( abs(zcc(ji,jj,jk)) - rdttrc(jk)*zcc(ji,jj,jk)**2  & 
274     &              /( e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk) ) ) &
275     &               *( zti(ji,jj,jk) - zti(ji,jj,jk - 1) )/  &
276     &                ( zti(ji,jj,jk) + zti(ji,jj,jk - 1) + rtrn )* rsc*( -1.)
277              END DO
278            END DO
279          END DO
280
281! 2.4 cross terms
282
283          IF (crosster) THEN
284              DO jk = 2,jpkm1
285                DO jj = 2,jpjm1
286                  DO ji = 2,jpim1
287                    zx(ji,jj,jk) = zx(ji,jj,jk) &
288     &                  - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,jk)*0.25* &
289     &                  (    (zbb(ji  ,jj - 1,jk  ) + zbb(ji + 1,jj - 1 &
290     &                  ,jk  ) + zbb(ji + 1,jj  ,jk  ) + zbb(ji  ,jj  &
291     &                  ,jk))* (zti(ji  ,jj + 1,jk  ) + zti(ji + 1,jj + &
292     &                  1,jk  ) - zti(ji + 1,jj - 1,jk  ) - zti(ji  ,jj &
293     &                  - 1,jk  ))/ (zti(ji  ,jj + 1,jk  ) + zti(ji + 1 &
294     &                  ,jj + 1,jk  ) + zti(ji + 1,jj - 1,jk  ) + zti(ji &
295     &                  ,jj - 1,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  ) + &
296     &                  zcc(ji + 1,jj  ,jk  ) + zcc(ji  ,jj  ,jk + 1) + &
297     &                  zcc(ji + 1,jj  ,jk + 1))* (zti(ji  ,jj  ,jk - 1) &
298     &                  + zti(ji + 1,jj  ,jk - 1) - zti(ji  ,jj  ,jk + 1 &
299     &                  )- zti(ji + 1,jj  ,jk + 1))/ (zti(ji  ,jj  ,jk - &
300     &                  1) + zti(ji + 1,jj  ,jk - 1) + zti(ji  ,jj  ,jk &
301     &                  +1) + zti(ji + 1,jj  ,jk + 1) + rtrn))/(e1u(ji &
302     &                  ,jj)*e2u(ji,jj)*fse3u(ji,jj,jk))*vmask(ji  ,jj - &
303     &                  1,jk  )*vmask(ji + 1,jj - 1,jk  )*vmask(ji + 1 &
304     &                  ,jj,jk)*vmask(ji  ,jj  ,jk  )*tmask(ji  ,jj  ,jk &
305     &                  )*tmask(ji + 1,jj  ,jk  )*tmask(ji  ,jj  ,jk + 1 &
306     &                  )*tmask(ji + 1,jj  ,jk + 1)
307
308                    zy(ji,jj,jk) = zy(ji,jj,jk)    &   
309     &                  - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,jk)*0.25*    &   
310     &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji - 1,jj + 1    &   
311     &                  ,jk  ) + zaa(ji  ,jj  ,jk  ) + zaa(ji  ,jj + 1    &   
312     &                  ,jk))* (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1,jj    &   
313     &                  ,jk  ) - zti(ji - 1,jj + 1,jk  ) - zti(ji - 1,jj    &   
314     &                  ,jk  ))/ (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1    &   
315     &                  ,jj  ,jk  ) + zti(ji - 1,jj + 1,jk  ) + zti(ji    &   
316     &                  - 1,jj  ,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  )    &   
317     &                  + zcc(ji  ,jj  ,jk + 1) + zcc(ji  ,jj + 1,jk  )    &   
318     &                  + zcc(ji  ,jj + 1,jk + 1))* (zti(ji  ,jj  ,jk -    &   
319     &                  1) + zti(ji  ,jj + 1,jk - 1) - zti(ji  ,jj  ,jk    &   
320     &                  +1) - zti(ji  ,jj + 1,jk + 1))/ (zti(ji  ,jj    &   
321     &                  ,jk- 1) + zti(ji  ,jj + 1,jk - 1) + zti(ji  ,jj    &   
322     &                  ,jk+ 1) + zti(ji  ,jj + 1,jk + 1) + rtrn))    &   
323     &                  /(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk))    &   
324     &                  *umask(ji - 1,jj,jk  )*umask(ji - 1,jj + 1,jk  )    &   
325     &                  *umask(ji  ,jj,jk  )*umask(ji  ,jj + 1,jk  )    &   
326     &                  *tmask(ji  ,jj,jk)*tmask(ji  ,jj  ,jk + 1)    &   
327     &                  *tmask(ji  ,jj + 1,jk)*tmask(ji  ,jj + 1,jk + 1)     
328
329                    zz(ji,jj,jk) = zz(ji,jj,jk)    &   
330     &                  - 0.5*rdttrc(jk)*rsc*zcc(ji,jj,jk)*0.25*    &   
331     &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji  ,jj  ,jk    &   
332     &                  ) + zaa(ji  ,jj  ,jk - 1) + zaa(ji - 1,jj  ,jk -    &   
333     &                  1))*(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1,jj    &   
334     &                  ,jk  ) - zti(ji - 1,jj  ,jk  ) - zti(ji - 1,jj    &   
335     &                  ,jk - 1))/(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1    &   
336     &                  ,jj,jk  ) + zti(ji - 1,jj  ,jk  ) + zti(ji - 1    &   
337     &                  ,jj,jk - 1) + rtrn) + (zbb(ji  ,jj - 1,jk  )    &   
338     &                  + zbb(ji  ,jj  ,jk  ) + zbb(ji  ,jj  ,jk - 1)    &   
339     &                  + zbb(ji  ,jj - 1,jk - 1))*(zti(ji  ,jj + 1,jk -    &   
340     &                  1) + zti(ji  ,jj + 1,jk  ) - zti(ji  ,jj - 1,jk    &   
341     &                  ) - zti(ji  ,jj - 1,jk - 1))/(zti(ji  ,jj + 1,jk    &   
342     &                  - 1) + zti(ji  ,jj + 1,jk  ) + zti(ji  ,jj - 1    &   
343     &                  ,jk  ) + zti(ji  ,jj - 1,jk - 1) + rtrn))    &   
344     &                  /(e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk))    &   
345     &                  *umask(ji - 1,jj,jk  )*umask(ji  ,jj  ,jk  )    &   
346     &                  *umask(ji  ,jj,jk- 1)*umask(ji - 1,jj  ,jk - 1)    &   
347     &                  *vmask(ji  ,jj- 1,jk)*vmask(ji  ,jj  ,jk  )    &   
348     &                  *vmask(ji  ,jj  ,jk-1)*vmask(ji  ,jj - 1,jk - 1)       
349                  END DO
350                END DO
351              END DO
352
353              DO jj = 2,jpjm1
354                DO ji = 2,jpim1
355                  zx(ji,jj,1) = zx(ji,jj,1)    &   
356     &                - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,1)*0.25*    &   
357     &                ( (zbb(ji  ,jj - 1,1  ) + zbb(ji + 1,jj - 1,1  )    &   
358     &                + zbb(ji + 1,jj  ,1  ) + zbb(ji  ,jj  ,1  ))    &   
359     &                *(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )    &   
360     &                - zti(ji + 1,jj - 1,1  ) - zti(ji  ,jj - 1,1  ))    &   
361     &                /(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )    &   
362     &                + zti(ji + 1,jj - 1,1  ) + zti(ji  ,jj - 1,1  ) +    &   
363     &                rtrn))/(e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,1))    &   
364     &                *vmask(ji  ,jj - 1,1  )*vmask(ji + 1,jj - 1,1  )    &   
365     &                *vmask(ji + 1,jj  ,1  )*vmask(ji  ,jj  ,1  )   
366
367                 zy(ji,jj,1) = zy(ji,jj,1)    &   
368     &                - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,1)*0.25*    &   
369     &                ( (zaa(ji-1  ,jj ,1  ) + zaa(ji - 1,jj + 1,1  )    &   
370     &                + zaa(ji ,jj  ,1  ) + zaa(ji  ,jj + 1  ,1  ))    &   
371     &                *(zti(ji + 1,jj + 1,1  ) + zti(ji + 1,jj ,1  )    &   
372     &                - zti(ji - 1,jj + 1,1  ) - zti(ji - 1,jj ,1  ))    &   
373     &                /(zti(ji + 1,jj + 1,1  ) + zti(ji + 1,jj ,1  )    &   
374     &                + zti(ji - 1,jj + 1,1  ) + zti(ji - 1,jj ,1  ) +    &   
375     &                rtrn))/(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,1))    &   
376     &                *umask(ji - 1,jj,1  )*umask(ji - 1,jj + 1,1  )    &   
377     &                *umask(ji    ,jj,1  )*umask(ji  ,jj + 1 ,1  )   
378
379                END DO
380              END DO
381          ENDIF
382
383          ! ... Lateral boundary conditions on z[xyz]
384          CALL lbc_lnk( zx, 'U', -1. )
385          CALL lbc_lnk( zy, 'V', -1. )
386          CALL lbc_lnk( zz, 'W',  1. )
387
388! 2.4 reinitialization
389
390          DO jk = 1,jpk
391            DO jj = 1,jpj
392              DO ji = 1,jpi
393                zaa(ji,jj,jk) = zx(ji,jj,jk)
394                zbb(ji,jj,jk) = zy(ji,jj,jk)
395                zcc(ji,jj,jk) = zz(ji,jj,jk)
396              END DO
397            END DO
398          END DO
399
400! 2.5 calcul of the final field:
401!    advection by antidiffusive mass fluxes and an upstream scheme
402
403          DO jk = 1,jpk
404             DO jj = 2,jpjm1
405                DO ji = 2,jpim1
406                   zfp_ui = 0.5 * ( zaa(ji,jj,jk) + ABS( zaa(ji,jj,jk) ) )
407                   zfp_vj = 0.5 * ( zbb(ji,jj,jk) + ABS( zbb(ji,jj,jk) ) )
408                   zfm_ui = 0.5 * ( zaa(ji,jj,jk) - ABS( zaa(ji,jj,jk) ) )
409                   zfm_vj = 0.5 * ( zbb(ji,jj,jk) - ABS( zbb(ji,jj,jk) ) )           
410                   zkx(ji,jj,jk) = zfp_ui * zti(ji,jj,jk) + zfm_ui * zti(ji+1,jj  ,jk) 
411                   zky(ji,jj,jk) = zfp_vj * zti(ji,jj,jk) + zfm_vj * zti(ji    ,jj+1,jk)             
412                END DO
413             END DO
414          END DO
415
416          DO jk = 2,jpk
417             DO jj = 1,jpj
418                DO ji = 1,jpi
419                   zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
420                   zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
421                   zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
422                END DO
423             END DO
424          END DO
425
426
427! ... Lateral boundary conditions on zk[xy]
428      CALL lbc_lnk( zkx, 'U', -1. )
429      CALL lbc_lnk( zky, 'V', -1. )
430
431
432! 2.6. calcul of after field using an upstream advection scheme
433
434          DO jk = 1,jpkm1
435            DO jj = 2,jpjm1
436              DO ji = 2,jpim1
437                zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
438                ztj(ji,jj,jk) = -zbtr*     & 
439     &              ( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk)    & 
440     &              + zky(ji,jj,jk) - zky(ji,jj - 1,jk)    & 
441     &              + zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
442#if defined key_trc_diatrd
443                trtrd(ji,jj,jk,jn,1) = trtrd(ji,jj,jk,jn,1) -    & 
444     &              zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )   
445
446                trtrd(ji,jj,jk,jn,2) = trtrd(ji,jj,jk,jn,2) -    & 
447     &              zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
448
449                trtrd(ji,jj,jk,jn,3) = trtrd(ji,jj,jk,jn,3) -    & 
450     &              zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
451#endif
452              END DO
453            END DO
454          END DO
455
456! 2.6 END of antidiffusive correction loop
457
458        END DO
459
460! 3. trend due to horizontal and vertical advection of tracer jn
461! --------------------------------------------------------------
462
463        DO jk = 1,jpk
464          DO jj = 2,jpjm1
465            DO ji = 2,jpim1
466              ztra = ( zbuf(ji,jj,jk) + ztj(ji,jj,jk) ) * tmask(ji,jj,jk)
467              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
468            END DO
469          END DO
470        END DO
471
472! 4.0 convert the transport trend into advection trend
473! ----------------------------------------------------
474
475#if defined key_trc_diatrd
476        DO jk = 1,jpk
477          DO jj = 2,jpjm1
478            DO  ji = 2,jpim1
479              zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
480              zgm = zbtr * trn(ji,jj,jk,jn) *     & 
481     &            ( zun(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk)    & 
482     &            -zun(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk))
483
484              zgz = zbtr * trn(ji,jj,jk,jn) *     & 
485     &            ( zvn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk)    & 
486     &            -zvn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk))
487
488              trtrd(ji,jj,jk,jn,1) = trtrd(ji,jj,jk,jn,1) + zgm
489              trtrd(ji,jj,jk,jn,2) = trtrd(ji,jj,jk,jn,2) + zgz
490              trtrd(ji,jj,jk,jn,3) = trtrd(ji,jj,jk,jn,3)    & 
491     &            - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
492            END DO
493          END DO
494        END DO
495
496        ! Lateral boundary conditions on trtrd:
497
498        CALL lbc_lnk( trtrd(1,1,1,jn,1), 'T', 1. )
499        CALL lbc_lnk( trtrd(1,1,1,jn,2), 'T', 1. )
500        CALL lbc_lnk( trtrd(1,1,1,jn,3), 'T', 1. )
501#endif
502
503        IF(l_ctl) THEN         ! print mean trends (used for debugging)
504           ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
505           WRITE(numout,*) ' trc/zad  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' smolar' 
506           tra_ctl(jn) = ztra 
507        ENDIF
508        ! END of tracer loop
509        ! ==================
510     ENDDO
511     
512  END SUBROUTINE trc_adv_smolar
513
514#else
515   !!----------------------------------------------------------------------
516   !!   Default option                                         Empty module
517   !!----------------------------------------------------------------------
518CONTAINS
519   SUBROUTINE trc_adv_smolar( kt ) 
520      INTEGER, INTENT(in) :: kt
521      WRITE(*,*) 'trc_adv_smolar: You should not have seen this print! error?', kt
522   END SUBROUTINE trc_adv_smolar
523#endif
524
525   !!======================================================================
526END MODULE trcadv_smolar
Note: See TracBrowser for help on using the repository browser.