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

Last change on this file since 197 was 186, checked in by opalod, 20 years ago

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