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