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

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

CT : UPDATE142 : Check the consistency between passive tracers transport modules (in TRP directory) and those used for the active tracers
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 22.0 KB
Line 
1MODULE trcadv_smolar
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_smolar  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_adv_smolar : update the passive tracer trend with the horizontal
9   !!                  and vertical advection trends using a Smolarkiewicz
10   !!                  FCT scheme
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce_trc             ! ocean dynamics and active tracers variables
14   USE trc                 ! ocean passive tracers variables
15   USE lbclnk              ! ocean lateral boundary conditions (or mpp link)
16   USE trcbbl              ! advective passive tracers in the BBL
17
18   IMPLICIT NONE
19   PRIVATE
20
21   !! * Accessibility
22   PUBLIC trc_adv_smolar    ! routine called by trcstp.F90
23
24   !! * Module variable
25   REAL(wp), DIMENSION(jpk) ::   &
26      rdttrc                     ! vertical profile of tracer time-step
27 
28   !! * Substitutions
29#  include "passivetrc_substitute.h90"
30   !!----------------------------------------------------------------------
31   !!   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
180         ! Surface value
181         IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN        ! free surface-constant volume
182            zkz(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn) * tmask(ji,jj,1)
183         ELSE                                                   ! rigid lid : flux set to zero
184            zkz(:,:, 1 ) = 0.e0 
185         ENDIF
186
187        DO jk = 2,jpk
188          DO jj = 1,jpj
189            DO ji = 1,jpi
190               zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
191               zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
192               zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
193            END DO
194          END DO
195        END DO
196
197! ... Lateral boundary conditions on zk[xy]
198      CALL lbc_lnk( zkx, 'U', -1. )
199      CALL lbc_lnk( zky, 'V', -1. )
200
201
202! 2. calcul of after field using an upstream advection scheme
203! -----------------------------------------------------------
204
205        DO jk = 1,jpkm1
206          DO jj = 2,jpjm1
207            DO ji = 2,jpim1
208              zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
209              ztj(ji,jj,jk) = -zbtr*    &
210     &            ( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk)  &
211     &            + zky(ji,jj,jk) - zky(ji,jj - 1,jk)  &
212     &            + zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
213#if defined key_trc_diatrd
214              trtrd(ji,jj,jk,jn,1) = trtrd(ji,jj,jk,jn,1) -  &
215     &                       zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )
216
217              trtrd(ji,jj,jk,jn,2) = trtrd(ji,jj,jk,jn,2) -  &
218     &            zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
219
220              trtrd(ji,jj,jk,jn,3) = trtrd(ji,jj,jk,jn,3) -  &
221     &            zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
222#endif
223            END DO
224          END DO
225        END DO
226
227! 2.1 start of antidiffusive correction loop
228
229        DO jt = 1,ncortrc
230
231! 2.2 calcul of intermediary field zti
232
233          DO jk = 1,jpkm1
234            DO jj = 2,jpjm1
235              DO ji = 2,jpim1
236                zti(ji,jj,jk) = zti(ji,jj,jk)+rdttrc(jk)*ztj(ji,jj,jk)
237                zbuf(ji,jj,jk) = zbuf(ji,jj,jk) + ztj(ji,jj,jk)
238              END DO
239            END DO
240          END DO
241
242! ... Lateral boundary conditions on zti
243      CALL lbc_lnk( zti, 'T', 1. )
244
245
246! 2.3 calcul of the antidiffusive flux
247
248          DO jk = 1,jpkm1
249            DO jj = 2,jpjm1
250              DO ji = 2,jpim1
251                zx(ji,jj,jk) = ( abs(zaa(ji,jj,jk)) - rdttrc(jk)       &
252     &              *zaa(ji,jj,jk)**2/                          &
253     &              (e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) )  &
254     &              *(zti(ji + 1,jj,jk) - zti( ji ,jj,jk))      &
255     &              /(zti( ji ,jj,jk) + zti(ji + 1,jj,jk) + rtrn)    &
256     &              * rsc
257
258                zy(ji,jj,jk) = ( abs(zbb(ji,jj,jk)) - rdttrc(jk)       &
259     &              *zbb(ji,jj,jk)**2/                          &
260     &              (e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) )  &
261     &              *(zti(ji,jj + 1,jk) - zti(ji, jj ,jk))      &
262     &              /(zti(ji, jj ,jk) + zti(ji,jj + 1,jk) + rtrn)    &
263     &              * rsc
264              END DO
265            END DO
266          END DO
267
268          DO jk = 2,jpkm1
269            DO jj = 2,jpjm1
270              DO ji = 2,jpim1
271                zz(ji,jj,jk) = ( abs(zcc(ji,jj,jk)) - rdttrc(jk)*zcc(ji,jj,jk)**2  & 
272     &              /( e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk) ) ) &
273     &               *( zti(ji,jj,jk) - zti(ji,jj,jk - 1) )/  &
274     &                ( zti(ji,jj,jk) + zti(ji,jj,jk - 1) + rtrn )* rsc*( -1.)
275              END DO
276            END DO
277          END DO
278
279! 2.4 cross terms
280
281          IF (crosster) THEN
282              DO jk = 2,jpkm1
283                DO jj = 2,jpjm1
284                  DO ji = 2,jpim1
285                    zx(ji,jj,jk) = zx(ji,jj,jk) &
286     &                  - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,jk)*0.25* &
287     &                  (    (zbb(ji  ,jj - 1,jk  ) + zbb(ji + 1,jj - 1 &
288     &                  ,jk  ) + zbb(ji + 1,jj  ,jk  ) + zbb(ji  ,jj  &
289     &                  ,jk))* (zti(ji  ,jj + 1,jk  ) + zti(ji + 1,jj + &
290     &                  1,jk  ) - zti(ji + 1,jj - 1,jk  ) - zti(ji  ,jj &
291     &                  - 1,jk  ))/ (zti(ji  ,jj + 1,jk  ) + zti(ji + 1 &
292     &                  ,jj + 1,jk  ) + zti(ji + 1,jj - 1,jk  ) + zti(ji &
293     &                  ,jj - 1,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  ) + &
294     &                  zcc(ji + 1,jj  ,jk  ) + zcc(ji  ,jj  ,jk + 1) + &
295     &                  zcc(ji + 1,jj  ,jk + 1))* (zti(ji  ,jj  ,jk - 1) &
296     &                  + zti(ji + 1,jj  ,jk - 1) - zti(ji  ,jj  ,jk + 1 &
297     &                  )- zti(ji + 1,jj  ,jk + 1))/ (zti(ji  ,jj  ,jk - &
298     &                  1) + zti(ji + 1,jj  ,jk - 1) + zti(ji  ,jj  ,jk &
299     &                  +1) + zti(ji + 1,jj  ,jk + 1) + rtrn))/(e1u(ji &
300     &                  ,jj)*e2u(ji,jj)*fse3u(ji,jj,jk))*vmask(ji  ,jj - &
301     &                  1,jk  )*vmask(ji + 1,jj - 1,jk  )*vmask(ji + 1 &
302     &                  ,jj,jk)*vmask(ji  ,jj  ,jk  )*tmask(ji  ,jj  ,jk &
303     &                  )*tmask(ji + 1,jj  ,jk  )*tmask(ji  ,jj  ,jk + 1 &
304     &                  )*tmask(ji + 1,jj  ,jk + 1)
305
306                    zy(ji,jj,jk) = zy(ji,jj,jk)    &   
307     &                  - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,jk)*0.25*    &   
308     &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji - 1,jj + 1    &   
309     &                  ,jk  ) + zaa(ji  ,jj  ,jk  ) + zaa(ji  ,jj + 1    &   
310     &                  ,jk))* (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1,jj    &   
311     &                  ,jk  ) - zti(ji - 1,jj + 1,jk  ) - zti(ji - 1,jj    &   
312     &                  ,jk  ))/ (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1    &   
313     &                  ,jj  ,jk  ) + zti(ji - 1,jj + 1,jk  ) + zti(ji    &   
314     &                  - 1,jj  ,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  )    &   
315     &                  + zcc(ji  ,jj  ,jk + 1) + zcc(ji  ,jj + 1,jk  )    &   
316     &                  + zcc(ji  ,jj + 1,jk + 1))* (zti(ji  ,jj  ,jk -    &   
317     &                  1) + zti(ji  ,jj + 1,jk - 1) - zti(ji  ,jj  ,jk    &   
318     &                  +1) - zti(ji  ,jj + 1,jk + 1))/ (zti(ji  ,jj    &   
319     &                  ,jk- 1) + zti(ji  ,jj + 1,jk - 1) + zti(ji  ,jj    &   
320     &                  ,jk+ 1) + zti(ji  ,jj + 1,jk + 1) + rtrn))    &   
321     &                  /(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk))    &   
322     &                  *umask(ji - 1,jj,jk  )*umask(ji - 1,jj + 1,jk  )    &   
323     &                  *umask(ji  ,jj,jk  )*umask(ji  ,jj + 1,jk  )    &   
324     &                  *tmask(ji  ,jj,jk)*tmask(ji  ,jj  ,jk + 1)    &   
325     &                  *tmask(ji  ,jj + 1,jk)*tmask(ji  ,jj + 1,jk + 1)     
326
327                    zz(ji,jj,jk) = zz(ji,jj,jk)    &   
328     &                  - 0.5*rdttrc(jk)*rsc*zcc(ji,jj,jk)*0.25*    &   
329     &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji  ,jj  ,jk    &   
330     &                  ) + zaa(ji  ,jj  ,jk - 1) + zaa(ji - 1,jj  ,jk -    &   
331     &                  1))*(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1,jj    &   
332     &                  ,jk  ) - zti(ji - 1,jj  ,jk  ) - zti(ji - 1,jj    &   
333     &                  ,jk - 1))/(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1    &   
334     &                  ,jj,jk  ) + zti(ji - 1,jj  ,jk  ) + zti(ji - 1    &   
335     &                  ,jj,jk - 1) + rtrn) + (zbb(ji  ,jj - 1,jk  )    &   
336     &                  + zbb(ji  ,jj  ,jk  ) + zbb(ji  ,jj  ,jk - 1)    &   
337     &                  + zbb(ji  ,jj - 1,jk - 1))*(zti(ji  ,jj + 1,jk -    &   
338     &                  1) + zti(ji  ,jj + 1,jk  ) - zti(ji  ,jj - 1,jk    &   
339     &                  ) - zti(ji  ,jj - 1,jk - 1))/(zti(ji  ,jj + 1,jk    &   
340     &                  - 1) + zti(ji  ,jj + 1,jk  ) + zti(ji  ,jj - 1    &   
341     &                  ,jk  ) + zti(ji  ,jj - 1,jk - 1) + rtrn))    &   
342     &                  /(e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk))    &   
343     &                  *umask(ji - 1,jj,jk  )*umask(ji  ,jj  ,jk  )    &   
344     &                  *umask(ji  ,jj,jk- 1)*umask(ji - 1,jj  ,jk - 1)    &   
345     &                  *vmask(ji  ,jj- 1,jk)*vmask(ji  ,jj  ,jk  )    &   
346     &                  *vmask(ji  ,jj  ,jk-1)*vmask(ji  ,jj - 1,jk - 1)       
347                  END DO
348                END DO
349              END DO
350
351              DO jj = 2,jpjm1
352                DO ji = 2,jpim1
353                  zx(ji,jj,1) = zx(ji,jj,1)    &   
354     &                - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,1)*0.25*    &   
355     &                ( (zbb(ji  ,jj - 1,1  ) + zbb(ji + 1,jj - 1,1  )    &   
356     &                + zbb(ji + 1,jj  ,1  ) + zbb(ji  ,jj  ,1  ))    &   
357     &                *(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )    &   
358     &                - zti(ji + 1,jj - 1,1  ) - zti(ji  ,jj - 1,1  ))    &   
359     &                /(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )    &   
360     &                + zti(ji + 1,jj - 1,1  ) + zti(ji  ,jj - 1,1  ) +    &   
361     &                rtrn))/(e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,1))    &   
362     &                *vmask(ji  ,jj - 1,1  )*vmask(ji + 1,jj - 1,1  )    &   
363     &                *vmask(ji + 1,jj  ,1  )*vmask(ji  ,jj  ,1  )   
364
365                 zy(ji,jj,1) = zy(ji,jj,1)    &   
366     &                - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,1)*0.25*    &   
367     &                ( (zaa(ji-1  ,jj ,1  ) + zaa(ji - 1,jj + 1,1  )    &   
368     &                + zaa(ji ,jj  ,1  ) + zaa(ji  ,jj + 1  ,1  ))    &   
369     &                *(zti(ji + 1,jj + 1,1  ) + zti(ji + 1,jj ,1  )    &   
370     &                - zti(ji - 1,jj + 1,1  ) - zti(ji - 1,jj ,1  ))    &   
371     &                /(zti(ji + 1,jj + 1,1  ) + zti(ji + 1,jj ,1  )    &   
372     &                + zti(ji - 1,jj + 1,1  ) + zti(ji - 1,jj ,1  ) +    &   
373     &                rtrn))/(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,1))    &   
374     &                *umask(ji - 1,jj,1  )*umask(ji - 1,jj + 1,1  )    &   
375     &                *umask(ji    ,jj,1  )*umask(ji  ,jj + 1 ,1  )   
376
377                END DO
378              END DO
379          ENDIF
380
381          ! ... Lateral boundary conditions on z[xyz]
382          CALL lbc_lnk( zx, 'U', -1. )
383          CALL lbc_lnk( zy, 'V', -1. )
384          CALL lbc_lnk( zz, 'W',  1. )
385
386! 2.4 reinitialization
387
388          DO jk = 1,jpk
389            DO jj = 1,jpj
390              DO ji = 1,jpi
391                zaa(ji,jj,jk) = zx(ji,jj,jk)
392                zbb(ji,jj,jk) = zy(ji,jj,jk)
393                zcc(ji,jj,jk) = zz(ji,jj,jk)
394              END DO
395            END DO
396          END DO
397
398! 2.5 calcul of the final field:
399!    advection by antidiffusive mass fluxes and an upstream scheme
400
401          DO jk = 1,jpk
402             DO jj = 2,jpjm1
403                DO ji = 2,jpim1
404                   zfp_ui = 0.5 * ( zaa(ji,jj,jk) + ABS( zaa(ji,jj,jk) ) )
405                   zfp_vj = 0.5 * ( zbb(ji,jj,jk) + ABS( zbb(ji,jj,jk) ) )
406                   zfm_ui = 0.5 * ( zaa(ji,jj,jk) - ABS( zaa(ji,jj,jk) ) )
407                   zfm_vj = 0.5 * ( zbb(ji,jj,jk) - ABS( zbb(ji,jj,jk) ) )           
408                   zkx(ji,jj,jk) = zfp_ui * zti(ji,jj,jk) + zfm_ui * zti(ji+1,jj  ,jk) 
409                   zky(ji,jj,jk) = zfp_vj * zti(ji,jj,jk) + zfm_vj * zti(ji    ,jj+1,jk)             
410                END DO
411             END DO
412          END DO
413
414          DO jk = 2,jpk
415             DO jj = 1,jpj
416                DO ji = 1,jpi
417                   zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
418                   zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
419                   zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
420                END DO
421             END DO
422          END DO
423
424
425! ... Lateral boundary conditions on zk[xy]
426      CALL lbc_lnk( zkx, 'U', -1. )
427      CALL lbc_lnk( zky, 'V', -1. )
428
429
430! 2.6. calcul of after field using an upstream advection scheme
431
432          DO jk = 1,jpkm1
433            DO jj = 2,jpjm1
434              DO ji = 2,jpim1
435                zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
436                ztj(ji,jj,jk) = -zbtr*     & 
437     &              ( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk)    & 
438     &              + zky(ji,jj,jk) - zky(ji,jj - 1,jk)    & 
439     &              + zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
440#if defined key_trc_diatrd
441                trtrd(ji,jj,jk,jn,1) = trtrd(ji,jj,jk,jn,1) -    & 
442     &              zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )   
443
444                trtrd(ji,jj,jk,jn,2) = trtrd(ji,jj,jk,jn,2) -    & 
445     &              zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
446
447                trtrd(ji,jj,jk,jn,3) = trtrd(ji,jj,jk,jn,3) -    & 
448     &              zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
449#endif
450              END DO
451            END DO
452          END DO
453
454! 2.6 END of antidiffusive correction loop
455
456        END DO
457
458! 3. trend due to horizontal and vertical advection of tracer jn
459! --------------------------------------------------------------
460
461        DO jk = 1,jpk
462          DO jj = 2,jpjm1
463            DO ji = 2,jpim1
464              ztra = ( zbuf(ji,jj,jk) + ztj(ji,jj,jk) ) * tmask(ji,jj,jk)
465              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
466            END DO
467          END DO
468        END DO
469
470! 4.0 convert the transport trend into advection trend
471! ----------------------------------------------------
472
473#if defined key_trc_diatrd
474        DO jk = 1,jpk
475          DO jj = 2,jpjm1
476            DO  ji = 2,jpim1
477              zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
478              zgm = zbtr * trn(ji,jj,jk,jn) *     & 
479     &            ( zun(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk)    & 
480     &            -zun(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk))
481
482              zgz = zbtr * trn(ji,jj,jk,jn) *     & 
483     &            ( zvn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk)    & 
484     &            -zvn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk))
485
486              trtrd(ji,jj,jk,jn,1) = trtrd(ji,jj,jk,jn,1) + zgm
487              trtrd(ji,jj,jk,jn,2) = trtrd(ji,jj,jk,jn,2) + zgz
488              trtrd(ji,jj,jk,jn,3) = trtrd(ji,jj,jk,jn,3)    & 
489     &            - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
490            END DO
491          END DO
492        END DO
493
494        ! Lateral boundary conditions on trtrd:
495
496        CALL lbc_lnk( trtrd(1,1,1,jn,1), 'T', 1. )
497        CALL lbc_lnk( trtrd(1,1,1,jn,2), 'T', 1. )
498        CALL lbc_lnk( trtrd(1,1,1,jn,3), 'T', 1. )
499#endif
500
501        IF(l_ctl) THEN         ! print mean trends (used for debugging)
502           ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
503           WRITE(numout,*) ' trc/zad  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' smolar' 
504           tra_ctl(jn) = ztra 
505        ENDIF
506        ! END of tracer loop
507        ! ==================
508     ENDDO
509     
510  END SUBROUTINE trc_adv_smolar
511
512#else
513   !!----------------------------------------------------------------------
514   !!   Default option                                         Empty module
515   !!----------------------------------------------------------------------
516CONTAINS
517   SUBROUTINE trc_adv_smolar( kt ) 
518      INTEGER, INTENT(in) :: kt
519      WRITE(*,*) 'trc_adv_smolar: You should not have seen this print! error?', kt
520   END SUBROUTINE trc_adv_smolar
521#endif
522
523   !!======================================================================
524END MODULE trcadv_smolar
Note: See TracBrowser for help on using the repository browser.