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_r13898_Tiling_Cleanup_MPI3/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/OCE/TRA/traadv_qck.F90 @ 13901

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

#2365: Merge in changes from dev_r13508_HPC-09_communications_cleanup up to [13701]
  • Property svn:keywords set to Id
File size: 21.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      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
94      ! NOTE: [tiling-comms-merge] These were changed to INTENT(inout) but they are not modified, so it is reverted
95      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
96      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
97      !!----------------------------------------------------------------------
98      !
99      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
100         IF( kt == kit000 )  THEN
101            IF(lwp) WRITE(numout,*)
102            IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
103            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
104            IF(lwp) WRITE(numout,*)
105         ENDIF
106         !
107         l_trd = .FALSE.
108         l_ptr = .FALSE.
109         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
110         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
111      ENDIF
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: [tiling] This can be A2D(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(A2D(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, nn_hls-1, nn_hls-1, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
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         IF (nn_hls.EQ.1) 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, nn_hls-1, 0, 1, jpkm1 )
158            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
159            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, nn_hls-1, 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         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
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, nn_hls-1, nn_hls-1, 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         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions
180
181         !
182         ! Tracer flux on the x-direction
183         DO_3D( 0, 0, 1, 0, 1, jpkm1 )
184            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
185            !--- If the second ustream point is a land point
186            !--- the flux is computed by the 1st order UPWIND scheme
187            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
188            zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
189            zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
190         END_3D
191         !
192         ! Computation of the trend
193         DO_3D( 0, 0, 0, 0, 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      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
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(A2D(nn_hls),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( nn_hls-1, nn_hls-1, 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         IF (nn_hls.EQ.1) 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( nn_hls-1, 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( nn_hls-1, 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         IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
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( nn_hls-1, nn_hls-1, 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         IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
273         !
274         ! Tracer flux on the x-direction
275         DO_3D( 1, 0, 0, 0, 1, jpkm1 )
276            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
277            !--- If the second ustream point is a land point
278            !--- the flux is computed by the 1st order UPWIND scheme
279            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
280            zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
281            zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
282         END_3D
283         !
284         ! Computation of the trend
285         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
286            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
287            ! horizontal advective trends
288            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
289            !--- add it to the general tracer trends
290            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
291         END_3D
292         !                                 ! trend diagnostics
293         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
294         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
295         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
296         !
297      END DO
298      !
299   END SUBROUTINE tra_adv_qck_j
300
301
302   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
303      !!----------------------------------------------------------------------
304      !!
305      !!----------------------------------------------------------------------
306      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
307      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
308      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
309      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
310      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
311      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
312      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
313      !
314      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
315      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwz   ! 3D workspace
316      !!----------------------------------------------------------------------
317      !
318      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
319      zwz(:,:,jpk) = 0._wp
320      !
321      !                                                          ! ===========
322      DO jn = 1, kjpt                                            ! tracer loop
323         !                                                       ! ===========
324         !
325         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
326            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)
327         END_3D
328         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
329            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
330               DO_2D( 0, 0, 0, 0 )
331                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
332               END_2D
333            ELSE                                   ! no ocean cavities (only ocean surface)
334               DO_2D( 0, 0, 0, 0 )
335                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
336               END_2D
337            ENDIF
338         ENDIF
339         !
340         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
341            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
342               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
343         END_3D
344         !                                 ! Send trends for diagnostic
345         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
346         !
347      END DO
348      !
349   END SUBROUTINE tra_adv_cen2_k
350
351
352   SUBROUTINE quickest( pfu, pfd, pfc, puc )
353      !!----------------------------------------------------------------------
354      !!
355      !! ** Purpose :  Computation of advective flux with Quickest scheme
356      !!
357      !! ** Method :   
358      !!----------------------------------------------------------------------
359      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfu   ! second upwind point
360      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfd   ! first douwning point
361      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
362      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
363      !!
364      INTEGER  ::  ji, jj, jk               ! dummy loop indices
365      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
366      REAL(wp) ::  zc, zcurv, zfho          !   -      -
367      !----------------------------------------------------------------------
368      !
369      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
370         zc     = puc(ji,jj,jk)                         ! Courant number
371         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
372         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
373         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
374         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
375         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
376         !
377         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
378         zcoef2 = ABS( zcoef1 )
379         zcoef3 = ABS( zcurv )
380         IF( zcoef3 >= zcoef2 ) THEN
381            zfho = pfc(ji,jj,jk) 
382         ELSE
383            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
384            IF( zcoef1 >= 0. ) THEN
385               zfho = MAX( pfc(ji,jj,jk), zfho ) 
386               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
387            ELSE
388               zfho = MIN( pfc(ji,jj,jk), zfho ) 
389               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
390            ENDIF
391         ENDIF
392         puc(ji,jj,jk) = zfho
393      END_3D
394      !
395   END SUBROUTINE quickest
396
397   !!======================================================================
398END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.