In computational mode 2, for choice #38, WERAMI outputs a file containing the  
N requested properties for the system of interest and its stable phases at 
the X-Y nodes of a two-dimensional grid. All data for a given node is written 
on a single line, thus the number of columns varies depending on the stable 
phases.

This option may not function for all 1-d computational modes 3 and 4. If you 
need the option in the 1-d modes please contact me and i'll check the code. 

The structure of this data file is illustrated by the following example in which 
density (kg/m3) and volume proportions (%) of a harzburgitic rock and its stable
phases have been calculated:

LINE   1: charz_pl                                                                                       
LINE   2:    10    10    2.40000        500.003        155555.    
LINE   3:    388.888    
LINE   4: P(bar)  T(K)    
LINE   5:  2
LINE   6:  4  3298.48      100.000     O(stx8)       3315.04      17.4623     Opx(stx8)     3315.04      17.4623     Cpx(stx8)     3315.04      17.4623     an            3315.04      17.4623    
LINE   7:  3  3893.99      100.000     Wad(stx8)     3866.66      3.38450     Gt(stx8)      3866.66      3.38450     stv           3866.66      3.38450    

...

LINE 105:  3  5515.36      100.000     Wus(stx7)     5207.98      75.2266     Ppv(stx8)     5207.98      75.2266     ca-pv         5207.98      75.2266   

The first five lines are the file header:

LINE 1 title
LINE 2 number of grid points in x (XPTS) & y (YPTS) direction; 
       minimum value of x variable (XMIN); 
       minimum value of y variable (YMIN); 
       increment of x variable (DX).
LINE 3 increment of y variable (DY).
LINE 4 x variable name, y variable name (8 characters each)
LINE 5 number of computed properties (N)

The next XPTS*YPTS lines correspond to the data at each node of the grid (one line 
per node). The first value (i2,1x) is a counter P of the number of phases stable at 
the nodal conditions. P may be zero if the optimization by VERTEX failed in the 
vicinity of the nodal coordinate.  For a grid on which N properties have been 
requested there are (N+1)*P - 1 additional entities (a12,1x). The first N entities 
are the bulk system properties; then for each phase N+1 entities are output, 
these are the phase name and the N phase properties. The properties are output 
in the order requested in WERAMI.


The nodal data is output with the row-wise, i.e., with the row index I
increasing from 1 to XPTS for each value of the column index J, which increases 
from 1 to YPTS. The X-Y coordinates at any node are X = XMIN + (I-1)*DX, 
Y = YMIN + (J-1)*DY. A fortran 77 program to read the data is (symbolic array 
indices and the code is untested):

      character file*100, title*100, phase(NP,NX,NY)*12, xname*8, yname*8,
     *          chars(N*NP)*12

      integer i, j, k, l, m, p, xpts, ypts, nprop, nphas(NX,NY)

      double precision xmin, ymin, dx, dy, sdata(N,NX,NY), pdata(N,NP,NX,NY)

      write (*,*) 'enter file name'
      read (*,*) file
      open (10,file)

      read (10,*) title
      read (10,*) xpts,ypts,xmin,ymin,dx,dy
      read (10,*) xname, yname
      read (10,*) nprop

      do i = 1, xpts
         do j = 1, ypts

            read (10,1000) p,(chars(k),k=1,p*(nprop+1)-1)

            do k = 1, nprop
               read (chars(k),*) sdata(k,i,j)
            end do
 
            m = nprop 

            do l = 1, p

               m = m + 1
               phase(l,i,j) = chars(m)

               do k = 1, nprop
                  m = m + 1
                  read (chars(m),*) pdata(l,k,i,j)
               end do 

            end do 

         end do 
      end do 

1000  format (i2,1x,200(a12,1x))

      end 

If the data are read into an arrays indexed by the nodal indices, the indices of 
the closest node to a given X,Y value can be computed as illustrated in the following 
fortran 77 subroutine. The logical variables DOWN and LEFT indicate the position of the 
X,Y value relative to node and can be used to determine the nodes to be used for an 
interpolation scheme. 

      subroutine xy2ij (X,Y,XMIN,YMIN,DX,DY,I,J,LEFT,DOWN)
c----------------------------------------------------------------------
c xy2ij - identifies the grid point associated with coordinate x-y 
c         and the sector of the x-y coordinate relative the nodal point        
c----------------------------------------------------------------------
      implicit none

      integer I, J

      logical LEFT, DOWN

      double precision RES, X, Y, XMIN, YMIN, DX, DY
c----------------------------------------------------------------------
      RES = (X-XMIN)/DX + 1d0
      I = int(RES)

      if (RES-dfloat(I).gt.0.5d0) then
         I = I + 1
         LEFT = .true.
      else 
         LEFT = .false.
      end if 

      RES = (Y-YMIN)/DY + 1d0
      J = int(RES) 

      if (RES-dfloat(J).gt.0.5d0) then
         J = J + 1
         DOWN = .true.
      else 
         DOWN = .false.
      end if 

      end