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_r12377_KERNEL-06_techene_e3/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_qck.F90 @ 12606

Last change on this file since 12606 was 12606, checked in by techene, 4 years ago

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