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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 11738

Last change on this file since 11738 was 7179, checked in by timgraham, 7 years ago

Manually merge in changes from v3.6_extra_CMIP6_diagnostics branch.
This change also includes a change of the domain_def.xml file so XIOS2 must be used from this revision onwards

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