source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 6140

Last change on this file since 6140 was 6140, checked in by timgraham, 5 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

  • Property svn:keywords set to Id
File size: 23.9 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 ) THEN 
352           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
353           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
354         ENDIF
355         !
356      END DO
357      !
358      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
359      !
360   END SUBROUTINE tra_adv_qck_j
361
362
363   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
364     &                                    ptn, pta, kjpt )
365      !!----------------------------------------------------------------------
366      !!
367      !!----------------------------------------------------------------------
368      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
369      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
370      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
371      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
372      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
373      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
374      !
375      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
376      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
377      !!----------------------------------------------------------------------
378      !
379      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
380      !
381      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
382      zwz(:,:,jpk) = 0._wp
383      !
384      !                                                          ! ===========
385      DO jn = 1, kjpt                                            ! tracer loop
386         !                                                       ! ===========
387         !
388         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
389            DO jj = 2, jpjm1
390               DO ji = fs_2, fs_jpim1   ! vector opt.
391                  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)
392               END DO
393            END DO
394         END DO
395         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
396            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
397               DO jj = 1, jpj
398                  DO ji = 1, jpi
399                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
400                  END DO
401               END DO   
402            ELSE                                   ! no ocean cavities (only ocean surface)
403               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
404            ENDIF
405         ENDIF
406         !
407         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
408            DO jj = 2, jpjm1
409               DO ji = fs_2, fs_jpim1   ! vector opt.
410                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
411                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
412               END DO
413            END DO
414         END DO
415         !                                 ! Send 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.