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

Last change on this file since 10806 was 10806, checked in by davestorkey, 21 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.

  • 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 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, ktlev, cdtype, p2dt, pu, pv, pw,      &
50      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_adv_qck  ***
53      !!
54      !! ** Purpose :   Compute the now trend due to the advection of tracers
55      !!      and add it to the general trend of passive tracer equations.
56      !!
57      !! ** Method :   The advection is evaluated by a third order scheme
58      !!             For a positive velocity u :              u(i)>0
59      !!                                          |--FU--|--FC--|--FD--|------|
60      !!                                             i-1    i      i+1   i+2
61      !!
62      !!             For a negative velocity u :              u(i)<0
63      !!                                          |------|--FD--|--FC--|--FU--|
64      !!                                             i-1    i      i+1   i+2
65      !!             where  FU is the second upwind point
66      !!                    FD is the first douwning point
67      !!                    FC is the central point (or the first upwind point)
68      !!
69      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
70      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
71      !!
72      !!         dt = 2*rdtra and the scalar values are tb and sb
73      !!
74      !!       On the vertical, the simple centered scheme used pt_lev2
75      !!
76      !!               The fluxes are bounded by the ULTIMATE limiter to
77      !!             guarantee the monotonicity of the solution and to
78      !!            prevent the appearance of spurious numerical oscillations
79      !!
80      !! ** Action : - update pt_rhs  with the now advective tracer trends
81      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
82      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
83      !!
84      !! ** Reference : Leonard (1979, 1991)
85      !!----------------------------------------------------------------------
86      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
87      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
88      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
89      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
90      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
91      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend
95      !!----------------------------------------------------------------------
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      l_ptr = .FALSE.
106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
107      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
108      !
109      !
110      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
111      CALL tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu, pt_lev1, pt_lev2, pt_rhs, kjpt ) 
112      CALL tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv, pt_lev1, pt_lev2, pt_rhs, kjpt ) 
113
114      !        ! vertical fluxes are computed with the 2nd order centered scheme
115      CALL tra_adv_cen2_k( kt, ktlev, cdtype, pw,         pt_lev2, pt_rhs, kjpt )
116      !
117   END SUBROUTINE tra_adv_qck
118
119
120   SUBROUTINE tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu,                  &
121      &                                        pt_lev1, pt_lev2, pt_rhs, kjpt   )
122      !!----------------------------------------------------------------------
123      !!
124      !!----------------------------------------------------------------------
125      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
126      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
127      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
128      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
129      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
130      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu        ! i-velocity components
131      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2   ! before and now tracer fields
132      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend
133      !!
134      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
135      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
136      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zfu, zfc, zfd
137      !----------------------------------------------------------------------
138      !
139      !                                                          ! ===========
140      DO jn = 1, kjpt                                            ! tracer loop
141         !                                                       ! ===========
142         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
143         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
144         !
145!!gm why not using a SHIFT instruction...
146         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
147            DO jj = 2, jpjm1
148               DO ji = fs_2, fs_jpim1   ! vector opt.
149                  zfc(ji,jj,jk) = pt_lev1(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
150                  zfd(ji,jj,jk) = pt_lev1(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
151               END DO
152            END DO
153         END DO
154         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
155         
156         !
157         ! Horizontal advective fluxes
158         ! ---------------------------
159         DO jk = 1, jpkm1                             
160            DO jj = 2, jpjm1
161               DO ji = fs_2, fs_jpim1   ! vector opt.         
162                  zdir = 0.5 + SIGN( 0.5, pu(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0
163                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
164               END DO
165            END DO
166         END DO
167         !
168         DO jk = 1, jpkm1 
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.   
171                  zdir = 0.5 + SIGN( 0.5, pu(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0
172                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,ktlev)
173                  zwx(ji,jj,jk)  = ABS( pu(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
174                  zfc(ji,jj,jk)  = zdir * pt_lev1(ji  ,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji+1,jj,jk,jn)  ! FC in the x-direction for T
175                  zfd(ji,jj,jk)  = zdir * pt_lev1(ji+1,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(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_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  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( 'traadv_qck', 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, pu(ji,jj,jk) )   ! if pu > 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) * pu(ji,jj,jk)
207               END DO
208            END DO
209         END DO
210         !
211         CALL lbc_lnk( 'traadv_qck', 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(ji,jj,jk,ktlev)
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                  pt_rhs(ji,jj,jk,jn) = pt_rhs(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, pu, pt_lev2(:,:,:,jn) )
227         !
228      END DO
229      !
230   END SUBROUTINE tra_adv_qck_i
231
232
233   SUBROUTINE tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv,                &
234      &                                        pt_lev1, pt_lev2, pt_rhs, kjpt )
235      !!----------------------------------------------------------------------
236      !!
237      !!----------------------------------------------------------------------
238      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
239      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
240      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
241      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
242      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
243      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pv        ! j-velocity components
244      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2   ! before and now tracer fields
245      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend
246      !!
247      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
248      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
249      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
250      !----------------------------------------------------------------------
251      !
252      !                                                          ! ===========
253      DO jn = 1, kjpt                                            ! tracer loop
254         !                                                       ! ===========
255         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
256         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
257         !                                                 
258         DO jk = 1, jpkm1                               
259            !                                             
260            !--- Computation of the ustream and downstream value of the tracer and the mask
261            DO jj = 2, jpjm1
262               DO ji = fs_2, fs_jpim1   ! vector opt.
263                  ! Upstream in the x-direction for the tracer
264                  zfc(ji,jj,jk) = pt_lev1(ji,jj-1,jk,jn)
265                  ! Downstream in the x-direction for the tracer
266                  zfd(ji,jj,jk) = pt_lev1(ji,jj+1,jk,jn)
267               END DO
268            END DO
269         END DO
270         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
271
272         
273         !
274         ! Horizontal advective fluxes
275         ! ---------------------------
276         !
277         DO jk = 1, jpkm1                             
278            DO jj = 2, jpjm1
279               DO ji = fs_2, fs_jpim1   ! vector opt.         
280                  zdir = 0.5 + SIGN( 0.5, pv(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0
281                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
282               END DO
283            END DO
284         END DO
285         !
286         DO jk = 1, jpkm1 
287            DO jj = 2, jpjm1
288               DO ji = fs_2, fs_jpim1   ! vector opt.   
289                  zdir = 0.5 + SIGN( 0.5, pv(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0
290                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,ktlev)
291                  zwy(ji,jj,jk)  = ABS( pv(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
292                  zfc(ji,jj,jk)  = zdir * pt_lev1(ji,jj  ,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj+1,jk,jn)  ! FC in the x-direction for T
293                  zfd(ji,jj,jk)  = zdir * pt_lev1(ji,jj+1,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj  ,jk,jn)  ! FD in the x-direction for T
294               END DO
295            END DO
296         END DO
297
298         !--- Lateral boundary conditions
299         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
300
301         !--- QUICKEST scheme
302         CALL quickest( zfu, zfd, zfc, zwy )
303         !
304         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
305         DO jk = 1, jpkm1 
306            DO jj = 2, jpjm1
307               DO ji = fs_2, fs_jpim1   ! vector opt.               
308                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
309               END DO
310            END DO
311         END DO
312         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions
313         !
314         ! Tracer flux on the x-direction
315         DO jk = 1, jpkm1 
316            !
317            DO jj = 2, jpjm1
318               DO ji = fs_2, fs_jpim1   ! vector opt.               
319                  zdir = 0.5 + SIGN( 0.5, pv(ji,jj,jk) )   ! if pu > 0 : zdir = 1 otherwise zdir = 0
320                  !--- If the second ustream point is a land point
321                  !--- the flux is computed by the 1st order UPWIND scheme
322                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
323                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
324                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pv(ji,jj,jk)
325               END DO
326            END DO
327         END DO
328         !
329         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
330         !
331         ! Computation of the trend
332         DO jk = 1, jpkm1 
333            DO jj = 2, jpjm1
334               DO ji = fs_2, fs_jpim1   ! vector opt. 
335                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
336                  ! horizontal advective trends
337                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
338                  !--- add it to the general tracer trends
339                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztra
340               END DO
341            END DO
342         END DO
343         !                                 ! trend diagnostics
344         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt_lev2(:,:,:,jn) )
345         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
346         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
347         !
348      END DO
349      !
350   END SUBROUTINE tra_adv_qck_j
351
352
353   SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pw,           &
354     &                                    pt_lev2, pt_rhs, kjpt )
355      !!----------------------------------------------------------------------
356      !!
357      !!----------------------------------------------------------------------
358      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
359      INTEGER                              , INTENT(in   ) ::   ktlev    ! time level index for source terms
360      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
361      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
362      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pw      ! vertical velocity
363      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev2      ! before and now tracer fields
364      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs      ! tracer trend
365      !
366      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
367      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
368      !!----------------------------------------------------------------------
369      !
370      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
371      zwz(:,:,jpk) = 0._wp
372      !
373      !                                                          ! ===========
374      DO jn = 1, kjpt                                            ! tracer loop
375         !                                                       ! ===========
376         !
377         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
378            DO jj = 2, jpjm1
379               DO ji = fs_2, fs_jpim1   ! vector opt.
380                  zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
381               END DO
382            END DO
383         END DO
384         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
385            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
386               DO jj = 1, jpj
387                  DO ji = 1, jpi
388                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface
389                  END DO
390               END DO   
391            ELSE                                   ! no ocean cavities (only ocean surface)
392               zwz(:,:,1) = pw(:,:,1) * pt_lev2(:,:,1,jn)
393            ENDIF
394         ENDIF
395         !
396         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
397            DO jj = 2, jpjm1
398               DO ji = fs_2, fs_jpim1   ! vector opt.
399                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
400                     &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
401               END DO
402            END DO
403         END DO
404         !                                 ! Send trends for diagnostic
405         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, pt_lev2(:,:,:,jn) )
406         !
407      END DO
408      !
409   END SUBROUTINE tra_adv_cen2_k
410
411
412   SUBROUTINE quickest( pfu, pfd, pfc, puc )
413      !!----------------------------------------------------------------------
414      !!
415      !! ** Purpose :  Computation of advective flux with Quickest scheme
416      !!
417      !! ** Method :   
418      !!----------------------------------------------------------------------
419      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
420      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
421      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
422      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
423      !!
424      INTEGER  ::  ji, jj, jk               ! dummy loop indices
425      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
426      REAL(wp) ::  zc, zcurv, zfho          !   -      -
427      !----------------------------------------------------------------------
428      !
429      DO jk = 1, jpkm1
430         DO jj = 1, jpj
431            DO ji = 1, jpi
432               zc     = puc(ji,jj,jk)                         ! Courant number
433               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
434               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
435               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
436               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
437               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
438               !
439               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
440               zcoef2 = ABS( zcoef1 )
441               zcoef3 = ABS( zcurv )
442               IF( zcoef3 >= zcoef2 ) THEN
443                  zfho = pfc(ji,jj,jk) 
444               ELSE
445                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
446                  IF( zcoef1 >= 0. ) THEN
447                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
448                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
449                  ELSE
450                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
451                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
452                  ENDIF
453               ENDIF
454               puc(ji,jj,jk) = zfho
455            END DO
456         END DO
457      END DO
458      !
459   END SUBROUTINE quickest
460
461   !!======================================================================
462END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.