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

source: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 7236

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

All changes related to diaptr (basin heat transports and transport components)
  • Property svn:keywords set to Id
File size: 23.8 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 wrk_nemo        ! Memory Allocation
28   USE timing          ! Timing
29   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_adv_qck   ! routine called by step.F90
35
36   LOGICAL  :: l_trd           ! flag to compute trends
37   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
49      &                                       ptb, ptn, pta, kjpt )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_adv_qck  ***
52      !!
53      !! ** Purpose :   Compute the now trend due to the advection of tracers
54      !!      and add it to the general trend of passive tracer equations.
55      !!
56      !! ** Method :   The advection is evaluated by a third order scheme
57      !!             For a positive velocity u :              u(i)>0
58      !!                                          |--FU--|--FC--|--FD--|------|
59      !!                                             i-1    i      i+1   i+2
60      !!
61      !!             For a negative velocity u :              u(i)<0
62      !!                                          |------|--FD--|--FC--|--FU--|
63      !!                                             i-1    i      i+1   i+2
64      !!             where  FU is the second upwind point
65      !!                    FD is the first douwning point
66      !!                    FC is the central point (or the first upwind point)
67      !!
68      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
69      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
70      !!
71      !!         dt = 2*rdtra and the scalar values are tb and sb
72      !!
73      !!       On the vertical, the simple centered scheme used ptn
74      !!
75      !!               The fluxes are bounded by the ULTIMATE limiter to
76      !!             guarantee the monotonicity of the solution and to
77      !!            prevent the appearance of spurious numerical oscillations
78      !!
79      !! ** Action : - update pta  with the now advective tracer trends
80      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
81      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
82      !!
83      !! ** Reference : Leonard (1979, 1991)
84      !!----------------------------------------------------------------------
85      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
86      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
87      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
88      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
89      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
90      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
91      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
92      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
93      !!----------------------------------------------------------------------
94      !
95      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
96      !
97      IF( kt == kit000 )  THEN
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
100         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
101         IF(lwp) WRITE(numout,*)
102      ENDIF
103      !
104      l_trd = .FALSE.
105      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
106      !
107      !        ! 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      !        ! 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)                             , INTENT(in   ) ::   p2dt       ! 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, 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._wp     ;   zfc(:,:,:) = 0._wp 
142         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
143         !
144!!gm why not using a SHIFT instruction...
145         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
149                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
150               END DO
151            END DO
152         END DO
153         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
154         
155         !
156         ! Horizontal advective fluxes
157         ! ---------------------------
158         DO jk = 1, jpkm1                             
159            DO jj = 2, jpjm1
160               DO ji = fs_2, fs_jpim1   ! vector opt.         
161                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
162                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
163               END DO
164            END DO
165         END DO
166         !
167         DO jk = 1, jpkm1 
168            DO jj = 2, jpjm1
169               DO ji = fs_2, fs_jpim1   ! vector opt.   
170                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
171                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
172                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
173                  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
174                  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
175               END DO
176            END DO
177         END DO 
178         !--- Lateral boundary conditions
179         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
180         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
181
182         !--- QUICKEST scheme
183         CALL quickest( zfu, zfd, zfc, zwx )
184         !
185         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
186         DO jk = 1, jpkm1 
187            DO jj = 2, jpjm1
188               DO ji = fs_2, fs_jpim1   ! vector opt.               
189                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
190               END DO
191            END DO
192         END DO
193         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
194
195         !
196         ! Tracer flux on the x-direction
197         DO jk = 1, jpkm1 
198            !
199            DO jj = 2, jpjm1
200               DO ji = fs_2, fs_jpim1   ! vector opt.               
201                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
202                  !--- If the second ustream point is a land point
203                  !--- the flux is computed by the 1st order UPWIND scheme
204                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
205                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
206                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
207               END DO
208            END DO
209         END DO
210         !
211         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
212         !
213         ! Computation of the trend
214         DO jk = 1, jpkm1 
215            DO jj = 2, jpjm1
216               DO ji = fs_2, fs_jpim1   ! vector opt. 
217                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
218                  ! horizontal advective trends
219                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
220                  !--- add it to the general tracer trends
221                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
222               END DO
223            END DO
224         END DO
225         !                                 ! trend diagnostics
226         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
227         !
228      END DO
229      !
230      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
231      !
232   END SUBROUTINE tra_adv_qck_i
233
234
235   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
236      &                                        ptb, ptn, pta, kjpt )
237      !!----------------------------------------------------------------------
238      !!
239      !!----------------------------------------------------------------------
240      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
241      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
242      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
243      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
244      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
245      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
246      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
247      !!
248      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
249      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
250      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
251      !----------------------------------------------------------------------
252      !
253      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
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         !--- Lateral boundary conditions
317         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
318         !
319         ! Tracer flux on the x-direction
320         DO jk = 1, jpkm1 
321            !
322            DO jj = 2, jpjm1
323               DO ji = fs_2, fs_jpim1   ! vector opt.               
324                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
325                  !--- If the second ustream point is a land point
326                  !--- the flux is computed by the 1st order UPWIND scheme
327                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
328                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
329                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
330               END DO
331            END DO
332         END DO
333         !
334         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
335         !
336         ! Computation of the trend
337         DO jk = 1, jpkm1 
338            DO jj = 2, jpjm1
339               DO ji = fs_2, fs_jpim1   ! vector opt. 
340                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
341                  ! horizontal advective trends
342                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
343                  !--- add it to the general tracer trends
344                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
345               END DO
346            END DO
347         END DO
348         !                                 ! trend diagnostics
349         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
350         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
351         IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:) )
352         !
353      END DO
354      !
355      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
356      !
357   END SUBROUTINE tra_adv_qck_j
358
359
360   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
361     &                                    ptn, pta, kjpt )
362      !!----------------------------------------------------------------------
363      !!
364      !!----------------------------------------------------------------------
365      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
366      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
367      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
368      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
369      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
370      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
371      !
372      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
373      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
374      !!----------------------------------------------------------------------
375      !
376      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
377      !
378      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
379      zwz(:,:,jpk) = 0._wp
380      !
381      !                                                          ! ===========
382      DO jn = 1, kjpt                                            ! tracer loop
383         !                                                       ! ===========
384         !
385         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
386            DO jj = 2, jpjm1
387               DO ji = fs_2, fs_jpim1   ! vector opt.
388                  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)
389               END DO
390            END DO
391         END DO
392         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
393            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
394               DO jj = 1, jpj
395                  DO ji = 1, jpi
396                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
397                  END DO
398               END DO   
399            ELSE                                   ! no ocean cavities (only ocean surface)
400               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
401            ENDIF
402         ENDIF
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                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
408                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
409               END DO
410            END DO
411         END DO
412         !                                 ! Send trends for diagnostic
413         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
414         !
415      END DO
416      !
417      CALL wrk_dealloc( jpi,jpj,jpk,   zwz )
418      !
419   END SUBROUTINE tra_adv_cen2_k
420
421
422   SUBROUTINE quickest( pfu, pfd, pfc, puc )
423      !!----------------------------------------------------------------------
424      !!
425      !! ** Purpose :  Computation of advective flux with Quickest scheme
426      !!
427      !! ** Method :   
428      !!----------------------------------------------------------------------
429      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
430      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
431      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
432      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
433      !!
434      INTEGER  ::  ji, jj, jk               ! dummy loop indices
435      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
436      REAL(wp) ::  zc, zcurv, zfho          !   -      -
437      !----------------------------------------------------------------------
438      !
439      IF( nn_timing == 1 )  CALL timing_start('quickest')
440      !
441      DO jk = 1, jpkm1
442         DO jj = 1, jpj
443            DO ji = 1, jpi
444               zc     = puc(ji,jj,jk)                         ! Courant number
445               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
446               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
447               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
448               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
449               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
450               !
451               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
452               zcoef2 = ABS( zcoef1 )
453               zcoef3 = ABS( zcurv )
454               IF( zcoef3 >= zcoef2 ) THEN
455                  zfho = pfc(ji,jj,jk) 
456               ELSE
457                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
458                  IF( zcoef1 >= 0. ) THEN
459                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
460                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
461                  ELSE
462                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
463                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
464                  ENDIF
465               ENDIF
466               puc(ji,jj,jk) = zfho
467            END DO
468         END DO
469      END DO
470      !
471      IF( nn_timing == 1 )  CALL timing_stop('quickest')
472      !
473   END SUBROUTINE quickest
474
475   !!======================================================================
476END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.