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.
traadv_qck.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 888

Last change on this file since 888 was 888, checked in by ctlod, 16 years ago

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.7 KB
Line 
1MODULE traadv_qck
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_qck  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   tra_adv_qck : update the tracer trend with the horizontal
9   !!                      advection trends using a 3st order
10   !!                      finite difference scheme
11   !!                      The vertical advection scheme is the 2nd centered scheme
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and active tracers
15   USE dom_oce         ! ocean space and time domain
16   USE dynspg_oce      !
17   USE trdmod_oce      ! ocean variables trends
18   USE trdmod          ! ocean active tracers trends
19   USE trabbl          ! advective term in the BBL
20   USE lib_mpp
21   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
22   USE in_out_manager  ! I/O manager
23   USE diaptr          ! poleward transport diagnostics
24   USE prtctl          ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Accessibility
30   PUBLIC tra_adv_qck    ! routine called by step.F90
31
32   !! * Module variables
33   REAL(wp), DIMENSION(jpi,jpj),     SAVE ::   &
34         zbtr2
35   REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   &
36         sl
37   REAL(wp) ::                                 &
38         cst1, cst2, dt, coef1                  ! temporary scalars
39   INTEGER  ::                                 & 
40         ji, jj, jk                             ! dummy loop indices
41   !!----------------------------------------------------------------------
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   OPA 9.0 , LOCEAN-IPSL (2005)
47   !! $Id$
48   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_adv_qck  ***
56      !!
57      !! ** Purpose :   Compute the now trend due to the advection of tracers
58      !!      and add it to the general trend of passive tracer equations.
59      !!
60      !! ** Method :   The advection is evaluated by a third order scheme
61      !!               For a positive velocity u :
62      !!
63      !!
64      !!                  i-1    i      i+1   i+2
65      !!
66      !!               |--FU--|--FC--|--FD--|------|
67      !!                           u(i)>0
68      !!
69      !!               For a negative velocity u :
70      !!
71      !!               |------|--FD--|--FC--|--FU--|
72      !!                           u(i)<0
73      !!
74      !!                FU is the second upwind point
75      !!                FD is the first douwning point
76      !!                FC is the central point (or the first upwind point)
77      !!
78      !!      Flux(i) = u(i) * {0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i)å?)/6)(FU+FD-2FC)}
79      !!                with C(i)=|u(i)|dx(i)/dt (Courant number)
80      !!
81      !!         dt = 2*rdtra and the scalar values are tb and sb
82      !!
83      !!       On the vertical, the simple centered scheme used tn and sn
84      !!
85      !!               The fluxes are bounded by the ULTIMATE limiter
86      !!               to guarantee the monotonicity of the solution and to
87      !!            prevent the appearance of spurious numerical oscillations
88      !!
89      !! ** Action : - update (ta,sa) with the now advective tracer trends
90      !!             - save the trends in (ttrdh,strdh) ('key_trdtra')
91      !!
92      !! ** Reference : Leonard (1979, 1991)
93      !! History :
94      !!  9.0     !  06-09  (G. Reffray)  Original code
95      !!----------------------------------------------------------------------
96      !! * Arguments
97      INTEGER, INTENT( in ) ::   kt             ! ocean time-step index
98      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun   ! effective ocean velocity, u_component
99      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pvn   ! effective ocean velocity, v_component
100      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pwn   ! effective ocean velocity, w_component
101      !!
102      REAL(wp) :: z2                                          ! temporary scalar
103      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrdt, ztrds       ! temporary 3D workspace
104
105      IF( kt == nit000 ) THEN
106         IF(lwp) WRITE(numout,*)
107         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme'
108         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
109         IF(lwp) WRITE(numout,*)
110
111         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
112         cst1 = 1./12.
113         cst2 = 2./3.
114         IF (l_trdtra ) THEN
115         CALL ctl_warn( ' Trends not yet implemented for PPM advection scheme ' )
116         ENDIF
117      ENDIF
118
119      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
120      ELSE                                        ;    z2 = 2.
121      ENDIF
122
123      ! Save ta and sa trends
124      IF( l_trdtra )   THEN      ! to be done
125         ztrdt(:,:,:) = ta(:,:,:) 
126         ztrds(:,:,:) = sa(:,:,:) 
127         l_adv = 'qst'
128      ENDIF
129
130      ! I. Slope estimation at the T-point for the limiter ULTIMATE
131      ! SL = Sum(1/C_out) with C_out : Courant number for the outflows
132      !---------------------------------------------------------------
133
134      sl(:,:,:) = 100.
135
136                                                      ! ===============
137      DO jk = 1, jpkm1                                ! Horizontal slab
138        !                                             ! ===============
139        dt  = z2 * rdttra(jk)
140        DO jj = 2, jpjm1
141          DO ji = 2, fs_jpim1   ! vector opt.
142            coef1 = 1.e-12
143            IF (pun(ji-1,jj  ,jk  ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj  ,jk  ))*dt/e1t(ji,jj)
144            IF (pun(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pun(ji  ,jj  ,jk  ))*dt/e1t(ji,jj)
145            IF (pvn(ji  ,jj-1,jk  ).LT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj-1,jk  ))*dt/e2t(ji,jj)
146            IF (pvn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj  ,jk  ))*dt/e2t(ji,jj)
147            IF (pwn(ji  ,jj  ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk+1))*dt/(fse3t(ji,jj,jk))
148            IF (pwn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk  ))*dt/(fse3t(ji,jj,jk))
149            sl(ji,jj,jk) = 1./coef1
150            sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk))
151            sl(ji,jj,jk) = MAX(1.  ,sl(ji,jj,jk))
152          ENDDO
153        ENDDO
154      ENDDO
155      !--- Lateral boundary conditions on the limiter slope
156      CALL lbc_lnk(   sl(:,:,:), 'T', 1. )
157
158      ! II. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
159      !---------------------------------------------------------------------------
160
161      CALL tra_adv_qck_hor( kt , pun, pvn, tb , ta , pht_adv , z2)
162      CALL tra_adv_qck_hor( kt , pun, pvn, sb , sa , pst_adv , z2) 
163
164      ! Save the horizontal advective trends for diagnostic
165      ! ---------------------------------------------------
166!      IF( l_trdtra )   THEN    ! to be done
167!         ! T/S ZONAL advection trends
168!      ENDIF
169
170      IF(ln_ctl)   THEN
171          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had  - Ta: ', mask1=tmask, &
172             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
173      ENDIF
174
175      ! III. The vertical fluxes are computed with the 2nd order centered scheme
176      !-------------------------------------------------------------------------
177
178      CALL tra_adv_qck_ver( pwn, tn , ta, z2 )
179      CALL tra_adv_qck_ver( pwn, sn , sa, z2 )
180
181      ! Save the vertical advective trends for diagnostic
182      ! -------------------------------------------------
183!      IF( l_trdtra )   THEN    ! to be done
184         ! Recompute the vertical advection zta & zsa trends computed
185         ! at the step 2. above in making the difference between the new
186         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
187         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
188!     ENDIF
189
190      IF(ln_ctl)   THEN
191          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad  - Ta: ', mask1=tmask, &
192             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
193      ENDIF
194
195   END SUBROUTINE tra_adv_qck
196
197   SUBROUTINE tra_adv_qck_hor ( kt , pun, pvn, tra , traa , phtra_adv ,z2 )
198      !!----------------------------------------------------------------------
199      !!
200      !!----------------------------------------------------------------------
201      !! * Arguments
202      INTEGER, INTENT( in ) ::   kt             ! ocean time-step index
203      REAL, INTENT( in )    ::   z2
204      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun, pvn            ! horizontal effective velocity
205
206      REAL(wp), INTENT ( out   ), DIMENSION(jpj)          ::     &
207         phtra_adv
208
209      REAL(wp), INTENT ( inout ), DIMENSION(jpi,jpj,jpk)  ::     &
210         tra, traa
211
212      REAL(wp) ::                                &
213         za, zbtr, e1, e2, c, dir, fu, fc, fd,   & ! temporary scalars
214         coef2, coef3, fho, mask, dx
215
216      REAL(wp), DIMENSION(jpi,jpj) ::            &
217         zee
218
219      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
220         zmask, zlap, dwst, lim 
221
222
223
224      !----------------------------------------------------------------------
225      ! 0. Initialization (should ot be needed on the whole array ???)
226      !----------------------------------------------------------------------
227
228      zmask = 0.0                                               
229      zlap  = 0.0                                                 
230      dwst  = 0.0                                                 
231      lim   = 0.0     
232                                           
233      !----------------------------------------------------------------------
234      ! I. Part 1 : x-direction
235      !----------------------------------------------------------------------
236
237                                                       ! ===============
238      DO jk = 1, jpkm1                                 ! Horizontal slab
239         !                                             ! ===============
240         ! Initialization of metric arrays (for z- or s-coordinates)
241         ! ---------------------------------------------------------
242         DO jj = 1, jpjm1
243            DO ji = 1, fs_jpim1   ! vector opt.
244#if defined key_zco
245               ! z-coordinates, no vertical scale factors
246               zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
247#else
248               ! vertical scale factor are used
249               zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
250#endif
251            END DO
252         END DO
253
254         ! Laplacian of tracers (at before time step)
255         ! ------------------------------------------
256         !--- First derivative (gradient)
257         DO jj = 1, jpjm1
258            DO ji = 1, fs_jpim1   ! vector opt.
259               zmask(ji,jj,jk) = zee(ji,jj) * ( tra(ji+1,jj  ,jk) - tra(ji,jj,jk) )
260            END DO
261         END DO
262         DO jj = 2, jpjm1
263            DO ji = fs_2, fs_jpim1   ! vector opt.
264#if defined key_zco
265               zee(ji,jj) = e1t(ji,jj) /  e2t(ji,jj)
266#else
267               zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk))
268#endif
269               zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji-1,jj,jk)  )
270            END DO
271         END DO
272         !--- Function lim=FU+SL*(FC-FU) used by the limiter
273         !--- Computation of the ustream and downstream lim at the T-points
274         DO jj = 2, jpjm1
275            DO ji = 2, fs_jpim1   ! vector opt.
276               ! Upstream in the x-direction for the tracer
277               zmask(ji,jj,jk)=tra(ji-1,jj,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji-1,jj,jk))
278               ! Downstream in the x-direction for the tracer
279               dwst (ji,jj,jk)=tra(ji+1,jj,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji+1,jj,jk))
280            ENDDO
281         ENDDO
282      END DO
283      !--- Lateral boundary conditions on the laplacian (unchanged sgn)
284      CALL lbc_lnk(   zlap(:,:,:), 'T', 1. ) 
285      !--- Lateral boundary conditions for the lim function
286      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. )
287                                                       ! ===============
288      DO jk = 1, jpkm1                                 ! Horizontal slab
289         !                                             ! ===============
290         ! --- lim at the U-points in function of the direction of the flow
291         ! ----------------------------------------------------------------
292         DO jj = 1, jpjm1
293            DO ji = 1, fs_jpim1   ! vector opt.
294               dir = 0.5 + sign(0.5,pun(ji,jj,jk))     ! if pun>0 : diru = 1 otherwise diru = 0
295               lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk)
296               ! Mask at the T-points in the x-direction (mask=0 or mask=1)
297               zmask(ji,jj,jk)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
298            END DO
299         END DO
300      END DO
301      !--- Lateral boundary conditions for the mask
302      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. ) 
303
304      ! Horizontal advective fluxes
305      ! ---------------------------
306                                                       ! ===============
307      DO jk = 1, jpkm1                                 ! Horizontal slab
308                                                       ! ===============
309         dt  = z2 * rdttra(jk)
310         !--- tracer flux at u and v-points
311         DO jj = 1, jpjm1
312            DO ji = 1, fs_jpim1   ! vector opt.
313#if defined key_zco
314               e2 = e2u(ji,jj)
315#else
316               e2 = e2u(ji,jj) * fse3u(ji,jj,jk)
317#endif
318               dir = 0.5 + sign(0.5,pun(ji,jj,jk))                ! if pun>0 : diru = 1 otherwise diru = 0
319
320               dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
321               c  = ABS(pun(ji,jj,jk))*dt/dx                      ! (0<cx<1 : Courant number on x-direction)
322
323               fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the x-direction for T
324               fc  = dir*tra(ji  ,jj,jk)+(1-dir)*tra(ji+1,jj,jk)  ! FC in the x-direction for T
325               fd  = dir*tra(ji+1,jj,jk)+(1-dir)*tra(ji  ,jj,jk)  ! FD in the x-direction for T
326
327               !--- QUICKEST scheme
328               ! Temperature on the x-direction
329               coef1 = 0.5*(fc+fd)
330               coef2 = 0.5*c*(fd-fc)
331               coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji+1,jj,jk) )
332               fho = coef1-coef2-coef3
333               fho = bound(fu,fd,fc,fho)
334               !--- If the second ustream point is a land point
335               !--- the flux is computed by the 1st order UPWIND scheme
336               mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji+1,jj,jk)
337               fho = mask*fho + (1-mask)*fc
338               dwst(ji,jj,jk)=e2*pun(ji,jj,jk)*fho
339            END DO
340         END DO
341
342         !--- Tracer flux divergence at t-point added to the general trend
343         DO jj = 2, jpjm1
344            DO ji = fs_2, fs_jpim1   ! vector opt.
345               !--- horizontal advective trends
346#if defined key_zco
347               zbtr = zbtr2(ji,jj)
348#else
349               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)             
350#endif
351               za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj  ,jk) )
352               !--- add it to the general tracer trends
353               traa(ji,jj,jk) = traa(ji,jj,jk) + za
354            END DO
355         END DO
356         !                                             ! ===============
357      END DO                                           !   End of slab
358      !                                                ! ===============
359      !----------------------------------------------------------------------
360      ! I. Part 2 : y-direction
361      !----------------------------------------------------------------------
362                                                       ! ==============
363      DO jk = 1, jpkm1                                 ! Horizontal slab
364         !                                             ! ===============
365         ! Initialization of metric arrays (for z- or s-coordinates)
366         ! ---------------------------------------------------------
367         DO jj = 1, jpjm1
368            DO ji = 1, fs_jpim1   ! vector opt.
369#if defined key_zco
370               ! z-coordinates, no vertical scale factors
371               zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
372#else
373               ! s-coordinates, vertical scale factor are used
374               zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
375#endif
376            END DO
377         END DO
378
379         ! Laplacian of tracers (at before time step)
380         ! ------------------------------------------
381         !--- First derivative (gradient)
382         DO jj = 1, jpjm1
383            DO ji = 1, fs_jpim1   ! vector opt.
384               zmask(ji,jj,jk) = zee(ji,jj) * ( tra(ji  ,jj+1,jk) - tra(ji,jj,jk) )
385            END DO
386         END DO
387         !--- Second derivative (divergence)
388         DO jj = 2, jpjm1
389            DO ji = fs_2, fs_jpim1   ! vector opt.
390#if defined key_zco
391               zee(ji,jj) = e2t(ji,jj) /  e1t(ji,jj)
392#else
393               zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk))
394#endif
395               zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji,jj-1,jk)  )
396            END DO
397         END DO
398         !--- Function lim=FU+SL*(FC-FU) used by the limiter
399         !--- Computation of the ustream and downstream lim at the T-points
400         DO jj = 2, jpjm1
401            DO ji = 2, fs_jpim1   ! vector opt.
402               ! Upstream in the y-direction for the tracer
403               zmask(ji,jj,jk)=tra(ji,jj-1,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji,jj-1,jk))
404               ! Downstream in the y-direction for the tracer
405               dwst (ji,jj,jk)=tra(ji,jj+1,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji,jj+1,jk))
406            ENDDO
407         ENDDO
408      END DO
409      !--- Lateral boundary conditions on the laplacian (unchanged sgn)
410      CALL lbc_lnk(   zlap(:,:,:), 'T', 1. )
411      !--- Lateral boundary conditions for the lim function
412      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. )
413
414      DO jk = 1, jpkm1                                 ! Horizontal slab
415         !                                             ! ===============
416         ! --- lim at the V-points in function of the direction of the flow
417         ! ----------------------------------------------------------------
418         DO jj = 1, jpjm1
419            DO ji = 1, fs_jpim1   ! vector opt.
420               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))      ! if pvn>0 : dirv = 1 otherwise dirv = 0
421               lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk)
422               ! Mask at the T-points in the y-direction (mask=0 or mask=1)
423               zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
424            END DO
425         END DO
426      END DO
427      !--- Lateral boundary conditions for the mask
428      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. ) 
429
430      ! Horizontal advective fluxes
431      ! -------------------------------
432                                                       ! ===============
433      DO jk = 1, jpkm1                                 ! Horizontal slab
434                                                       ! ===============
435         dt  = z2 * rdttra(jk)
436         !--- tracer flux at u and v-points
437         DO jj = 1, jpjm1
438            DO ji = 1, fs_jpim1   ! vector opt.
439#if defined key_zco
440               e1 = e1v(ji,jj)
441#else
442               e1 = e1v(ji,jj) * fse3v(ji,jj,jk)
443#endif
444               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))                ! if pvn>0 : dirv = 1 otherwise dirv = 0
445
446               dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
447               c  = ABS(pvn(ji,jj,jk))*dt/dx                      ! (0<cy<1 : Courant number on y-direction)
448
449               fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the y-direction for T
450               fc  = dir*tra(ji,jj  ,jk)+(1-dir)*tra(ji,jj+1,jk)  ! FC in the y-direction for T
451               fd  = dir*tra(ji,jj+1,jk)+(1-dir)*tra(ji,jj  ,jk)  ! FD in the y-direction for T
452
453               !--- QUICKEST scheme
454               ! Temperature on the y-direction
455               coef1 = 0.5*(fc+fd)
456               coef2 = 0.5*c*(fd-fc)
457               coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji,jj+1,jk) )
458               fho = coef1-coef2-coef3
459               fho = bound(fu,fd,fc,fho)
460               !--- If the second ustream point is a land point
461               !--- the flux is computed by the 1st order UPWIND scheme
462               mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji,jj+1,jk)
463               fho = mask*fho + (1-mask)*fc
464               dwst(ji,jj,jk)=e1*pvn(ji,jj,jk)*fho
465            END DO
466         END DO
467
468         !--- Tracer flux divergence at t-point added to the general trend
469         DO jj = 2, jpjm1
470            DO ji = fs_2, fs_jpim1   ! vector opt.
471               !--- horizontal advective trends
472#if defined key_zco
473               zbtr = zbtr2(ji,jj)
474#else
475               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
476#endif
477               za = - zbtr * (  dwst(ji,jj,jk) - dwst(ji  ,jj-1,jk) )
478               !--- add it to the general tracer trends
479               traa(ji,jj,jk) = traa(ji,jj,jk) + za
480            END DO
481         END DO
482         !                                             ! ===============
483      END DO                                           !   End of slab
484      !                                                ! ===============
485
486      ! "zonal" mean advective heat and salt transport
487      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
488#if defined key_zco
489         DO jk = 1, jpkm1
490            DO jj = 2, jpjm1
491               DO ji = fs_2, fs_jpim1   ! vector opt.
492                  dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk)
493               END DO
494            END DO
495         END DO
496         phtra_adv(:) = ptr_vj( dwst(:,:,:) )
497#else
498         phtra_adv(:) = ptr_vj( dwst(:,:,:) )
499# endif
500      ENDIF
501
502   END SUBROUTINE tra_adv_qck_hor
503
504   SUBROUTINE tra_adv_qck_ver ( pwn, tra , traa, z2  )     
505      !!----------------------------------------------------------------------
506      !!
507      !!----------------------------------------------------------------------
508      !! * Arguments
509
510      REAL(wp), INTENT ( in ) :: z2
511      REAL(wp), INTENT ( in ), DIMENSION(jpi,jpj,jpk)  ::     &
512         pwn
513      REAL(wp), INTENT ( inout ), DIMENSION(jpi,jpj,jpk)  ::     &
514         tra, traa
515
516      REAL(wp) ::                             &
517         za, ze3tr, dt, dir, fc, fd             ! temporary scalars
518
519      ! Vertical advection
520      ! ------------------
521
522      ! 1. Vertical advective fluxes
523      ! ----------------------------
524
525      !Bottom value : flux set to zero
526      sl(:,:,jpk) = 0.e0
527
528      ! Surface value
529      IF( lk_dynspg_rl .OR. lk_vvl ) THEN
530         ! rigid lid : flux set to zero
531         sl(:,:, 1 ) = 0.e0
532      ELSE
533         ! free surface-constant volume
534         sl(:,:, 1 ) = pwn(:,:,1) * tra(:,:,1)
535      ENDIF
536
537      ! Second order centered tracer flux at w-point
538
539      DO jk = 2, jpkm1
540         dt  = z2 *  rdttra(jk)
541         DO jj = 2, jpjm1
542            DO ji = fs_2, fs_jpim1   ! vector opt.
543               dir = 0.5 + sign(0.5,pwn(ji,jj,jk))                                         ! if pwn>0 : dirw = 1 otherwise dirw = 0
544               fc = dir*tra(ji,jj,jk  )*fse3t(ji,jj,jk-1)+(1-dir)*tra(ji,jj,jk-1)*fse3t(ji,jj,jk  )   ! FC in the z-direction for T
545               fd = dir*tra(ji,jj,jk-1)*fse3t(ji,jj,jk  )+(1-dir)*tra(ji,jj,jk  )*fse3t(ji,jj,jk-1)   ! FD in the z-direction for T
546               !--- Second order centered scheme
547               sl(ji,jj,jk)=pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk))
548            END DO
549         END DO
550      END DO
551
552      ! 2. Tracer flux divergence at t-point added to the general trend
553      ! ---------------------------------------------------------------
554
555      DO jk = 1, jpkm1
556         DO jj = 2, jpjm1
557            DO ji = fs_2, fs_jpim1   ! vector opt.
558               ze3tr = 1. / fse3t(ji,jj,jk)
559               ! vertical advective trends
560               za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) )
561               ! add it to the general tracer trends
562               traa(ji,jj,jk) =  traa(ji,jj,jk) + za
563            END DO
564         END DO
565      END DO
566
567   END SUBROUTINE tra_adv_qck_ver
568
569   REAL FUNCTION bound(fu,fd,fc,fho)
570      real     ::  fu,fd,fc,fho,fref1,fref2
571      fref1 = fu
572      fref2 = MAX(MIN(fc,fd),MIN(MAX(fc,fd),fref1))
573      bound = MAX(MIN(fho,fc),MIN(MAX(fho,fc),fref2))
574   END FUNCTION
575
576   !!======================================================================
577END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.