source: NEMO/trunk/src/OCE/TRA/traadv_qck.F90 @ 11993

Last change on this file since 11993 was 11993, checked in by cetlod, 3 months ago

trunk : undo bad commit. Oups

  • Property svn:keywords set to Id
File size: 23.2 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   !
24   USE in_out_manager  ! I/O manager
25   USE lib_mpp         ! distribued memory computing
26   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
27   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_qck   ! routine called by step.F90
33
34   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
35
36   LOGICAL  ::   l_trd   ! flag to compute trends
37   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
38
39
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
50      &                                       ptb, ptn, pta, kjpt )
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 ptn
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 pta  with the now advective tracer trends
81      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
82      !!             - htr_adv, str_adv : 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   ) ::   kit000          ! first time step index
88      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
89      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
90      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
91      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
92      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
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. ln_diaptr )                                               l_ptr = .TRUE. 
107      !
108      !
109      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
110      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
111      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
112
113      !        ! vertical fluxes are computed with the 2nd order centered scheme
114      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
115      !
116   END SUBROUTINE tra_adv_qck
117
118
119   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
120      &                                        ptb, ptn, pta, kjpt   )
121      !!----------------------------------------------------------------------
122      !!
123      !!----------------------------------------------------------------------
124      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
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   ) ::   pun        ! i-velocity components
129      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
130      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
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 jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
145            DO jj = 2, jpjm1
146               DO ji = fs_2, fs_jpim1   ! vector opt.
147                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
148                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
149               END DO
150            END DO
151         END DO
152         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
153         
154         !
155         ! Horizontal advective fluxes
156         ! ---------------------------
157         DO jk = 1, jpkm1                             
158            DO jj = 2, jpjm1
159               DO ji = fs_2, fs_jpim1   ! vector opt.         
160                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
161                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
162               END DO
163            END DO
164         END DO
165         !
166         DO jk = 1, jpkm1 
167            DO jj = 2, jpjm1
168               DO ji = fs_2, fs_jpim1   ! vector opt.   
169                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
170                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
171                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
172                  zfc(ji,jj,jk)  = zdir * ptb(ji  ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn)  ! FC in the x-direction for T
173                  zfd(ji,jj,jk)  = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji  ,jj,jk,jn)  ! FD in the x-direction for T
174               END DO
175            END DO
176         END DO 
177         !--- Lateral boundary conditions
178         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. )
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 jk = 1, jpkm1 
185            DO jj = 2, jpjm1
186               DO ji = fs_2, fs_jpim1   ! vector opt.               
187                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
188               END DO
189            END DO
190         END DO
191         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
192
193         !
194         ! Tracer flux on the x-direction
195         DO jk = 1, jpkm1 
196            !
197            DO jj = 2, jpjm1
198               DO ji = fs_2, fs_jpim1   ! vector opt.               
199                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
200                  !--- If the second ustream point is a land point
201                  !--- the flux is computed by the 1st order UPWIND scheme
202                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
203                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
204                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
205               END DO
206            END DO
207         END DO
208         !
209         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
210         !
211         ! Computation of the trend
212         DO jk = 1, jpkm1 
213            DO jj = 2, jpjm1
214               DO ji = fs_2, fs_jpim1   ! vector opt. 
215                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
216                  ! horizontal advective trends
217                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
218                  !--- add it to the general tracer trends
219                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
220               END DO
221            END DO
222         END DO
223         !                                 ! trend diagnostics
224         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
225         !
226      END DO
227      !
228   END SUBROUTINE tra_adv_qck_i
229
230
231   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
232      &                                        ptb, ptn, pta, kjpt )
233      !!----------------------------------------------------------------------
234      !!
235      !!----------------------------------------------------------------------
236      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
237      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
238      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
239      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
240      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
241      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
242      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
243      !!
244      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
245      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
246      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
247      !----------------------------------------------------------------------
248      !
249      !                                                          ! ===========
250      DO jn = 1, kjpt                                            ! tracer loop
251         !                                                       ! ===========
252         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
253         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
254         !                                                 
255         DO jk = 1, jpkm1                               
256            !                                             
257            !--- Computation of the ustream and downstream value of the tracer and the mask
258            DO jj = 2, jpjm1
259               DO ji = fs_2, fs_jpim1   ! vector opt.
260                  ! Upstream in the x-direction for the tracer
261                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
262                  ! Downstream in the x-direction for the tracer
263                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
264               END DO
265            END DO
266         END DO
267         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
268
269         
270         !
271         ! Horizontal advective fluxes
272         ! ---------------------------
273         !
274         DO jk = 1, jpkm1                             
275            DO jj = 2, jpjm1
276               DO ji = fs_2, fs_jpim1   ! vector opt.         
277                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
278                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
279               END DO
280            END DO
281         END DO
282         !
283         DO jk = 1, jpkm1 
284            DO jj = 2, jpjm1
285               DO ji = fs_2, fs_jpim1   ! vector opt.   
286                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
287                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
288                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
289                  zfc(ji,jj,jk)  = zdir * ptb(ji,jj  ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn)  ! FC in the x-direction for T
290                  zfd(ji,jj,jk)  = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj  ,jk,jn)  ! FD in the x-direction for T
291               END DO
292            END DO
293         END DO
294
295         !--- Lateral boundary conditions
296         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
297
298         !--- QUICKEST scheme
299         CALL quickest( zfu, zfd, zfc, zwy )
300         !
301         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
302         DO jk = 1, jpkm1 
303            DO jj = 2, jpjm1
304               DO ji = fs_2, fs_jpim1   ! vector opt.               
305                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
306               END DO
307            END DO
308         END DO
309         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions
310         !
311         ! Tracer flux on the x-direction
312         DO jk = 1, jpkm1 
313            !
314            DO jj = 2, jpjm1
315               DO ji = fs_2, fs_jpim1   ! vector opt.               
316                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
317                  !--- If the second ustream point is a land point
318                  !--- the flux is computed by the 1st order UPWIND scheme
319                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
320                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
321                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
322               END DO
323            END DO
324         END DO
325         !
326         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
327         !
328         ! Computation of the trend
329         DO jk = 1, jpkm1 
330            DO jj = 2, jpjm1
331               DO ji = fs_2, fs_jpim1   ! vector opt. 
332                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
333                  ! horizontal advective trends
334                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
335                  !--- add it to the general tracer trends
336                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
337               END DO
338            END DO
339         END DO
340         !                                 ! trend diagnostics
341         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
342         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
343         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
344         !
345      END DO
346      !
347   END SUBROUTINE tra_adv_qck_j
348
349
350   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
351     &                                    ptn, pta, kjpt )
352      !!----------------------------------------------------------------------
353      !!
354      !!----------------------------------------------------------------------
355      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
356      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
357      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
358      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
359      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
360      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
361      !
362      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
363      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
364      !!----------------------------------------------------------------------
365      !
366      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
367      zwz(:,:,jpk) = 0._wp
368      !
369      !                                                          ! ===========
370      DO jn = 1, kjpt                                            ! tracer loop
371         !                                                       ! ===========
372         !
373         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
374            DO jj = 2, jpjm1
375               DO ji = fs_2, fs_jpim1   ! vector opt.
376                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
377               END DO
378            END DO
379         END DO
380         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
381            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
382               DO jj = 1, jpj
383                  DO ji = 1, jpi
384                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
385                  END DO
386               END DO   
387            ELSE                                   ! no ocean cavities (only ocean surface)
388               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
389            ENDIF
390         ENDIF
391         !
392         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
393            DO jj = 2, jpjm1
394               DO ji = fs_2, fs_jpim1   ! vector opt.
395                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
396                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
397               END DO
398            END DO
399         END DO
400         !                                 ! Send trends for diagnostic
401         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
402         !
403      END DO
404      !
405   END SUBROUTINE tra_adv_cen2_k
406
407
408   SUBROUTINE quickest( pfu, pfd, pfc, puc )
409      !!----------------------------------------------------------------------
410      !!
411      !! ** Purpose :  Computation of advective flux with Quickest scheme
412      !!
413      !! ** Method :   
414      !!----------------------------------------------------------------------
415      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
416      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
417      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
418      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
419      !!
420      INTEGER  ::  ji, jj, jk               ! dummy loop indices
421      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
422      REAL(wp) ::  zc, zcurv, zfho          !   -      -
423      !----------------------------------------------------------------------
424      !
425      DO jk = 1, jpkm1
426         DO jj = 1, jpj
427            DO ji = 1, jpi
428               zc     = puc(ji,jj,jk)                         ! Courant number
429               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
430               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
431               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
432               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
433               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
434               !
435               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
436               zcoef2 = ABS( zcoef1 )
437               zcoef3 = ABS( zcurv )
438               IF( zcoef3 >= zcoef2 ) THEN
439                  zfho = pfc(ji,jj,jk) 
440               ELSE
441                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
442                  IF( zcoef1 >= 0. ) THEN
443                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
444                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
445                  ELSE
446                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
447                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
448                  ENDIF
449               ENDIF
450               puc(ji,jj,jk) = zfho
451            END DO
452         END DO
453      END DO
454      !
455   END SUBROUTINE quickest
456
457   !!======================================================================
458END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.