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

Last change on this file since 13295 was 13295, checked in by acc, 9 months ago

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision

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