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 @ 13906

Last change on this file since 13906 was 13906, checked in by mocavero, 3 years ago

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