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

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

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