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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_qck.F90 @ 12205

Last change on this file since 12205 was 12193, checked in by davestorkey, 4 years ago

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