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/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/TRA/traadv_qck.F90 @ 14801

Last change on this file since 14801 was 14801, checked in by francesca, 3 years ago

add loop fusion to DYN and TRA modules - ticket #2607

  • 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#if defined key_loop_fusion
30   USE traadv_qck_lf   ! QCK    scheme            (tra_adv_qck  routine - loop fusion version)
31#endif
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_adv_qck   ! routine called by step.F90
37
38   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
39
40   LOGICAL  ::   l_trd   ! flag to compute trends
41   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
42
43
44   !! * Substitutions
45#  include "do_loop_substitute.h90"
46#  include "domzgr_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE tra_adv_qck  ***
57      !!
58      !! ** Purpose :   Compute the now trend due to the advection of tracers
59      !!      and add it to the general trend of passive tracer equations.
60      !!
61      !! ** Method :   The advection is evaluated by a third order scheme
62      !!             For a positive velocity u :              u(i)>0
63      !!                                          |--FU--|--FC--|--FD--|------|
64      !!                                             i-1    i      i+1   i+2
65      !!
66      !!             For a negative velocity u :              u(i)<0
67      !!                                          |------|--FD--|--FC--|--FU--|
68      !!                                             i-1    i      i+1   i+2
69      !!             where  FU is the second upwind point
70      !!                    FD is the first douwning point
71      !!                    FC is the central point (or the first upwind point)
72      !!
73      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
74      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
75      !!
76      !!         dt = 2*rdtra and the scalar values are tb and sb
77      !!
78      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm)
79      !!
80      !!               The fluxes are bounded by the ULTIMATE limiter to
81      !!             guarantee the monotonicity of the solution and to
82      !!            prevent the appearance of spurious numerical oscillations
83      !!
84      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
85      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
86      !!             - poleward advective heat and salt transport (ln_diaptr=T)
87      !!
88      !! ** Reference : Leonard (1979, 1991)
89      !!----------------------------------------------------------------------
90      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
91      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
92      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
93      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
94      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
95      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
96      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
97      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
98      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
99      !!----------------------------------------------------------------------
100      !
101#if defined key_loop_fusion
102      CALL tra_adv_qck_lf ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
103#else
104      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
105         IF( kt == kit000 )  THEN
106            IF(lwp) WRITE(numout,*)
107            IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
108            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
109            IF(lwp) WRITE(numout,*)
110         ENDIF
111         !
112         l_trd = .FALSE.
113         l_ptr = .FALSE.
114         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
115         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
116      ENDIF
117      !
118      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
119      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
120      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
121
122      !        ! vertical fluxes are computed with the 2nd order centered scheme
123      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
124      !
125#endif
126   END SUBROUTINE tra_adv_qck
127
128
129   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
130      !!----------------------------------------------------------------------
131      !!
132      !!----------------------------------------------------------------------
133      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
134      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
135      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
136      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
137      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
138      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
139      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components
140      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
141      !!
142      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
143      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
144      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwx, zfu, zfc, zfd
145      !----------------------------------------------------------------------
146      !
147      !                                                          ! ===========
148      DO jn = 1, kjpt                                            ! tracer loop
149         !                                                       ! ===========
150         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp
151         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp
152         !
153!!gm why not using a SHIFT instruction...
154         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
155            zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
156            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
157         END_3D
158         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary conditions
159
160         !
161         ! Horizontal advective fluxes
162         ! ---------------------------
163         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
164            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
165            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
166         END_3D
167         !
168         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
169            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
170            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
171            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
172            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
173            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
174         END_3D
175         !--- Lateral boundary conditions
176         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
177
178         !--- QUICKEST scheme
179         CALL quickest( zfu, zfd, zfc, zwx )
180         !
181         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
182         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )
183            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
184         END_3D
185         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )      ! Lateral boundary conditions
186
187         !
188         ! Tracer flux on the x-direction
189         DO_3D( 1, 0, 0, 0, 1, jpkm1 )
190            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
191            !--- If the second ustream point is a land point
192            !--- the flux is computed by the 1st order UPWIND scheme
193            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
194            zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
195            zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
196         END_3D
197         !
198         ! Computation of the trend
199         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
200            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
201            ! horizontal advective trends
202            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
203            !--- add it to the general tracer trends
204            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
205         END_3D
206         !                                 ! trend diagnostics
207         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
208         !
209      END DO
210      !
211   END SUBROUTINE tra_adv_qck_i
212
213
214   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
215      !!----------------------------------------------------------------------
216      !!
217      !!----------------------------------------------------------------------
218      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
219      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
220      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
221      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
222      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
223      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
224      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
225      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
226      !!
227      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
228      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
229      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
230      !----------------------------------------------------------------------
231      !
232      !                                                          ! ===========
233      DO jn = 1, kjpt                                            ! tracer loop
234         !                                                       ! ===========
235         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0
236         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0
237         !
238         DO jk = 1, jpkm1
239            !
240            !--- Computation of the ustream and downstream value of the tracer and the mask
241            DO_2D( 0, 0, nn_hls-1, nn_hls-1 )
242               ! Upstream in the x-direction for the tracer
243               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
244               ! Downstream in the x-direction for the tracer
245               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
246            END_2D
247         END DO
248         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary conditions
249
250         !
251         ! Horizontal advective fluxes
252         ! ---------------------------
253         !
254         DO_3D( 0, 0, nn_hls-1, 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, nn_hls-1, 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         IF (nn_hls==1) CALL lbc_lnk( '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, nn_hls-1, nn_hls-1, 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         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )    !--- Lateral boundary conditions
278         !
279         ! Tracer flux on the x-direction
280         DO_3D( 0, 0, 1, 0, 1, jpkm1 )
281            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
282            !--- If the second ustream point is a land point
283            !--- the flux is computed by the 1st order UPWIND scheme
284            zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
285            zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
286            zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
287         END_3D
288         !
289         ! Computation of the trend
290         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
291            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
292            ! horizontal advective trends
293            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
294            !--- add it to the general tracer trends
295            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
296         END_3D
297         !                                 ! trend diagnostics
298         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
299         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
300         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
301         !
302      END DO
303      !
304   END SUBROUTINE tra_adv_qck_j
305
306
307   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
308      !!----------------------------------------------------------------------
309      !!
310      !!----------------------------------------------------------------------
311      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
312      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
313      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
314      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
315      ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
316      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
317      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
318      !
319      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
320      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwz   ! 3D workspace
321      !!----------------------------------------------------------------------
322      !
323      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
324      zwz(:,:,jpk) = 0._wp
325      !
326      !                                                          ! ===========
327      DO jn = 1, kjpt                                            ! tracer loop
328         !                                                       ! ===========
329         !
330         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
331            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)
332         END_3D
333         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
334            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
335               DO_2D( 0, 0, 0, 0 )
336                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
337               END_2D
338            ELSE                                   ! no ocean cavities (only ocean surface)
339               DO_2D( 0, 0, 0, 0 )
340                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
341               END_2D
342            ENDIF
343         ENDIF
344         !
345         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
346            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
347               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
348         END_3D
349         !                                 ! Send trends for diagnostic
350         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
351         !
352      END DO
353      !
354   END SUBROUTINE tra_adv_cen2_k
355
356
357   SUBROUTINE quickest( pfu, pfd, pfc, puc )
358      !!----------------------------------------------------------------------
359      !!
360      !! ** Purpose :  Computation of advective flux with Quickest scheme
361      !!
362      !! ** Method :
363      !!----------------------------------------------------------------------
364      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfu   ! second upwind point
365      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfd   ! first douwning point
366      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
367      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
368      !!
369      INTEGER  ::  ji, jj, jk               ! dummy loop indices
370      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars
371      REAL(wp) ::  zc, zcurv, zfho          !   -      -
372      !----------------------------------------------------------------------
373      !
374      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
375         zc     = puc(ji,jj,jk)                         ! Courant number
376         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
377         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
378         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
379         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
380         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
381         !
382         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
383         zcoef2 = ABS( zcoef1 )
384         zcoef3 = ABS( zcurv )
385         IF( zcoef3 >= zcoef2 ) THEN
386            zfho = pfc(ji,jj,jk)
387         ELSE
388            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
389            IF( zcoef1 >= 0. ) THEN
390               zfho = MAX( pfc(ji,jj,jk), zfho )
391               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
392            ELSE
393               zfho = MIN( pfc(ji,jj,jk), zfho )
394               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
395            ENDIF
396         ENDIF
397         puc(ji,jj,jk) = zfho
398      END_3D
399      !
400   END SUBROUTINE quickest
401
402   !!======================================================================
403END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.