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