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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3150

Last change on this file since 3150 was 3116, checked in by cetlod, 13 years ago

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