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

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

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