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

Last change on this file since 1119 was 1119, checked in by cetlod, 16 years ago

style of all top namelist has been modified ; update modules to take it into account, see ticket:196
  • Property svn:executable set to *
File size: 22.6 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 trp_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 trctrp_lec          ! passive tracers transport
19   USE prtctl_trc          ! Print control for debbuging
20
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC trc_adv_smolar    ! routine called by trcstp.F90
25
26   REAL(wp), DIMENSION(jpk) ::   rdttrc       ! vertical profile of tracer time-step
27 
28   !! * Substitutions
29#  include "top_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#if defined key_trc_diatrd
99      REAL(wp) :: zgm, zgz
100#endif
101
102      REAL(wp) :: zbtr, ztra
103      REAL(wp) :: zfp_ui, zfp_vj, zfm_ui, zfm_vj, zfp_w, zfm_w
104      CHARACTER (len=22) :: charout
105      !!----------------------------------------------------------------------
106
107
108      IF( kt == nittrc000  .AND. lwp ) THEN
109         WRITE(numout,*)
110         WRITE(numout,*) 'trc_adv_smolar : SMOLARKIEWICZ advection scheme'
111         WRITE(numout,*) '~~~~~~~~~~~~~~~'
112         rdttrc(:) = rdttra(:) * FLOAT(ndttrc)
113      ENDIF
114
115
116#if defined key_trcbbl_adv       
117      ! Advective bottom boundary layer
118      ! -------------------------------
119      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
120      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
121      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl( :,:,:)
122#endif
123
124      ! tracer loop parallelized (macrotasking)
125      ! =======================================
126     
127      DO jn = 1, jptra
128         
129         ! 1. tracer flux in the 3 directions
130         ! ----------------------------------
131         
132         ! 1.1 mass flux at u v and t-points and initialization
133
134        DO jk = 1,jpk
135
136           DO jj = 1,jpj
137              DO ji = 1,jpi
138                 zaa(ji,jj,jk) = e2u(ji,jj)*fse3u(ji,jj,jk) * zun(ji,jj,jk)
139                 zbb(ji,jj,jk) = e1v(ji,jj)*fse3v(ji,jj,jk) * zvn(ji,jj,jk)
140                 zcc(ji,jj,jk) = e1t(ji,jj)*e2t(ji,jj)      * zwn(ji,jj,jk)
141                 zbuf(ji,jj,jk) = 0.
142                 ztj(ji,jj,jk) = 0.
143                 zx(ji,jj,jk) = 0.
144                 zy(ji,jj,jk) = 0.
145                 zz(ji,jj,jk) = 0.
146                 zti(ji,jj,jk) = trn(ji,jj,jk,jn)
147#if defined key_trc_diatrd
148                 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = 0.
149                 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = 0.
150                 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = 0.
151#endif
152              END DO
153           END DO
154           
155           ! 1.2 calcul of intermediate field with an upstream advection scheme
156           !     and mass fluxes calculated above
157           
158           ! calcul of tracer flux in the i and j direction
159           
160           DO jj=1,jpj
161              zkx(  1,jj,jk)=0.
162              zkx(jpi,jj,jk)=0.
163           END DO
164           
165           DO ji=1,jpi
166              zky(ji,  1,jk)=0.
167              zky(ji,jpj,jk)=0.
168           END DO
169           
170           DO jj = 2,jpjm1
171              DO ji = 2,jpim1
172                 zfp_ui = 0.5 * ( zaa(ji,jj,jk) + ABS( zaa(ji,jj,jk) ) )
173                 zfp_vj = 0.5 * ( zbb(ji,jj,jk) + ABS( zbb(ji,jj,jk) ) )
174                 zfm_ui = 0.5 * ( zaa(ji,jj,jk) - ABS( zaa(ji,jj,jk) ) )
175                 zfm_vj = 0.5 * ( zbb(ji,jj,jk) - ABS( zbb(ji,jj,jk) ) )           
176                 zkx(ji,jj,jk) = zfp_ui * zti(ji,jj,jk) + zfm_ui * zti(ji+1,jj  ,jk) 
177                 zky(ji,jj,jk) = zfp_vj * zti(ji,jj,jk) + zfm_vj * zti(ji    ,jj+1,jk)             
178              END DO
179           END DO
180
181        END DO
182
183         ! II. Vertical advection
184         ! ----------------------
185
186         ! Surface value
187         IF( lk_dynspg_rl ) THEN        ! rigid lid : flux set to zero
188            zkz(:,:, 1 ) = 0.e0 
189         ELSE                           ! free surface
190            zkz(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn) * tmask(ji,jj,1)
191         ENDIF
192
193        DO jk = 2,jpk
194          DO jj = 1,jpj
195            DO ji = 1,jpi
196               zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
197               zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
198               zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
199            END DO
200          END DO
201        END DO
202
203! ... Lateral boundary conditions on zk[xy]
204      CALL lbc_lnk( zkx, 'U', -1. )
205      CALL lbc_lnk( zky, 'V', -1. )
206
207
208! 2. calcul of after field using an upstream advection scheme
209! -----------------------------------------------------------
210
211        DO jk = 1,jpkm1
212          DO jj = 2,jpjm1
213            DO ji = 2,jpim1
214              zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
215              ztj(ji,jj,jk) = -zbtr*    &
216     &            ( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk)  &
217     &            + zky(ji,jj,jk) - zky(ji,jj - 1,jk)  &
218     &            + zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
219#if defined key_trc_diatrd
220              IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = trtrd(ji,jj,jk,ikeep(jn),1) -  &
221     &                       zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )
222
223              IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = trtrd(ji,jj,jk,ikeep(jn),2) -  &
224     &            zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
225
226              IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = trtrd(ji,jj,jk,ikeep(jn),3) -  &
227     &            zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
228#endif
229            END DO
230          END DO
231        END DO
232
233! 2.1 start of antidiffusive correction loop
234
235        DO jt = 1,ncortrc
236
237! 2.2 calcul of intermediary field zti
238
239          DO jk = 1,jpkm1
240            DO jj = 2,jpjm1
241              DO ji = 2,jpim1
242                zti(ji,jj,jk) = zti(ji,jj,jk)+rdttrc(jk)*ztj(ji,jj,jk)
243                zbuf(ji,jj,jk) = zbuf(ji,jj,jk) + ztj(ji,jj,jk)
244              END DO
245            END DO
246          END DO
247
248! ... Lateral boundary conditions on zti
249      CALL lbc_lnk( zti, 'T', 1. )
250
251
252! 2.3 calcul of the antidiffusive flux
253
254          DO jk = 1,jpkm1
255            DO jj = 2,jpjm1
256              DO ji = 2,jpim1
257                zx(ji,jj,jk) = ( abs(zaa(ji,jj,jk)) - rdttrc(jk)       &
258     &              *zaa(ji,jj,jk)**2/                          &
259     &              (e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) )  &
260     &              *(zti(ji + 1,jj,jk) - zti( ji ,jj,jk))      &
261     &              /(zti( ji ,jj,jk) + zti(ji + 1,jj,jk) + rtrn)    &
262     &              * rsc
263
264                zy(ji,jj,jk) = ( abs(zbb(ji,jj,jk)) - rdttrc(jk)       &
265     &              *zbb(ji,jj,jk)**2/                          &
266     &              (e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) )  &
267     &              *(zti(ji,jj + 1,jk) - zti(ji, jj ,jk))      &
268     &              /(zti(ji, jj ,jk) + zti(ji,jj + 1,jk) + rtrn)    &
269     &              * rsc
270              END DO
271            END DO
272          END DO
273
274          DO jk = 2,jpkm1
275            DO jj = 2,jpjm1
276              DO ji = 2,jpim1
277                zz(ji,jj,jk) = ( abs(zcc(ji,jj,jk)) - rdttrc(jk)*zcc(ji,jj,jk)**2  & 
278     &              /( e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk) ) ) &
279     &               *( zti(ji,jj,jk) - zti(ji,jj,jk - 1) )/  &
280     &                ( zti(ji,jj,jk) + zti(ji,jj,jk - 1) + rtrn )* rsc*( -1.)
281              END DO
282            END DO
283          END DO
284
285! 2.4 cross terms
286
287          IF (crosster) THEN
288              DO jk = 2,jpkm1
289                DO jj = 2,jpjm1
290                  DO ji = 2,jpim1
291                    zx(ji,jj,jk) = zx(ji,jj,jk) &
292     &                  - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,jk)*0.25* &
293     &                  (    (zbb(ji  ,jj - 1,jk  ) + zbb(ji + 1,jj - 1 &
294     &                  ,jk  ) + zbb(ji + 1,jj  ,jk  ) + zbb(ji  ,jj  &
295     &                  ,jk))* (zti(ji  ,jj + 1,jk  ) + zti(ji + 1,jj + &
296     &                  1,jk  ) - zti(ji + 1,jj - 1,jk  ) - zti(ji  ,jj &
297     &                  - 1,jk  ))/ (zti(ji  ,jj + 1,jk  ) + zti(ji + 1 &
298     &                  ,jj + 1,jk  ) + zti(ji + 1,jj - 1,jk  ) + zti(ji &
299     &                  ,jj - 1,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  ) + &
300     &                  zcc(ji + 1,jj  ,jk  ) + zcc(ji  ,jj  ,jk + 1) + &
301     &                  zcc(ji + 1,jj  ,jk + 1))* (zti(ji  ,jj  ,jk - 1) &
302     &                  + zti(ji + 1,jj  ,jk - 1) - zti(ji  ,jj  ,jk + 1 &
303     &                  )- zti(ji + 1,jj  ,jk + 1))/ (zti(ji  ,jj  ,jk - &
304     &                  1) + zti(ji + 1,jj  ,jk - 1) + zti(ji  ,jj  ,jk &
305     &                  +1) + zti(ji + 1,jj  ,jk + 1) + rtrn))/(e1u(ji &
306     &                  ,jj)*e2u(ji,jj)*fse3u(ji,jj,jk))*vmask(ji  ,jj - &
307     &                  1,jk  )*vmask(ji + 1,jj - 1,jk  )*vmask(ji + 1 &
308     &                  ,jj,jk)*vmask(ji  ,jj  ,jk  )*tmask(ji  ,jj  ,jk &
309     &                  )*tmask(ji + 1,jj  ,jk  )*tmask(ji  ,jj  ,jk + 1 &
310     &                  )*tmask(ji + 1,jj  ,jk + 1)
311
312                    zy(ji,jj,jk) = zy(ji,jj,jk)    &   
313     &                  - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,jk)*0.25*    &   
314     &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji - 1,jj + 1    &   
315     &                  ,jk  ) + zaa(ji  ,jj  ,jk  ) + zaa(ji  ,jj + 1    &   
316     &                  ,jk))* (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1,jj    &   
317     &                  ,jk  ) - zti(ji - 1,jj + 1,jk  ) - zti(ji - 1,jj    &   
318     &                  ,jk  ))/ (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1    &   
319     &                  ,jj  ,jk  ) + zti(ji - 1,jj + 1,jk  ) + zti(ji    &   
320     &                  - 1,jj  ,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  )    &   
321     &                  + zcc(ji  ,jj  ,jk + 1) + zcc(ji  ,jj + 1,jk  )    &   
322     &                  + zcc(ji  ,jj + 1,jk + 1))* (zti(ji  ,jj  ,jk -    &   
323     &                  1) + zti(ji  ,jj + 1,jk - 1) - zti(ji  ,jj  ,jk    &   
324     &                  +1) - zti(ji  ,jj + 1,jk + 1))/ (zti(ji  ,jj    &   
325     &                  ,jk- 1) + zti(ji  ,jj + 1,jk - 1) + zti(ji  ,jj    &   
326     &                  ,jk+ 1) + zti(ji  ,jj + 1,jk + 1) + rtrn))    &   
327     &                  /(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk))    &   
328     &                  *umask(ji - 1,jj,jk  )*umask(ji - 1,jj + 1,jk  )    &   
329     &                  *umask(ji  ,jj,jk  )*umask(ji  ,jj + 1,jk  )    &   
330     &                  *tmask(ji  ,jj,jk)*tmask(ji  ,jj  ,jk + 1)    &   
331     &                  *tmask(ji  ,jj + 1,jk)*tmask(ji  ,jj + 1,jk + 1)     
332
333                    zz(ji,jj,jk) = zz(ji,jj,jk)    &   
334     &                  - 0.5*rdttrc(jk)*rsc*zcc(ji,jj,jk)*0.25*    &   
335     &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji  ,jj  ,jk    &   
336     &                  ) + zaa(ji  ,jj  ,jk - 1) + zaa(ji - 1,jj  ,jk -    &   
337     &                  1))*(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1,jj    &   
338     &                  ,jk  ) - zti(ji - 1,jj  ,jk  ) - zti(ji - 1,jj    &   
339     &                  ,jk - 1))/(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1    &   
340     &                  ,jj,jk  ) + zti(ji - 1,jj  ,jk  ) + zti(ji - 1    &   
341     &                  ,jj,jk - 1) + rtrn) + (zbb(ji  ,jj - 1,jk  )    &   
342     &                  + zbb(ji  ,jj  ,jk  ) + zbb(ji  ,jj  ,jk - 1)    &   
343     &                  + zbb(ji  ,jj - 1,jk - 1))*(zti(ji  ,jj + 1,jk -    &   
344     &                  1) + zti(ji  ,jj + 1,jk  ) - zti(ji  ,jj - 1,jk    &   
345     &                  ) - zti(ji  ,jj - 1,jk - 1))/(zti(ji  ,jj + 1,jk    &   
346     &                  - 1) + zti(ji  ,jj + 1,jk  ) + zti(ji  ,jj - 1    &   
347     &                  ,jk  ) + zti(ji  ,jj - 1,jk - 1) + rtrn))    &   
348     &                  /(e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk))    &   
349     &                  *umask(ji - 1,jj,jk  )*umask(ji  ,jj  ,jk  )    &   
350     &                  *umask(ji  ,jj,jk- 1)*umask(ji - 1,jj  ,jk - 1)    &   
351     &                  *vmask(ji  ,jj- 1,jk)*vmask(ji  ,jj  ,jk  )    &   
352     &                  *vmask(ji  ,jj  ,jk-1)*vmask(ji  ,jj - 1,jk - 1)       
353                  END DO
354                END DO
355              END DO
356
357              DO jj = 2,jpjm1
358                DO ji = 2,jpim1
359                  zx(ji,jj,1) = zx(ji,jj,1)    &   
360     &                - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,1)*0.25*    &   
361     &                ( (zbb(ji  ,jj - 1,1  ) + zbb(ji + 1,jj - 1,1  )    &   
362     &                + zbb(ji + 1,jj  ,1  ) + zbb(ji  ,jj  ,1  ))    &   
363     &                *(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )    &   
364     &                - zti(ji + 1,jj - 1,1  ) - zti(ji  ,jj - 1,1  ))    &   
365     &                /(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )    &   
366     &                + zti(ji + 1,jj - 1,1  ) + zti(ji  ,jj - 1,1  ) +    &   
367     &                rtrn))/(e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,1))    &   
368     &                *vmask(ji  ,jj - 1,1  )*vmask(ji + 1,jj - 1,1  )    &   
369     &                *vmask(ji + 1,jj  ,1  )*vmask(ji  ,jj  ,1  )   
370
371                 zy(ji,jj,1) = zy(ji,jj,1)    &   
372     &                - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,1)*0.25*    &   
373     &                ( (zaa(ji-1  ,jj ,1  ) + zaa(ji - 1,jj + 1,1  )    &   
374     &                + zaa(ji ,jj  ,1  ) + zaa(ji  ,jj + 1  ,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     &                + zti(ji - 1,jj + 1,1  ) + zti(ji - 1,jj ,1  ) +    &   
379     &                rtrn))/(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,1))    &   
380     &                *umask(ji - 1,jj,1  )*umask(ji - 1,jj + 1,1  )    &   
381     &                *umask(ji    ,jj,1  )*umask(ji  ,jj + 1 ,1  )   
382
383                END DO
384              END DO
385          ENDIF
386
387          ! ... Lateral boundary conditions on z[xyz]
388          CALL lbc_lnk( zx, 'U', -1. )
389          CALL lbc_lnk( zy, 'V', -1. )
390          CALL lbc_lnk( zz, 'W',  1. )
391
392! 2.4 reinitialization
393
394          DO jk = 1,jpk
395            DO jj = 1,jpj
396              DO ji = 1,jpi
397                zaa(ji,jj,jk) = zx(ji,jj,jk)
398                zbb(ji,jj,jk) = zy(ji,jj,jk)
399                zcc(ji,jj,jk) = zz(ji,jj,jk)
400              END DO
401            END DO
402          END DO
403
404! 2.5 calcul of the final field:
405!    advection by antidiffusive mass fluxes and an upstream scheme
406
407          DO jk = 1,jpk
408             DO jj = 2,jpjm1
409                DO ji = 2,jpim1
410                   zfp_ui = 0.5 * ( zaa(ji,jj,jk) + ABS( zaa(ji,jj,jk) ) )
411                   zfp_vj = 0.5 * ( zbb(ji,jj,jk) + ABS( zbb(ji,jj,jk) ) )
412                   zfm_ui = 0.5 * ( zaa(ji,jj,jk) - ABS( zaa(ji,jj,jk) ) )
413                   zfm_vj = 0.5 * ( zbb(ji,jj,jk) - ABS( zbb(ji,jj,jk) ) )           
414                   zkx(ji,jj,jk) = zfp_ui * zti(ji,jj,jk) + zfm_ui * zti(ji+1,jj  ,jk) 
415                   zky(ji,jj,jk) = zfp_vj * zti(ji,jj,jk) + zfm_vj * zti(ji    ,jj+1,jk)             
416                END DO
417             END DO
418          END DO
419
420          DO jk = 2,jpk
421             DO jj = 1,jpj
422                DO ji = 1,jpi
423                   zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
424                   zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
425                   zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
426                END DO
427             END DO
428          END DO
429
430
431! ... Lateral boundary conditions on zk[xy]
432      CALL lbc_lnk( zkx, 'U', -1. )
433      CALL lbc_lnk( zky, 'V', -1. )
434
435
436! 2.6. calcul of after field using an upstream advection scheme
437
438          DO jk = 1,jpkm1
439            DO jj = 2,jpjm1
440              DO ji = 2,jpim1
441                zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
442                ztj(ji,jj,jk) = -zbtr*     & 
443     &              ( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk)    & 
444     &              + zky(ji,jj,jk) - zky(ji,jj - 1,jk)    & 
445     &              + zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
446#if defined key_trc_diatrd
447                IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = trtrd(ji,jj,jk,ikeep(jn),1) -    & 
448     &              zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )   
449
450                IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = trtrd(ji,jj,jk,ikeep(jn),2) -    & 
451     &              zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
452
453                IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = trtrd(ji,jj,jk,ikeep(jn),3) -    & 
454     &              zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
455#endif
456              END DO
457            END DO
458          END DO
459
460! 2.6 END of antidiffusive correction loop
461
462        END DO
463
464! 3. trend due to horizontal and vertical advection of tracer jn
465! --------------------------------------------------------------
466
467        DO jk = 1,jpk
468          DO jj = 2,jpjm1
469            DO ji = 2,jpim1
470              ztra = ( zbuf(ji,jj,jk) + ztj(ji,jj,jk) ) * tmask(ji,jj,jk)
471              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
472            END DO
473          END DO
474        END DO
475
476! 4.0 convert the transport trend into advection trend
477! ----------------------------------------------------
478
479#if defined key_trc_diatrd
480        DO jk = 1,jpk
481          DO jj = 2,jpjm1
482            DO  ji = 2,jpim1
483              zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
484              zgm = zbtr * trn(ji,jj,jk,jn) *     & 
485     &            ( zun(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk)    & 
486     &            -zun(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk))
487
488              zgz = zbtr * trn(ji,jj,jk,jn) *     & 
489     &            ( zvn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk)    & 
490     &            -zvn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk))
491
492              IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = trtrd(ji,jj,jk,ikeep(jn),1) + zgm
493              IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = trtrd(ji,jj,jk,ikeep(jn),2) + zgz
494              IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = trtrd(ji,jj,jk,ikeep(jn),3)    & 
495     &            - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
496            END DO
497          END DO
498        END DO
499
500        ! Lateral boundary conditions on trtrd:
501
502        IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),1), 'T', 1. )
503        IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),2), 'T', 1. )
504        IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),3), 'T', 1. )
505#endif
506
507 
508        ! END of tracer loop
509        ! ==================
510     ENDDO
511
512      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
513         WRITE(charout, FMT="('smolar - adv')")
514         CALL prt_ctl_trc_info(charout)
515         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
516      ENDIF
517     
518  END SUBROUTINE trc_adv_smolar
519
520#else
521   !!----------------------------------------------------------------------
522   !!   Default option                                         Empty module
523   !!----------------------------------------------------------------------
524CONTAINS
525   SUBROUTINE trc_adv_smolar( kt ) 
526      INTEGER, INTENT(in) :: kt
527      WRITE(*,*) 'trc_adv_smolar: You should not have seen this print! error?', kt
528   END SUBROUTINE trc_adv_smolar
529#endif
530
531   !!======================================================================
532END MODULE trcadv_smolar
Note: See TracBrowser for help on using the repository browser.