source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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