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 tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 23.1 KB
RevLine 
[1878]1MODULE traadv_qck
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_qck  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  3.0  !  2008-07  (G. Reffray)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_adv_qck      : update the tracer trend with the horizontal advection
11   !!                      trends using a 3rd order finite difference scheme
12   !!   tra_adv_qck_i  :
13   !!   tra_adv_qck_j  :
14   !!   tra_adv_cen2_k : 2nd centered scheme for the vertical advection
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and active tracers
17   USE dom_oce         ! ocean space and time domain
18   USE trdmod          ! ocean active tracers trends
19   USE trdmod_oce      ! ocean variables trends
20   USE trabbl          ! advective term in the BBL
21   USE lib_mpp         ! distribued memory computing
22   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
23   USE dynspg_oce      ! surface pressure gradient variables
24   USE in_out_manager  ! I/O manager
25   USE diaptr          ! poleward transport diagnostics
26   USE prtctl          ! Print control
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   tra_adv_qck   ! routine called by step.F90
32
33   REAL(wp), DIMENSION(jpi,jpj) ::   btr2
34   REAL(wp)                     ::   r1_6
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
41   !! $Id$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_adv_qck  ***
50      !!
51      !! ** Purpose :   Compute the now trend due to the advection of tracers
52      !!      and add it to the general trend of passive tracer equations.
53      !!
54      !! ** Method :   The advection is evaluated by a third order scheme
55      !!             For a positive velocity u :              u(i)>0
56      !!                                          |--FU--|--FC--|--FD--|------|
57      !!                                             i-1    i      i+1   i+2
58      !!
59      !!             For a negative velocity u :              u(i)<0
60      !!                                          |------|--FD--|--FC--|--FU--|
61      !!                                             i-1    i      i+1   i+2
62      !!             where  FU is the second upwind point
63      !!                    FD is the first douwning point
64      !!                    FC is the central point (or the first upwind point)
65      !!
66      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
67      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
68      !!
69      !!         dt = 2*rdtra and the scalar values are tb and sb
70      !!
71      !!       On the vertical, the simple centered scheme used tn and sn
72      !!
73      !!               The fluxes are bounded by the ULTIMATE limiter to
74      !!             guarantee the monotonicity of the solution and to
75      !!            prevent the appearance of spurious numerical oscillations
76      !!
77      !! ** Action : - update (ta,sa) with the now advective tracer trends
78      !!             - save the trends ('key_trdtra')
79      !!
80      !! ** Reference : Leonard (1979, 1991)
81      !!----------------------------------------------------------------------
82      USE oce, ONLY :   ztrdt => ua   ! use ua as workspace
83      USE oce, ONLY :   ztrds => va   ! use va as workspace
84      !!
85      INTEGER , INTENT(in)                         ::  kt  ! ocean time-step index
86      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun ! effective ocean velocity, u_component
87      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pvn ! effective ocean velocity, v_component
88      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pwn ! effective ocean velocity, w_component
89      !!
90      INTEGER  ::   ji, jj, jk                          ! dummy loop indices
91      REAL(wp) ::   z_hdivn_x, z_hdivn_y, z_hdivn       ! temporary scalars
92      REAL(wp) ::   zbtr, z2                            !    "         "
93      !!----------------------------------------------------------------------
94
95      IF( kt == nit000 ) THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme'
98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
99         IF(lwp) WRITE(numout,*)
100         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
101         r1_6      = 1. / 6.
102      ENDIF
103
104      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
105      ELSE                                        ;    z2 = 2.
106      ENDIF
107
108      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
109      !---------------------------------------------------------------------------
110
111      CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2)
112      CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2)
113
114      IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
115
116      CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2)
117      CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2) 
118
119      IF( l_trdtra ) THEN
120         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
121         !
122         ztrdt(:,:,:) = ta(:,:,:)    ! Save the horizontal up-to-date ta/sa trends
123         ztrds(:,:,:) = sa(:,:,:)
124      END IF
125
126      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' qck had  - Ta: ', mask1=tmask, &
127         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
128
129      ! II. The vertical fluxes are computed with the 2nd order centered scheme
130      !-------------------------------------------------------------------------
131      !
132      CALL tra_adv_cen2_k( pwn, tn, ta )
133      CALL tra_adv_cen2_k( pwn, sn, sa )
134      !
135      !Save the vertical advective trends for diagnostic
136      ! ----------------------------------------------------
137      IF( l_trdtra )   THEN
138         ! Recompute the vertical advection zta & zsa trends computed
139         ! at the step 2. above in making the difference between the new
140         ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
141         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
142
143         DO jk = 1, jpkm1
144            DO jj = 2, jpjm1
145               DO ji = fs_2, fs_jpim1   ! vector opt.
146#if defined key_zco
147                  zbtr      = btr2(ji,jj)
148                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
149                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
150#else
151                  zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk)
152                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
153                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
154#endif
155                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
156                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
157                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
158               END DO
159            END DO
160         END DO
161         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
162      ENDIF
163
164      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' qck zad  - Ta: ', mask1=tmask, &
165         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
166      !
167   END SUBROUTINE tra_adv_qck
168
169
170   SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 )
171      !!----------------------------------------------------------------------
172      !!
173      !!----------------------------------------------------------------------
174      REAL,     INTENT(in)                            :: z2
175      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran  ! horizontal effective velocity
176      REAL(wp), INTENT(out)  , DIMENSION(jpi,jpj,jpk) :: ztrdtra
177      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
178      !
179      INTEGER  :: ji, jj, jk
180      REAL(wp) :: za, zbtr, dir, dx, dt     ! temporary scalars
181      REAL(wp) :: z_hdivn_x
182      REAL(wp), DIMENSION(jpi,jpj)     ::  zmask, zupst, zdwst, zc_cfl
183      REAL(wp), DIMENSION(jpi,jpj)     ::  zfu, zfc, zfd, zfho, zmskl, zsc_e 
184      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zflux
185      !----------------------------------------------------------------------
186
187      zfu   (:,jpj) = 0.e0   ;   zfc   (:,jpj) = 0.e0
188      zfd   (:,jpj) = 0.e0   ;   zc_cfl(:,jpj) = 0.e0
189      zsc_e (:,jpj) = 0.e0   ;   zmskl (:,jpj) = 0.e0
190      zfho  (:,jpj) = 0.e0
191                                                       ! ===============
192      DO jk = 1, jpkm1                                 ! Horizontal slab
193         !                                             ! ===============
194         !--- Computation of the ustream and downstream value of the tracer and the mask
195         DO jj = 2, jpjm1
196            DO ji = 2, fs_jpim1   ! vector opt.
197               ! Upstream in the x-direction for the tracer
198               zupst(ji,jj)=tra(ji-1,jj,jk)
199               ! Downstream in the x-direction for the tracer
200               zdwst(ji,jj)=tra(ji+1,jj,jk)
201               ! Mask at the T-points in the x-direction (mask=0 or mask=1)
202               zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
203            END DO
204         END DO
205         !
206         !--- Lateral boundary conditions
207         CALL lbc_lnk( zupst(:,:), 'T', 1. )
208         CALL lbc_lnk( zdwst(:,:), 'T', 1. ) 
209         CALL lbc_lnk( zmask(:,:), 'T', 1. ) 
210         !
211         ! Horizontal advective fluxes
212         ! ---------------------------
213         !
214         dt = z2 * rdttra(jk)
215         !--- tracer flux at u-points
216         DO jj = 1, jpjm1
217            DO ji = 1, jpi
218#if defined key_zco
219               zsc_e(ji,jj) = e2u(ji,jj)
220#else
221               zsc_e(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk)
222#endif
223               dir = 0.5 + sign(0.5,pun(ji,jj,jk))                             ! if pun>0 : dir = 1 otherwise dir = 0
224               dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
225               zc_cfl (ji,jj) = ABS(pun(ji,jj,jk))*dt/dx                       ! (0<zc_cfl<1 : Courant number on x-direction)
226
227               zfu(ji,jj)   = dir*zupst(ji  ,jj   )+(1-dir)*zdwst(ji+1,jj   )  ! FU in the x-direction for T
228               zfc(ji,jj)   = dir*tra  (ji  ,jj,jk)+(1-dir)*tra  (ji+1,jj,jk)  ! FC in the x-direction for T
229               zfd(ji,jj)   = dir*tra  (ji+1,jj,jk)+(1-dir)*tra  (ji  ,jj,jk)  ! FD in the x-direction for T
230               zmskl(ji,jj) = dir*zmask(ji  ,jj)   +(1-dir)*zmask(ji+1,jj)
231           END DO
232         END DO
233         !
234         !--- QUICKEST scheme
235         ! Tracer flux on the x-direction
236         CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
237         !--- If the second ustream point is a land point
238         !--- the flux is computed by the 1st order UPWIND scheme
239         zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:)
240         !--- Computation of fluxes
241         zflux(:,:,jk) = zsc_e(:,:)*pun(:,:,jk)*zfho(:,:)
242         !
243         !--- Tracer flux divergence at t-point added to the general trend
244         DO jj = 2, jpjm1
245            DO ji = fs_2, fs_jpim1   ! vector opt.
246               !--- horizontal advective trends
247#if defined key_zco
248               zbtr = btr2(ji,jj)
249#else
250               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)             
251#endif
252               za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) )
253               !--- add it to the general tracer trends
254               traa(ji,jj,jk) = traa(ji,jj,jk) + za
255            END DO
256         END DO
257         !                                             ! ===============
258      END DO                                           !   End of slab
259      !                                                ! ===============
260      !
261      !  Save the horizontal advective trends for diagnostic
262      ! -----------------------------------------------------
263      IF( l_trdtra ) THEN
264         ! T/S ZONAL advection trends
265         ztrdtra(:,:,:) = 0.e0
266         !
267         DO jk = 1, jpkm1
268            DO jj = 2, jpjm1
269               DO ji = fs_2, fs_jpim1   ! vector opt.
270                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
271                  !   N.B. This computation is not valid along OBCs (if any)
272#if defined key_zco
273                  zbtr      = btr2(ji,jj)
274                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
275                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
276#else
277                  zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk)
278                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
279                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
280#endif
281                  ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) + tran(ji,jj,jk) * z_hdivn_x
282               END DO
283            END DO
284         END DO
285      END IF
286
287   END SUBROUTINE tra_adv_qck_i
288
289
290   SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 )
291      !!----------------------------------------------------------------------
292      !!
293      !!----------------------------------------------------------------------
294      INTEGER,  INTENT(in)                            :: kt              ! ocean time-step index
295      REAL,     INTENT(in)                            :: z2
296      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pvn, tra, tran  ! horizontal effective velocity
297      REAL(wp), INTENT(out)  , DIMENSION(jpj)         :: trd_adv
298      REAL(wp), INTENT(out)  , DIMENSION(jpi,jpj,jpk) :: ztrdtra
299      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
300      !!
301      INTEGER  :: ji, jj, jk
302      REAL(wp) :: za, zbtr, dir, dx, dt     ! temporary scalars
303      REAL(wp) :: z_hdivn_y
304      REAL(wp), DIMENSION(jpi,jpj)     ::  zmask, zupst, zdwst, zc_cfl
305      REAL(wp), DIMENSION(jpi,jpj)     ::  zfu, zfc, zfd, zfho, zmskl, zsc_e
306      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zflux
307      !----------------------------------------------------------------------
308      ! II. Part 2 : y-direction
309      !----------------------------------------------------------------------
310
311      zfu   (:,jpj) = 0.e0   ;   zfc   (:,jpj) = 0.e0
312      zfd   (:,jpj) = 0.e0   ;   zc_cfl(:,jpj) = 0.e0
313      zsc_e (:,jpj) = 0.e0   ;   zmskl (:,jpj) = 0.e0
314      zfho  (:,jpj) = 0.e0
315
316                                                       ! ===============
317      DO jk = 1, jpkm1                                 ! Horizontal slab
318         !                                             ! ===============
319         !--- Computation of the ustream and downstream value of the tracer and the mask
320         DO jj = 2, jpjm1
321            DO ji = 2, fs_jpim1   ! vector opt.
322               ! Upstream in the x-direction for the tracer
323               zupst(ji,jj)=tra(ji,jj-1,jk)
324               ! Downstream in the x-direction for the tracer
325               zdwst(ji,jj)=tra(ji,jj+1,jk)
326               ! Mask at the T-points in the x-direction (mask=0 or mask=1)
327               zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
328            END DO
329         END DO
330         !
331         !--- Lateral boundary conditions
332         CALL lbc_lnk( zupst(:,:), 'T', 1. )
333         CALL lbc_lnk( zdwst(:,:), 'T', 1. )
334         CALL lbc_lnk( zmask(:,:), 'T', 1. )
335         !
336         ! Horizontal advective fluxes
337         ! ---------------------------
338         !
339         dt = z2 * rdttra(jk)
340         !--- tracer flux at v-points
341         DO jj = 1, jpjm1
342            DO ji = 1, jpi
343#if defined key_zco
344               zsc_e(ji,jj) = e1v(ji,jj)
345#else
346               zsc_e(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk)
347#endif
348               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))                             ! if pvn>0 : dir = 1 otherwise dir = 0
349               dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
350               zc_cfl(ji,jj) = ABS(pvn(ji,jj,jk))*dt/dx                        ! (0<zc_cfl<1 : Courant number on y-direction)
351
352               zfu(ji,jj)   = dir*zupst(ji,jj     )+(1-dir)*zdwst(ji,jj+1   )  ! FU in the y-direction for T
353               zfc(ji,jj)   = dir*tra  (ji,jj  ,jk)+(1-dir)*tra  (ji,jj+1,jk)  ! FC in the y-direction for T
354               zfd(ji,jj)   = dir*tra  (ji,jj+1,jk)+(1-dir)*tra  (ji,jj  ,jk)  ! FD in the y-direction for T
355               zmskl(ji,jj) = dir*zmask(ji,jj     )+(1-dir)*zmask(ji,jj+1)
356            END DO
357         END DO
358         !
359         !--- QUICKEST scheme
360         ! Tracer flux on the y-direction
361         CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
362         !--- If the second ustream point is a land point
363         !--- the flux is computed by the 1st order UPWIND scheme
364         zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:)
365         !--- Computation of fluxes
366         zflux(:,:,jk) = zsc_e(:,:)*pvn(:,:,jk)*zfho(:,:)
367         !
368         !--- Tracer flux divergence at t-point added to the general trend
369         DO jj = 2, jpjm1
370            DO ji = fs_2, fs_jpim1   ! vector opt.
371               !--- horizontal advective trends
372#if defined key_zco
373               zbtr = btr2(ji,jj)
374#else
375               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
376#endif
377               za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) )
378               !--- add it to the general tracer trends
379               traa(ji,jj,jk) = traa(ji,jj,jk) + za
380            END DO
381         END DO
382         !                                             ! ===============
383      END DO                                           !   End of slab
384      !                                                ! ===============
385      !
386      !  Save the horizontal advective trends for diagnostic
387      ! -----------------------------------------------------
388      IF( l_trdtra ) THEN
389         ! T/S MERIDIONAL advection trends
390         DO jk = 1, jpkm1
391            DO jj = 2, jpjm1
392               DO ji = fs_2, fs_jpim1   ! vector opt.
393                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
394                  !   N.B. This computation is not valid along OBCs (if any)
395#if defined key_zco
396                  zbtr      = btr2(ji,jj)
397                  z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              &
398                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
399#else
400                  zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk)
401                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
402                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
403#endif
404                  ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) + tran(ji,jj,jk) * z_hdivn_y
405               END DO
406            END DO
407         END DO
408      END IF
409
410      ! "zonal" mean advective heat and salt transport
411      ! ----------------------------------------------
412
413      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
414         IF( lk_zco ) THEN
415            DO jk = 1, jpkm1
416               DO jj = 2, jpjm1
417                  DO ji = fs_2, fs_jpim1   ! vector opt.
418                    zflux(ji,jj,jk) = zflux(ji,jj,jk) * fse3v(ji,jj,jk)
419                  END DO
420               END DO
421            END DO
422         ENDIF
423         trd_adv(:) = ptr_vj( zflux(:,:,:) )
424      ENDIF
425
426   END SUBROUTINE tra_adv_qck_j
427
428
429   SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta )
430      !!----------------------------------------------------------------------
431      !!
432      !!----------------------------------------------------------------------
433      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)  :: pwn   ! vertical effective velocity
434      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)  :: ptn   ! now tracer
435      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)  :: pta   ! tracer general trend
436      !!
437      INTEGER  ::   ji, jj, jk   ! dummy loop indices
438      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zflux   ! 3D workspace
439      !!----------------------------------------------------------------------
440      !
441      !                         !==  Vertical advective fluxes  ==!
442      zflux(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
443      !
444      !                                 ! Surface value
445      IF( lk_vvl ) THEN   ;   zflux(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
446      ELSE                ;   zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! Constant volume : advective flux through the surface
447      ENDIF
448      !
449      DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
450         DO jj = 2, jpjm1
451            DO ji = fs_2, fs_jpim1   ! vector opt.
452               zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) )
453            END DO
454         END DO
455      END DO
456      !
457      DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
458         DO jj = 2, jpjm1
459            DO ji = fs_2, fs_jpim1   ! vector opt.
460               pta(ji,jj,jk) =  pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) )   &
461                  &                           /   fse3t(ji,jj,jk)
462            END DO
463         END DO
464      END DO
465      !
466   END SUBROUTINE tra_adv_cen2_k
467
468
469   SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl )
470      !!----------------------------------------------------------------------
471      !!
472      !!----------------------------------------------------------------------
473      REAL(wp), INTENT(in)  , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl 
474      REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho
475      REAL(wp)              , DIMENSION(jpi,jpj) :: zcurv, zcoef1, zcoef2, zcoef3  ! temporary scalars
476      !
477      zcurv (:,:) = fd(:,:) + fu(:,:) - 2.*fc(:,:)
478      zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) )
479      zcoef2(:,:) = 0.5*fc_cfl(:,:)*( fd(:,:) - fc(:,:) )
480      zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*r1_6 )*zcurv(:,:)
481      fho   (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:)                         ! phi_f QUICKEST
482      !
483      zcoef1(:,:) = fd(:,:) - fu(:,:)                                               ! DEL
484      zcoef2(:,:) = ABS( zcoef1(:,:) )                                              ! ABS(DEL)
485      zcoef3(:,:) = ABS( zcurv(:,:) )                                               ! ABS(CURV)
486      !
487      WHERE ( zcoef3(:,:) >= zcoef2(:,:) ) 
488        fho(:,:) = fc(:,:)
489      ELSEWHERE
490        zcoef3(:,:) = fu(:,:) + ( ( fc(:,:) - fu(:,:) )/MAX(fc_cfl(:,:),1.e-9) )    ! phi_REF
491          WHERE ( zcoef1(:,:) >= 0.e0 )
492            fho(:,:) = MAX(fc(:,:),fho(:,:))
493            fho(:,:) = MIN(fho(:,:),MIN(zcoef3(:,:),fd(:,:)))
494          ELSEWHERE
495            fho(:,:) = MIN(fc(:,:),fho(:,:))
496            fho(:,:) = MAX(fho(:,:),MAX(zcoef3(:,:),fd(:,:)))
497          ENDWHERE
498      ENDWHERE
499      !
500   END SUBROUTINE quickest
501
502   !!======================================================================
503END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.