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 branches/dev_001_GM/NEMO/TOP_SRC/TRP – NEMO

source: branches/dev_001_GM/NEMO/TOP_SRC/TRP/trcadv_smolar.F90 @ 772

Last change on this file since 772 was 772, checked in by gm, 16 years ago

dev_001_GM - change the name of cpp key to key_top, key_lobster, key_pisces, key_kriest and the corresponding lk_

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