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 branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA – NEMO

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2083

Last change on this file since 2083 was 2083, checked in by cetlod, 14 years ago

correction of minor bugs
File size: 23.6 KB
Line 
1MODULE traadv_qck
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_qck  ***
4   !! Ocean tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  3.0  !  2008-07  (G. Reffray)  Original code
7   !!            3.3  !  2010-05  (C.Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_adv_qck      : update the tracer trend with the horizontal advection
12   !!                      trends using a 3rd order finite difference scheme
13   !!   tra_adv_qck_i  :
14   !!   tra_adv_qck_j  :
15   !!   tra_adv_cen2_k : 2nd centered scheme for the vertical advection
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and active tracers
18   USE dom_oce         ! ocean space and time domain
19   USE trdmod_oce         ! ocean space and time domain
20   USE trdtra      ! ocean tracers trends
21   USE trabbl          ! advective term in the BBL
22   USE lib_mpp         ! distribued memory computing
23   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
24   USE dynspg_oce      ! surface pressure gradient variables
25   USE in_out_manager  ! I/O manager
26   USE diaptr          ! poleward transport diagnostics
27   USE trc_oce         ! share passive tracers/Ocean variables
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_qck   ! routine called by step.F90
33
34   REAL(wp)  :: r1_6 = 1./ 6.
35   LOGICAL   :: l_trd    ! flag to compute trends
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
42   !! $Id: traadv_qck.F90 2024 2010-07-29 10:57:35Z cetlod $
43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn, &
49      &                                       ptb, ptn, pta, kjpt   )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_adv_qck  ***
52      !!
53      !! ** Purpose :   Compute the now trend due to the advection of tracers
54      !!      and add it to the general trend of passive tracer equations.
55      !!
56      !! ** Method :   The advection is evaluated by a third order scheme
57      !!             For a positive velocity u :              u(i)>0
58      !!                                          |--FU--|--FC--|--FD--|------|
59      !!                                             i-1    i      i+1   i+2
60      !!
61      !!             For a negative velocity u :              u(i)<0
62      !!                                          |------|--FD--|--FC--|--FU--|
63      !!                                             i-1    i      i+1   i+2
64      !!             where  FU is the second upwind point
65      !!                    FD is the first douwning point
66      !!                    FC is the central point (or the first upwind point)
67      !!
68      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
69      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
70      !!
71      !!         dt = 2*rdtra and the scalar values are tb and sb
72      !!
73      !!       On the vertical, the simple centered scheme used ptn
74      !!
75      !!               The fluxes are bounded by the ULTIMATE limiter to
76      !!             guarantee the monotonicity of the solution and to
77      !!            prevent the appearance of spurious numerical oscillations
78      !!
79      !! ** Action : - update (pta) with the now advective tracer trends
80      !!             - save the trends
81      !!
82      !! ** Reference : Leonard (1979, 1991)
83      !!----------------------------------------------------------------------
84      !!
85      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
86      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
87      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
88      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
89      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
90      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
91      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
92      !!----------------------------------------------------------------------
93
94      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
97         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
98         IF(lwp) WRITE(numout,*)
99         !
100         l_trd = .FALSE.
101         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
102      ENDIF
103
104      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
105      !---------------------------------------------------------------------------
106
107      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
108      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
109
110      ! II. The vertical fluxes are computed with the 2nd order centered scheme
111      !-------------------------------------------------------------------------
112      !
113      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
114      !
115   END SUBROUTINE tra_adv_qck
116
117   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,    &
118      &                                        ptb, ptn, pta, kjpt   )
119      !!----------------------------------------------------------------------
120      !!
121      !!----------------------------------------------------------------------
122      USE oce         , zwx => ua   ! use ua as workspace
123      !!
124      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
125      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
126      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
127      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
128      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun             ! zonal velocity component
129      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
130      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
131      !!
132      INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
133      REAL(wp) :: ztra, zbtr               ! temporary scalars
134      REAL(wp) :: zdir, zdx, zdt, zmsk     ! temporary scalars
135      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfu, zfc, zfd
136      !----------------------------------------------------------------------
137
138     
139      DO jn = 1, kjpt                                            ! tracer loop
140         !                                                       ! ===========
141         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
142         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
143         !                                                 
144         DO jk = 1, jpkm1                               
145            !                                             
146            !--- Computation of the ustream and downstream value of the tracer and the mask
147            DO jj = 2, jpjm1
148               DO ji = fs_2, fs_jpim1   ! vector opt.
149                  ! Upstream in the x-direction for the tracer
150                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
151                  ! Downstream in the x-direction for the tracer
152                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
153               END DO
154            END DO
155         END DO
156         !
157         !--- Lateral boundary conditions
158         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 
159         
160         !
161         ! Horizontal advective fluxes
162         ! ---------------------------
163         !
164         DO jk = 1, jpkm1                             
165            DO jj = 2, jpjm1
166               DO ji = fs_2, fs_jpim1   ! vector opt.         
167                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
168                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
169               END DO
170            END DO
171         END DO
172         !
173         DO jk = 1, jpkm1 
174            zdt =  p2dt(jk)
175            DO jj = 2, jpjm1
176               DO ji = fs_2, fs_jpim1   ! vector opt.   
177                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
178                  zdx  = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
179                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
180                  zfc(ji,jj,jk)  = zdir * ptb(ji  ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn)  ! FC in the x-direction for T
181                  zfd(ji,jj,jk)  = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji  ,jj,jk,jn)  ! FD in the x-direction for T
182               END DO
183            END DO
184         END DO      !
185
186         !--- Lateral boundary conditions
187         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
188         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
189
190         !--- QUICKEST scheme
191         CALL quickest( zfu, zfd, zfc, zwx )
192         !
193         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
194         DO jk = 1, jpkm1 
195            DO jj = 2, jpjm1
196               DO ji = fs_2, fs_jpim1   ! vector opt.               
197                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
198               ENDDO
199            END DO
200         END DO
201         !--- Lateral boundary conditions
202         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
203         !
204         ! Tracer flux on the x-direction
205         DO jk = 1, jpkm1 
206            !
207            DO jj = 2, jpjm1
208               DO ji = fs_2, fs_jpim1   ! vector opt.               
209                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
210                  !--- If the second ustream point is a land point
211                  !--- the flux is computed by the 1st order UPWIND scheme
212                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
213                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
214                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
215               END DO
216            END DO
217            !
218            ! Computation of the trend
219            DO jj = 2, jpjm1
220               DO ji = fs_2, fs_jpim1   ! vector opt. 
221                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
222                  ! horizontal advective trends
223                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
224                  !--- add it to the general tracer trends
225                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
226               END DO
227            END DO
228            !
229         END DO
230         !                                 ! trend diagnostics (contribution of upstream fluxes)
231         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
232         !
233      END DO
234      !
235   END SUBROUTINE tra_adv_qck_i
236
237   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,   &
238      &                                        ptb, ptn, pta, kjpt   )
239      !!----------------------------------------------------------------------
240      !!
241      !!----------------------------------------------------------------------
242      !!
243      USE oce         , zwy => ua   ! use ua as workspace
244      !!
245      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
246      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
247      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
248      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
249      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pvn             ! meridional velocity component
250      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
251      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
252      !!
253      INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
254      REAL(wp) :: ztra, zbtr               ! temporary scalars
255      REAL(wp) :: zdir, zdx, zdt, zmsk     ! temporary scalars
256      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfu, zfc, zfd
257      !----------------------------------------------------------------------
258
259      DO jn = 1, kjpt                                            ! tracer loop
260         !                                                       ! ===========
261         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
262         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
263         !                                                 
264         DO jk = 1, jpkm1                               
265            !                                             
266            !--- Computation of the ustream and downstream value of the tracer and the mask
267            DO jj = 2, jpjm1
268               DO ji = fs_2, fs_jpim1   ! vector opt.
269                  ! Upstream in the x-direction for the tracer
270                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
271                  ! Downstream in the x-direction for the tracer
272                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
273               END DO
274            END DO
275         END DO
276         !
277         !--- Lateral boundary conditions
278         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 
279         
280         !
281         ! Horizontal advective fluxes
282         ! ---------------------------
283         !
284         DO jk = 1, jpkm1                             
285            DO jj = 2, jpjm1
286               DO ji = fs_2, fs_jpim1   ! vector opt.         
287                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
288                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
289               END DO
290            END DO
291         END DO
292         !
293         DO jk = 1, jpkm1 
294            zdt =  p2dt(jk)
295            DO jj = 2, jpjm1
296               DO ji = fs_2, fs_jpim1   ! vector opt.   
297                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
298                  zdx  = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
299                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
300                  zfc(ji,jj,jk)  = zdir * ptb(ji,jj  ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn)  ! FC in the x-direction for T
301                  zfd(ji,jj,jk)  = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj  ,jk,jn)  ! FD in the x-direction for T
302               END DO
303            END DO
304         END DO      !
305
306         !--- Lateral boundary conditions
307         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
308         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
309
310         !--- QUICKEST scheme
311         CALL quickest( zfu, zfd, zfc, zwy )
312         !
313         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
314         DO jk = 1, jpkm1 
315            DO jj = 2, jpjm1
316               DO ji = fs_2, fs_jpim1   ! vector opt.               
317                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
318               ENDDO
319            END DO
320         END DO
321         !--- Lateral boundary conditions
322         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
323         !
324         ! Tracer flux on the x-direction
325         DO jk = 1, jpkm1 
326            !
327            DO jj = 2, jpjm1
328               DO ji = fs_2, fs_jpim1   ! vector opt.               
329                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
330                  !--- If the second ustream point is a land point
331                  !--- the flux is computed by the 1st order UPWIND scheme
332                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
333                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
334                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
335               END DO
336            END DO
337            !
338            ! Computation of the trend
339            DO jj = 2, jpjm1
340               DO ji = fs_2, fs_jpim1   ! vector opt. 
341                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
342                  ! horizontal advective trends
343                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
344                  !--- add it to the general tracer trends
345                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
346               END DO
347            END DO
348            !
349         END DO
350         !                                 ! trend diagnostics (contribution of upstream fluxes)
351         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
352         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
353         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
354           IF( jn == jp_tem )  pht_adv(:) = ptr_vj( zwy(:,:,:) )
355           IF( jn == jp_sal )  pst_adv(:) = ptr_vj( zwy(:,:,:) )
356         ENDIF
357         !
358      END DO
359
360   END SUBROUTINE tra_adv_qck_j
361
362   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,   &
363     &                                    ptn, pta, kjpt )
364      !!----------------------------------------------------------------------
365      !!
366      !!----------------------------------------------------------------------
367      !!
368      USE oce         , zwz => ua   ! use ua as workspace
369      !!
370      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
371      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
372      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
373      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pwn             ! vertical velocity component
374      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn           ! now tracer field
375      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
376      !!
377      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
378      REAL(wp) ::   zbtr , ztra      ! temporary scalars
379      !!----------------------------------------------------------------------
380
381      !
382      DO jn = 1, kjpt                                            ! tracer loop
383         !                                                       ! ===========
384         ! 1. Bottom value : flux set to zero
385         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
386         !
387         !                                 ! Surface value
388         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
389         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
390         ENDIF
391         !
392         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
393            DO jj = 2, jpjm1
394               DO ji = fs_2, fs_jpim1   ! vector opt.
395                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
396               END DO
397            END DO
398         END DO
399         !
400         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
401            DO jj = 2, jpjm1
402               DO ji = fs_2, fs_jpim1   ! vector opt.
403                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
404                  ! k- vertical advective trends
405                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
406                  ! added to the general tracer trends
407                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
408               END DO
409            END DO
410         END DO
411         !                                 ! Save the vertical advective trends for diagnostic
412         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
413         !
414      END DO
415      !
416   END SUBROUTINE tra_adv_cen2_k
417
418
419   SUBROUTINE quickest( pfu, pfd, pfc, puc )
420      !!----------------------------------------------------------------------
421      !!
422      !! ** Purpose :  Computation of advective flux with Quickest scheme
423      !!
424      !! ** Method :   
425      !!----------------------------------------------------------------------
426      REAL(wp), INTENT(in)    , DIMENSION(jpi,jpj,jpk) :: pfu   ! second upwind point
427      REAL(wp), INTENT(in)    , DIMENSION(jpi,jpj,jpk) :: pfd   ! first douwning point
428      REAL(wp), INTENT(in)    , DIMENSION(jpi,jpj,jpk) :: pfc   ! the central point (or the first upwind point)
429      REAL(wp), INTENT(inout) , DIMENSION(jpi,jpj,jpk) :: puc   ! input as Courant number ; output as flux
430      !!
431      INTEGER  ::  ji, jj, jk               ! dummy loop indices
432      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! temporary scalars         
433      REAL(wp) ::  zc, zcurv, zfho          !
434      !----------------------------------------------------------------------
435
436      DO jk = 1, jpkm1
437         DO jj = 1, jpj
438            DO ji = 1, jpi
439               zc     = puc(ji,jj,jk)                         ! Courant number
440               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
441               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
442               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
443               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
444               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
445               !
446               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
447               zcoef2 = ABS( zcoef1 )
448               zcoef3 = ABS( zcurv )
449               IF( zcoef3 >= zcoef2 ) THEN
450                  zfho = pfc(ji,jj,jk) 
451               ELSE
452                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
453                  IF( zcoef1 >= 0. ) THEN
454                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
455                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
456                  ELSE
457                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
458                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
459                  ENDIF
460               ENDIF
461               puc(ji,jj,jk) = zfho
462            ENDDO
463         ENDDO
464      ENDDO
465      !
466   END SUBROUTINE quickest
467
468   !!======================================================================
469END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.