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

Last change on this file since 1528 was 1528, checked in by rblod, 15 years ago

Suppress rigid-lid option, see ticket #486
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 30.5 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_vvl ) THEN
169            ! variable volume: flux set to zero
170            zkz(:,:, 1 ) = 0.e0 
171         ELSE
172            ! free surface-constant volume
173            zkz(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn) * tmask(ji,jj,1)
174         ENDIF
175         
176         DO jk = 2, jpk     ! Vector opt. ???
177            DO jj = 1, jpj
178               DO ji = 1, jpi
179                  zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
180                  zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
181                  zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
182               END DO
183            END DO
184         END DO
185         
186         ! 2. Compute after field using an upstream advection scheme
187         ! ---------------------------------------------------------
188         
189         DO jk = 1, jpkm1
190            DO jj = 2, jpjm1
191               DO ji = fs_2, fs_jpim1  ! Vector opt.
192                  zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
193                  ztj(ji,jj,jk) = - zbtr *                                &
194                       &         ( zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )    &
195                       &         + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )    &
196                       &         + zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1) )
197#if defined key_trc_diatrd
198                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = trtrd(ji,jj,jk,ikeep(jn),1) -  &
199                       &                                        zbtr*( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )
200
201                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = trtrd(ji,jj,jk,ikeep(jn),2) -  &
202                       &                                        zbtr*( zky(ji,jj,jk) - zky(ji,jj-1,jk) )
203
204                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = trtrd(ji,jj,jk,ikeep(jn),3) -  &
205                       &                                        zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
206#endif
207
208               END DO
209            END DO
210         END DO
211         
212         ! 3. Diagnose passive tracer trends (part 1/3)
213         ! --------------------------------------------
214
215         IF( l_trdtrc ) THEN
216            DO jk = 1, jpkm1
217               DO jj = 2, jpjm1
218                  DO ji = fs_2, fs_jpim1   ! Vector opt.
219                     zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
220                     ztrtrd(ji,jj,jk,1) = ztrtrd(ji,jj,jk,1) - zbtr*( zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  ) )
221                     ztrtrd(ji,jj,jk,2) = ztrtrd(ji,jj,jk,2) - zbtr*( zky(ji,jj,jk) - zky(ji  ,jj-1,jk  ) )
222                     ztrtrd(ji,jj,jk,3) = ztrtrd(ji,jj,jk,3) - zbtr*( zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1) )
223                  END DO
224               END DO
225            END DO
226         ENDIF
227         
228         ! 4. Antidiffusive correction loop
229         ! --------------------------------
230         !                                      ! -----------------------------
231         DO jt = 1, ncortrc                     ! antidiffusive correction loop
232            !                                   ! -----------------------------
233           
234            !--4.1 Compute intermediate field zti
235            DO jk = 1, jpkm1
236               zti (:,:,jk) = zti (:,:,jk) + rdttrc(jk) * ztj(:,:,jk)
237            END DO
238            zbuf(:,:,:) = zbuf(:,:,:) + ztj(:,:,:)
239           
240            CALL lbc_lnk( zti, 'T', 1. )    ! lateral boundary
241           
242            !--4.2 Compute the antidiffusive fluxes
243            DO jk = 1, jpkm1
244               DO jj = 2, jpjm1
245                  DO ji = fs_2, fs_jpim1    ! Vector opt.
246                     zx(ji,jj,jk) = ( abs(zaa(ji,jj,jk)) - rdttrc(jk)                     &
247                          &              *zaa(ji,jj,jk)**2/                               &
248                          &              (e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) )       &
249                          &              *(zti(ji + 1,jj,jk) - zti( ji ,jj,jk))           &
250                          &              /(zti( ji ,jj,jk) + zti(ji + 1,jj,jk) + rtrn)    &
251                          &              * rsc
252                     
253                     zy(ji,jj,jk) = ( abs(zbb(ji,jj,jk)) - rdttrc(jk)                     &
254                          &              *zbb(ji,jj,jk)**2/                               &
255                          &              (e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) )       &
256                          &              *(zti(ji,jj + 1,jk) - zti(ji, jj ,jk))           &
257                          &              /(zti(ji, jj ,jk) + zti(ji,jj + 1,jk) + rtrn)    &
258                          &              * rsc
259                  END DO
260               END DO
261            END DO
262           
263            DO jk = 2, jpkm1
264               DO jj = 2, jpjm1
265                  DO ji = fs_2, fs_jpim1    ! Vector opt.
266                     zz(ji,jj,jk) = ( abs(zcc(ji,jj,jk)) - rdttrc(jk)*zcc(ji,jj,jk)**2    & 
267                          &              /( e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk) ) )     &
268                          &               *( zti(ji,jj,jk) - zti(ji,jj,jk - 1) )/         &
269                          &                ( zti(ji,jj,jk) + zti(ji,jj,jk - 1) + rtrn )* rsc*( -1.)
270                  END DO
271               END DO
272            END DO
273           
274            !--4.3 Cross terms
275            CROSSTERMS: IF( crosster ) THEN
276               !
277               DO jk = 2, jpkm1
278                  DO jj = 2, jpjm1
279                     DO ji = fs_2, fs_jpim1    ! Vector opt.
280                        zx(ji,jj,jk) = zx(ji,jj,jk)                                               &
281                             &                  - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,jk)*0.25*          &
282                             &                  (    (zbb(ji  ,jj - 1,jk  ) + zbb(ji + 1,jj - 1   &
283                             &                  ,jk  ) + zbb(ji + 1,jj  ,jk  ) + zbb(ji  ,jj      &
284                             &                  ,jk))* (zti(ji  ,jj + 1,jk  ) + zti(ji + 1,jj +   &
285                             &                  1,jk  ) - zti(ji + 1,jj - 1,jk  ) - zti(ji  ,jj   &
286                             &                  - 1,jk  ))/ (zti(ji  ,jj + 1,jk  ) + zti(ji + 1   &
287                             &                  ,jj + 1,jk  ) + zti(ji + 1,jj - 1,jk  ) + zti(ji  &
288                             &                  ,jj - 1,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  ) +    &
289                             &                  zcc(ji + 1,jj  ,jk  ) + zcc(ji  ,jj  ,jk + 1) +   &
290                             &                  zcc(ji + 1,jj  ,jk + 1))* (zti(ji  ,jj  ,jk - 1)  &
291                             &                  + zti(ji + 1,jj  ,jk - 1) - zti(ji  ,jj  ,jk + 1  &
292                             &                  )- zti(ji + 1,jj  ,jk + 1))/ (zti(ji  ,jj  ,jk -  &
293                             &                  1) + zti(ji + 1,jj  ,jk - 1) + zti(ji  ,jj  ,jk   &
294                             &                  +1) + zti(ji + 1,jj  ,jk + 1) + rtrn))/(e1u(ji    &
295                             &                  ,jj)*e2u(ji,jj)*fse3u(ji,jj,jk))*vmask(ji  ,jj -  &
296                             &                  1,jk  )*vmask(ji + 1,jj - 1,jk  )*vmask(ji + 1    &
297                             &                  ,jj,jk)*vmask(ji  ,jj  ,jk  )*tmask(ji  ,jj  ,jk  &
298                             &                  )*tmask(ji + 1,jj  ,jk  )*tmask(ji  ,jj  ,jk + 1  &
299                             &                  )*tmask(ji + 1,jj  ,jk + 1)
300                        zy(ji,jj,jk) = zy(ji,jj,jk)                                               &   
301                             &                  - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,jk)*0.25*          &   
302                             &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji - 1,jj + 1   &   
303                             &                  ,jk  ) + zaa(ji  ,jj  ,jk  ) + zaa(ji  ,jj + 1    &   
304                             &                  ,jk))* (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1,jj   &   
305                             &                  ,jk  ) - zti(ji - 1,jj + 1,jk  ) - zti(ji - 1,jj  &   
306                             &                  ,jk  ))/ (zti(ji + 1,jj + 1,jk  ) + zti(ji + 1    &   
307                             &                  ,jj  ,jk  ) + zti(ji - 1,jj + 1,jk  ) + zti(ji    &   
308                             &                  - 1,jj  ,jk  ) + rtrn) + (zcc(ji  ,jj  ,jk  )     &   
309                             &                  + zcc(ji  ,jj  ,jk + 1) + zcc(ji  ,jj + 1,jk  )   &   
310                             &                  + zcc(ji  ,jj + 1,jk + 1))* (zti(ji  ,jj  ,jk -   &   
311                             &                  1) + zti(ji  ,jj + 1,jk - 1) - zti(ji  ,jj  ,jk   &   
312                             &                  +1) - zti(ji  ,jj + 1,jk + 1))/ (zti(ji  ,jj      &   
313                             &                  ,jk- 1) + zti(ji  ,jj + 1,jk - 1) + zti(ji  ,jj   &   
314                             &                  ,jk+ 1) + zti(ji  ,jj + 1,jk + 1) + rtrn))        &   
315                             &                  /(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk))          &   
316                             &                  *umask(ji - 1,jj,jk  )*umask(ji - 1,jj + 1,jk  )  &   
317                             &                  *umask(ji  ,jj,jk  )*umask(ji  ,jj + 1,jk  )      &   
318                             &                  *tmask(ji  ,jj,jk)*tmask(ji  ,jj  ,jk + 1)        &   
319                             &                  *tmask(ji  ,jj + 1,jk)*tmask(ji  ,jj + 1,jk + 1)     
320                        zz(ji,jj,jk) = zz(ji,jj,jk)                                               &   
321                             &                  - 0.5*rdttrc(jk)*rsc*zcc(ji,jj,jk)*0.25*          &   
322                             &                  (    (zaa(ji - 1,jj  ,jk  ) + zaa(ji  ,jj  ,jk    &   
323                             &                  ) + zaa(ji  ,jj  ,jk - 1) + zaa(ji - 1,jj  ,jk -  &   
324                             &                  1))*(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1,jj      &   
325                             &                  ,jk  ) - zti(ji - 1,jj  ,jk  ) - zti(ji - 1,jj    &   
326                             &                  ,jk - 1))/(zti(ji + 1,jj  ,jk - 1) + zti(ji + 1   &   
327                             &                  ,jj,jk  ) + zti(ji - 1,jj  ,jk  ) + zti(ji - 1    &   
328                             &                  ,jj,jk - 1) + rtrn) + (zbb(ji  ,jj - 1,jk  )      &   
329                             &                  + zbb(ji  ,jj  ,jk  ) + zbb(ji  ,jj  ,jk - 1)     &   
330                             &                  + zbb(ji  ,jj - 1,jk - 1))*(zti(ji  ,jj + 1,jk -  &   
331                             &                  1) + zti(ji  ,jj + 1,jk  ) - zti(ji  ,jj - 1,jk   &   
332                             &                  ) - zti(ji  ,jj - 1,jk - 1))/(zti(ji  ,jj + 1,jk  &   
333                             &                  - 1) + zti(ji  ,jj + 1,jk  ) + zti(ji  ,jj - 1    &   
334                             &                  ,jk  ) + zti(ji  ,jj - 1,jk - 1) + rtrn))         &   
335                             &                  /(e1t(ji,jj)*e2t(ji,jj)*fse3w(ji,jj,jk))          &   
336                             &                  *umask(ji - 1,jj,jk  )*umask(ji  ,jj  ,jk  )      &   
337                             &                  *umask(ji  ,jj,jk- 1)*umask(ji - 1,jj  ,jk - 1)   &   
338                             &                  *vmask(ji  ,jj- 1,jk)*vmask(ji  ,jj  ,jk  )       &   
339                             &                  *vmask(ji  ,jj  ,jk-1)*vmask(ji  ,jj - 1,jk - 1)       
340                     END DO
341                  END DO
342               END DO
343               
344               DO jj = 2,jpjm1
345                     DO ji = fs_2, fs_jpim1    ! Vector opt.
346                     zx(ji,jj,1) = zx(ji,jj,1)                                                    &   
347                          &                - 0.5*rdttrc(jk)*rsc*zaa(ji,jj,1)*0.25*                &   
348                          &                ( (zbb(ji  ,jj - 1,1  ) + zbb(ji + 1,jj - 1,1  )       &   
349                          &                + zbb(ji + 1,jj  ,1  ) + zbb(ji  ,jj  ,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                          &                /(zti(ji  ,jj + 1,1  ) + zti(ji + 1,jj + 1,1  )        &   
353                          &                + zti(ji + 1,jj - 1,1  ) + zti(ji  ,jj - 1,1  ) +      &   
354                          &                rtrn))/(e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,1))          &   
355                          &                *vmask(ji  ,jj - 1,1  )*vmask(ji + 1,jj - 1,1  )       &   
356                          &                *vmask(ji + 1,jj  ,1  )*vmask(ji  ,jj  ,1  )   
357                     zy(ji,jj,1) = zy(ji,jj,1)                                                    &   
358                          &                - 0.5*rdttrc(jk)*rsc*zbb(ji,jj,1)*0.25*                &   
359                          &                ( (zaa(ji-1  ,jj ,1  ) + zaa(ji - 1,jj + 1,1  )        &   
360                          &                + zaa(ji ,jj  ,1  ) + zaa(ji  ,jj + 1  ,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                          &                + zti(ji - 1,jj + 1,1  ) + zti(ji - 1,jj ,1  ) +       &   
365                          &                rtrn))/(e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,1))          &   
366                          &                *umask(ji - 1,jj,1  )*umask(ji - 1,jj + 1,1  )         &   
367                          &                *umask(ji    ,jj,1  )*umask(ji  ,jj + 1 ,1  )   
368                  END DO
369               END DO
370               !
371            ENDIF CROSSTERMS
372           
373            ! ... Lateral boundary conditions on z[xyz]
374            CALL lbc_lnk( zx, 'U', -1. )    ;    CALL lbc_lnk( zy, 'V', -1. )
375            CALL lbc_lnk( zz, 'W',  1. )
376
377            !--4.4 Reinitialization
378            zaa(:,:,:) = zx(:,:,:)
379            zbb(:,:,:) = zy(:,:,:)
380            zcc(:,:,:) = zz(:,:,:)
381
382            ! 5. Advection by antidiffusive mass fluxes & upstream scheme
383            ! -----------------------------------------------------------
384           
385            ! ... Horizontal
386            DO jk = 1, jpk
387               DO jj = 2, jpjm1
388                  DO ji = fs_2, fs_jpim1    ! Vector opt.
389                     zfp_ui = 0.5 * ( zaa(ji,jj,jk) + ABS( zaa(ji,jj,jk) ) )
390                     zfp_vj = 0.5 * ( zbb(ji,jj,jk) + ABS( zbb(ji,jj,jk) ) )
391                     zfm_ui = 0.5 * ( zaa(ji,jj,jk) - ABS( zaa(ji,jj,jk) ) )
392                     zfm_vj = 0.5 * ( zbb(ji,jj,jk) - ABS( zbb(ji,jj,jk) ) )           
393                     zkx(ji,jj,jk) = zfp_ui * zti(ji,jj,jk) + zfm_ui * zti(ji+1,jj  ,jk) 
394                     zky(ji,jj,jk) = zfp_vj * zti(ji,jj,jk) + zfm_vj * zti(ji  ,jj+1,jk)             
395                  END DO
396               END DO
397            END DO
398
399            ! ... Lateral boundary conditions on zk[xy]
400            CALL lbc_lnk( zkx, 'U', -1. )
401            CALL lbc_lnk( zky, 'V', -1. )
402           
403            ! ... Vertical
404            DO jk = 2, jpk
405               DO jj = 1, jpj
406                  DO ji = fs_2, fs_jpim1    ! Vector opt.
407                     zfp_w = 0.5 * ( zcc(ji,jj,jk) + ABS( zcc(ji,jj,jk) ) )
408                     zfm_w = 0.5 * ( zcc(ji,jj,jk) - ABS( zcc(ji,jj,jk) ) )       
409                     zkz(ji,jj,jk) = zfp_w * zti(ji,jj,jk) + zfm_w * zti(ji,jj,jk-1)     
410                  END DO
411               END DO
412            END DO
413
414            ! ... Compute ztj
415            DO jk = 1,jpkm1
416               DO jj = 2,jpjm1
417                  DO ji = fs_2, fs_jpim1    ! Vector opt.
418                     zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
419                     ztj(ji,jj,jk) = - zbtr *                                     & 
420                          &              ( zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )    & 
421                          &              + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )    & 
422                          &              + zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1) )
423#if defined key_trc_diatrd
424                     IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = trtrd(ji,jj,jk,ikeep(jn),1) -    &
425                          &                                        zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )
426
427                     IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = trtrd(ji,jj,jk,ikeep(jn),2) -    &
428                          &                                        zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
429
430                     IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = trtrd(ji,jj,jk,ikeep(jn),3) -    &
431                          &                                        zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
432#endif
433
434                  END DO
435               END DO
436            END DO
437
438            ! 6. Diagnose passive tracer trends (part 2/3)
439            ! --------------------------------------------
440            IF( l_trdtrc ) THEN
441               DO jk = 1, jpkm1
442                  DO jj = 2, jpjm1
443                     DO ji = fs_2, fs_jpim1    ! Vector opt.
444                        zbtr = 1./(e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk))
445                        ztrtrd(ji,jj,jk,1) = ztrtrd(ji,jj,jk,1) - zbtr*( zkx(ji,jj,jk) - zkx(ji - 1,jj,jk) )   
446                        ztrtrd(ji,jj,jk,2) = ztrtrd(ji,jj,jk,2) - zbtr*( zky(ji,jj,jk) - zky(ji,jj - 1,jk) )
447                        ztrtrd(ji,jj,jk,3) = ztrtrd(ji,jj,jk,3) - zbtr*( zkz(ji,jj,jk) - zkz(ji,jj,jk + 1) )
448                     END DO
449                  END DO
450               END DO
451            ENDIF
452            !                            ! ------------------------------------
453         END DO                          ! End of antidiffusive correction loop
454         !                               ! ------------------------------------
455
456         ! 7. Trend due to horizontal and vertical advection of tracer jn
457         ! --------------------------------------------------------------
458         
459         DO jk = 1, jpk
460            DO jj = 2, jpjm1
461               DO ji = fs_2, fs_jpim1     ! Vector opt.
462                  ztra = ( zbuf(ji,jj,jk) + ztj(ji,jj,jk) ) * tmask(ji,jj,jk)
463                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
464               END DO
465            END DO
466         END DO
467
468
469         ! 8. Convert the transport trend into advection trend (part 3/3)
470         ! --------------------------------------------------------------
471
472         IF( l_trdtrc ) THEN
473            ! ... Update the trends array
474            DO jk = 1, jpk
475               DO jj = 2, jpjm1
476                  DO  ji = fs_2, fs_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                     ztrtrd(ji,jj,jk,1) = ztrtrd(ji,jj,jk,1) + zgm
487                     ztrtrd(ji,jj,jk,2) = ztrtrd(ji,jj,jk,2) + zgz
488                     ztrtrd(ji,jj,jk,3) = ztrtrd(ji,jj,jk,3) - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
489                  END DO
490               END DO
491            END DO
492 
493            ! ... Lateral boundary conditions on trtrd:
494            CALL lbc_lnk( ztrtrd(:,:,:,1), 'T', 1. )
495            CALL lbc_lnk( ztrtrd(:,:,:,2), 'T', 1. )
496            CALL lbc_lnk( ztrtrd(:,:,:,3), 'T', 1. )
497 
498            ! ... Miscellaneous trends diagnostics
499            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd(:,:,:,1), jn, jptrc_trd_xad, kt )
500            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd(:,:,:,2), jn, jptrc_trd_yad, kt )
501            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd(:,:,:,3), jn, jptrc_trd_zad, kt )
502         ENDIF
503
504         !  Convert the transport trend into advection trend
505         ! ---------------------------------------------------
506
507#if defined key_trc_diatrd
508        DO jk = 1,jpk
509          DO jj = 2,jpjm1
510             DO ji = fs_2, fs_jpim1     ! Vector opt.
511                zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
512                zgm = zbtr * trn(ji,jj,jk,jn) *                                       &
513                     &    (  zun(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk)       &
514                     &     - zun(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk)    )
515
516                zgz = zbtr * trn(ji,jj,jk,jn) *                                       &
517                     &    (  zvn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk)       &
518                     &     - zvn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk)   )
519
520                IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = trtrd(ji,jj,jk,ikeep(jn),1) + zgm
521                IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = trtrd(ji,jj,jk,ikeep(jn),2) + zgz
522                IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = trtrd(ji,jj,jk,ikeep(jn),3)         &
523                     &                                      - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
524             END DO
525          END DO
526       END DO
527
528       ! Lateral boundary conditions on trtrd
529       IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),1), 'T', 1. )
530       IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),2), 'T', 1. )
531       IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),3), 'T', 1. )
532#endif
533
534       
535       !                                                 ! ==================
536    END DO                                               ! END of tracer loop
537    !                                                    ! ==================
538
539    IF( l_trdtrc ) DEALLOCATE( ztrtrd )
540 
541    IF( ln_ctl ) THEN      ! print mean trends (used for debugging)
542       WRITE(charout, FMT="('smolar - adv')")
543       CALL prt_ctl_trc_info(charout)
544       CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
545    ENDIF
546   
547  END SUBROUTINE trc_adv_smolar
548
549#else
550   !!----------------------------------------------------------------------
551   !!   Default option                                         Empty module
552   !!----------------------------------------------------------------------
553CONTAINS
554   SUBROUTINE trc_adv_smolar( kt ) 
555      INTEGER, INTENT(in) :: kt
556      WRITE(*,*) 'trc_adv_smolar: You should not have seen this print! error?', kt
557   END SUBROUTINE trc_adv_smolar
558#endif
559
560   !!======================================================================
561END MODULE trcadv_smolar
Note: See TracBrowser for help on using the repository browser.