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/branches/2020/dev_r13508_HPC-09_communications_cleanup/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13508_HPC-09_communications_cleanup/src/OCE/TRA/traadv_qck.F90 @ 13701

Last change on this file since 13701 was 13701, checked in by francesca, 4 years ago

end of cleanup of the TRA modules - ticket #2367
  • Property svn:keywords set to Id
File size: 21.3 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      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(inout) ::   pU, pV, pW      ! 3 ocean volume transport components
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
95      !!----------------------------------------------------------------------
96      !
97      IF( kt == kit000 )  THEN
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
100         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
101         IF(lwp) WRITE(numout,*)
102      ENDIF
103      !
104      l_trd = .FALSE.
105      l_ptr = .FALSE.
106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
107      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 
108      !
109      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
110      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
111      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
112
113      !        ! vertical fluxes are computed with the 2nd order centered scheme
114      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
115      !
116   END SUBROUTINE tra_adv_qck
117
118
119   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
120      !!----------------------------------------------------------------------
121      !!
122      !!----------------------------------------------------------------------
123      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
124      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
125      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
126      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
127      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
128      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components
129      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
130      !!
131      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
132      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
133      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zfu, zfc, zfd
134      !----------------------------------------------------------------------
135      !
136      !                                                          ! ===========
137      DO jn = 1, kjpt                                            ! tracer loop
138         !                                                       ! ===========
139         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
140         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
141         !
142!!gm why not using a SHIFT instruction...
143         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
144            zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
145            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
146         END_3D
147         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
148         
149         !
150         ! Horizontal advective fluxes
151         ! ---------------------------
152         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
153            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
154            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
155         END_3D
156         !
157         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
158            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
159            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
160            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
161            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
162            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
163         END_3D
164         !--- Lateral boundary conditions
165         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
166
167         !--- QUICKEST scheme
168         CALL quickest( zfu, zfd, zfc, zwx )
169         !
170         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
171         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )
172            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
173         END_3D
174         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions
175
176         !
177         ! Tracer flux on the x-direction
178         DO_3D( 0, 0, 1, 0, 1, jpkm1 )
179            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
180            !--- If the second ustream point is a land point
181            !--- the flux is computed by the 1st order UPWIND scheme
182            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
183            zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
184            zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
185         END_3D
186         !
187         ! Computation of the trend
188         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
189            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
190            ! horizontal advective trends
191            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
192            !--- add it to the general tracer trends
193            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
194         END_3D
195         !                                 ! trend diagnostics
196         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
197         !
198      END DO
199      !
200   END SUBROUTINE tra_adv_qck_i
201
202
203   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
204      !!----------------------------------------------------------------------
205      !!
206      !!----------------------------------------------------------------------
207      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
208      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
209      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
210      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
211      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
212      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
213      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
214      !!
215      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
216      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
217      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
218      !----------------------------------------------------------------------
219      !
220      !                                                          ! ===========
221      DO jn = 1, kjpt                                            ! tracer loop
222         !                                                       ! ===========
223         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
224         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
225         !                                                 
226         DO jk = 1, jpkm1                               
227            !                                             
228            !--- Computation of the ustream and downstream value of the tracer and the mask
229            DO_2D( nn_hls-1, nn_hls-1, 0, 0 )
230               ! Upstream in the x-direction for the tracer
231               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
232               ! Downstream in the x-direction for the tracer
233               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
234            END_2D
235         END DO
236         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
237
238         
239         !
240         ! Horizontal advective fluxes
241         ! ---------------------------
242         !
243         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
244            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
245            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
246         END_3D
247         !
248         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
249            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
250            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
251            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
252            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
253            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
254         END_3D
255
256         !--- Lateral boundary conditions
257         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
258
259         !--- QUICKEST scheme
260         CALL quickest( zfu, zfd, zfc, zwy )
261         !
262         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
263         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )
264            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
265         END_3D
266         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
267         !
268         ! Tracer flux on the x-direction
269         DO_3D( 1, 0, 0, 0, 1, jpkm1 )
270            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
271            !--- If the second ustream point is a land point
272            !--- the flux is computed by the 1st order UPWIND scheme
273            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
274            zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
275            zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
276         END_3D
277         !
278         ! Computation of the trend
279         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
280            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
281            ! horizontal advective trends
282            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
283            !--- add it to the general tracer trends
284            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
285         END_3D
286         !                                 ! trend diagnostics
287         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
288         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
289         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
290         !
291      END DO
292      !
293   END SUBROUTINE tra_adv_qck_j
294
295
296   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
297      !!----------------------------------------------------------------------
298      !!
299      !!----------------------------------------------------------------------
300      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
301      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
302      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
303      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
304      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
305      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
306      !
307      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
308      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
309      !!----------------------------------------------------------------------
310      !
311      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
312      zwz(:,:,jpk) = 0._wp
313      !
314      !                                                          ! ===========
315      DO jn = 1, kjpt                                            ! tracer loop
316         !                                                       ! ===========
317         !
318         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
319            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)
320         END_3D
321         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
322            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
323               DO_2D( 0, 0, 0, 0 )
324                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
325               END_2D
326            ELSE                                   ! no ocean cavities (only ocean surface)
327               DO_2D( 0, 0, 0, 0 )
328                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
329               END_2D
330            ENDIF
331         ENDIF
332         !
333         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
334            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
335               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
336         END_3D
337         !                                 ! Send trends for diagnostic
338         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
339         !
340      END DO
341      !
342   END SUBROUTINE tra_adv_cen2_k
343
344
345   SUBROUTINE quickest( pfu, pfd, pfc, puc )
346      !!----------------------------------------------------------------------
347      !!
348      !! ** Purpose :  Computation of advective flux with Quickest scheme
349      !!
350      !! ** Method :   
351      !!----------------------------------------------------------------------
352      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
353      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
354      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
355      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
356      !!
357      INTEGER  ::  ji, jj, jk               ! dummy loop indices
358      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
359      REAL(wp) ::  zc, zcurv, zfho          !   -      -
360      !----------------------------------------------------------------------
361      !
362      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
363         zc     = puc(ji,jj,jk)                         ! Courant number
364         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
365         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
366         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
367         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
368         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
369         !
370         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
371         zcoef2 = ABS( zcoef1 )
372         zcoef3 = ABS( zcurv )
373         IF( zcoef3 >= zcoef2 ) THEN
374            zfho = pfc(ji,jj,jk) 
375         ELSE
376            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
377            IF( zcoef1 >= 0. ) THEN
378               zfho = MAX( pfc(ji,jj,jk), zfho ) 
379               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
380            ELSE
381               zfho = MIN( pfc(ji,jj,jk), zfho ) 
382               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
383            ENDIF
384         ENDIF
385         puc(ji,jj,jk) = zfho
386      END_3D
387      !
388   END SUBROUTINE quickest
389
390   !!======================================================================
391END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.