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_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA/traadv_qck.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

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