New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
traadv_qck.F90 in NEMO/trunk/src/OCE/TRA – NEMO

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

Last change on this file since 13226 was 13226, checked in by orioltp, 4 years ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk
  • Property svn:keywords set to Id
File size: 21.1 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   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id$
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_adv_qck  ***
53      !!
54      !! ** Purpose :   Compute the now trend due to the advection of tracers
55      !!      and add it to the general trend of passive tracer equations.
56      !!
57      !! ** Method :   The advection is evaluated by a third order scheme
58      !!             For a positive velocity u :              u(i)>0
59      !!                                          |--FU--|--FC--|--FD--|------|
60      !!                                             i-1    i      i+1   i+2
61      !!
62      !!             For a negative velocity u :              u(i)<0
63      !!                                          |------|--FD--|--FC--|--FU--|
64      !!                                             i-1    i      i+1   i+2
65      !!             where  FU is the second upwind point
66      !!                    FD is the first douwning point
67      !!                    FC is the central point (or the first upwind point)
68      !!
69      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
70      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
71      !!
72      !!         dt = 2*rdtra and the scalar values are tb and sb
73      !!
74      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm)
75      !!
76      !!               The fluxes are bounded by the ULTIMATE limiter to
77      !!             guarantee the monotonicity of the solution and to
78      !!            prevent the appearance of spurious numerical oscillations
79      !!
80      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
81      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
82      !!             - poleward advective heat and salt transport (ln_diaptr=T)
83      !!
84      !! ** Reference : Leonard (1979, 1991)
85      !!----------------------------------------------------------------------
86      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
87      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
88      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
89      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
90      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
91      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
92      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
94      !!----------------------------------------------------------------------
95      !
96      IF( kt == kit000 )  THEN
97         IF(lwp) WRITE(numout,*)
98         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
99         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
100         IF(lwp) WRITE(numout,*)
101      ENDIF
102      !
103      l_trd = .FALSE.
104      l_ptr = .FALSE.
105      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
106      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 
107      !
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_00_00( 1, jpkm1 )
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         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_00_00( 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_00_00( 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         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_00_00( 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         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_00_00
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         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions
191         !
192         ! Computation of the trend
193         DO_3D_00_00( 1, jpkm1 )
194            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
195            ! horizontal advective trends
196            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
197            !--- add it to the general tracer trends
198            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
199         END_3D
200         !                                 ! trend diagnostics
201         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
202         !
203      END DO
204      !
205   END SUBROUTINE tra_adv_qck_i
206
207
208   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
209      !!----------------------------------------------------------------------
210      !!
211      !!----------------------------------------------------------------------
212      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
213      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
214      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
215      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
216      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
217      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
218      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
219      !!
220      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
221      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
222      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
223      !----------------------------------------------------------------------
224      !
225      !                                                          ! ===========
226      DO jn = 1, kjpt                                            ! tracer loop
227         !                                                       ! ===========
228         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
229         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
230         !                                                 
231         DO jk = 1, jpkm1                               
232            !                                             
233            !--- Computation of the ustream and downstream value of the tracer and the mask
234            DO_2D_00_00
235               ! Upstream in the x-direction for the tracer
236               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
237               ! Downstream in the x-direction for the tracer
238               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
239            END_2D
240         END DO
241         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
242
243         
244         !
245         ! Horizontal advective fluxes
246         ! ---------------------------
247         !
248         DO_3D_00_00( 1, jpkm1 )
249            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
250            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
251         END_3D
252         !
253         DO_3D_00_00( 1, jpkm1 )
254            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
255            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
256            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
257            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
258            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
259         END_3D
260
261         !--- Lateral boundary conditions
262         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 )
263
264         !--- QUICKEST scheme
265         CALL quickest( zfu, zfd, zfc, zwy )
266         !
267         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
268         DO_3D_00_00( 1, jpkm1 )
269            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
270         END_3D
271         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
272         !
273         ! Tracer flux on the x-direction
274         DO jk = 1, jpkm1 
275            !
276            DO_2D_00_00
277               zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
278               !--- If the second ustream point is a land point
279               !--- the flux is computed by the 1st order UPWIND scheme
280               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
281               zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
282               zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
283            END_2D
284         END DO
285         !
286         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions
287         !
288         ! Computation of the trend
289         DO_3D_00_00( 1, jpkm1 )
290            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
291            ! horizontal advective trends
292            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
293            !--- add it to the general tracer trends
294            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
295         END_3D
296         !                                 ! trend diagnostics
297         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
298         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
299         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
300         !
301      END DO
302      !
303   END SUBROUTINE tra_adv_qck_j
304
305
306   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
307      !!----------------------------------------------------------------------
308      !!
309      !!----------------------------------------------------------------------
310      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
311      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
312      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
313      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
314      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
315      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
316      !
317      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
318      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
319      !!----------------------------------------------------------------------
320      !
321      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
322      zwz(:,:,jpk) = 0._wp
323      !
324      !                                                          ! ===========
325      DO jn = 1, kjpt                                            ! tracer loop
326         !                                                       ! ===========
327         !
328         DO_3D_00_00( 2, jpkm1 )
329            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)
330         END_3D
331         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
332            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
333               DO_2D_11_11
334                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
335               END_2D
336            ELSE                                   ! no ocean cavities (only ocean surface)
337               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm)
338            ENDIF
339         ENDIF
340         !
341         DO_3D_00_00( 1, jpkm1 )
342            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
343               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
344         END_3D
345         !                                 ! Send trends for diagnostic
346         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
347         !
348      END DO
349      !
350   END SUBROUTINE tra_adv_cen2_k
351
352
353   SUBROUTINE quickest( pfu, pfd, pfc, puc )
354      !!----------------------------------------------------------------------
355      !!
356      !! ** Purpose :  Computation of advective flux with Quickest scheme
357      !!
358      !! ** Method :   
359      !!----------------------------------------------------------------------
360      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
361      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
362      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
363      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
364      !!
365      INTEGER  ::  ji, jj, jk               ! dummy loop indices
366      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
367      REAL(wp) ::  zc, zcurv, zfho          !   -      -
368      !----------------------------------------------------------------------
369      !
370      DO_3D_11_11( 1, jpkm1 )
371         zc     = puc(ji,jj,jk)                         ! Courant number
372         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
373         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
374         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
375         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
376         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
377         !
378         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
379         zcoef2 = ABS( zcoef1 )
380         zcoef3 = ABS( zcurv )
381         IF( zcoef3 >= zcoef2 ) THEN
382            zfho = pfc(ji,jj,jk) 
383         ELSE
384            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
385            IF( zcoef1 >= 0. ) THEN
386               zfho = MAX( pfc(ji,jj,jk), zfho ) 
387               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
388            ELSE
389               zfho = MIN( pfc(ji,jj,jk), zfho ) 
390               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
391            ENDIF
392         ENDIF
393         puc(ji,jj,jk) = zfho
394      END_3D
395      !
396   END SUBROUTINE quickest
397
398   !!======================================================================
399END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.