source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

  • Property svn:keywords set to Id
File size: 24.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!!gm   USE trabbl          ! advective term in the BBL
23   USE dynspg_oce      ! surface pressure gradient variables
24   USE diaptr          ! poleward transport diagnostics
25   !
26   USE lib_mpp         ! distribued memory computing
27   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
28   USE in_out_manager  ! I/O manager
29   USE wrk_nemo        ! Memory Allocation
30   USE timing          ! Timing
31   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_adv_qck   ! routine called by step.F90
37
38   LOGICAL  :: l_trd           ! flag to compute trends
39   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
52      &                                       ptb, ptn, pta, kjpt )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_adv_qck  ***
55      !!
56      !! ** Purpose :   Compute the now trend due to the advection of tracers
57      !!      and add it to the general trend of passive tracer equations.
58      !!
59      !! ** Method :   The advection is evaluated by a third order scheme
60      !!             For a positive velocity u :              u(i)>0
61      !!                                          |--FU--|--FC--|--FD--|------|
62      !!                                             i-1    i      i+1   i+2
63      !!
64      !!             For a negative velocity u :              u(i)<0
65      !!                                          |------|--FD--|--FC--|--FU--|
66      !!                                             i-1    i      i+1   i+2
67      !!             where  FU is the second upwind point
68      !!                    FD is the first douwning point
69      !!                    FC is the central point (or the first upwind point)
70      !!
71      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
72      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
73      !!
74      !!         dt = 2*rdtra and the scalar values are tb and sb
75      !!
76      !!       On the vertical, the simple centered scheme used ptn
77      !!
78      !!               The fluxes are bounded by the ULTIMATE limiter to
79      !!             guarantee the monotonicity of the solution and to
80      !!            prevent the appearance of spurious numerical oscillations
81      !!
82      !! ** Action : - update (pta) with the now advective tracer trends
83      !!             - save the trends
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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of 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      !
98      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
99      !
100      IF( kt == kit000 )  THEN
101         IF(lwp) WRITE(numout,*)
102         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
103         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
104         IF(lwp) WRITE(numout,*)
105         !
106         l_trd = .FALSE.
107         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
108      ENDIF
109
110      ! I. The 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      ! II. The vertical fluxes are computed with the 2nd order centered scheme
115      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
116      !
117      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
118      !
119   END SUBROUTINE tra_adv_qck
120
121
122   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
123      &                                        ptb, ptn, pta, kjpt   )
124      !!----------------------------------------------------------------------
125      !!
126      !!----------------------------------------------------------------------
127      USE oce     , ONLY:   zwx => ua       ! ua used as workspace
128      !
129      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
130      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
131      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
132      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
133      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
134      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
135      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
136      !!
137      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
138      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
139      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
140      !----------------------------------------------------------------------
141      !
142      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
143      !                                                          ! ===========
144      DO jn = 1, kjpt                                            ! tracer loop
145         !                                                       ! ===========
146         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
147         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
148         !                                                 
149         DO jk = 1, jpkm1                               
150            !                                             
151            !--- Computation of the ustream and downstream value of the tracer and the mask
152            DO jj = 2, jpjm1
153               DO ji = fs_2, fs_jpim1   ! vector opt.
154                  ! Upstream in the x-direction for the tracer
155                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
156                  ! Downstream in the x-direction for the tracer
157                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
158               END DO
159            END DO
160         END DO
161         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
162         
163         !
164         ! Horizontal advective fluxes
165         ! ---------------------------
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                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
172               END DO
173            END DO
174         END DO
175         !
176         DO jk = 1, jpkm1 
177            zdt =  p2dt(jk)
178            DO jj = 2, jpjm1
179               DO ji = fs_2, fs_jpim1   ! vector opt.   
180                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
181                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
182                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
183                  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
184                  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
185               END DO
186            END DO
187         END DO 
188         !--- Lateral boundary conditions
189         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
190         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
191
192         !--- QUICKEST scheme
193         CALL quickest( zfu, zfd, zfc, zwx )
194         !
195         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
196         DO jk = 1, jpkm1 
197            DO jj = 2, jpjm1
198               DO ji = fs_2, fs_jpim1   ! vector opt.               
199                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
200               END DO
201            END DO
202         END DO
203         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
204
205         !
206         ! Tracer flux on the x-direction
207         DO jk = 1, jpkm1 
208            !
209            DO jj = 2, jpjm1
210               DO ji = fs_2, fs_jpim1   ! vector opt.               
211                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
212                  !--- If the second ustream point is a land point
213                  !--- the flux is computed by the 1st order UPWIND scheme
214                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
215                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
216                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
217               END DO
218            END DO
219         END DO
220         !
221         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
222         !
223         ! Computation of the trend
224         DO jk = 1, jpkm1 
225            DO jj = 2, jpjm1
226               DO ji = fs_2, fs_jpim1   ! vector opt. 
227                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
228                  ! horizontal advective trends
229                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
230                  !--- add it to the general tracer trends
231                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
232               END DO
233            END DO
234         END DO
235         !                                 ! trend diagnostics (contribution of upstream fluxes)
236         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
237         !
238      END DO
239      !
240      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
241      !
242   END SUBROUTINE tra_adv_qck_i
243
244
245   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
246      &                                        ptb, ptn, pta, kjpt )
247      !!----------------------------------------------------------------------
248      !!
249      !!----------------------------------------------------------------------
250      USE oce     , ONLY:   zwy => ua       ! ua used as workspace
251      !
252      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
253      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
254      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
255      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
256      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
257      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
258      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
259      !!
260      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
261      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
262      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
263      !----------------------------------------------------------------------
264      !
265      CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
266      !
267      !                                                          ! ===========
268      DO jn = 1, kjpt                                            ! tracer loop
269         !                                                       ! ===========
270         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
271         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
272         !                                                 
273         DO jk = 1, jpkm1                               
274            !                                             
275            !--- Computation of the ustream and downstream value of the tracer and the mask
276            DO jj = 2, jpjm1
277               DO ji = fs_2, fs_jpim1   ! vector opt.
278                  ! Upstream in the x-direction for the tracer
279                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
280                  ! Downstream in the x-direction for the tracer
281                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
282               END DO
283            END DO
284         END DO
285         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
286
287         
288         !
289         ! Horizontal advective fluxes
290         ! ---------------------------
291         !
292         DO jk = 1, jpkm1                             
293            DO jj = 2, jpjm1
294               DO ji = fs_2, fs_jpim1   ! vector opt.         
295                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
296                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
297               END DO
298            END DO
299         END DO
300         !
301         DO jk = 1, jpkm1 
302            zdt =  p2dt(jk)
303            DO jj = 2, jpjm1
304               DO ji = fs_2, fs_jpim1   ! vector opt.   
305                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
306                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
307                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
308                  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
309                  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
310               END DO
311            END DO
312         END DO
313
314         !--- Lateral boundary conditions
315         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
316         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
317
318         !--- QUICKEST scheme
319         CALL quickest( zfu, zfd, zfc, zwy )
320         !
321         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
322         DO jk = 1, jpkm1 
323            DO jj = 2, jpjm1
324               DO ji = fs_2, fs_jpim1   ! vector opt.               
325                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
326               END DO
327            END DO
328         END DO
329         !--- Lateral boundary conditions
330         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
331         !
332         ! Tracer flux on the x-direction
333         DO jk = 1, jpkm1 
334            !
335            DO jj = 2, jpjm1
336               DO ji = fs_2, fs_jpim1   ! vector opt.               
337                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
338                  !--- If the second ustream point is a land point
339                  !--- the flux is computed by the 1st order UPWIND scheme
340                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
341                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
342                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
343               END DO
344            END DO
345         END DO
346         !
347         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
348         !
349         ! Computation of the trend
350         DO jk = 1, jpkm1 
351            DO jj = 2, jpjm1
352               DO ji = fs_2, fs_jpim1   ! vector opt. 
353                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
354                  ! horizontal advective trends
355                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
356                  !--- add it to the general tracer trends
357                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
358               END DO
359            END DO
360         END DO
361         !                                 ! trend diagnostics (contribution of upstream fluxes)
362         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
363         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
364         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
365           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
366           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
367         ENDIF
368         !
369      END DO
370      !
371      CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
372      !
373   END SUBROUTINE tra_adv_qck_j
374
375
376   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
377     &                                    ptn, pta, kjpt )
378      !!----------------------------------------------------------------------
379      !!
380      !!----------------------------------------------------------------------
381      USE oce, ONLY:   zwz => ua   ! ua used as workspace
382      !
383      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
384      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
385      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
386      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
387      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
388      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
389      !
390      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
391      REAL(wp) ::   zbtr , ztra      ! local scalars
392      !!----------------------------------------------------------------------
393
394      !                                                          ! ===========
395      DO jn = 1, kjpt                                            ! tracer loop
396         !                                                       ! ===========
397         ! 1. Bottom value : flux set to zero
398         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
399         !
400         !                                 ! Surface value
401         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
402         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
403         ENDIF
404         !
405         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
406            DO jj = 2, jpjm1
407               DO ji = fs_2, fs_jpim1   ! vector opt.
408                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
409               END DO
410            END DO
411         END DO
412         !
413         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
414            DO jj = 2, jpjm1
415               DO ji = fs_2, fs_jpim1   ! vector opt.
416                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
417                  ! k- vertical advective trends
418                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
419                  ! added to the general tracer trends
420                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
421               END DO
422            END DO
423         END DO
424         !                                 ! Save the vertical advective trends for diagnostic
425         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
426         !
427      END DO
428      !
429   END SUBROUTINE tra_adv_cen2_k
430
431
432   SUBROUTINE quickest( pfu, pfd, pfc, puc )
433      !!----------------------------------------------------------------------
434      !!
435      !! ** Purpose :  Computation of advective flux with Quickest scheme
436      !!
437      !! ** Method :   
438      !!----------------------------------------------------------------------
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
440      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
441      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
442      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
443      !!
444      INTEGER  ::  ji, jj, jk               ! dummy loop indices
445      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
446      REAL(wp) ::  zc, zcurv, zfho          !   -      -
447      !----------------------------------------------------------------------
448      !
449      IF( nn_timing == 1 )  CALL timing_start('quickest')
450      !
451      DO jk = 1, jpkm1
452         DO jj = 1, jpj
453            DO ji = 1, jpi
454               zc     = puc(ji,jj,jk)                         ! Courant number
455               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
456               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
457               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
458               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
459               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
460               !
461               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
462               zcoef2 = ABS( zcoef1 )
463               zcoef3 = ABS( zcurv )
464               IF( zcoef3 >= zcoef2 ) THEN
465                  zfho = pfc(ji,jj,jk) 
466               ELSE
467                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
468                  IF( zcoef1 >= 0. ) THEN
469                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
470                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
471                  ELSE
472                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
473                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
474                  ENDIF
475               ENDIF
476               puc(ji,jj,jk) = zfho
477            END DO
478         END DO
479      END DO
480      !
481      IF( nn_timing == 1 )  CALL timing_stop('quickest')
482      !
483   END SUBROUTINE quickest
484
485   !!======================================================================
486END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.