program loops
c
C-- Data Partitioning Using Weighted Inertial Bisection
c
c--Sweep over an unstructured mesh
c
       integer nedge, nnode
       parameter (nedge=3136,nnode=545)


      integer end_pts(3136,2)
      integer wt(545)
      real*8  flux(545), econs
      real*8  w(545,4),p(545) 
      real*8  cc(545)
      real*8  ftemp, flux_sum
      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, wt
CHPF$    ALIGN (:) WITH reg  ::x, y, z


      filen = 'm6.545.part'
      call read_mesh2(end_pts,myvals,X,Y,Z,3136,545,filen, wt)

C  Now, obtain a new data distribution using
C  Weighted Inertial Bisection (WIBS) partitioner
C  and redistribute data
C
CUMD$    CONSTRUCT G (545, LOAD(wt), GEOMETRY(3,X,Y,Z))
CUMD$    SET distfmt BY PARTITIONING G USING WIBS
CUMD$    REDISTRIBUTE reg(distfmt)
      
      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


      econs = 1.0
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