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 branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2598

Last change on this file since 2598 was 2590, checked in by trackstand2, 13 years ago

Merge branch 'dynamic_memory' into master-svn-dyn
  • Property svn:keywords set to Id
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 trdmod_oce      ! ocean space and time domain
20   USE trdtra          ! ocean tracers trends
21   USE trabbl          ! advective term in the BBL
22   USE lib_mpp         ! distribued memory computing
23   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
24   USE dynspg_oce      ! surface pressure gradient variables
25   USE in_out_manager  ! I/O manager
26   USE diaptr          ! poleward transport diagnostics
27   USE trc_oce         ! share passive tracers/Ocean variables
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_qck   ! routine called by step.F90
33
34   LOGICAL  :: l_trd           ! flag to compute trends
35   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn,      &
48      &                                       ptb, ptn, pta, kjpt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE tra_adv_qck  ***
51      !!
52      !! ** Purpose :   Compute the now trend due to the advection of tracers
53      !!      and add it to the general trend of passive tracer equations.
54      !!
55      !! ** Method :   The advection is evaluated by a third order scheme
56      !!             For a positive velocity u :              u(i)>0
57      !!                                          |--FU--|--FC--|--FD--|------|
58      !!                                             i-1    i      i+1   i+2
59      !!
60      !!             For a negative velocity u :              u(i)<0
61      !!                                          |------|--FD--|--FC--|--FU--|
62      !!                                             i-1    i      i+1   i+2
63      !!             where  FU is the second upwind point
64      !!                    FD is the first douwning point
65      !!                    FC is the central point (or the first upwind point)
66      !!
67      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
68      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
69      !!
70      !!         dt = 2*rdtra and the scalar values are tb and sb
71      !!
72      !!       On the vertical, the simple centered scheme used ptn
73      !!
74      !!               The fluxes are bounded by the ULTIMATE limiter to
75      !!             guarantee the monotonicity of the solution and to
76      !!            prevent the appearance of spurious numerical oscillations
77      !!
78      !! ** Action : - update (pta) with the now advective tracer trends
79      !!             - save the trends
80      !!
81      !! ** Reference : Leonard (1979, 1991)
82      !!----------------------------------------------------------------------
83      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
84      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
85      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
86      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
87      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
88      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
89      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
90      !!----------------------------------------------------------------------
91
92      IF( kt == nit000 )  THEN
93         IF(lwp) WRITE(numout,*)
94         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
95         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
96         IF(lwp) WRITE(numout,*)
97         !
98         l_trd = .FALSE.
99         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
100      ENDIF
101
102      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
103      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
104      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
105
106      ! II. The vertical fluxes are computed with the 2nd order centered scheme
107      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
108      !
109   END SUBROUTINE tra_adv_qck
110
111
112   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
113      &                                        ptb, ptn, pta, kjpt   )
114      !!----------------------------------------------------------------------
115      !!
116      !!----------------------------------------------------------------------
117      USE oce         , zwx => ua   ! use ua as workspace
118      USE wrk_nemo, ONLY: wrk_use, wrk_release
119      USE wrk_nemo, ONLY: zfu => wrk_3d_1, zfc => wrk_3d_2, zfd => wrk_3d_3
120      !!
121      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
122      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
123      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
124      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
125      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun             ! i-velocity components
126      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
127      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
128      !!
129      INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
130      REAL(wp) :: ztra, zbtr               ! local scalars
131      REAL(wp) :: zdir, zdx, zdt, zmsk     ! local scalars
132      !----------------------------------------------------------------------
133      !
134      IF(.not. wrk_use(3, 1,2,3))THEN
135         CALL ctl_stop('tra_adv_qck_i: ERROR: requested workspace arrays unavailable')
136         RETURN
137      END IF
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               ENDDO
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            !
215            ! Computation of the trend
216            DO jj = 2, jpjm1
217               DO ji = fs_2, fs_jpim1   ! vector opt. 
218                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
219                  ! horizontal advective trends
220                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
221                  !--- add it to the general tracer trends
222                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
223               END DO
224            END DO
225            !
226         END DO
227         !                                 ! trend diagnostics (contribution of upstream fluxes)
228         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
229         !
230      END DO
231      !
232      IF(.not. wrk_release(3, 1,2,3))THEN
233         CALL ctl_stop('tra_adv_qck_i: ERROR: failed to release workspace arrays')
234      END IF
235      !
236   END SUBROUTINE tra_adv_qck_i
237
238
239   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
240      &                                        ptb, ptn, pta, kjpt )
241      !!----------------------------------------------------------------------
242      !!
243      !!----------------------------------------------------------------------
244      USE oce         , zwy => ua   ! use ua as workspace
245      USE wrk_nemo, ONLY: wrk_use, wrk_release
246      USE wrk_nemo, ONLY: zfu => wrk_3d_1, zfc => wrk_3d_2, zfd => wrk_3d_3
247      !!
248      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
249      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
250      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
251      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
252      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn             ! j-velocity components
253      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
254      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
255      !!
256      INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
257      REAL(wp) :: ztra, zbtr               ! local scalars
258      REAL(wp) :: zdir, zdx, zdt, zmsk     ! local scalars
259      !----------------------------------------------------------------------
260      !
261      IF(.not. wrk_use(3, 1,2,3))THEN
262         CALL ctl_stop('tra_adv_qck_j: ERROR: requested workspace arrays unavailable')
263         RETURN
264      END IF
265      !                                                          ! ===========
266      DO jn = 1, kjpt                                            ! tracer loop
267         !                                                       ! ===========
268         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
269         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
270         !                                                 
271         DO jk = 1, jpkm1                               
272            !                                             
273            !--- Computation of the ustream and downstream value of the tracer and the mask
274            DO jj = 2, jpjm1
275               DO ji = fs_2, fs_jpim1   ! vector opt.
276                  ! Upstream in the x-direction for the tracer
277                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
278                  ! Downstream in the x-direction for the tracer
279                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
280               END DO
281            END DO
282         END DO
283         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
284
285         
286         !
287         ! Horizontal advective fluxes
288         ! ---------------------------
289         !
290         DO jk = 1, jpkm1                             
291            DO jj = 2, jpjm1
292               DO ji = fs_2, fs_jpim1   ! vector opt.         
293                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
294                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
295               END DO
296            END DO
297         END DO
298         !
299         DO jk = 1, jpkm1 
300            zdt =  p2dt(jk)
301            DO jj = 2, jpjm1
302               DO ji = fs_2, fs_jpim1   ! vector opt.   
303                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
304                  zdx  = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
305                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
306                  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
307                  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
308               END DO
309            END DO
310         END DO
311
312         !--- Lateral boundary conditions
313         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
314         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
315
316         !--- QUICKEST scheme
317         CALL quickest( zfu, zfd, zfc, zwy )
318         !
319         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
320         DO jk = 1, jpkm1 
321            DO jj = 2, jpjm1
322               DO ji = fs_2, fs_jpim1   ! vector opt.               
323                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
324               END DO
325            END DO
326         END DO
327         !--- Lateral boundary conditions
328         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
329         !
330         ! Tracer flux on the x-direction
331         DO jk = 1, jpkm1 
332            !
333            DO jj = 2, jpjm1
334               DO ji = fs_2, fs_jpim1   ! vector opt.               
335                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
336                  !--- If the second ustream point is a land point
337                  !--- the flux is computed by the 1st order UPWIND scheme
338                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
339                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
340                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
341               END DO
342            END DO
343            !
344            ! Computation of the trend
345            DO jj = 2, jpjm1
346               DO ji = fs_2, fs_jpim1   ! vector opt. 
347                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
348                  ! horizontal advective trends
349                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
350                  !--- add it to the general tracer trends
351                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
352               END DO
353            END DO
354            !
355         END DO
356         !                                 ! trend diagnostics (contribution of upstream fluxes)
357         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
358         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
359         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
360           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
361           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
362         ENDIF
363         !
364      END DO
365      !
366      IF(.not. wrk_release(3, 1,2,3))THEN
367         CALL ctl_stop('tra_adv_qck_j: ERROR: failed to release workspace arrays')
368      END IF
369      !
370   END SUBROUTINE tra_adv_qck_j
371
372
373   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
374     &                                    ptn, pta, kjpt )
375      !!----------------------------------------------------------------------
376      !!
377      !!----------------------------------------------------------------------
378      USE oce         , zwz => ua   ! use ua as workspace
379      !!
380      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
381      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
382      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
383      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn             ! vertical velocity
384      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn             ! before and now tracer fields
385      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
386      !!
387      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
388      REAL(wp) ::   zbtr , ztra      ! temporary scalars
389      !!----------------------------------------------------------------------
390
391      !                                                          ! ===========
392      DO jn = 1, kjpt                                            ! tracer loop
393         !                                                       ! ===========
394         ! 1. Bottom value : flux set to zero
395         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
396         !
397         !                                 ! Surface value
398         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
399         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
400         ENDIF
401         !
402         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
403            DO jj = 2, jpjm1
404               DO ji = fs_2, fs_jpim1   ! vector opt.
405                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
406               END DO
407            END DO
408         END DO
409         !
410         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
411            DO jj = 2, jpjm1
412               DO ji = fs_2, fs_jpim1   ! vector opt.
413                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
414                  ! k- vertical advective trends
415                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
416                  ! added to the general tracer trends
417                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
418               END DO
419            END DO
420         END DO
421         !                                 ! Save the vertical advective trends for diagnostic
422         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
423         !
424      END DO
425      !
426   END SUBROUTINE tra_adv_cen2_k
427
428
429   SUBROUTINE quickest( pfu, pfd, pfc, puc )
430      !!----------------------------------------------------------------------
431      !!
432      !! ** Purpose :  Computation of advective flux with Quickest scheme
433      !!
434      !! ** Method :   
435      !!----------------------------------------------------------------------
436      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
437      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
438      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
440      !!
441      INTEGER  ::  ji, jj, jk               ! dummy loop indices
442      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
443      REAL(wp) ::  zc, zcurv, zfho          !   -      -
444      !----------------------------------------------------------------------
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   END SUBROUTINE quickest
477
478   !!======================================================================
479END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.