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 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation
  • Property svn:keywords set to Id
File size: 23.9 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 wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
118      USE oce     , ONLY:   zwx => ua       ! ua used as workspace
119      USE wrk_nemo, ONLY:   zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3   ! 3D workspace
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, zdir, zdx, zdt, zmsk   ! local scalars
131      !----------------------------------------------------------------------
132      !
133      IF( wrk_in_use(3, 1,2,3) ) THEN
134         CALL ctl_stop('tra_adv_qck_i: requested workspace arrays unavailable')   ;   RETURN
135      ENDIF
136      !                                                          ! ===========
137      DO jn = 1, kjpt                                            ! tracer loop
138         !                                                       ! ===========
139         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
140         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
141         !                                                 
142         DO jk = 1, jpkm1                               
143            !                                             
144            !--- Computation of the ustream and downstream value of the tracer and the mask
145            DO jj = 2, jpjm1
146               DO ji = fs_2, fs_jpim1   ! vector opt.
147                  ! Upstream in the x-direction for the tracer
148                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
149                  ! Downstream in the x-direction for the tracer
150                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
151               END DO
152            END DO
153         END DO
154         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
155         
156         !
157         ! Horizontal advective fluxes
158         ! ---------------------------
159         !
160         DO jk = 1, jpkm1                             
161            DO jj = 2, jpjm1
162               DO ji = fs_2, fs_jpim1   ! vector opt.         
163                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
164                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
165               END DO
166            END DO
167         END DO
168         !
169         DO jk = 1, jpkm1 
170            zdt =  p2dt(jk)
171            DO jj = 2, jpjm1
172               DO ji = fs_2, fs_jpim1   ! vector opt.   
173                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
174                  zdx  = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
175                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
176                  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
177                  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
178               END DO
179            END DO
180         END DO 
181         !--- Lateral boundary conditions
182         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
183         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
184
185         !--- QUICKEST scheme
186         CALL quickest( zfu, zfd, zfc, zwx )
187         !
188         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
189         DO jk = 1, jpkm1 
190            DO jj = 2, jpjm1
191               DO ji = fs_2, fs_jpim1   ! vector opt.               
192                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
193               END DO
194            END DO
195         END DO
196         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
197
198         !
199         ! Tracer flux on the x-direction
200         DO jk = 1, jpkm1 
201            !
202            DO jj = 2, jpjm1
203               DO ji = fs_2, fs_jpim1   ! vector opt.               
204                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
205                  !--- If the second ustream point is a land point
206                  !--- the flux is computed by the 1st order UPWIND scheme
207                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
208                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
209                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
210               END DO
211            END DO
212            !
213            ! Computation of the trend
214            DO jj = 2, jpjm1
215               DO ji = fs_2, fs_jpim1   ! vector opt. 
216                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(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            !
224         END DO
225         !                                 ! trend diagnostics (contribution of upstream fluxes)
226         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
227         !
228      END DO
229      !
230      IF( wrk_not_released(3, 1,2,3) )   CALL ctl_stop('tra_adv_qck_i: failed to release workspace arrays')
231      !
232   END SUBROUTINE tra_adv_qck_i
233
234
235   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
236      &                                        ptb, ptn, pta, kjpt )
237      !!----------------------------------------------------------------------
238      !!
239      !!----------------------------------------------------------------------
240      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
241      USE oce     , ONLY:   zwy => ua       ! ua used as workspace
242      USE wrk_nemo, ONLY:   zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3   ! 3D workspace
243      !
244      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
245      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
246      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
247      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
248      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
249      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
250      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
251      !!
252      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
253      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
254      !----------------------------------------------------------------------
255      !
256      IF(wrk_in_use(3, 1,2,3))THEN
257         CALL ctl_stop('tra_adv_qck_j: ERROR: requested workspace arrays unavailable')
258         RETURN
259      END IF
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            !
339            ! Computation of the trend
340            DO jj = 2, jpjm1
341               DO ji = fs_2, fs_jpim1   ! vector opt. 
342                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
343                  ! horizontal advective trends
344                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
345                  !--- add it to the general tracer trends
346                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
347               END DO
348            END DO
349            !
350         END DO
351         !                                 ! trend diagnostics (contribution of upstream fluxes)
352         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
353         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
354         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
355           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
356           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
357         ENDIF
358         !
359      END DO
360      !
361      IF( wrk_not_released(3, 1,2,3) )   CALL ctl_stop('tra_adv_qck_j: failed to release workspace arrays')
362      !
363   END SUBROUTINE tra_adv_qck_j
364
365
366   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
367     &                                    ptn, pta, kjpt )
368      !!----------------------------------------------------------------------
369      !!
370      !!----------------------------------------------------------------------
371      USE oce, ONLY:   zwz => ua   ! ua used as workspace
372      !
373      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
374      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
375      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
376      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
377      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
378      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
379      !
380      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
381      REAL(wp) ::   zbtr , ztra      ! local scalars
382      !!----------------------------------------------------------------------
383
384      !                                                          ! ===========
385      DO jn = 1, kjpt                                            ! tracer loop
386         !                                                       ! ===========
387         ! 1. Bottom value : flux set to zero
388         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
389         !
390         !                                 ! Surface value
391         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
392         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
393         ENDIF
394         !
395         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
396            DO jj = 2, jpjm1
397               DO ji = fs_2, fs_jpim1   ! vector opt.
398                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
399               END DO
400            END DO
401         END DO
402         !
403         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
404            DO jj = 2, jpjm1
405               DO ji = fs_2, fs_jpim1   ! vector opt.
406                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
407                  ! k- vertical advective trends
408                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
409                  ! added to the general tracer trends
410                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
411               END DO
412            END DO
413         END DO
414         !                                 ! Save the vertical advective trends for diagnostic
415         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
416         !
417      END DO
418      !
419   END SUBROUTINE tra_adv_cen2_k
420
421
422   SUBROUTINE quickest( pfu, pfd, pfc, puc )
423      !!----------------------------------------------------------------------
424      !!
425      !! ** Purpose :  Computation of advective flux with Quickest scheme
426      !!
427      !! ** Method :   
428      !!----------------------------------------------------------------------
429      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
430      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
431      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
432      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
433      !!
434      INTEGER  ::  ji, jj, jk               ! dummy loop indices
435      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
436      REAL(wp) ::  zc, zcurv, zfho          !   -      -
437      !----------------------------------------------------------------------
438
439      DO jk = 1, jpkm1
440         DO jj = 1, jpj
441            DO ji = 1, jpi
442               zc     = puc(ji,jj,jk)                         ! Courant number
443               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
444               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
445               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
446               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
447               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
448               !
449               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
450               zcoef2 = ABS( zcoef1 )
451               zcoef3 = ABS( zcurv )
452               IF( zcoef3 >= zcoef2 ) THEN
453                  zfho = pfc(ji,jj,jk) 
454               ELSE
455                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
456                  IF( zcoef1 >= 0. ) THEN
457                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
458                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
459                  ELSE
460                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
461                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
462                  ENDIF
463               ENDIF
464               puc(ji,jj,jk) = zfho
465            END DO
466         END DO
467      END DO
468      !
469   END SUBROUTINE quickest
470
471   !!======================================================================
472END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.