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

Last change on this file since 21 was 21, checked in by cholod, 12 years ago

Update of the trunk with the v3.4_r3306
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   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: traadv_qck.F90 3301 2012-02-08 16:56:27Z cbricaud $
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) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
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         END DO
218         !
219         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
220         !
221         ! Computation of the trend
222         DO jk = 1, jpkm1 
223            DO jj = 2, jpjm1
224               DO ji = fs_2, fs_jpim1   ! vector opt. 
225                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
226                  ! horizontal advective trends
227                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
228                  !--- add it to the general tracer trends
229                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
230               END DO
231            END DO
232         END DO
233         !                                 ! trend diagnostics (contribution of upstream fluxes)
234         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
235         !
236      END DO
237      !
238      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
239      !
240   END SUBROUTINE tra_adv_qck_i
241
242
243   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
244      &                                        ptb, ptn, pta, kjpt )
245      !!----------------------------------------------------------------------
246      !!
247      !!----------------------------------------------------------------------
248      USE oce     , ONLY:   zwy => ua       ! ua used as workspace
249      !
250      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
251      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
252      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
253      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
254      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
255      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
256      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
257      !!
258      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
259      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
260      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
261      !----------------------------------------------------------------------
262      !
263      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
264      !
265      !                                                          ! ===========
266      DO jn = 1, kjpt                                            ! tracer loop
267         !                                                       ! ===========
268         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
269         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
270         !                                                 
271         DO jk = 1, jpkm1                               
272            !                                             
273            !--- Computation of the ustream and downstream value of the tracer and the mask
274            DO jj = 2, jpjm1
275               DO ji = fs_2, fs_jpim1   ! vector opt.
276                  ! Upstream in the x-direction for the tracer
277                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
278                  ! Downstream in the x-direction for the tracer
279                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
280               END DO
281            END DO
282         END DO
283         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
284
285         
286         !
287         ! Horizontal advective fluxes
288         ! ---------------------------
289         !
290         DO jk = 1, jpkm1                             
291            DO jj = 2, jpjm1
292               DO ji = fs_2, fs_jpim1   ! vector opt.         
293                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
294                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
295               END DO
296            END DO
297         END DO
298         !
299         DO jk = 1, jpkm1 
300            zdt =  p2dt(jk)
301            DO jj = 2, jpjm1
302               DO ji = fs_2, fs_jpim1   ! vector opt.   
303                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
304                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
305                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
306                  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
307                  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
308               END DO
309            END DO
310         END DO
311
312         !--- Lateral boundary conditions
313         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
314         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
315
316         !--- QUICKEST scheme
317         CALL quickest( zfu, zfd, zfc, zwy )
318         !
319         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
320         DO jk = 1, jpkm1 
321            DO jj = 2, jpjm1
322               DO ji = fs_2, fs_jpim1   ! vector opt.               
323                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
324               END DO
325            END DO
326         END DO
327         !--- Lateral boundary conditions
328         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
329         !
330         ! Tracer flux on the x-direction
331         DO jk = 1, jpkm1 
332            !
333            DO jj = 2, jpjm1
334               DO ji = fs_2, fs_jpim1   ! vector opt.               
335                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
336                  !--- If the second ustream point is a land point
337                  !--- the flux is computed by the 1st order UPWIND scheme
338                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
339                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
340                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
341               END DO
342            END DO
343         END DO
344         !
345         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
346         !
347         ! Computation of the trend
348         DO jk = 1, jpkm1 
349            DO jj = 2, jpjm1
350               DO ji = fs_2, fs_jpim1   ! vector opt. 
351                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
352                  ! horizontal advective trends
353                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
354                  !--- add it to the general tracer trends
355                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
356               END DO
357            END DO
358         END DO
359         !                                 ! trend diagnostics (contribution of upstream fluxes)
360         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
361         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
362         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
363           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
364           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
365         ENDIF
366         !
367      END DO
368      !
369      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
370      !
371   END SUBROUTINE tra_adv_qck_j
372
373
374   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
375     &                                    ptn, pta, kjpt )
376      !!----------------------------------------------------------------------
377      !!
378      !!----------------------------------------------------------------------
379      USE oce, ONLY:   zwz => ua   ! ua used as workspace
380      !
381      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
382      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
383      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
384      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
385      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
386      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
387      !
388      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
389      REAL(wp) ::   zbtr , ztra      ! local scalars
390      !!----------------------------------------------------------------------
391
392      !                                                          ! ===========
393      DO jn = 1, kjpt                                            ! tracer loop
394         !                                                       ! ===========
395         ! 1. Bottom value : flux set to zero
396         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
397         !
398         !                                 ! Surface value
399         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
400         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
401         ENDIF
402         !
403         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
404            DO jj = 2, jpjm1
405               DO ji = fs_2, fs_jpim1   ! vector opt.
406                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
407               END DO
408            END DO
409         END DO
410         !
411         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
412            DO jj = 2, jpjm1
413               DO ji = fs_2, fs_jpim1   ! vector opt.
414                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
415                  ! k- vertical advective trends
416                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
417                  ! added to the general tracer trends
418                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
419               END DO
420            END DO
421         END DO
422         !                                 ! Save the vertical advective trends for diagnostic
423         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
424         !
425      END DO
426      !
427   END SUBROUTINE tra_adv_cen2_k
428
429
430   SUBROUTINE quickest( pfu, pfd, pfc, puc )
431      !!----------------------------------------------------------------------
432      !!
433      !! ** Purpose :  Computation of advective flux with Quickest scheme
434      !!
435      !! ** Method :   
436      !!----------------------------------------------------------------------
437      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
438      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
440      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
441      !!
442      INTEGER  ::  ji, jj, jk               ! dummy loop indices
443      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
444      REAL(wp) ::  zc, zcurv, zfho          !   -      -
445      !----------------------------------------------------------------------
446      !
447      IF( nn_timing == 1 )  CALL timing_start('quickest')
448      !
449      DO jk = 1, jpkm1
450         DO jj = 1, jpj
451            DO ji = 1, jpi
452               zc     = puc(ji,jj,jk)                         ! Courant number
453               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
454               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
455               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
456               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
457               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
458               !
459               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
460               zcoef2 = ABS( zcoef1 )
461               zcoef3 = ABS( zcurv )
462               IF( zcoef3 >= zcoef2 ) THEN
463                  zfho = pfc(ji,jj,jk) 
464               ELSE
465                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
466                  IF( zcoef1 >= 0. ) THEN
467                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
468                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
469                  ELSE
470                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
471                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
472                  ENDIF
473               ENDIF
474               puc(ji,jj,jk) = zfho
475            END DO
476         END DO
477      END DO
478      !
479      IF( nn_timing == 1 )  CALL timing_stop('quickest')
480      !
481   END SUBROUTINE quickest
482
483   !!======================================================================
484END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.