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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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