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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90 @ 10880

Last change on this file since 10880 was 10880, checked in by davestorkey, 3 years ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps:

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