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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3299

Last change on this file since 3299 was 3299, checked in by cbricaud, 12 years ago

fix bug for ticket 924

  • Property svn:keywords set to Id
File size: 23.8 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  : apply QUICK scheme in i-direction
14   !!   tra_adv_qck_j  : apply QUICK scheme in j-direction
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   USE wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_adv_qck   ! routine called by step.F90
35
36   LOGICAL  :: l_trd           ! flag to compute trends
37   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
50      &                                       ptb, ptn, pta, kjpt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_adv_qck  ***
53      !!
54      !! ** Purpose :   Compute the now trend due to the advection of tracers
55      !!      and add it to the general trend of passive tracer equations.
56      !!
57      !! ** Method :   The advection is evaluated by a third order scheme
58      !!             For a positive velocity u :              u(i)>0
59      !!                                          |--FU--|--FC--|--FD--|------|
60      !!                                             i-1    i      i+1   i+2
61      !!
62      !!             For a negative velocity u :              u(i)<0
63      !!                                          |------|--FD--|--FC--|--FU--|
64      !!                                             i-1    i      i+1   i+2
65      !!             where  FU is the second upwind point
66      !!                    FD is the first douwning point
67      !!                    FC is the central point (or the first upwind point)
68      !!
69      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
70      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
71      !!
72      !!         dt = 2*rdtra and the scalar values are tb and sb
73      !!
74      !!       On the vertical, the simple centered scheme used ptn
75      !!
76      !!               The fluxes are bounded by the ULTIMATE limiter to
77      !!             guarantee the monotonicity of the solution and to
78      !!            prevent the appearance of spurious numerical oscillations
79      !!
80      !! ** Action : - update (pta) with the now advective tracer trends
81      !!             - save the trends
82      !!
83      !! ** Reference : Leonard (1979, 1991)
84      !!----------------------------------------------------------------------
85      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
86      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
87      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
88      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
89      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
90      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
91      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
92      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
93      !!----------------------------------------------------------------------
94
95      !
96      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
97      !
98      IF( kt == kit000 )  THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
102         IF(lwp) WRITE(numout,*)
103         !
104         l_trd = .FALSE.
105         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
106      ENDIF
107
108      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
109      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
110      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
111
112      ! II. The vertical fluxes are computed with the 2nd order centered scheme
113      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
114      !
115      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
116      !
117   END SUBROUTINE tra_adv_qck
118
119
120   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
121      &                                        ptb, ptn, pta, kjpt   )
122      !!----------------------------------------------------------------------
123      !!
124      !!----------------------------------------------------------------------
125      USE oce     , ONLY:   zwx => ua       ! ua used as workspace
126      !
127      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
128      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
129      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
130      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
131      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
132      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
133      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
134      !!
135      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
136      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
137      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
138      !----------------------------------------------------------------------
139      !
140      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
141      !                                                          ! ===========
142      DO jn = 1, kjpt                                            ! tracer loop
143         !                                                       ! ===========
144         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
145         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
146         !                                                 
147         DO jk = 1, jpkm1                               
148            !                                             
149            !--- Computation of the ustream and downstream value of the tracer and the mask
150            DO jj = 2, jpjm1
151               DO ji = fs_2, fs_jpim1   ! vector opt.
152                  ! Upstream in the x-direction for the tracer
153                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
154                  ! Downstream in the x-direction for the tracer
155                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
156               END DO
157            END DO
158         END DO
159         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
160         
161         !
162         ! Horizontal advective fluxes
163         ! ---------------------------
164         !
165         DO jk = 1, jpkm1                             
166            DO jj = 2, jpjm1
167               DO ji = fs_2, fs_jpim1   ! vector opt.         
168                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
169                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
170               END DO
171            END DO
172         END DO
173         !
174         DO jk = 1, jpkm1 
175            zdt =  p2dt(jk)
176            DO jj = 2, jpjm1
177               DO ji = fs_2, fs_jpim1   ! vector opt.   
178                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
179                  zdx  = zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj)
180                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
181                  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
182                  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
183               END DO
184            END DO
185         END DO 
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               END DO
199            END DO
200         END DO
201         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
202
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      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
236      !
237   END SUBROUTINE tra_adv_qck_i
238
239
240   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
241      &                                        ptb, ptn, pta, kjpt )
242      !!----------------------------------------------------------------------
243      !!
244      !!----------------------------------------------------------------------
245      USE oce     , ONLY:   zwy => ua       ! ua used as workspace
246      !
247      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
248      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
249      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
250      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
251      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
252      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
253      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
254      !!
255      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
256      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
257      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
258      !----------------------------------------------------------------------
259      !
260      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
261      !
262      !                                                          ! ===========
263      DO jn = 1, kjpt                                            ! tracer loop
264         !                                                       ! ===========
265         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
266         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
267         !                                                 
268         DO jk = 1, jpkm1                               
269            !                                             
270            !--- Computation of the ustream and downstream value of the tracer and the mask
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1   ! vector opt.
273                  ! Upstream in the x-direction for the tracer
274                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
275                  ! Downstream in the x-direction for the tracer
276                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
277               END DO
278            END DO
279         END DO
280         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
281
282         
283         !
284         ! Horizontal advective fluxes
285         ! ---------------------------
286         !
287         DO jk = 1, jpkm1                             
288            DO jj = 2, jpjm1
289               DO ji = fs_2, fs_jpim1   ! vector opt.         
290                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
291                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
292               END DO
293            END DO
294         END DO
295         !
296         DO jk = 1, jpkm1 
297            zdt =  p2dt(jk)
298            DO jj = 2, jpjm1
299               DO ji = fs_2, fs_jpim1   ! vector opt.   
300                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
301                  zdx  = zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1)
302                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
303                  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
304                  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
305               END DO
306            END DO
307         END DO
308
309         !--- Lateral boundary conditions
310         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
311         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
312
313         !--- QUICKEST scheme
314         CALL quickest( zfu, zfd, zfc, zwy )
315         !
316         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
317         DO jk = 1, jpkm1 
318            DO jj = 2, jpjm1
319               DO ji = fs_2, fs_jpim1   ! vector opt.               
320                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
321               END DO
322            END DO
323         END DO
324         !--- Lateral boundary conditions
325         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
326         !
327         ! Tracer flux on the x-direction
328         DO jk = 1, jpkm1 
329            !
330            DO jj = 2, jpjm1
331               DO ji = fs_2, fs_jpim1   ! vector opt.               
332                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
333                  !--- If the second ustream point is a land point
334                  !--- the flux is computed by the 1st order UPWIND scheme
335                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
336                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
337                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
338               END DO
339            END DO
340            !
341            ! Computation of the trend
342            DO jj = 2, jpjm1
343               DO ji = fs_2, fs_jpim1   ! vector opt. 
344                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
345                  ! horizontal advective trends
346                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
347                  !--- add it to the general tracer trends
348                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
349               END DO
350            END DO
351            !
352         END DO
353         !                                 ! trend diagnostics (contribution of upstream fluxes)
354         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
355         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
356         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
357           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
358           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
359         ENDIF
360         !
361      END DO
362      !
363      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
364      !
365   END SUBROUTINE tra_adv_qck_j
366
367
368   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
369     &                                    ptn, pta, kjpt )
370      !!----------------------------------------------------------------------
371      !!
372      !!----------------------------------------------------------------------
373      USE oce, ONLY:   zwz => ua   ! ua used as workspace
374      !
375      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
376      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
377      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
378      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
379      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
380      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
381      !
382      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
383      REAL(wp) ::   zbtr , ztra      ! local scalars
384      !!----------------------------------------------------------------------
385
386      !                                                          ! ===========
387      DO jn = 1, kjpt                                            ! tracer loop
388         !                                                       ! ===========
389         ! 1. Bottom value : flux set to zero
390         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
391         !
392         !                                 ! Surface value
393         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
394         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
395         ENDIF
396         !
397         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
398            DO jj = 2, jpjm1
399               DO ji = fs_2, fs_jpim1   ! vector opt.
400                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
401               END DO
402            END DO
403         END DO
404         !
405         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
406            DO jj = 2, jpjm1
407               DO ji = fs_2, fs_jpim1   ! vector opt.
408                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
409                  ! k- vertical advective trends
410                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
411                  ! added to the general tracer trends
412                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
413               END DO
414            END DO
415         END DO
416         !                                 ! Save the vertical advective trends for diagnostic
417         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
418         !
419      END DO
420      !
421   END SUBROUTINE tra_adv_cen2_k
422
423
424   SUBROUTINE quickest( pfu, pfd, pfc, puc )
425      !!----------------------------------------------------------------------
426      !!
427      !! ** Purpose :  Computation of advective flux with Quickest scheme
428      !!
429      !! ** Method :   
430      !!----------------------------------------------------------------------
431      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
432      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
433      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
434      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
435      !!
436      INTEGER  ::  ji, jj, jk               ! dummy loop indices
437      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
438      REAL(wp) ::  zc, zcurv, zfho          !   -      -
439      !----------------------------------------------------------------------
440      !
441      IF( nn_timing == 1 )  CALL timing_start('quickest')
442      !
443      DO jk = 1, jpkm1
444         DO jj = 1, jpj
445            DO ji = 1, jpi
446               zc     = puc(ji,jj,jk)                         ! Courant number
447               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
448               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
449               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
450               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
451               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
452               !
453               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
454               zcoef2 = ABS( zcoef1 )
455               zcoef3 = ABS( zcurv )
456               IF( zcoef3 >= zcoef2 ) THEN
457                  zfho = pfc(ji,jj,jk) 
458               ELSE
459                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
460                  IF( zcoef1 >= 0. ) THEN
461                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
462                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
463                  ELSE
464                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
465                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
466                  ENDIF
467               ENDIF
468               puc(ji,jj,jk) = zfho
469            END DO
470         END DO
471      END DO
472      !
473      IF( nn_timing == 1 )  CALL timing_stop('quickest')
474      !
475   END SUBROUTINE quickest
476
477   !!======================================================================
478END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.