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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

Ediag branche: #927 split TRA/DYN trd computation
  • Property svn:keywords set to Id
File size: 24.0 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 trd_oce         ! trends: ocean variables
20   USE trdtra          ! trends manager: tracers
21   USE lib_mpp         ! distribued memory computing
22   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
23   USE dynspg_oce      ! surface pressure gradient variables
24   USE in_out_manager  ! I/O manager
25   USE diaptr          ! poleward transport diagnostics
26   USE trc_oce         ! share passive tracers/Ocean variables
27   USE wrk_nemo        ! Memory Allocation
28   USE timing          ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_adv_qck   ! routine called by step.F90
34
35   LOGICAL  :: l_trd           ! flag to compute trends
36   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE tra_adv_qck ( kt, kit000, 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      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
85      INTEGER                              , INTENT(in   ) ::   kit000          ! first 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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
89      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
90      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
91      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
92      !!----------------------------------------------------------------------
93
94      !
95      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
96      !
97      IF( kt == kit000 )  THEN
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
100         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
101         IF(lwp) WRITE(numout,*)
102         !
103         l_trd = .FALSE.
104         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
105      ENDIF
106
107      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
108      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
109      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
110
111      ! II. The vertical fluxes are computed with the 2nd order centered scheme
112      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
113      !
114      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
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 oce     , ONLY:   zwx => ua       ! ua used as workspace
125      !
126      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
127      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
128      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
129      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
130      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
131      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
132      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
133      !!
134      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
135      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
136      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
137      !----------------------------------------------------------------------
138      !
139      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
140      !                                                          ! ===========
141      DO jn = 1, kjpt                                            ! tracer loop
142         !                                                       ! ===========
143         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
144         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
145         !                                                 
146         DO jk = 1, jpkm1                               
147            !                                             
148            !--- Computation of the ustream and downstream value of the tracer and the mask
149            DO jj = 2, jpjm1
150               DO ji = fs_2, fs_jpim1   ! vector opt.
151                  ! Upstream in the x-direction for the tracer
152                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
153                  ! Downstream in the x-direction for the tracer
154                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
155               END DO
156            END DO
157         END DO
158         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
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         !--- Lateral boundary conditions
186         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
187         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
188
189         !--- QUICKEST scheme
190         CALL quickest( zfu, zfd, zfc, zwx )
191         !
192         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
193         DO jk = 1, jpkm1 
194            DO jj = 2, jpjm1
195               DO ji = fs_2, fs_jpim1   ! vector opt.               
196                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
197               END DO
198            END DO
199         END DO
200         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
201
202         !
203         ! Tracer flux on the x-direction
204         DO jk = 1, jpkm1 
205            !
206            DO jj = 2, jpjm1
207               DO ji = fs_2, fs_jpim1   ! vector opt.               
208                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
209                  !--- If the second ustream point is a land point
210                  !--- the flux is computed by the 1st order UPWIND scheme
211                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
212                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
213                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
214               END DO
215            END DO
216         END DO
217         !
218         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
219         !
220         ! Computation of the trend
221         DO jk = 1, jpkm1 
222            DO jj = 2, jpjm1
223               DO ji = fs_2, fs_jpim1   ! vector opt. 
224                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
225                  ! horizontal advective trends
226                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
227                  !--- add it to the general tracer trends
228                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
229               END DO
230            END DO
231         END DO
232         !                                 ! trend diagnostics (contribution of upstream fluxes)
233         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
234         !
235      END DO
236      !
237      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
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 oce     , ONLY:   zwy => ua       ! ua used as workspace
248      !
249      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
250      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
251      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
252      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
253      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
254      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
255      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
256      !!
257      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
258      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
259      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
260      !----------------------------------------------------------------------
261      !
262      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
263      !
264      !                                                          ! ===========
265      DO jn = 1, kjpt                                            ! tracer loop
266         !                                                       ! ===========
267         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
268         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
269         !                                                 
270         DO jk = 1, jpkm1                               
271            !                                             
272            !--- Computation of the ustream and downstream value of the tracer and the mask
273            DO jj = 2, jpjm1
274               DO ji = fs_2, fs_jpim1   ! vector opt.
275                  ! Upstream in the x-direction for the tracer
276                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
277                  ! Downstream in the x-direction for the tracer
278                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
279               END DO
280            END DO
281         END DO
282         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
283
284         
285         !
286         ! Horizontal advective fluxes
287         ! ---------------------------
288         !
289         DO jk = 1, jpkm1                             
290            DO jj = 2, jpjm1
291               DO ji = fs_2, fs_jpim1   ! vector opt.         
292                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
293                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
294               END DO
295            END DO
296         END DO
297         !
298         DO jk = 1, jpkm1 
299            zdt =  p2dt(jk)
300            DO jj = 2, jpjm1
301               DO ji = fs_2, fs_jpim1   ! vector opt.   
302                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
303                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
304                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
305                  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
306                  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
307               END DO
308            END DO
309         END DO
310
311         !--- Lateral boundary conditions
312         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
313         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
314
315         !--- QUICKEST scheme
316         CALL quickest( zfu, zfd, zfc, zwy )
317         !
318         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
319         DO jk = 1, jpkm1 
320            DO jj = 2, jpjm1
321               DO ji = fs_2, fs_jpim1   ! vector opt.               
322                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
323               END DO
324            END DO
325         END DO
326         !--- Lateral boundary conditions
327         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
328         !
329         ! Tracer flux on the x-direction
330         DO jk = 1, jpkm1 
331            !
332            DO jj = 2, jpjm1
333               DO ji = fs_2, fs_jpim1   ! vector opt.               
334                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
335                  !--- If the second ustream point is a land point
336                  !--- the flux is computed by the 1st order UPWIND scheme
337                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
338                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
339                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
340               END DO
341            END DO
342         END DO
343         !
344         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
345         !
346         ! Computation of the trend
347         DO jk = 1, jpkm1 
348            DO jj = 2, jpjm1
349               DO ji = fs_2, fs_jpim1   ! vector opt. 
350                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
351                  ! horizontal advective trends
352                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
353                  !--- add it to the general tracer trends
354                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
355               END DO
356            END DO
357         END DO
358         !                                 ! trend diagnostics (contribution of upstream fluxes)
359         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_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      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
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_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      IF( nn_timing == 1 )  CALL timing_start('quickest')
447      !
448      DO jk = 1, jpkm1
449         DO jj = 1, jpj
450            DO ji = 1, jpi
451               zc     = puc(ji,jj,jk)                         ! Courant number
452               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
453               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
454               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
455               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
456               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
457               !
458               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
459               zcoef2 = ABS( zcoef1 )
460               zcoef3 = ABS( zcurv )
461               IF( zcoef3 >= zcoef2 ) THEN
462                  zfho = pfc(ji,jj,jk) 
463               ELSE
464                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
465                  IF( zcoef1 >= 0. ) THEN
466                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
467                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
468                  ELSE
469                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
470                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
471                  ENDIF
472               ENDIF
473               puc(ji,jj,jk) = zfho
474            END DO
475         END DO
476      END DO
477      !
478      IF( nn_timing == 1 )  CALL timing_stop('quickest')
479      !
480   END SUBROUTINE quickest
481
482   !!======================================================================
483END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.