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_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_qck.F90 @ 13551

Last change on this file since 13551 was 13551, checked in by hadcv, 4 years ago

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