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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch

  • Property svn:keywords set to Id
File size: 23.6 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 diaptr          ! poleward transport diagnostics
23   !
24   USE lib_mpp         ! distribued memory computing
25   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
26   USE in_out_manager  ! I/O manager
27   USE timing          ! Timing
28   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_adv_qck   ! routine called by step.F90
34
35   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
36
37   LOGICAL  ::   l_trd   ! flag to compute trends
38   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
39
40
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
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      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
83      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
84      !!
85      !! ** Reference : Leonard (1979, 1991)
86      !!----------------------------------------------------------------------
87      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
88      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
89      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
90      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
91      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
95      !!----------------------------------------------------------------------
96      !
97      IF( ln_timing )   CALL timing_start('tra_adv_qck')
98      !
99      IF( kt == kit000 )  THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
103         IF(lwp) WRITE(numout,*)
104      ENDIF
105      !
106      l_trd = .FALSE.
107      l_ptr = .FALSE.
108      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
109      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
110      !
111      !
112      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
113      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
114      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
115
116      !        ! vertical fluxes are computed with the 2nd order centered scheme
117      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
118      !
119      IF( ln_timing )   CALL timing_stop('tra_adv_qck')
120      !
121   END SUBROUTINE tra_adv_qck
122
123
124   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
125      &                                        ptb, ptn, pta, kjpt   )
126      !!----------------------------------------------------------------------
127      !!
128      !!----------------------------------------------------------------------
129      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
130      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
131      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
132      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
133      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
134      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
135      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
136      !!
137      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
138      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
139      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zfu, zfc, zfd
140      !----------------------------------------------------------------------
141      !
142      !                                                          ! ===========
143      DO jn = 1, kjpt                                            ! tracer loop
144         !                                                       ! ===========
145         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
146         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
147         !
148!!gm why not using a SHIFT instruction...
149         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
150            DO jj = 2, jpjm1
151               DO ji = fs_2, fs_jpim1   ! vector opt.
152                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
153                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
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         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            DO jj = 2, jpjm1
173               DO ji = fs_2, fs_jpim1   ! vector opt.   
174                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
175                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
176                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
177                  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
178                  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
179               END DO
180            END DO
181         END DO 
182         !--- Lateral boundary conditions
183         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
184         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
185
186         !--- QUICKEST scheme
187         CALL quickest( zfu, zfd, zfc, zwx )
188         !
189         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
190         DO jk = 1, jpkm1 
191            DO jj = 2, jpjm1
192               DO ji = fs_2, fs_jpim1   ! vector opt.               
193                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
194               END DO
195            END DO
196         END DO
197         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
198
199         !
200         ! Tracer flux on the x-direction
201         DO jk = 1, jpkm1 
202            !
203            DO jj = 2, jpjm1
204               DO ji = fs_2, fs_jpim1   ! vector opt.               
205                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
206                  !--- If the second ustream point is a land point
207                  !--- the flux is computed by the 1st order UPWIND scheme
208                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
209                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
210                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
211               END DO
212            END DO
213         END DO
214         !
215         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
216         !
217         ! Computation of the trend
218         DO jk = 1, jpkm1 
219            DO jj = 2, jpjm1
220               DO ji = fs_2, fs_jpim1   ! vector opt. 
221                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
222                  ! horizontal advective trends
223                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
224                  !--- add it to the general tracer trends
225                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
226               END DO
227            END DO
228         END DO
229         !                                 ! trend diagnostics
230         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
231         !
232      END DO
233      !
234   END SUBROUTINE tra_adv_qck_i
235
236
237   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
238      &                                        ptb, ptn, pta, kjpt )
239      !!----------------------------------------------------------------------
240      !!
241      !!----------------------------------------------------------------------
242      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
243      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
244      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
245      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
246      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
247      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
248      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
249      !!
250      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
251      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
252      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
253      !----------------------------------------------------------------------
254      !
255      !                                                          ! ===========
256      DO jn = 1, kjpt                                            ! tracer loop
257         !                                                       ! ===========
258         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
259         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
260         !                                                 
261         DO jk = 1, jpkm1                               
262            !                                             
263            !--- Computation of the ustream and downstream value of the tracer and the mask
264            DO jj = 2, jpjm1
265               DO ji = fs_2, fs_jpim1   ! vector opt.
266                  ! Upstream in the x-direction for the tracer
267                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
268                  ! Downstream in the x-direction for the tracer
269                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
270               END DO
271            END DO
272         END DO
273         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
274
275         
276         !
277         ! Horizontal advective fluxes
278         ! ---------------------------
279         !
280         DO jk = 1, jpkm1                             
281            DO jj = 2, jpjm1
282               DO ji = fs_2, fs_jpim1   ! vector opt.         
283                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
284                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
285               END DO
286            END DO
287         END DO
288         !
289         DO jk = 1, jpkm1 
290            DO jj = 2, jpjm1
291               DO ji = fs_2, fs_jpim1   ! vector opt.   
292                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
293                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
294                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
295                  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
296                  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
297               END DO
298            END DO
299         END DO
300
301         !--- Lateral boundary conditions
302         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
303         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
304
305         !--- QUICKEST scheme
306         CALL quickest( zfu, zfd, zfc, zwy )
307         !
308         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
309         DO jk = 1, jpkm1 
310            DO jj = 2, jpjm1
311               DO ji = fs_2, fs_jpim1   ! vector opt.               
312                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
313               END DO
314            END DO
315         END DO
316         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions
317         !
318         ! Tracer flux on the x-direction
319         DO jk = 1, jpkm1 
320            !
321            DO jj = 2, jpjm1
322               DO ji = fs_2, fs_jpim1   ! vector opt.               
323                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
324                  !--- If the second ustream point is a land point
325                  !--- the flux is computed by the 1st order UPWIND scheme
326                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
327                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
328                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
329               END DO
330            END DO
331         END DO
332         !
333         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
334         !
335         ! Computation of the trend
336         DO jk = 1, jpkm1 
337            DO jj = 2, jpjm1
338               DO ji = fs_2, fs_jpim1   ! vector opt. 
339                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
340                  ! horizontal advective trends
341                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
342                  !--- add it to the general tracer trends
343                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
344               END DO
345            END DO
346         END DO
347         !                                 ! trend diagnostics
348         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
349         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
350         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
351         !
352      END DO
353      !
354   END SUBROUTINE tra_adv_qck_j
355
356
357   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
358     &                                    ptn, pta, kjpt )
359      !!----------------------------------------------------------------------
360      !!
361      !!----------------------------------------------------------------------
362      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
363      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
364      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
365      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
366      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
367      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
368      !
369      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
370      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
371      !!----------------------------------------------------------------------
372      !
373      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
374      zwz(:,:,jpk) = 0._wp
375      !
376      !                                                          ! ===========
377      DO jn = 1, kjpt                                            ! tracer loop
378         !                                                       ! ===========
379         !
380         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
381            DO jj = 2, jpjm1
382               DO ji = fs_2, fs_jpim1   ! vector opt.
383                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
384               END DO
385            END DO
386         END DO
387         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
388            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
389               DO jj = 1, jpj
390                  DO ji = 1, jpi
391                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
392                  END DO
393               END DO   
394            ELSE                                   ! no ocean cavities (only ocean surface)
395               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
396            ENDIF
397         ENDIF
398         !
399         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
400            DO jj = 2, jpjm1
401               DO ji = fs_2, fs_jpim1   ! vector opt.
402                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
403                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
404               END DO
405            END DO
406         END DO
407         !                                 ! Send trends for diagnostic
408         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
409         !
410      END DO
411      !
412   END SUBROUTINE tra_adv_cen2_k
413
414
415   SUBROUTINE quickest( pfu, pfd, pfc, puc )
416      !!----------------------------------------------------------------------
417      !!
418      !! ** Purpose :  Computation of advective flux with Quickest scheme
419      !!
420      !! ** Method :   
421      !!----------------------------------------------------------------------
422      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
423      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
424      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
425      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
426      !!
427      INTEGER  ::  ji, jj, jk               ! dummy loop indices
428      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
429      REAL(wp) ::  zc, zcurv, zfho          !   -      -
430      !----------------------------------------------------------------------
431      !
432      IF( ln_timing )   CALL timing_start('quickest')
433      !
434      DO jk = 1, jpkm1
435         DO jj = 1, jpj
436            DO ji = 1, jpi
437               zc     = puc(ji,jj,jk)                         ! Courant number
438               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
439               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
440               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
441               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
442               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
443               !
444               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
445               zcoef2 = ABS( zcoef1 )
446               zcoef3 = ABS( zcurv )
447               IF( zcoef3 >= zcoef2 ) THEN
448                  zfho = pfc(ji,jj,jk) 
449               ELSE
450                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
451                  IF( zcoef1 >= 0. ) THEN
452                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
453                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
454                  ELSE
455                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
456                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
457                  ENDIF
458               ENDIF
459               puc(ji,jj,jk) = zfho
460            END DO
461         END DO
462      END DO
463      !
464      IF( ln_timing )   CALL timing_stop('quickest')
465      !
466   END SUBROUTINE quickest
467
468   !!======================================================================
469END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.