!---------------------------------------
!             Read tables of load balance
!             ANalyze the results
Program Read_Balance
Integer, Parameter :: ndiv_x=5 ,ndiv_y=5 ,ndiv_z=5,  &
                      ng   = 256,                    &
                      n_nodes = ndiv_x*ndiv_y*ndiv_z, &
                      iStepB = 400
Integer            :: l_divX( ndiv_x +1), l_divY(ndiv_x,ndiv_y+1), &
                      l_divZ(ndiv_x,ndiv_y,ndiv_z +1),             &
                      ldx(ndiv_x+1)
Real               :: w_node(ndiv_x,ndiv_y,ndiv_z),wx(ndiv_x),     &
                      wy(ndiv_x,ndiv_y),                           &
                      ww(ng),wf(ng)
Character          :: Line*80,str(4)*10


Open(1,file='load_balance.log')
   Read(1,'(a)')      Line
   write(*,'(2a)') 'Line=',Line
 20   read(1,*,err=100,end=100) iStep,str(1),str(2),Tmax,str(3),Tmin,str(4),Tav, &
                w_node
      read(1,*)
      read(1,*)str(1),  (l_divX(i),i=2,ndiv_x)
      read(1,*)str(2),                                &
                       ((l_divY(i,j),j=2,ndiv_y),i=1,ndiv_x)
      write(*,'(120("-"))')
      write(*,'(i5,2x,a,3f9.2)')iStep,str(2),Tmax,Tmin,Tav
      write(*,'((5(5f5.2,3x)))') w_node
      write(*,'(a,10i4)')str(1),  (l_divX(i),i=2,ndiv_x)
      write(*,'(a,5(4i4,2x))')str(2),                                &
                       ((l_divY(i,j),j=2,ndiv_y),i=1,ndiv_x)
      l_divX(1) =1 
      l_divX(ndiv_x+1) = ng+1

      Do i=1,ndiv_x
        wx(i) = sum(w_node(i,:,:))
         iLeft =l_divX(i)    ! find width of i-nodes
         iRight=l_divX(i+1)
         aver = wx(i)/MAX(iRight-iLeft,1) ! average weight
                                                   ! of a cell in i-th strip
         Do j=iLeft ,iRight-1           ! set weights for ng elements
            ww(j) =aver
         EndDo 
         write(*,'(3i4,2g12.4)') i,iLeft,iRight,aver,wx(i)
      EndDo
      Do i=2,ng-1
         wf(i) = (ww(i-1)+ww(i)+ww(i+1))/3.
      EndDo 
      wf(1) = (ww(ng)+ww(1)+ww(2))/3.
      wf(ng) = (ww(ng)+ww(1)+ww(ng-1))/3.
      Do i=1,ng
         ww(i) = wf(i)
      EndDo 
      wideal = sum(wx)/ndiv_x                  ! x- direction: set l_divx
      write(*,*) ' ideal =',wideal
      wc          = 0.
      wcc        = 0.
      node_c  = 1
      node_w = 0
      Do i=1,ng
         wc          = wc          +ww(i)
        wcc         = wcc         +ww(i)
         node_w = node_w +1
         If(wc.gt.wideal.or.wc.gt.wideal-ww(i+1)/2.)Then 
            node_c = node_c +1
            ldx(node_c) = i
!               write(*,'(" X:   weight=",F7.2," ideal=",F7.2,
!     &                      " left=",F7.2)')wc,wideal,wwx-wcc
            wc       = 0
            node_w = 0
         EndIf 
      EndDo
!                       Ldx is ideal position of boundaries
!                       Move l_divX toward Ldx, but limit
!                       displacements by iStepB         
      write(*,'(a,100i10)') 'OldivX  =',l_divX
      Do i=2,ndiv_x
         idispl = MIN(iStepB,abs(l_divX(i)-ldx(i)))
         l_divX(i) =l_divX(i) + SIGN(idispl,ldx(i)-l_divX(i))
      EndDo 
      l_divX(1) = 1                          ! set boundary values for l_div'x
      l_divX(ndiv_x+1) = ng+1
      write(*,'(a,10i10)') 'New divX=',l_divX

      

      Do j=1,ndiv_y
      Do i=1,ndiv_x
        wy(i,j) = sum(w_node(i,j,:))
      EndDo
      EndDo
      Do i=1,ndiv_x
         write(*,'(5x,a,f9.2,a,100f9.2)')'sum X=',wx(i),' Y=',wy(i,:)
      EndDo
      goto 20



100 Stop
End Program Read_Balance