!--------------------------------------------------------------------------
   !---------------------------- advect_horiz ----------------------------------
   ! evaluate tracer field at point cc using piecewise linear reconstruction
   ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
   ! sign*hfluxt>0 iff outgoing
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         IF(ne_right*hfluxt(ij+u_right,l)>0.) THEN
            qe = qi(ij,l)
            qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
            qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
            qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
         ELSE
            qe = qi(ij+t_right,l)
            qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij+t_right,1))*gradq3d(ij+t_right,l,1)
            qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij+t_right,2))*gradq3d(ij+t_right,l,2)
            qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij+t_right,3))*gradq3d(ij+t_right,l,3)
         END IF
         qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
         IF(diagflux_on) qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
         IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
            qe = qi(ij,l)
            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
         ELSE
            qe = qi(ij+t_lup,l)
            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij+t_lup,1))*gradq3d(ij+t_lup,l,1)
            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij+t_lup,2))*gradq3d(ij+t_lup,l,2)
            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij+t_lup,3))*gradq3d(ij+t_lup,l,3)
         END IF
         qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
         IF(diagflux_on) qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
         IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
            qe = qi(ij,l)
            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
         ELSE
            qe = qi(ij+t_ldown,l)
            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij+t_ldown,1))*gradq3d(ij+t_ldown,l,1)
            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij+t_ldown,2))*gradq3d(ij+t_ldown,l,2)
            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij+t_ldown,3))*gradq3d(ij+t_ldown,l,3)
         END IF
         qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
         IF(diagflux_on) qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
      END DO
   END DO
   ! update q and, if update_mass, update mass
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         dmass=0.
         dq=0.
         dmass = dmass + ne_rup*hfluxt(ij+u_rup,l)
         dq = dq + ne_rup*qflux(ij+u_rup,l)
         dmass = dmass + ne_lup*hfluxt(ij+u_lup,l)
         dq = dq + ne_lup*qflux(ij+u_lup,l)
         dmass = dmass + ne_left*hfluxt(ij+u_left,l)
         dq = dq + ne_left*qflux(ij+u_left,l)
         dmass = dmass + ne_ldown*hfluxt(ij+u_ldown,l)
         dq = dq + ne_ldown*qflux(ij+u_ldown,l)
         dmass = dmass + ne_rdown*hfluxt(ij+u_rdown,l)
         dq = dq + ne_rdown*qflux(ij+u_rdown,l)
         dmass = dmass + ne_right*hfluxt(ij+u_right,l)
         dq = dq + ne_right*qflux(ij+u_right,l)
         newmass = mass(ij,l) - dmass/Ai(ij)
         qi(ij,l) = (qi(ij,l)*mass(ij,l)-dq/Ai(ij)) / newmass
         IF(update_mass) mass(ij,l)=newmass
      END DO
   END DO
   !---------------------------- advect_horiz ----------------------------------
   !--------------------------------------------------------------------------