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 NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traadv_qck.F90 @ 14433

Last change on this file since 14433 was 14433, checked in by smasson, 3 years ago

trunk: merge dev_r14312_MPI_Interface into the trunk, #2598
  • Property svn:keywords set to Id
File size: 21.6 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 trc_oce         ! share passive tracers/Ocean variables
20   USE trd_oce         ! trends: ocean variables
21   USE trdtra          ! trends manager: tracers
22   USE diaptr          ! poleward transport diagnostics
23   USE iom
24   !
25   USE in_out_manager  ! I/O manager
26   USE lib_mpp         ! distribued memory computing
27   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
28   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_adv_qck   ! routine called by step.F90
34
35   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
36
37   LOGICAL  ::   l_trd   ! flag to compute trends
38   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
39
40
41   !! * Substitutions
42#  include "do_loop_substitute.h90"
43#  include "domzgr_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
46   !! $Id$
47   !! Software governed by the CeCILL license (see ./LICENSE)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_adv_qck  ***
54      !!
55      !! ** Purpose :   Compute the now trend due to the advection of tracers
56      !!      and add it to the general trend of passive tracer equations.
57      !!
58      !! ** Method :   The advection is evaluated by a third order scheme
59      !!             For a positive velocity u :              u(i)>0
60      !!                                          |--FU--|--FC--|--FD--|------|
61      !!                                             i-1    i      i+1   i+2
62      !!
63      !!             For a negative velocity u :              u(i)<0
64      !!                                          |------|--FD--|--FC--|--FU--|
65      !!                                             i-1    i      i+1   i+2
66      !!             where  FU is the second upwind point
67      !!                    FD is the first douwning point
68      !!                    FC is the central point (or the first upwind point)
69      !!
70      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
71      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
72      !!
73      !!         dt = 2*rdtra and the scalar values are tb and sb
74      !!
75      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm)
76      !!
77      !!               The fluxes are bounded by the ULTIMATE limiter to
78      !!             guarantee the monotonicity of the solution and to
79      !!            prevent the appearance of spurious numerical oscillations
80      !!
81      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
82      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
83      !!             - poleward advective heat and salt transport (ln_diaptr=T)
84      !!
85      !! ** Reference : Leonard (1979, 1991)
86      !!----------------------------------------------------------------------
87      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
88      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
89      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
90      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
91      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
92      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
93      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
94      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
95      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
96      !!----------------------------------------------------------------------
97      !
98      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
99         IF( kt == kit000 )  THEN
100            IF(lwp) WRITE(numout,*)
101            IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
102            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
103            IF(lwp) WRITE(numout,*)
104         ENDIF
105         !
106         l_trd = .FALSE.
107         l_ptr = .FALSE.
108         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
109         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
110      ENDIF
111      !
112      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
113      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
114      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
115
116      !        ! vertical fluxes are computed with the 2nd order centered scheme
117      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
118      !
119   END SUBROUTINE tra_adv_qck
120
121
122   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
123      !!----------------------------------------------------------------------
124      !!
125      !!----------------------------------------------------------------------
126      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
127      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
128      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
129      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
130      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
131      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
132      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components
133      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
134      !!
135      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
136      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
137      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwx, zfu, zfc, zfd
138      !----------------------------------------------------------------------
139      !
140      !                                                          ! ===========
141      DO jn = 1, kjpt                                            ! tracer loop
142         !                                                       ! ===========
143         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp
144         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp
145         !
146!!gm why not using a SHIFT instruction...
147         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
148            zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
149            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
150         END_3D
151         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
152
153         !
154         ! Horizontal advective fluxes
155         ! ---------------------------
156         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
157            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
158            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
159         END_3D
160         !
161         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
162            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
163            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
164            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
165            zfc(ji,jj,jk)  = zdir * pt(ji  ,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji+1,jj,jk,jn,Kbb)  ! FC in the x-direction for T
166            zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T
167         END_3D
168         !--- Lateral boundary conditions
169         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
170
171         !--- QUICKEST scheme
172         CALL quickest( zfu, zfd, zfc, zwx )
173         !
174         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
175         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )
176            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
177         END_3D
178         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions
179
180         !
181         ! Tracer flux on the x-direction
182         DO_3D( 1, 0, 0, 0, 1, jpkm1 )
183            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
184            !--- If the second ustream point is a land point
185            !--- the flux is computed by the 1st order UPWIND scheme
186            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
187            zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
188            zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
189         END_3D
190         !
191         ! Computation of the trend
192         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
193            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
194            ! horizontal advective trends
195            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
196            !--- add it to the general tracer trends
197            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
198         END_3D
199         !                                 ! trend diagnostics
200         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
201         !
202      END DO
203      !
204   END SUBROUTINE tra_adv_qck_i
205
206
207   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
208      !!----------------------------------------------------------------------
209      !!
210      !!----------------------------------------------------------------------
211      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
212      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
213      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
214      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
215      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
216      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
217      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
218      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
219      !!
220      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
221      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
222      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
223      !----------------------------------------------------------------------
224      !
225      !                                                          ! ===========
226      DO jn = 1, kjpt                                            ! tracer loop
227         !                                                       ! ===========
228         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0
229         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0
230         !
231         DO jk = 1, jpkm1
232            !
233            !--- Computation of the ustream and downstream value of the tracer and the mask
234            DO_2D( 0, 0, nn_hls-1, nn_hls-1 )
235               ! Upstream in the x-direction for the tracer
236               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
237               ! Downstream in the x-direction for the tracer
238               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
239            END_2D
240         END DO
241         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
242
243         !
244         ! Horizontal advective fluxes
245         ! ---------------------------
246         !
247         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
248            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
249            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
250         END_3D
251         !
252         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
253            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
254            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
255            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
256            zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T
257            zfd(ji,jj,jk)  = zdir * pt(ji,jj+1,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj  ,jk,jn,Kbb)  ! FD in the x-direction for T
258         END_3D
259
260         !--- Lateral boundary conditions
261         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
262
263         !--- QUICKEST scheme
264         CALL quickest( zfu, zfd, zfc, zwy )
265         !
266         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
267         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )
268            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
269         END_3D
270         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
271         !
272         ! Tracer flux on the x-direction
273         DO_3D( 0, 0, 1, 0, 1, jpkm1 )
274            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
275            !--- If the second ustream point is a land point
276            !--- the flux is computed by the 1st order UPWIND scheme
277            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
278            zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
279            zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
280         END_3D
281         !
282         ! Computation of the trend
283         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
284            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
285            ! horizontal advective trends
286            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
287            !--- add it to the general tracer trends
288            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
289         END_3D
290         !                                 ! trend diagnostics
291         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
292         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
293         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
294         !
295      END DO
296      !
297   END SUBROUTINE tra_adv_qck_j
298
299
300   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
301      !!----------------------------------------------------------------------
302      !!
303      !!----------------------------------------------------------------------
304      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
305      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
306      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
307      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
308      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
309      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
310      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
311      !
312      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
313      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwz   ! 3D workspace
314      !!----------------------------------------------------------------------
315      !
316      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
317      zwz(:,:,jpk) = 0._wp
318      !
319      !                                                          ! ===========
320      DO jn = 1, kjpt                                            ! tracer loop
321         !                                                       ! ===========
322         !
323         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
324            zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
325         END_3D
326         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
327            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
328               DO_2D( 0, 0, 0, 0 )
329                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
330               END_2D
331            ELSE                                   ! no ocean cavities (only ocean surface)
332               DO_2D( 0, 0, 0, 0 )
333                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
334               END_2D
335            ENDIF
336         ENDIF
337         !
338         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
339            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
340               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
341         END_3D
342         !                                 ! Send trends for diagnostic
343         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
344         !
345      END DO
346      !
347   END SUBROUTINE tra_adv_cen2_k
348
349
350   SUBROUTINE quickest( pfu, pfd, pfc, puc )
351      !!----------------------------------------------------------------------
352      !!
353      !! ** Purpose :  Computation of advective flux with Quickest scheme
354      !!
355      !! ** Method :
356      !!----------------------------------------------------------------------
357      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfu   ! second upwind point
358      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfd   ! first douwning point
359      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
360      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
361      !!
362      INTEGER  ::  ji, jj, jk               ! dummy loop indices
363      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars
364      REAL(wp) ::  zc, zcurv, zfho          !   -      -
365      !----------------------------------------------------------------------
366      !
367      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
368         zc     = puc(ji,jj,jk)                         ! Courant number
369         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
370         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
371         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
372         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
373         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
374         !
375         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
376         zcoef2 = ABS( zcoef1 )
377         zcoef3 = ABS( zcurv )
378         IF( zcoef3 >= zcoef2 ) THEN
379            zfho = pfc(ji,jj,jk)
380         ELSE
381            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
382            IF( zcoef1 >= 0. ) THEN
383               zfho = MAX( pfc(ji,jj,jk), zfho )
384               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
385            ELSE
386               zfho = MIN( pfc(ji,jj,jk), zfho )
387               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
388            ENDIF
389         ENDIF
390         puc(ji,jj,jk) = zfho
391      END_3D
392      !
393   END SUBROUTINE quickest
394
395   !!======================================================================
396END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.