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 @ 13660

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

communications cleanup of the TRA modules - ticket #2367
  • Property svn:keywords set to Id
File size: 21.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  : 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 jk = 1, jpkm1 
179            !
180            DO_2D( 0, 0, 1, 0 )
181               zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
182               !--- If the second ustream point is a land point
183               !--- the flux is computed by the 1st order UPWIND scheme
184               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
185               zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
186               zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
187            END_2D
188         END DO
189         !
190         ! Computation of the trend
191         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
192            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
193            ! horizontal advective trends
194            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
195            !--- add it to the general tracer trends
196            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
197         END_3D
198         !                                 ! trend diagnostics
199         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
200         !
201      END DO
202      !
203   END SUBROUTINE tra_adv_qck_i
204
205
206   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
207      !!----------------------------------------------------------------------
208      !!
209      !!----------------------------------------------------------------------
210      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
211      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
212      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
213      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
214      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
215      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
216      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
217      !!
218      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
219      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
220      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
221      !----------------------------------------------------------------------
222      !
223      !                                                          ! ===========
224      DO jn = 1, kjpt                                            ! tracer loop
225         !                                                       ! ===========
226         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
227         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
228         !                                                 
229         DO jk = 1, jpkm1                               
230            !                                             
231            !--- Computation of the ustream and downstream value of the tracer and the mask
232            DO_2D( nn_hls-1, nn_hls-1, 0, 0 )
233               ! Upstream in the x-direction for the tracer
234               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
235               ! Downstream in the x-direction for the tracer
236               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
237            END_2D
238         END DO
239         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
240
241         
242         !
243         ! Horizontal advective fluxes
244         ! ---------------------------
245         !
246         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
247            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
248            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
249         END_3D
250         !
251         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
252            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
253            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
254            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
255            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
256            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
257         END_3D
258
259         !--- Lateral boundary conditions
260         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 )
261
262         !--- QUICKEST scheme
263         CALL quickest( zfu, zfd, zfc, zwy )
264         !
265         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
266         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )
267            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
268         END_3D
269         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
270         !
271         ! Tracer flux on the x-direction
272         DO jk = 1, jpkm1 
273            !
274            DO_2D( 1, 0, 0, 0 )
275               zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
276               !--- If the second ustream point is a land point
277               !--- the flux is computed by the 1st order UPWIND scheme
278               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
279               zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
280               zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
281            END_2D
282         END DO
283         !
284         ! Computation of the trend
285         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
286            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
287            ! horizontal advective trends
288            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
289            !--- add it to the general tracer trends
290            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
291         END_3D
292         !                                 ! trend diagnostics
293         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
294         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
295         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
296         !
297      END DO
298      !
299   END SUBROUTINE tra_adv_qck_j
300
301
302   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
303      !!----------------------------------------------------------------------
304      !!
305      !!----------------------------------------------------------------------
306      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
307      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
308      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
309      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
310      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
311      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
312      !
313      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
314      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
315      !!----------------------------------------------------------------------
316      !
317      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
318      zwz(:,:,jpk) = 0._wp
319      !
320      !                                                          ! ===========
321      DO jn = 1, kjpt                                            ! tracer loop
322         !                                                       ! ===========
323         !
324         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
325            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)
326         END_3D
327         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
328            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
329               DO_2D( 0, 0, 0, 0 )
330                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
331               END_2D
332            ELSE                                   ! no ocean cavities (only ocean surface)
333               DO_2D( 0, 0, 0, 0 )
334                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
335               END_2D
336            ENDIF
337         ENDIF
338         !
339         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
340            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
341               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
342         END_3D
343         !                                 ! Send trends for diagnostic
344         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
345         !
346      END DO
347      !
348   END SUBROUTINE tra_adv_cen2_k
349
350
351   SUBROUTINE quickest( pfu, pfd, pfc, puc )
352      !!----------------------------------------------------------------------
353      !!
354      !! ** Purpose :  Computation of advective flux with Quickest scheme
355      !!
356      !! ** Method :   
357      !!----------------------------------------------------------------------
358      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
359      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
360      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
361      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
362      !!
363      INTEGER  ::  ji, jj, jk               ! dummy loop indices
364      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
365      REAL(wp) ::  zc, zcurv, zfho          !   -      -
366      !----------------------------------------------------------------------
367      !
368      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
369         zc     = puc(ji,jj,jk)                         ! Courant number
370         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
371         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
372         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
373         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
374         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
375         !
376         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
377         zcoef2 = ABS( zcoef1 )
378         zcoef3 = ABS( zcurv )
379         IF( zcoef3 >= zcoef2 ) THEN
380            zfho = pfc(ji,jj,jk) 
381         ELSE
382            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
383            IF( zcoef1 >= 0. ) THEN
384               zfho = MAX( pfc(ji,jj,jk), zfho ) 
385               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
386            ELSE
387               zfho = MIN( pfc(ji,jj,jk), zfho ) 
388               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
389            ENDIF
390         ENDIF
391         puc(ji,jj,jk) = zfho
392      END_3D
393      !
394   END SUBROUTINE quickest
395
396   !!======================================================================
397END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.