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

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

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