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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

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