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

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

Last change on this file since 11989 was 11989, checked in by cetlod, 4 years ago

cmip6 diag : finalisation of CMIP6 diags implementation
  • 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   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, pun, pvn, pwn,      &
51      &                                       ptb, ptn, pta, kjpt )
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 ptn
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 pta  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   ) ::   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   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
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, pun, ptb, ptn, pta, kjpt ) 
112      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
113
114      !        ! vertical fluxes are computed with the 2nd order centered scheme
115      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
116      !
117   END SUBROUTINE tra_adv_qck
118
119
120   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
121      &                                        ptb, ptn, pta, kjpt   )
122      !!----------------------------------------------------------------------
123      !!
124      !!----------------------------------------------------------------------
125      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
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   ) ::   pun        ! i-velocity components
130      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
131      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
132      !!
133      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
134      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
135      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zfu, zfc, zfd
136      !----------------------------------------------------------------------
137      !
138      !                                                          ! ===========
139      DO jn = 1, kjpt                                            ! tracer loop
140         !                                                       ! ===========
141         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
142         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
143         !
144!!gm why not using a SHIFT instruction...
145         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
149                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
150               END DO
151            END DO
152         END DO
153         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
154         
155         !
156         ! Horizontal advective fluxes
157         ! ---------------------------
158         DO jk = 1, jpkm1                             
159            DO jj = 2, jpjm1
160               DO ji = fs_2, fs_jpim1   ! vector opt.         
161                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
162                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
163               END DO
164            END DO
165         END DO
166         !
167         DO jk = 1, jpkm1 
168            DO jj = 2, jpjm1
169               DO ji = fs_2, fs_jpim1   ! vector opt.   
170                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
171                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
172                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
173                  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
174                  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
175               END DO
176            END DO
177         END DO 
178         !--- Lateral boundary conditions
179         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. )
180
181         !--- QUICKEST scheme
182         CALL quickest( zfu, zfd, zfc, zwx )
183         !
184         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
185         DO jk = 1, jpkm1 
186            DO jj = 2, jpjm1
187               DO ji = fs_2, fs_jpim1   ! vector opt.               
188                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
189               END DO
190            END DO
191         END DO
192         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
193
194         !
195         ! Tracer flux on the x-direction
196         DO jk = 1, jpkm1 
197            !
198            DO jj = 2, jpjm1
199               DO ji = fs_2, fs_jpim1   ! vector opt.               
200                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
201                  !--- If the second ustream point is a land point
202                  !--- the flux is computed by the 1st order UPWIND scheme
203                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
204                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
205                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
206               END DO
207            END DO
208         END DO
209         !
210         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
211         !
212         ! Computation of the trend
213         DO jk = 1, jpkm1 
214            DO jj = 2, jpjm1
215               DO ji = fs_2, fs_jpim1   ! vector opt. 
216                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
217                  ! horizontal advective trends
218                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
219                  !--- add it to the general tracer trends
220                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
221               END DO
222            END DO
223         END DO
224         !                                 ! trend diagnostics
225         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
226         !
227      END DO
228      !
229   END SUBROUTINE tra_adv_qck_i
230
231
232   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
233      &                                        ptb, ptn, pta, kjpt )
234      !!----------------------------------------------------------------------
235      !!
236      !!----------------------------------------------------------------------
237      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
238      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
239      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
240      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
241      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
242      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
243      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
244      !!
245      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
246      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
247      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
248      !----------------------------------------------------------------------
249      !
250      !                                                          ! ===========
251      DO jn = 1, kjpt                                            ! tracer loop
252         !                                                       ! ===========
253         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
254         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
255         !                                                 
256         DO jk = 1, jpkm1                               
257            !                                             
258            !--- Computation of the ustream and downstream value of the tracer and the mask
259            DO jj = 2, jpjm1
260               DO ji = fs_2, fs_jpim1   ! vector opt.
261                  ! Upstream in the x-direction for the tracer
262                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
263                  ! Downstream in the x-direction for the tracer
264                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
265               END DO
266            END DO
267         END DO
268         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
269
270         
271         !
272         ! Horizontal advective fluxes
273         ! ---------------------------
274         !
275         DO jk = 1, jpkm1                             
276            DO jj = 2, jpjm1
277               DO ji = fs_2, fs_jpim1   ! vector opt.         
278                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
279                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
280               END DO
281            END DO
282         END DO
283         !
284         DO jk = 1, jpkm1 
285            DO jj = 2, jpjm1
286               DO ji = fs_2, fs_jpim1   ! vector opt.   
287                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
288                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
289                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
290                  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
291                  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
292               END DO
293            END DO
294         END DO
295
296         !--- Lateral boundary conditions
297         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
298
299         !--- QUICKEST scheme
300         CALL quickest( zfu, zfd, zfc, zwy )
301         !
302         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
303         DO jk = 1, jpkm1 
304            DO jj = 2, jpjm1
305               DO ji = fs_2, fs_jpim1   ! vector opt.               
306                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
307               END DO
308            END DO
309         END DO
310         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions
311         !
312         ! Tracer flux on the x-direction
313         DO jk = 1, jpkm1 
314            !
315            DO jj = 2, jpjm1
316               DO ji = fs_2, fs_jpim1   ! vector opt.               
317                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
318                  !--- If the second ustream point is a land point
319                  !--- the flux is computed by the 1st order UPWIND scheme
320                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
321                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
322                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
323               END DO
324            END DO
325         END DO
326         !
327         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
328         !
329         ! Computation of the trend
330         DO jk = 1, jpkm1 
331            DO jj = 2, jpjm1
332               DO ji = fs_2, fs_jpim1   ! vector opt. 
333                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
334                  ! horizontal advective trends
335                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
336                  !--- add it to the general tracer trends
337                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
338               END DO
339            END DO
340         END DO
341         !                                 ! trend diagnostics
342         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
343         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
344         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
345         !
346      END DO
347      !
348   END SUBROUTINE tra_adv_qck_j
349
350
351   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
352     &                                    ptn, pta, kjpt )
353      !!----------------------------------------------------------------------
354      !!
355      !!----------------------------------------------------------------------
356      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
357      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
358      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
359      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
360      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
361      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
362      !
363      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
364      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
365      !!----------------------------------------------------------------------
366      !
367      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
368      zwz(:,:,jpk) = 0._wp
369      !
370      !                                                          ! ===========
371      DO jn = 1, kjpt                                            ! tracer loop
372         !                                                       ! ===========
373         !
374         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
375            DO jj = 2, jpjm1
376               DO ji = fs_2, fs_jpim1   ! vector opt.
377                  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)
378               END DO
379            END DO
380         END DO
381         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
382            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
383               DO jj = 1, jpj
384                  DO ji = 1, jpi
385                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
386                  END DO
387               END DO   
388            ELSE                                   ! no ocean cavities (only ocean surface)
389               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
390            ENDIF
391         ENDIF
392         !
393         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
394            DO jj = 2, jpjm1
395               DO ji = fs_2, fs_jpim1   ! vector opt.
396                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
397                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
398               END DO
399            END DO
400         END DO
401         !                                 ! Send trends for diagnostic
402         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
403         !
404      END DO
405      !
406   END SUBROUTINE tra_adv_cen2_k
407
408
409   SUBROUTINE quickest( pfu, pfd, pfc, puc )
410      !!----------------------------------------------------------------------
411      !!
412      !! ** Purpose :  Computation of advective flux with Quickest scheme
413      !!
414      !! ** Method :   
415      !!----------------------------------------------------------------------
416      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
417      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
418      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
419      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
420      !!
421      INTEGER  ::  ji, jj, jk               ! dummy loop indices
422      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
423      REAL(wp) ::  zc, zcurv, zfho          !   -      -
424      !----------------------------------------------------------------------
425      !
426      DO jk = 1, jpkm1
427         DO jj = 1, jpj
428            DO ji = 1, jpi
429               zc     = puc(ji,jj,jk)                         ! Courant number
430               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
431               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
432               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
433               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
434               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
435               !
436               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
437               zcoef2 = ABS( zcoef1 )
438               zcoef3 = ABS( zcurv )
439               IF( zcoef3 >= zcoef2 ) THEN
440                  zfho = pfc(ji,jj,jk) 
441               ELSE
442                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
443                  IF( zcoef1 >= 0. ) THEN
444                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
445                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
446                  ELSE
447                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
448                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
449                  ENDIF
450               ENDIF
451               puc(ji,jj,jk) = zfho
452            END DO
453         END DO
454      END DO
455      !
456   END SUBROUTINE quickest
457
458   !!======================================================================
459END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.