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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2034

Last change on this file since 2034 was 2034, checked in by cetlod, 14 years ago

cosmetic changes
File size: 23.5 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  :
14   !!   tra_adv_qck_j  :
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
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   tra_adv_qck   ! routine called by step.F90
32
33   REAL(wp)  :: r1_6 = 1./ 6.
34   LOGICAL   :: l_trd    ! flag to compute trends
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
41   !! $Id: traadv_qck.F90 2024 2010-07-29 10:57:35Z cetlod $
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE tra_adv_qck ( kt, cdtype, 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      !!
84      INTEGER         , INTENT(in   )                               ::   kt              ! ocean 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)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
88      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
89      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
90      !!
91      REAL(wp) ::   z2                            !    temporary scalar
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'
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      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
105      ELSE                                        ;    z2 = 2.
106      ENDIF
107
108      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
109      !---------------------------------------------------------------------------
110
111      CALL tra_adv_qck_i( kt, cdtype, pun, z2, ptb, ptn, pta, kjpt ) 
112      CALL tra_adv_qck_j( kt, cdtype, pvn, z2, ptb, ptn, pta, kjpt ) 
113
114      ! II. The vertical fluxes are computed with the 2nd order centered scheme
115      !-------------------------------------------------------------------------
116      !
117      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
118      !
119   END SUBROUTINE tra_adv_qck
120
121   SUBROUTINE tra_adv_qck_i( kt, cdtype, pun, pz2,   &
122      &                                  ptb, ptn, pta, kjpt   )
123      !!----------------------------------------------------------------------
124      !!
125      !!----------------------------------------------------------------------
126      USE oce         , zwx => ua   ! use ua 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)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun             ! zonal velocity component
132      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
133      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
134      REAL(wp)        , INTENT(in   )                               ::   pz2
135      !!
136      INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
137      REAL(wp) :: ztra, zbtr               ! temporary scalars
138      REAL(wp) :: zdir, zdx, zdt, zmsk     ! temporary scalars
139      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfu, zfc, zfd
140      !----------------------------------------------------------------------
141
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         !
161         !--- Lateral boundary conditions
162         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 
163         
164         !
165         ! Horizontal advective fluxes
166         ! ---------------------------
167         !
168         DO jk = 1, jpkm1                             
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.         
171                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
172                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
173               END DO
174            END DO
175         END DO
176         !
177         DO jk = 1, jpkm1 
178            zdt =  pz2 * rdttra(jk)
179            DO jj = 2, jpjm1
180               DO ji = fs_2, fs_jpim1   ! vector opt.   
181                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
182                  zdx  = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
183                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
184                  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
185                  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
186               END DO
187            END DO
188         END DO      !
189
190         !--- Lateral boundary conditions
191         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
192         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
193
194         !--- QUICKEST scheme
195         CALL quickest( zfu, zfd, zfc, zwx )
196         !
197         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
198         DO jk = 1, jpkm1 
199            DO jj = 2, jpjm1
200               DO ji = fs_2, fs_jpim1   ! vector opt.               
201                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
202               ENDDO
203            END DO
204         END DO
205         !--- Lateral boundary conditions
206         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
207         !
208         ! Tracer flux on the x-direction
209         DO jk = 1, jpkm1 
210            !
211            DO jj = 2, jpjm1
212               DO ji = fs_2, fs_jpim1   ! vector opt.               
213                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
214                  !--- If the second ustream point is a land point
215                  !--- the flux is computed by the 1st order UPWIND scheme
216                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
217                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
218                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
219               END DO
220            END DO
221            !
222            ! Computation of the trend
223            DO jj = 2, jpjm1
224               DO ji = fs_2, fs_jpim1   ! vector opt. 
225                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
226                  ! horizontal advective trends
227                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
228                  !--- add it to the general tracer trends
229                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
230               END DO
231            END DO
232            !
233         END DO
234         !                                 ! trend diagnostics (contribution of upstream fluxes)
235         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jpt_trd_xad, zwx, pun, ptn(:,:,:,jn) )
236         !
237      END DO
238      !
239   END SUBROUTINE tra_adv_qck_i
240
241   SUBROUTINE tra_adv_qck_j( kt, cdtype, pvn, pz2,   &
242      &                                  ptb, ptn, pta, kjpt   )
243      !!----------------------------------------------------------------------
244      !!
245      !!----------------------------------------------------------------------
246      !!
247      USE oce         , zwy => ua   ! use ua 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)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pvn             ! meridional velocity component
253      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
254      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
255      REAL(wp)        , INTENT(in   )                               ::   pz2
256      !!
257      INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
258      REAL(wp) :: ztra, zbtr               ! temporary scalars
259      REAL(wp) :: zdir, zdx, zdt, zmsk     ! temporary scalars
260      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfu, zfc, zfd
261      !----------------------------------------------------------------------
262
263      DO jn = 1, kjpt                                            ! tracer loop
264         !                                                       ! ===========
265         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
266         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
267         !                                                 
268         DO jk = 1, jpkm1                               
269            !                                             
270            !--- Computation of the ustream and downstream value of the tracer and the mask
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1   ! vector opt.
273                  ! Upstream in the x-direction for the tracer
274                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
275                  ! Downstream in the x-direction for the tracer
276                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
277               END DO
278            END DO
279         END DO
280         !
281         !--- Lateral boundary conditions
282         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 
283         
284         !
285         ! Horizontal advective fluxes
286         ! ---------------------------
287         !
288         DO jk = 1, jpkm1                             
289            DO jj = 2, jpjm1
290               DO ji = fs_2, fs_jpim1   ! vector opt.         
291                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
292                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
293               END DO
294            END DO
295         END DO
296         !
297         DO jk = 1, jpkm1 
298            zdt =  pz2 * rdttra(jk)
299            DO jj = 2, jpjm1
300               DO ji = fs_2, fs_jpim1   ! vector opt.   
301                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
302                  zdx  = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
303                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
304                  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
305                  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
306               END DO
307            END DO
308         END DO      !
309
310         !--- Lateral boundary conditions
311         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
312         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
313
314         !--- QUICKEST scheme
315         CALL quickest( zfu, zfd, zfc, zwy )
316         !
317         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
318         DO jk = 1, jpkm1 
319            DO jj = 2, jpjm1
320               DO ji = fs_2, fs_jpim1   ! vector opt.               
321                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
322               ENDDO
323            END DO
324         END DO
325         !--- Lateral boundary conditions
326         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
327         !
328         ! Tracer flux on the x-direction
329         DO jk = 1, jpkm1 
330            !
331            DO jj = 2, jpjm1
332               DO ji = fs_2, fs_jpim1   ! vector opt.               
333                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
334                  !--- If the second ustream point is a land point
335                  !--- the flux is computed by the 1st order UPWIND scheme
336                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
337                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
338                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
339               END DO
340            END DO
341            !
342            ! Computation of the trend
343            DO jj = 2, jpjm1
344               DO ji = fs_2, fs_jpim1   ! vector opt. 
345                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
346                  ! horizontal advective trends
347                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
348                  !--- add it to the general tracer trends
349                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
350               END DO
351            END DO
352            !
353         END DO
354         !                                 ! trend diagnostics (contribution of upstream fluxes)
355         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jpt_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
356         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
357         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
358           IF( jn == jp_tem )  pht_adv(:) = ptr_vj( zwy(:,:,:) )
359           IF( jn == jp_sal )  pst_adv(:) = ptr_vj( zwy(:,:,:) )
360         ENDIF
361         !
362      END DO
363
364   END SUBROUTINE tra_adv_qck_j
365
366   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,   &
367     &                                    ptn, pta, kjpt )
368      !!----------------------------------------------------------------------
369      !!
370      !!----------------------------------------------------------------------
371      !!
372      USE oce         , zwz => ua   ! use ua 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)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pwn             ! vertical velocity component
378      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn           ! now tracer field
379      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
380      !!
381      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
382      REAL(wp) ::   zbtr , ztra      ! temporary 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, jpt_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), INTENT(in)    , DIMENSION(jpi,jpj,jpk) :: pfu   ! second upwind point
431      REAL(wp), INTENT(in)    , DIMENSION(jpi,jpj,jpk) :: pfd   ! first douwning point
432      REAL(wp), INTENT(in)    , DIMENSION(jpi,jpj,jpk) :: pfc   ! the central point (or the first upwind point)
433      REAL(wp), INTENT(inout) , DIMENSION(jpi,jpj,jpk) :: puc   ! input as Courant number ; output as flux
434      !!
435      INTEGER  ::  ji, jj, jk               ! dummy loop indices
436      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! temporary 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            ENDDO
467         ENDDO
468      ENDDO
469      !
470   END SUBROUTINE quickest
471
472   !!======================================================================
473END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.