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

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3625

Last change on this file since 3625 was 3625, checked in by acc, 11 years ago

Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF
  • Property svn:keywords set to Id
File size: 24.2 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   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
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, kit000, 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      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
88      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
89      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
90      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
91      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
92      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
94      !!----------------------------------------------------------------------
95
96      !
97      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
98      !
99      IF( kt == kit000 )  THEN
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      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
117      !
118   END SUBROUTINE tra_adv_qck
119
120
121   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
122      &                                        ptb, ptn, pta, kjpt   )
123      !!----------------------------------------------------------------------
124      !!
125      !!----------------------------------------------------------------------
126      USE oce     , ONLY:   zwx => ua       ! ua used as 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      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
139      !----------------------------------------------------------------------
140      !
141      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
142      !                                                          ! ===========
143      DO jn = 1, kjpt                                            ! tracer loop
144         !                                                       ! ===========
145         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
146         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
147         !                                                 
148         DO jk = 1, jpkm1                               
149            !                                             
150            !--- Computation of the ustream and downstream value of the tracer and the mask
151            DO jj = 2, jpjm1
152               DO ji = fs_2, fs_jpim1   ! vector opt.
153                  ! Upstream in the x-direction for the tracer
154                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
155                  ! Downstream in the x-direction for the tracer
156                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
157               END DO
158            END DO
159         END DO
160         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
161         
162         !
163         ! Horizontal advective fluxes
164         ! ---------------------------
165         !
166         DO jk = 1, jpkm1                             
167            DO jj = 2, jpjm1
168               DO ji = fs_2, fs_jpim1   ! vector opt.         
169                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
170                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
171               END DO
172            END DO
173         END DO
174         !
175         DO jk = 1, jpkm1 
176            zdt =  p2dt(jk)
177            DO jj = 2, jpjm1
178               DO ji = fs_2, fs_jpim1   ! vector opt.   
179                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
180                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
181                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
182                  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
183                  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
184               END DO
185            END DO
186         END DO 
187         !--- Lateral boundary conditions
188         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
189         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
190
191         !--- QUICKEST scheme
192         CALL quickest( zfu, zfd, zfc, zwx )
193         !
194         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
195         DO jk = 1, jpkm1 
196            DO jj = 2, jpjm1
197               DO ji = fs_2, fs_jpim1   ! vector opt.               
198                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
199               END DO
200            END DO
201         END DO
202         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
203
204         !
205         ! Tracer flux on the x-direction
206         DO jk = 1, jpkm1 
207            !
208            DO jj = 2, jpjm1
209               DO ji = fs_2, fs_jpim1   ! vector opt.               
210                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
211                  !--- If the second ustream point is a land point
212                  !--- the flux is computed by the 1st order UPWIND scheme
213                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
214                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
215                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
216               END DO
217            END DO
218         END DO
219         !
220         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
221         !
222         ! Computation of the trend
223         DO jk = 1, jpkm1 
224            DO jj = 2, jpjm1
225               DO ji = fs_2, fs_jpim1   ! vector opt. 
226                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
227                  ! horizontal advective trends
228                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
229                  !--- add it to the general tracer trends
230                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
231               END DO
232            END DO
233         END DO
234         !                                 ! trend diagnostics (contribution of upstream fluxes)
235         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
236         !
237      END DO
238      !
239      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
240      !
241   END SUBROUTINE tra_adv_qck_i
242
243
244   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
245      &                                        ptb, ptn, pta, kjpt )
246      !!----------------------------------------------------------------------
247      !!
248      !!----------------------------------------------------------------------
249      USE oce     , ONLY:   zwy => ua       ! ua used as 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      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
262      !----------------------------------------------------------------------
263      !
264      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
265      !
266      !                                                          ! ===========
267      DO jn = 1, kjpt                                            ! tracer loop
268         !                                                       ! ===========
269         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
270         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
271         !                                                 
272         DO jk = 1, jpkm1                               
273            !                                             
274            !--- Computation of the ustream and downstream value of the tracer and the mask
275            DO jj = 2, jpjm1
276               DO ji = fs_2, fs_jpim1   ! vector opt.
277                  ! Upstream in the x-direction for the tracer
278                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
279                  ! Downstream in the x-direction for the tracer
280                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
281               END DO
282            END DO
283         END DO
284         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
285
286         
287         !
288         ! Horizontal advective fluxes
289         ! ---------------------------
290         !
291         DO jk = 1, jpkm1                             
292            DO jj = 2, jpjm1
293               DO ji = fs_2, fs_jpim1   ! vector opt.         
294                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
295                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
296               END DO
297            END DO
298         END DO
299         !
300         DO jk = 1, jpkm1 
301            zdt =  p2dt(jk)
302            DO jj = 2, jpjm1
303               DO ji = fs_2, fs_jpim1   ! vector opt.   
304                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
305                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
306                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
307                  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
308                  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
309               END DO
310            END DO
311         END DO
312
313         !--- Lateral boundary conditions
314         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
315         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
316
317         !--- QUICKEST scheme
318         CALL quickest( zfu, zfd, zfc, zwy )
319         !
320         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
321         DO jk = 1, jpkm1 
322            DO jj = 2, jpjm1
323               DO ji = fs_2, fs_jpim1   ! vector opt.               
324                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
325               END DO
326            END DO
327         END DO
328         !--- Lateral boundary conditions
329         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
330         !
331         ! Tracer flux on the x-direction
332         DO jk = 1, jpkm1 
333            !
334            DO jj = 2, jpjm1
335               DO ji = fs_2, fs_jpim1   ! vector opt.               
336                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
337                  !--- If the second ustream point is a land point
338                  !--- the flux is computed by the 1st order UPWIND scheme
339                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
340                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
341                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
342               END DO
343            END DO
344         END DO
345         !
346         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
347         !
348         ! Computation of the trend
349         DO jk = 1, jpkm1 
350            DO jj = 2, jpjm1
351               DO ji = fs_2, fs_jpim1   ! vector opt. 
352                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
353                  ! horizontal advective trends
354                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
355                  !--- add it to the general tracer trends
356                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
357               END DO
358            END DO
359         END DO
360         !                                 ! trend diagnostics (contribution of upstream fluxes)
361         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
362         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
363         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
364           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
365           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
366         ENDIF
367         !
368      END DO
369      !
370      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
371      !
372   END SUBROUTINE tra_adv_qck_j
373
374
375   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
376     &                                    ptn, pta, kjpt )
377      !!----------------------------------------------------------------------
378      !!
379      !!----------------------------------------------------------------------
380      USE oce, ONLY:   zwz => ua   ! ua used as workspace
381      !
382      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
383      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
384      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
385      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
386      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
387      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
388      !
389      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
390      REAL(wp) ::   zbtr , ztra      ! local scalars
391      !!----------------------------------------------------------------------
392
393      !                                                          ! ===========
394      DO jn = 1, kjpt                                            ! tracer loop
395         !                                                       ! ===========
396         ! 1. Bottom value : flux set to zero
397         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
398         !
399         !                                 ! Surface value
400         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
401         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
402         ENDIF
403         !
404         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
405            DO jj = 2, jpjm1
406               DO ji = fs_2, fs_jpim1   ! vector opt.
407                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
408               END DO
409            END DO
410         END DO
411         !
412         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
413            DO jj = 2, jpjm1
414               DO ji = fs_2, fs_jpim1   ! vector opt.
415                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
416                  ! k- vertical advective trends
417                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
418                  ! added to the general tracer trends
419                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
420               END DO
421            END DO
422         END DO
423         !                                 ! Save the vertical advective trends for diagnostic
424         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
425         !
426      END DO
427      !
428   END SUBROUTINE tra_adv_cen2_k
429
430
431   SUBROUTINE quickest( pfu, pfd, pfc, puc )
432      !!----------------------------------------------------------------------
433      !!
434      !! ** Purpose :  Computation of advective flux with Quickest scheme
435      !!
436      !! ** Method :   
437      !!----------------------------------------------------------------------
438      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
440      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
441      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
442      !!
443      INTEGER  ::  ji, jj, jk               ! dummy loop indices
444      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
445      REAL(wp) ::  zc, zcurv, zfho          !   -      -
446      !----------------------------------------------------------------------
447      !
448      IF( nn_timing == 1 )  CALL timing_start('quickest')
449      !
450      DO jk = 1, jpkm1
451         DO jj = 1, jpj
452            DO ji = 1, jpi
453               zc     = puc(ji,jj,jk)                         ! Courant number
454               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
455               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
456               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
457               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
458               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
459               !
460               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
461               zcoef2 = ABS( zcoef1 )
462               zcoef3 = ABS( zcurv )
463               IF( zcoef3 >= zcoef2 ) THEN
464                  zfho = pfc(ji,jj,jk) 
465               ELSE
466                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
467                  IF( zcoef1 >= 0. ) THEN
468                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
469                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
470                  ELSE
471                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
472                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
473                  ENDIF
474               ENDIF
475               puc(ji,jj,jk) = zfho
476            END DO
477         END DO
478      END DO
479      !
480      IF( nn_timing == 1 )  CALL timing_stop('quickest')
481      !
482   END SUBROUTINE quickest
483
484   !!======================================================================
485END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.