plot file data
c-------------------------------------
          write (n4,1000) loopx, loopy, jinc(1) 
c                                 output graphics data
         do i = 1, loopx

            if (i.ne.1.and.igrd(i,1).eq.0) igrd(i,1) = igrd(i-1,1)

            kst = 1

20          jst = kst
            if (i.ne.1.and.igrd(i,jst).eq.0) igrd(i,jst) = igrd(i-1,jst)
            kd = igrd(i,jst)
            ltic = -1

            do j = jst, loopy

               if (i.ne.1.and.igrd(i,j).eq.0) igrd(i,j) = igrd(i-1,j)

               if (igrd(i,j).eq.0.or.igrd(i,j).eq.kd) then
                  ltic = ltic + 1
                  if (j.eq.loopy) then 
                     write (n4,1000) ltic,kd
                     kst = 1
                  end if 
               else 
                  write (n4,1000) ltic,kd
                  kst = j
                  goto 20 
               end if 
            end do 
         end do     

c                                 write assemblage list
      write (n4,*) iasct

      do i = 1, iasct
         write (n4,*) iavar(1,i),iavar(2,i),iavar(3,i)
         write (n4,*) (idasls(j,i), j = 1, iavar(3,i))
      end do 

c----------------------------------

loopx, loopy, jinc = 2**(jlev-1)

number of nodes associated with a result - 1, pointer to the result (in bpl) [read with j changing fastest]

followed by 

iasct = number of asssemblages 

and for each assemblage i

iavar(1,i) = number of solution phases, iavar(2,i) = stoichiometric phases, iavar(3,i) = sum

followed by a pointer for each phase in the assemblage (either to a solution phase or a stoichiometric phase)



========================================================================================

bplot file data write statements

c-----------------------------------
c                                graphics output  
      write (n5,1000) ico,jco,iap(ibulk)
c                                phase molar amounts
      write (n5,1010) (b(i),i=1,np+ncpd)
c                                solution phase compositions
      do i = 1, np
         ids = kkp(i)
         write (n5,1010) ((x3(i,j,k),k=1,ispg(ids,j)),j=1,istg(ids))
      end do 

1000  format (10(i8,1x))
1010  format (6(g16.8,1x))
c-----------------------------------

this i/o is repeasted for every optimization (the numbering of the optimizations is sequential, i.e., the first data is for optimization 1 = ibulk)

ico, jco = the i,j coordinate at which the optimization was done, iap(ibulk) is the pointer to the phase assemblage (ibulk = 1..iasct)

b(i) i = 1...iavar(iavar(3,iap(ibulk))

then 
   for each of the solution phases in the assemblage (iavar(1,iap(ibulk)) read the following data

     for solution ids = idasls(j,iap(ibulk)) 
        for site j = 1...istg(ids)
           read the compositions of each of the k = 1..ispg(ids,j) species