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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago

#1613: vvl by default, step III: Merge with the trunk (free surface simplification) (see wiki)
  • Property svn:keywords set to Id
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 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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of 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), 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._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            zdt =  p2dt(jk)
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.   
171                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
172                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
173                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
174                  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
175                  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
176               END DO
177            END DO
178         END DO 
179         !--- Lateral boundary conditions
180         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
181         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
182
183         !--- QUICKEST scheme
184         CALL quickest( zfu, zfd, zfc, zwx )
185         !
186         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
187         DO jk = 1, jpkm1 
188            DO jj = 2, jpjm1
189               DO ji = fs_2, fs_jpim1   ! vector opt.               
190                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
191               END DO
192            END DO
193         END DO
194         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
195
196         !
197         ! Tracer flux on the x-direction
198         DO jk = 1, jpkm1 
199            !
200            DO jj = 2, jpjm1
201               DO ji = fs_2, fs_jpim1   ! vector opt.               
202                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
203                  !--- If the second ustream point is a land point
204                  !--- the flux is computed by the 1st order UPWIND scheme
205                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
206                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
207                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
208               END DO
209            END DO
210         END DO
211         !
212         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
213         !
214         ! Computation of the trend
215         DO jk = 1, jpkm1 
216            DO jj = 2, jpjm1
217               DO ji = fs_2, fs_jpim1   ! vector opt. 
218                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
219                  ! horizontal advective trends
220                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
221                  !--- add it to the general tracer trends
222                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
223               END DO
224            END DO
225         END DO
226         !                                 ! trend diagnostics
227         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
228         !
229      END DO
230      !
231      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
232      !
233   END SUBROUTINE tra_adv_qck_i
234
235
236   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
237      &                                        ptb, ptn, pta, kjpt )
238      !!----------------------------------------------------------------------
239      !!
240      !!----------------------------------------------------------------------
241      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
242      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
243      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
244      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
245      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
246      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
247      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
248      !!
249      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
250      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
251      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
252      !----------------------------------------------------------------------
253      !
254      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
255      !
256      !                                                          ! ===========
257      DO jn = 1, kjpt                                            ! tracer loop
258         !                                                       ! ===========
259         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
260         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
261         !                                                 
262         DO jk = 1, jpkm1                               
263            !                                             
264            !--- Computation of the ustream and downstream value of the tracer and the mask
265            DO jj = 2, jpjm1
266               DO ji = fs_2, fs_jpim1   ! vector opt.
267                  ! Upstream in the x-direction for the tracer
268                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
269                  ! Downstream in the x-direction for the tracer
270                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
271               END DO
272            END DO
273         END DO
274         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
275
276         
277         !
278         ! Horizontal advective fluxes
279         ! ---------------------------
280         !
281         DO jk = 1, jpkm1                             
282            DO jj = 2, jpjm1
283               DO ji = fs_2, fs_jpim1   ! vector opt.         
284                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
285                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
286               END DO
287            END DO
288         END DO
289         !
290         DO jk = 1, jpkm1 
291            zdt =  p2dt(jk)
292            DO jj = 2, jpjm1
293               DO ji = fs_2, fs_jpim1   ! vector opt.   
294                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
295                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
296                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
297                  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
298                  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
299               END DO
300            END DO
301         END DO
302
303         !--- Lateral boundary conditions
304         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
305         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
306
307         !--- QUICKEST scheme
308         CALL quickest( zfu, zfd, zfc, zwy )
309         !
310         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
311         DO jk = 1, jpkm1 
312            DO jj = 2, jpjm1
313               DO ji = fs_2, fs_jpim1   ! vector opt.               
314                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
315               END DO
316            END DO
317         END DO
318         !--- Lateral boundary conditions
319         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
320         !
321         ! Tracer flux on the x-direction
322         DO jk = 1, jpkm1 
323            !
324            DO jj = 2, jpjm1
325               DO ji = fs_2, fs_jpim1   ! vector opt.               
326                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
327                  !--- If the second ustream point is a land point
328                  !--- the flux is computed by the 1st order UPWIND scheme
329                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
330                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
331                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
332               END DO
333            END DO
334         END DO
335         !
336         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
337         !
338         ! Computation of the trend
339         DO jk = 1, jpkm1 
340            DO jj = 2, jpjm1
341               DO ji = fs_2, fs_jpim1   ! vector opt. 
342                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
343                  ! horizontal advective trends
344                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
345                  !--- add it to the general tracer trends
346                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
347               END DO
348            END DO
349         END DO
350         !                                 ! trend diagnostics
351         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
352         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
353         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
354           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
355           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
356         ENDIF
357         !
358      END DO
359      !
360      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
361      !
362   END SUBROUTINE tra_adv_qck_j
363
364
365   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
366     &                                    ptn, pta, kjpt )
367      !!----------------------------------------------------------------------
368      !!
369      !!----------------------------------------------------------------------
370      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
371      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
372      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
373      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
374      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
375      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
376      !
377      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
378      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
379      !!----------------------------------------------------------------------
380      !
381      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
382      !
383      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
384      zwz(:,:,jpk) = 0._wp
385      !
386      !                                                          ! ===========
387      DO jn = 1, kjpt                                            ! tracer loop
388         !                                                       ! ===========
389         !
390         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
391            DO jj = 2, jpjm1
392               DO ji = fs_2, fs_jpim1   ! vector opt.
393                  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)
394               END DO
395            END DO
396         END DO
397         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
398            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
399               DO jj = 1, jpj
400                  DO ji = 1, jpi
401                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
402                  END DO
403               END DO   
404            ELSE                                   ! no ocean cavities (only ocean surface)
405               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
406            ENDIF
407         ENDIF
408         !
409         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
410            DO jj = 2, jpjm1
411               DO ji = fs_2, fs_jpim1   ! vector opt.
412                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
413                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
414               END DO
415            END DO
416         END DO
417         !                                 ! Send trends for diagnostic
418         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
419         !
420      END DO
421      !
422      CALL wrk_dealloc( jpi,jpj,jpk,   zwz )
423      !
424   END SUBROUTINE tra_adv_cen2_k
425
426
427   SUBROUTINE quickest( pfu, pfd, pfc, puc )
428      !!----------------------------------------------------------------------
429      !!
430      !! ** Purpose :  Computation of advective flux with Quickest scheme
431      !!
432      !! ** Method :   
433      !!----------------------------------------------------------------------
434      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
435      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
436      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
437      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
438      !!
439      INTEGER  ::  ji, jj, jk               ! dummy loop indices
440      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
441      REAL(wp) ::  zc, zcurv, zfho          !   -      -
442      !----------------------------------------------------------------------
443      !
444      IF( nn_timing == 1 )  CALL timing_start('quickest')
445      !
446      DO jk = 1, jpkm1
447         DO jj = 1, jpj
448            DO ji = 1, jpi
449               zc     = puc(ji,jj,jk)                         ! Courant number
450               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
451               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
452               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
453               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
454               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
455               !
456               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
457               zcoef2 = ABS( zcoef1 )
458               zcoef3 = ABS( zcurv )
459               IF( zcoef3 >= zcoef2 ) THEN
460                  zfho = pfc(ji,jj,jk) 
461               ELSE
462                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
463                  IF( zcoef1 >= 0. ) THEN
464                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
465                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
466                  ELSE
467                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
468                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
469                  ENDIF
470               ENDIF
471               puc(ji,jj,jk) = zfho
472            END DO
473         END DO
474      END DO
475      !
476      IF( nn_timing == 1 )  CALL timing_stop('quickest')
477      !
478   END SUBROUTINE quickest
479
480   !!======================================================================
481END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.