!--------------------------------------------------------------------------
   !---------------------------- advect_horiz ----------------------------------

   !$acc data create(qflux(:,:)) present(qi(:,:), cc(:,:,:), gradq3d(:,:,:), mass(:,:), hfluxt(:,:), Ai(:), xyz_i(:,:)) async

   ! 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
   !$acc parallel loop collapse(2) async
   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
            ij_tmp = ij
         ELSE
            ij_tmp = ij+t_right
         END IF

         qe = qi(ij_tmp,l)
         qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
         qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
         qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)

         qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe

         IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
            ij_tmp = ij
         ELSE
            ij_tmp = ij+t_lup
         END IF

         qe = qi(ij_tmp,l)
         qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
         qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
         qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)

         qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe

         IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
            ij_tmp = ij
         ELSE
            ij_tmp = ij+t_ldown
         END IF

         qe = qi(ij_tmp,l)
         qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
         qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
         qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)

         qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
      END DO
   END DO

   IF(diagflux_on) THEN
     !$acc parallel loop collapse(2) copy(qfluxt(:,:)) async
     DO l = ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin_ext, ij_end_ext
           qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
           qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
           qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
        END DO
     END DO
   END IF

   ! update q and, if update_mass, update mass
   !$acc parallel loop collapse(2) async
   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
   
   !$acc end data

   !---------------------------- advect_horiz ----------------------------------
   !--------------------------------------------------------------------------