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/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

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