MODULE traadv_qck
   !!==============================================================================
   !!                       ***  MODULE  traadv_qck  ***
   !! Ocean active tracers:  horizontal & vertical advective trend
   !!==============================================================================
   !! History :  3.0  !  2008-07  (G. Reffray)  Original code
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   tra_adv_qck      : update the tracer trend with the horizontal advection
   !!                      trends using a 3rd order finite difference scheme
   !!   tra_adv_qck_i  : 
   !!   tra_adv_qck_j  : 
   !!   tra_adv_cen2_k : 2nd centered scheme for the vertical advection
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and active tracers
   USE dom_oce         ! ocean space and time domain
   USE trdmod          ! ocean active tracers trends 
   USE trdmod_oce      ! ocean variables trends
   USE trabbl          ! advective term in the BBL
   USE lib_mpp         ! distribued memory computing
   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
   USE dynspg_oce      ! surface pressure gradient variables
   USE in_out_manager  ! I/O manager
   USE diaptr          ! poleward transport diagnostics
   USE prtctl          ! Print control

   IMPLICIT NONE
   PRIVATE

   PUBLIC   tra_adv_qck   ! routine called by step.F90

   REAL(wp), DIMENSION(jpi,jpj) ::   btr2
   REAL(wp)                     ::   r1_6

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
   !! $Id$
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE tra_adv_qck  ***
      !!
      !! ** Purpose :   Compute the now trend due to the advection of tracers
      !!      and add it to the general trend of passive tracer equations.
      !!
      !! ** Method :   The advection is evaluated by a third order scheme
      !!             For a positive velocity u :              u(i)>0
      !!                                          |--FU--|--FC--|--FD--|------|
      !!                                             i-1    i      i+1   i+2
      !!
      !!             For a negative velocity u :              u(i)<0
      !!                                          |------|--FD--|--FC--|--FU--|
      !!                                             i-1    i      i+1   i+2
      !!             where  FU is the second upwind point
      !!                    FD is the first douwning point
      !!                    FC is the central point (or the first upwind point)
      !!
      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
      !!
      !!         dt = 2*rdtra and the scalar values are tb and sb
      !!
      !!       On the vertical, the simple centered scheme used tn and sn
      !!
      !!               The fluxes are bounded by the ULTIMATE limiter to
      !!             guarantee the monotonicity of the solution and to
      !!            prevent the appearance of spurious numerical oscillations
      !!
      !! ** Action : - update (ta,sa) with the now advective tracer trends
      !!             - save the trends ('key_trdtra')
      !!
      !! ** Reference : Leonard (1979, 1991)
      !!----------------------------------------------------------------------
      USE oce, ONLY :   ztrdt => ua   ! use ua as workspace
      USE oce, ONLY :   ztrds => va   ! use va as workspace
      !!
      INTEGER , INTENT(in)                         ::  kt  ! ocean time-step index
      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun ! effective ocean velocity, u_component
      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pvn ! effective ocean velocity, v_component
      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pwn ! effective ocean velocity, w_component
      !!
      INTEGER  ::   ji, jj, jk                          ! dummy loop indices
      REAL(wp) ::   z_hdivn_x, z_hdivn_y, z_hdivn       ! temporary scalars
      REAL(wp) ::   zbtr, z2                            !    "         "
      !!----------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme'
         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
         IF(lwp) WRITE(numout,*)
         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
         r1_6      = 1. / 6.
      ENDIF

      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
      ELSE                                        ;    z2 = 2.
      ENDIF

      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
      !---------------------------------------------------------------------------

      CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2)
      CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2)

      IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)

      CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2)
      CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2)  

      IF( l_trdtra ) THEN
         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
         !
         ztrdt(:,:,:) = ta(:,:,:)    ! Save the horizontal up-to-date ta/sa trends
         ztrds(:,:,:) = sa(:,:,:)
      END IF

      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' qck had  - Ta: ', mask1=tmask, &
         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )

      ! II. The vertical fluxes are computed with the 2nd order centered scheme
      !-------------------------------------------------------------------------
      !
      CALL tra_adv_cen2_k( pwn, tn, ta )
      CALL tra_adv_cen2_k( pwn, sn, sa )
      !
      !Save the vertical advective trends for diagnostic
      ! ----------------------------------------------------
      IF( l_trdtra )   THEN
         ! Recompute the vertical advection zta & zsa trends computed
         ! at the step 2. above in making the difference between the new
         ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends

         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
#if defined key_zco
                  zbtr      = btr2(ji,jj)
                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
#else
                  zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk)
                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
#endif
                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
               END DO
            END DO
         END DO
         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
      ENDIF

      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' qck zad  - Ta: ', mask1=tmask, &
         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
      !
   END SUBROUTINE tra_adv_qck


   SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 )
      !!----------------------------------------------------------------------
      !!
      !!----------------------------------------------------------------------
      REAL,     INTENT(in)                            :: z2
      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran  ! horizontal effective velocity
      REAL(wp), INTENT(out)  , DIMENSION(jpi,jpj,jpk) :: ztrdtra
      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
      !
      INTEGER  :: ji, jj, jk
      REAL(wp) :: za, zbtr, dir, dx, dt     ! temporary scalars
      REAL(wp) :: z_hdivn_x
      REAL(wp), DIMENSION(jpi,jpj)     ::  zmask, zupst, zdwst, zc_cfl
      REAL(wp), DIMENSION(jpi,jpj)     ::  zfu, zfc, zfd, zfho, zmskl, zsc_e 
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zflux
      !----------------------------------------------------------------------

      zfu   (:,jpj) = 0.e0   ;   zfc   (:,jpj) = 0.e0
      zfd   (:,jpj) = 0.e0   ;   zc_cfl(:,jpj) = 0.e0
      zsc_e (:,jpj) = 0.e0   ;   zmskl (:,jpj) = 0.e0
      zfho  (:,jpj) = 0.e0
                                                       ! ===============
      DO jk = 1, jpkm1                                 ! Horizontal slab
         !                                             ! ===============
         !--- Computation of the ustream and downstream value of the tracer and the mask
         DO jj = 2, jpjm1
            DO ji = 2, fs_jpim1   ! vector opt.
               ! Upstream in the x-direction for the tracer
               zupst(ji,jj)=tra(ji-1,jj,jk)
               ! Downstream in the x-direction for the tracer
               zdwst(ji,jj)=tra(ji+1,jj,jk)
               ! Mask at the T-points in the x-direction (mask=0 or mask=1)
               zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
            END DO
         END DO
         !
         !--- Lateral boundary conditions 
         CALL lbc_lnk( zupst(:,:), 'T', 1. )
         CALL lbc_lnk( zdwst(:,:), 'T', 1. ) 
         CALL lbc_lnk( zmask(:,:), 'T', 1. ) 
         !
         ! Horizontal advective fluxes
         ! ---------------------------
         !
         dt = z2 * rdttra(jk)
         !--- tracer flux at u-points
         DO jj = 1, jpjm1
            DO ji = 1, jpi
#if defined key_zco
               zsc_e(ji,jj) = e2u(ji,jj)
#else
               zsc_e(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk)
#endif
               dir = 0.5 + sign(0.5,pun(ji,jj,jk))                             ! if pun>0 : dir = 1 otherwise dir = 0
               dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
               zc_cfl (ji,jj) = ABS(pun(ji,jj,jk))*dt/dx                       ! (0<zc_cfl<1 : Courant number on x-direction)

               zfu(ji,jj)   = dir*zupst(ji  ,jj   )+(1-dir)*zdwst(ji+1,jj   )  ! FU in the x-direction for T
               zfc(ji,jj)   = dir*tra  (ji  ,jj,jk)+(1-dir)*tra  (ji+1,jj,jk)  ! FC in the x-direction for T
               zfd(ji,jj)   = dir*tra  (ji+1,jj,jk)+(1-dir)*tra  (ji  ,jj,jk)  ! FD in the x-direction for T
               zmskl(ji,jj) = dir*zmask(ji  ,jj)   +(1-dir)*zmask(ji+1,jj)
           END DO
         END DO
         !
         !--- QUICKEST scheme
         ! Tracer flux on the x-direction
         CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
         !--- If the second ustream point is a land point
         !--- the flux is computed by the 1st order UPWIND scheme
         zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:)
         !--- Computation of fluxes
         zflux(:,:,jk) = zsc_e(:,:)*pun(:,:,jk)*zfho(:,:)
         !
         !--- Tracer flux divergence at t-point added to the general trend
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               !--- horizontal advective trends
#if defined key_zco
               zbtr = btr2(ji,jj)
#else
               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)              
#endif
               za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) )
               !--- add it to the general tracer trends
               traa(ji,jj,jk) = traa(ji,jj,jk) + za
            END DO
         END DO
         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============
      !
      !  Save the horizontal advective trends for diagnostic
      ! -----------------------------------------------------
      IF( l_trdtra ) THEN
         ! T/S ZONAL advection trends
         ztrdtra(:,:,:) = 0.e0
         !
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
                  !   N.B. This computation is not valid along OBCs (if any)
#if defined key_zco
                  zbtr      = btr2(ji,jj)
                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
#else
                  zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk)
                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
#endif
                  ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) + tran(ji,jj,jk) * z_hdivn_x
               END DO
            END DO
         END DO
      END IF

   END SUBROUTINE tra_adv_qck_i


   SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 )
      !!----------------------------------------------------------------------
      !!
      !!----------------------------------------------------------------------
      INTEGER,  INTENT(in)                            :: kt              ! ocean time-step index
      REAL,     INTENT(in)                            :: z2
      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pvn, tra, tran  ! horizontal effective velocity
      REAL(wp), INTENT(out)  , DIMENSION(jpj)         :: trd_adv
      REAL(wp), INTENT(out)  , DIMENSION(jpi,jpj,jpk) :: ztrdtra
      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
      !!
      INTEGER  :: ji, jj, jk
      REAL(wp) :: za, zbtr, dir, dx, dt     ! temporary scalars
      REAL(wp) :: z_hdivn_y
      REAL(wp), DIMENSION(jpi,jpj)     ::  zmask, zupst, zdwst, zc_cfl
      REAL(wp), DIMENSION(jpi,jpj)     ::  zfu, zfc, zfd, zfho, zmskl, zsc_e
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zflux
      !----------------------------------------------------------------------
      ! II. Part 2 : y-direction
      !----------------------------------------------------------------------

      zfu   (:,jpj) = 0.e0   ;   zfc   (:,jpj) = 0.e0
      zfd   (:,jpj) = 0.e0   ;   zc_cfl(:,jpj) = 0.e0
      zsc_e (:,jpj) = 0.e0   ;   zmskl (:,jpj) = 0.e0
      zfho  (:,jpj) = 0.e0

                                                       ! ===============
      DO jk = 1, jpkm1                                 ! Horizontal slab
         !                                             ! ===============
         !--- Computation of the ustream and downstream value of the tracer and the mask
         DO jj = 2, jpjm1
            DO ji = 2, fs_jpim1   ! vector opt.
               ! Upstream in the x-direction for the tracer
               zupst(ji,jj)=tra(ji,jj-1,jk)
               ! Downstream in the x-direction for the tracer
               zdwst(ji,jj)=tra(ji,jj+1,jk)
               ! Mask at the T-points in the x-direction (mask=0 or mask=1)
               zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
            END DO
         END DO
         !
         !--- Lateral boundary conditions
         CALL lbc_lnk( zupst(:,:), 'T', 1. )
         CALL lbc_lnk( zdwst(:,:), 'T', 1. )
         CALL lbc_lnk( zmask(:,:), 'T', 1. )
         !
         ! Horizontal advective fluxes
         ! ---------------------------
         !
         dt = z2 * rdttra(jk)
         !--- tracer flux at v-points
         DO jj = 1, jpjm1
            DO ji = 1, jpi
#if defined key_zco
               zsc_e(ji,jj) = e1v(ji,jj)
#else
               zsc_e(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk)
#endif
               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))                             ! if pvn>0 : dir = 1 otherwise dir = 0
               dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
               zc_cfl(ji,jj) = ABS(pvn(ji,jj,jk))*dt/dx                        ! (0<zc_cfl<1 : Courant number on y-direction)

               zfu(ji,jj)   = dir*zupst(ji,jj     )+(1-dir)*zdwst(ji,jj+1   )  ! FU in the y-direction for T
               zfc(ji,jj)   = dir*tra  (ji,jj  ,jk)+(1-dir)*tra  (ji,jj+1,jk)  ! FC in the y-direction for T
               zfd(ji,jj)   = dir*tra  (ji,jj+1,jk)+(1-dir)*tra  (ji,jj  ,jk)  ! FD in the y-direction for T
               zmskl(ji,jj) = dir*zmask(ji,jj     )+(1-dir)*zmask(ji,jj+1)
            END DO
         END DO
         !
         !--- QUICKEST scheme
         ! Tracer flux on the y-direction
         CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
         !--- If the second ustream point is a land point
         !--- the flux is computed by the 1st order UPWIND scheme
         zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:)
         !--- Computation of fluxes
         zflux(:,:,jk) = zsc_e(:,:)*pvn(:,:,jk)*zfho(:,:)
         !
         !--- Tracer flux divergence at t-point added to the general trend
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               !--- horizontal advective trends
#if defined key_zco
               zbtr = btr2(ji,jj)
#else
               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
#endif
               za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) )
               !--- add it to the general tracer trends
               traa(ji,jj,jk) = traa(ji,jj,jk) + za
            END DO
         END DO
         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============
      !
      !  Save the horizontal advective trends for diagnostic
      ! -----------------------------------------------------
      IF( l_trdtra ) THEN
         ! T/S MERIDIONAL advection trends
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
                  !   N.B. This computation is not valid along OBCs (if any)
#if defined key_zco
                  zbtr      = btr2(ji,jj)
                  z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              &
                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
#else
                  zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk)
                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
#endif
                  ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) + tran(ji,jj,jk) * z_hdivn_y
               END DO
            END DO
         END DO
      END IF

      ! "zonal" mean advective heat and salt transport
      ! ----------------------------------------------

      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
         IF( lk_zco ) THEN
            DO jk = 1, jpkm1
               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1   ! vector opt.
                    zflux(ji,jj,jk) = zflux(ji,jj,jk) * fse3v(ji,jj,jk)
                  END DO
               END DO
            END DO
         ENDIF
         trd_adv(:) = ptr_vj( zflux(:,:,:) )
      ENDIF

   END SUBROUTINE tra_adv_qck_j


   SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta )
      !!----------------------------------------------------------------------
      !!
      !!----------------------------------------------------------------------
      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)  :: pwn   ! vertical effective velocity
      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)  :: ptn   ! now tracer
      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)  :: pta   ! tracer general trend
      !!
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zflux   ! 3D workspace
      !!----------------------------------------------------------------------
      !
      !                         !==  Vertical advective fluxes  ==!
      zflux(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
      !
      !                                 ! Surface value
      IF( lk_vvl ) THEN   ;   zflux(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
      ELSE                ;   zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! Constant volume : advective flux through the surface
      ENDIF
      !
      DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) )
            END DO
         END DO
      END DO
      !
      DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               pta(ji,jj,jk) =  pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) )   &
                  &                           /   fse3t(ji,jj,jk)
            END DO
         END DO
      END DO
      !
   END SUBROUTINE tra_adv_cen2_k


   SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl )
      !!----------------------------------------------------------------------
      !!
      !!----------------------------------------------------------------------
      REAL(wp), INTENT(in)  , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl  
      REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho
      REAL(wp)              , DIMENSION(jpi,jpj) :: zcurv, zcoef1, zcoef2, zcoef3  ! temporary scalars
      !
      zcurv (:,:) = fd(:,:) + fu(:,:) - 2.*fc(:,:)
      zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) )
      zcoef2(:,:) = 0.5*fc_cfl(:,:)*( fd(:,:) - fc(:,:) )
      zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*r1_6 )*zcurv(:,:)
      fho   (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:)                         ! phi_f QUICKEST 
      !
      zcoef1(:,:) = fd(:,:) - fu(:,:)                                               ! DEL
      zcoef2(:,:) = ABS( zcoef1(:,:) )                                              ! ABS(DEL)
      zcoef3(:,:) = ABS( zcurv(:,:) )                                               ! ABS(CURV)
      !
      WHERE ( zcoef3(:,:) >= zcoef2(:,:) ) 
        fho(:,:) = fc(:,:)
      ELSEWHERE
        zcoef3(:,:) = fu(:,:) + ( ( fc(:,:) - fu(:,:) )/MAX(fc_cfl(:,:),1.e-9) )    ! phi_REF
          WHERE ( zcoef1(:,:) >= 0.e0 )
            fho(:,:) = MAX(fc(:,:),fho(:,:))
            fho(:,:) = MIN(fho(:,:),MIN(zcoef3(:,:),fd(:,:)))
          ELSEWHERE
            fho(:,:) = MIN(fc(:,:),fho(:,:))
            fho(:,:) = MAX(fho(:,:),MAX(zcoef3(:,:),fd(:,:)))
          ENDWHERE
      ENDWHERE
      !
   END SUBROUTINE quickest

   !!======================================================================
END MODULE traadv_qck