program loops
c--Sweep over an unstructured mesh--partitioning by Spectral Bisection
      integer nedge, nnode
      parameter (nedge=3136,nnode=545)
      integer end_pts(3136,2)
      real*8  flux(545), econs
      real*8  w(545,4),p(545) 
      real*8  cc(545)
      real*8  x(545), y(545), z(545)
      integer i, j
      character*80 filen

CHPF$  PROCESSORS P(4)
C start with HPF BLOCK distribution
CHPF$    TEMPLATE reg(545), reg2(3136)
CHPF$    DISTRIBUTE reg(MBLOCK), reg2(MBLOCK)
CHPF$    ALIGN (:,*) WITH reg2 ::end_pts
CHPF$    ALIGN (:,*) WITH reg  ::w
CHPF$    ALIGN (:) WITH reg  :: p, cc, flux
CHPF$    ALIGN (:) WITH reg  ::x, y, z
      filen = 'm6.545.part'
      call read_mesh(end_pts,myvals,X,Y,Z,3136,545,filen)

C  Recursive Spectral  Bisection (RSB)
CUMD$    CONSTRUCT G (545,LINK(3136,end_pts(1,1),end_pts(1,2)))
CUMD$    SET distfmt BY PARTITIONING G USING RSB
CUMD$    REDISTRIBUTE reg(distfmt)

      econs = 1.0
      forall(i=1:545) flux(i) = 0.0
      forall(i=1:545) w(i,1)  = 1.0
      forall(i=1:545) w(i,2)  = 1.0
      forall(i=1:545) w(i,3)  = 1.0
      forall(i=1:545) w(i,4)  = 1.0
      forall(i=1:545) cc(i)   = 1.0
      forall(i=1:545) p(i)    = 1.0

      do 100 j = 1, 5
C
C Sweep over edges
C
      mforall(i=1:3136:1)
      flux(end_pts(i,1)) = flux(end_pts(i,1))
     .                + ((econs*w(end_pts(i,1),2) +
     .                   econs*w(end_pts(i,1),3) +
     .                   econs*w(end_pts(i,1),4)) / w(end_pts(i,1),1) +
     .                   (econs*w(end_pts(i,2),2) +
     .                   econs*w(end_pts(i,2),3)  +
     .                   econs*w(end_pts(i,2),4)) / w(end_pts(i,2),1))/2
     .                + cc(end_pts(i,1))*econs

       flux(end_pts(i,2)) = flux(end_pts(i,2))
     .                + ((econs*w(end_pts(i,1),2) +
     .                   econs*w(end_pts(i,1),3) +
     .                   econs*w(end_pts(i,1),4)) / w(end_pts(i,1),1) +
     .                   (econs*w(end_pts(i,2),2) +
     .                   econs*w(end_pts(i,2),3)  +
     .                   econs*w(end_pts(i,2),4)) / w(end_pts(i,2),1))/2
     .                + cc(end_pts(i,2))*econs
       endforall

  100 end do

      end