* program MINTAB-XPT *---------------------------------------------------------------------------- * Version 1.0 (beta) 1/99 * * For compilation with Fortran90 on PC or UNIX platform * (modified for not using EOF function of Fortran95) * * This program reads plot files from CHIM-XPT, GASWORKS, and SOLVEQ and * generates an ASCII file with data in tabular (colums) format for input * into commercial plotting packages. * * N.Spycher, 12/98 * Developed from program MINPLOT originally written by: * Nicolas Spycher U. of O. 1982-1983 * Anthony Michaels, Mark Reed, and Jeff DeBraal U. of O. 1989-1990 * * The type of data to output in tabular format (i.e. temperature, pressure * etc.) can be specified interactively or (faster) by a list of command * strings entered in an ASCII file (MINTAB.PAR). This file can created * with a text editor or by running MINTAB. The default name for this file * is MINTAB.PAR but any other file name can be specified. * ASCII character strings in MINTAB.PAR must be entered one per record, without * heading blanks. These strings identify the parameters to tabulate. String * formats are described in subroutine ReadParam. * * Labels (column headings) of general parameters (i.e. temperature, pressure, * enthalpy, porosity etc.) are assigned to array parlab0 in the BLOCK DATA * program unit. Depending on which version of CHIM-XPT, GASWORKS, or SOLVEQ * you are using, these can be easily reassigned in the BLOCK DATA. * * Maximum allowed array dimensions are defined, and can be modified, * in file MINTABDFx.INC. Most common blocks are set in file MINTABCMx.INC. * * Note: with MINTAB, the number of datasets in the CHIM-XPT, GASWORKS, or * SOLVEQ files is unlimited (no maximum; limited only by the amount of * disk storage). * * 2006-2008 modified to read output from SOLVEQ-XPT and then CHIM-XPT. J.Palandri * 2009 Oct modified for new SOLVEQ-XPT with gas fugacities and their sum grouped * above mineral log Q/K's * *-------------------------------------------------------------------------------- * * Definitions of variables: * * Array maximum dimensions * max1 = max. number of general parameters * max2 = max. number of species (component + derived) * max3 = max. number of gases * max4 = max. number of component species * max5 = max. number of variables (data columns) to output * max6 = unused * max7 = max. number of minerals in simulation * * Number of components in simulation * ibstot component species * nst all aqueous species * ngas gases * nmin minerals * * Type of input plot files * ipar_typ =1 CHIM-XPT * =2 GASWORKS * =3 SOLVEQ-XPT * * General parameters (such as those on the first record of data sets * in the CHIM-XPT or GASWORKS plot file). The index of these parameters * are in the same order as for array PARLAB0(i,ipar_typ) in the BLOCK * DATA program unit: * par(1) (e.g. temperature) * par(2) (e.g. Pfluid) * par(3) ...and so on... * par(4) * par(5) * par(6) * par(7) * par(8) * * Other variables read in the input plot files * rmtry = molalities of individual species * rmtot = total moles of component species * aqm = aqueous moles of component species * (total molality of component species for Solveq) * gasnam = the name of the gas * fuga = fugacity of gas * rname = name of saturated mineral or gas * rmol = saturated mineral/gas moles * gra = saturated mineral/gas grams * solf = saturated mineral/gas mole frac. in solid solution * per = saturated mineral/gas weight percent * qklog = log(Q/K) * * Related to input plot files * nset number of datasets read (i.e. # of temperature steps) * ipar_typ flag for CHIM-XPT (1), GASWORKS(2) or SOLVEQ(3) data * * Number of pointers, and pointer values * npar, ipar number and index of par's to save * nmtry,imtry number and index of rmtry's to save * nmtot,imtot number and index of rmtot's to save * naqm,iaqm number and index of aqm's to save * nfuga,ifuga number and index of fuga's to save * nmol,imol number and index of mol's to save * ngra,igra number and index of gra's to save * nsol,isol number and index of solf's to save * nper,iper number and index of per's to save * nqk,iqk number and index of qk's to save * * Labels for tabulated output * parlab labels assigned to selected par values * fugalab labels assigned to selected fuga values * mtrylab labels assigned to selected rmtry values * mtotlab labels assigned to selected rmtot values * aqmlab labels assigned to selected array aqm * mollab labels assigned to selected rmol values * gralab labels assigned to selected gra value * solflab labels assigned to selected solf values * perlab labels assigned to selected per values * qklab labels assigned to selected qk values * * Other labels read in input plot files: * specnam aq. species names * rname mineral names * gasname gas names * * * I/O device numbers (defined in mintabdfx.inc) * * inp1 main plot input file (i.e. chillplt.dat) * inp2 parameter input file (inp2 can be output if created) * inp3 label file (CHIM-XPT/gasworks data only) * iout1 output file (tabulated data) * iunit1 chillout file * ************************ block data ************************ implicit real*4 (a-h,o-z) implicit integer (i-n) character*12 parlab0(8,3) !default labels for general parameters ! (1=CHIM-XPT, 2=solvgas and 3=solveq) character*8 blank, blank8 character*12 blank12 character*20 blank20 character*1 pos data blank8 /' '/ data blank12 /' '/ data blank20 /' '/ data pos /'y'/ * Note: the labels below must correspond, in same order, to the variables * output in plot files by CHIM-XPT, gasworks, and solveq. You can change * these labels to match plot files of various program versions (that's * all you need to do to, then recompile). * Note: if parlab0 for pH is 'pH', input a(H+) is converted automatically * to pH. Otherwise, a(H)+ is output. * parlab0(i,1) i=1,8 describes the data in par(i) for CHIM-XPT plot files * parlab0(i,2) i=1,8 describes the data in par(i) for GASWORKS plot files * parlab0(i,3) i=1,8 describes the data in par(i) for SOLVEQ output files data parlab0/'Temp ','Pfluid ','pH ','Mix.frac', !CHIM-XPT + 'Tempc ','Enthalpy','Porosity','Gas_wt.%', !CHIM-XPT + 'Temp ','Pfluid ','f(O2) ','Mix.frac', !gasworks + ' ',' ',' ',' ', !gasworks + 'Temp ','Pfluid ','Alkalin.','pH ', !solveq + ' ',' ',' ',' '/ !solveq common /parlab/parlab0 common /strings/pos,blank,blank8,blank12,blank20 end ******************************** * main program ******************************** program mintab_xpt include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations character*1 ans character*8 dummy character*20 name character*8 title1(10), title2(10) character*80 fil character*12 fildef(3) logical exists fildef(1)='chillplt.dat' fildef(2)='workplot.dat' fildef(3)='solout.dat' * Version number write(*,"(///5x,'*** MINTAB-XPT version 1.0 ***'///)") * Initialization do i=1,max7 divmin(i)=1.0 !mineral division factor end do * Define the type of input plot file write(*,"(/5x,'Type of data file to read:' & //10x,'1 chim-xpt plot file'/10x,'2 gasworks plot file' & /10x,'3 solveq output file' & //5x,'Enter selection (default=1) :> '$)") read(*,"(i1)") ipar_typ if(ipar_typ.le.1.or.ipar_typ.gt.3) ipar_typ=1 * Open files write(*,"(/5x,'Enter name of plot data input file (default=' & a12,')',/5x,':> '$)") fildef(ipar_typ) fil=fildef(ipar_typ) call open_old(inp1,fil) * Default parameter file name paramfil='mintab.par' * Default output file name filout = fil do k = 1,80 if (fil(k:k) .eq. '.') then filout(k+1:k+3) = 'tab' exit end if end do * Default label file name for chim-xpt, gasworks labelfil = fil do k = 1,80 if (fil(k:k) .eq. '.') then labelfil(k+1:k+3) = 'lab' exit end if end do * Look to see if label file exists (chim-xpt/gasworks files only) if(ipar_typ.lt.3) then inquire (file = labelfil, exist = exists) if(.not.exists) call GenLabels(ipar_typ) end if write(*,"(/5x,'Enter name of tabulated output data file' & /5x,'(default = ',a12,') :> '$)") filout call open_new(iout1,filout) *---------------------------------------------------------------- * Reads the plot file (first time through) * For CHIM-XPT/GASWORKS plot file only if(ipar_typ.lt.3) then * Read in the numbers ibstot, nst, and ngas in the input file * (correspond to the number of components, complexes, * and the number of gases respectively). read(inp1, *) ibstot, nst, ngas write(*, "( //2x, 'Ibstot = ',i3, 2x,'Nst = ',i3, & 2x,'Ngas = ', i3,//)") ibstot, nst, ngas if(ibstot.gt.max4) then write(*,"(/5x,'Number of component species (ibstot) exceeded.' & /,5x,'Maximum allowed is: ',i3)") max4 stop elseif(nst.gt.max2) then write(*,"(/5x,'Number of derived species (nst) exceeded.'/ & ,5x,'Maximum allowed is: ',i3)") max2 stop elseif(ngas.gt.max3) then write(*,"(/5x,'Number of gases (ngas) exceeded.'/ & ,5x,'Maximum allowed is: ',i3)") max3 stop endif read (inp1, "( 10a8/10a8)") title1, title2 write (*, "( //2x, 10a8/2x, 10a8/5x, & 'Reading plot file to check and order data...'/)")title1,title2 end if *--------------------------------------------------------------------- * Now we need to read the entire plot file to get the mineral and gas * names and sort them. We will read the file again later, set by * set, to tabulate the data. This way we are not limited by the * plot file size (i.e. do not need to store all data in memory). * For SOLVEQ data, we also get species names here nsets=0 !index of data set iend=0 !flag for end-of-file when reading top of set do k = 1,max7 !initialize mineral name variables rname(k)=blank20 end do write(*,*) do if(ipar_typ.lt.3) then !CHIM-XPT/GASWORKS plot file call ReadPlotC(iend,ipar_typ) else !SOLVEQ output file call ReadPlotS(iend) end if if(iend.eq.1) exit !end of file when reading top of set nsets=nsets+1 write (*,"(5x,'Done reading data set number : ',i4, 8x, & ' Number of mins : ', i3)") nsets, nmin end do * Rewind the plot file and advance to the top of the first data set rewind (inp1) if(ipar_typ.lt.3) read (inp1,"(//)") !skip first 3 records *----------------------------------------------------------------------- * Now we select the data to tabulate write(*,"(/5x,'Use an existing plot parameter' & 'file (default=yes) ? :> '$)") read(*,"(a1)") ans if(ans.eq.'N'.or.ans.eq.'n') then * Select data to tabulate and generate the parameter file do call Select(ipar_typ) write (*, '( ///5x, ''Do you want to tabulate more '', + ''data ? (default=yes) :> '', $)') read ( * , '( a1)') ans if(ans.eq.' ') ans = pos if ((ans .ne. 'y') .and. (ans .ne. 'Y')) exit end do call GenParam else * Reads the Parameter file (parameters to tabulate) call ReadParam(ipar_typ) end if *----------------------------------------------------------------- * * Writes labels at top of the output file *----------------------------------------------------------------- * General parameters do i=1,npar write(iout1,"(1x,a12,$)") parlab(i) end do write (iout1,'(a12,$)') ' ' * Gas fugacities do i=1,nfuga write(iout1,"(1x,a24,$)") fugalab(i) !fugacity end do * Selected aq. species do i=1, nmtry write (iout1,"(1x,a16,$)") mtrylab(i) !molality end do do i=1, nmtot write (iout1,"(1x,a16,$)") mtotlab(i) !total moles end do do i=1, naqm write (iout1,"(1x,a16,$)") aqmlab(i) !total aq. molality end do * Selected minerals do i=1,nmol write (iout1,"(1x,a24,$)") mollab(i) !moles end do do i=1,ngra write (iout1,"(1x,a24,$)") gralab(i) !grams end do do i=1,nsol write (iout1,"(1x,a24,$)") solflab(i) !solid sol. fraction end do do i=1,nper write (iout1,"(1x,a24,$)") perlab(i) !weight percent end do do i=1,nqk write (iout1,"(1x,a24,$)") qklab(i) !log q/k end do write (iout1,"(' ')") !carriage control (next line) *----------------------------------------------------------------- * * Reads each plot data set and output the data in the output file * --------------------------------------------------------------- * do while (.not.eof(inp1)) !F95 compilers only - use iend instead iend=0 !eof flag do do k = 1,max7 !reset mineral variables for current step rmol(k)=0. per(k)=0. solf(k)=0. gra(k)=0. end do if(ipar_typ.lt.3) then call ReadPlotC(iend,ipar_typ) !read one CHIM-XPT/GASWORKS data set else call ReadPlotS(iend) !read one SOLVEQ data set end if if(iend.eq.1) exit !end of file *Output selected general parameters do i=1,npar write(iout1,"(1x,g12.5,$)") par(ipar(i)) end do *Output selected gas fugacities do i=1,nfuga write(iout1,"(15x,e10.4,$)") fuga(ifuga(i)) !fugacity end do *Output selected data for aq. species do i=1, nmtry write (iout1,"(7x,e10.4,$)") rmtry(imtry(i)) !molality end do do i=1, nmtot write (iout1,"(7x,e10.4,$)") rmtot(imtot(i)) !total moles end do do i=1, naqm write (iout1,"(7x,e10.4,$)") aqm(iaqm(i)) !total aq. molality end do *Output selected data for minerals do i=1,nmol write (iout1,"(15x,e10.4,$)") rmol(imol(i))/divmin(imol(i)) !mineral moles end do do i=1,ngra write (iout1,"(15x,e10.4,$)") gra(igra(i)) !mineral grams end do do i=1,nsol write (iout1,"(15x,e10.4,$)") solf(isol(i)) !mineral solid sol. fraction end do do i=1,nper write (iout1,"(15x,e10.4,$)") per(iper(i)) !mineral weight percent end do do i=1,nqk write (iout1,"(15x,f10.4,$)") qklog(iqk(i)) !mineral log q/k end do write (iout1,"(1x)") !carriage control (next line) end do *--- Finished output------------------------------ *Get here when eof reached close (unit = inp1) close (unit = iout1) write(*, "( //2x, ' .....Done!'/)") stop end subroutine Select(ipar_typ) *********************************************************** * Allows user to select the plot data to be tabulate. * ipar_typ: CHIM-XPT = 1, GASWORKS = 2, SOLVEQ = 3 include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations integer iplt(max5) character*1 ans character*8 dummy character*20 name character*150 lineinp data icall/0/ save icall !flag: nb of time routine is called icall=icall+1 if(icall.gt.1) then write(*,"(/5x,'Re-enter general parameters ? (default=no) :> ' + ,$)") read(*,"(a1)") ans if(ans.eq.'Y') ans='y' end if if(icall.eq.1.or.ans.eq.'y') then !general parameters * Queries the user to select general parameters to tabulate (one or more) if(icall.eq.1) then !first time do i=1, max5 ipar(i)=0 parlab(i)=' ' end do end if write (*,"(//5x,'General parameters to tabulate:'/)") do i=1,4 write (*, "(12x,i1,' = ',A8,10x,i1,' = ',A8)") + i,parlab0(i,ipar_typ),i+4, parlab0(i+4,ipar_typ) end do write (*, "(/5x,'Select a list of indexes:> '$)") read(*,"(a150)") lineinp read (lineinp,*,end=111) (ipar(i), i=1,8) 111 continue end if * General parameters, assign labels do i = 1, max5 parlab(i)=parlab0(ipar(i),ipar_typ) if(ipar(i).eq.0) then npar = i-1 exit end if end do * Queries the user for other data to tabulate (one at a time) write (*, '( //5x, ''Other parameters:'')') if (ipar_typ.eq.1) then write (*, '( /12x, ''1 = Individual species molalities'', 1 /12x, ''2 = Total moles of components'', 2 /12x, ''3 = Total aqueous moles'', 3 '' (total molalities of component species)'', 4 /12x, ''4 = Log Fugacities of gases'', 5 /12x, ''5 = Minerals''//)') elseif (ipar_typ.eq.2) then write (*,'( /12x,''1 = Individual species mole numbers'', 1 /12x, ''2 = Total moles of components'', 2 /12x, ''3 = Total gas moles'', 4 /12x, ''4 = '', 5 /12x, ''5 = Minerals''//)') elseif (ipar_typ.eq.3) then write (*, '( /12x, ''1 = Individual species molalities'', 1 /12x, ''2 = Total moles of components'', 2 /12x, ''3 = Total aqueous moles'', 3 '' (total molalities of component species)'', 4 /12x,''4 = Log(Q/K) of Minerals or Gases''//)') endif write (*, '(5x,''Enter one selection (return for none):> ''$)') read (*,"(a150)") lineinp read (lineinp ,*) indy if(indy.eq.0) return !done selecting * If the user had chosen to plot individual species molalities * then call the subprogram to get (or generate) the label file. * From this list of species, the user will be able to choose * the required species. if (indy .eq. 1) then !molalities write(*, '( //)') key = -1 call GetLabels(key,imtry,nmtry,mtrylab,ipar_typ) *c do k=1,nmtry *c l=len_trim(mtotlab(k)) *c mtotlab(k)(l+1:l+4)='_mly' *c end do *cjp do k=1,nmtry *cjp l=len_trim(mtrylab(k)) *cjp mtrylab(k)(l+1:l+4)='_mly' *cjp end do *------------------------------------------------------- elseif (indy .eq. 2) then !total moles key = ibstot call GetLabels(key,imtot,nmtot,mtotlab,ipar_typ) do k=1,nmtot l=len_trim(mtotlab(k)) mtotlab(k)(l+1:l+4)='_tot' end do elseif (indy .eq. 3) then !total aqueous molalities key = ibstot call GetLabels(key,iaqm,naqm,aqmlab,ipar_typ) do k=1,naqm l=len_trim(aqmlab(k)) aqmlab(k)(l+1:l+4)='_aqm' end do end if * CHIM-XPT/GASWORKS data - minerals and gases if(ipar_typ.lt.3) then if (indy .eq. 4) then do i=1,max5 fugalab(i)=' ' ifuga(i)=0 end do do k = 1, ngas write (*, '( 5x, i2, 5x, a8)') k, gasnam(k) end do write (*, "( //5x,'Enter list of gas indexes:> '$)") read(*,"(a150)") lineinp read (lineinp,*,end=222) (ifuga(k), k = 1, max5) 222 continue do k = 1, max5 if(ifuga(k).eq.0) then nfuga=k-1 exit else fugalab(k)= gasnam(ifuga(k)) endif end do elseif (indy .eq. 5) then do i=1,max5 iplt(i)=0 end do write (*, '( //5x, ''Available mineral data: ''/, 1 /12x, ''1 = Moles'', /12x, ''2 = Grams'', 2 /12x, ''3 = %'', /12x ''4 = Mole Fraction''/)') write (*,"( 5x,'Enter one selection (return if none):> '$)") read(*,"(a150)") lineinp read (lineinp,*) inmy minf = 1 *jp do n = 1, 13 *jp n1 = n + 13 *jp n2 = n + 26 *jp write (*, '( 3( 2x, i3, 2x, a8))') *jp 1 n, rname(n), n1, rname(n1), n2, rname(n2) *jp end do do n = 1, 20 n1 = n + 20 n2 = n + 40 write (*, '( 4( 2x, i3, 2x, a20))') 1 n, rname(n), n1, rname(n1), n2, rname(n2) end do write (*,"( /5x,'Enter list of mineral indexes:> '$)") read(*,"(a150)") lineinp read (lineinp,*,end=333) (iplt(i), i = 1, max5) 333 continue if (inmy .eq. 1) then !moles minerals write (*, "( //5x, 'Which minerals do you want to', 1 ' divide, if any?'/, 5x,'Input mineral', 2 ' index and dividing factor (e.g. 6,5.0).', 3 /5X,'Press [enter] to skip this step.')") do read(*,"(a150)") lineinp if(lineinp(1:1).eq.'') exit read (lineinp,*,end=444) k, divmin(k) 444 continue * if (k .eq. 0) exit write (*, '( /5x, ''Mineral : '', a8, 1 '' divided by '', f8.3)') rname(k), divmin(k) end do do i=1,max5 mollab(i)=rname(iplt(i)) l=len_trim(mollab(i)) mollab(i)(l+1:l+4)='_mol' imol(i)=iplt(i) if(iplt(i).eq.0) then nmol=i-1 exit end if end do elseif (inmy .eq. 2) then !grams minerals do i=1,max5 gralab(i)=rname(iplt(i)) l=len_trim(gralab(i)) gralab(i)(l+1:l+4)='_gra' igra(i)=iplt(i) if(iplt(i).eq.0) then ngra=i-1 exit end if end do elseif (inmy .eq. 3) then !percent minerals do i=1,max5 perlab(i)=rname(iplt(i)) l=len_trim(perlab(i)) perlab(i)(l+1:l+4)='_per' iper(i)=iplt(i) if(iplt(i).eq.0) then nper=i-1 exit end if end do elseif (inmy .eq. 4) then !mole fractions do i=1,max5 solflab(i)=rname(iplt(i)) l=len_trim(solflab(i)) solflab(i)(l+1:l+4)='_frc' isol(i)=iplt(i) if(iplt(i).eq.0) then nsol=i-1 exit end if end do end if end if * SOLVEQ data - minerals log Q/K else if (indy .eq. 4) then *ptx nloop=nmin/65 *ptx if(mod(nmin,65).ne.0) nloop=nloop+1 *ptx istart=loop*65-64 *ptx do n = istart, istart+12 *ptx n1 = n + 13 *ptx n2 = n + 26 *ptx n3 = n + 39 *ptx n4 = n + 52 nloop=nmin/48 if(mod(nmin,48).ne.0) nloop=nloop+1 nqk=0 do loop=1,nloop istart=loop*48-47 do n = istart, istart+15 n1 = n + 16 n2 = n + 32 write (*, '( 5( 2x, i3, 2x, a20))') 1 n, rname(n), n1, rname(n1), n2, rname(n2) end do write (*,"( /5x, 'Enter list of mineral indexes:> '$)") read(*,"(a150)") lineinp read (lineinp,*,end=555) (iplt(i), i = nqk+1, max5) 555 continue do i=nqk+1,max5 if(iplt(i).eq.0) then nqk=i-1 exit end if qklab(i)=rname(iplt(i)) l=len_trim(qklab(i)) qklab(i)(l+1:l+4)='_lqk' iqk(i)=iplt(i) end do end do end if end if return end subroutine open_old (iunit,default) *********************************** * Opens an old file (By N.S. 5/98) implicit integer (i-n) implicit real*4 (a-h,o-z) logical exists character*80 filenam,default ifirst=1 do if(ifirst.eq.0) write(*,"(/5x,'Enter another file name' + ' (or q to quit):> '$)") ifirst=0 read(*,"(a40)") filenam if(filenam.ne.' ') default=filenam if((filenam(1:1).eq.'Q'.or.filenam(1:1).eq.'q').and. + filenam(2:2).eq.' ') stop inquire (file = default, exist = exists) if(exists) then open (unit=iunit,file=default,status='old',err=100) return endif write(*,"(5x,'Warning! File does not exist!'/)") end do 100 write(*,"(/'Error while opening the file - stop'/)") stop end subroutine open_new (iunit,default) *********************************** * Opens a new file (By N.S. 5/98) implicit integer (i-n) implicit real*4 (a-h,o-z) logical exists character*80 filenam,default character*1 ans ifirst=1 10 if(ifirst.eq.0) write(*,"(/5x,'Enter another file name:> '$)") ifirst=0 read(*,"(a40)") filenam if(filenam.ne.' ') default=filenam inquire (file = default, exist = exists) if(exists) then write(*,"(/5x,'File already exists. Replace ? (y/n) :> '$)") read(*,"(a1)") ans if(ans.ne.'Y'.and.ans.ne.'y') goto 10 end if open (unit = iunit, file = default, status = 'unknown',err=100) return 100 write(*,"(/'Error while opening the file - stop'/)") stop end subroutine GetLabels(key,iplt,index,label,ipar_typ) ********************************************************** * This subprogram displays the species on the screen. * The user can select any of these species for graphing * by selecting the corresponding keys. The selected * list of species are then returned to min in order to * extract the selected data corresponding to each species. * GetLabels is called by the main program MINPLOT * iplt array of selected species indeces * nplt total number of species selected include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations character*16 label(max5),dummy character*20 fmt1, fmt2 character*150 lineinp integer iplt(max5) data fmt1(1:14) /'(1x,i3,1x,a12)'/ data fmt2(1:17) /'(6(1x,i3,1x,a12))'/ icount=nst *-------------------------------------------------------------------- * CHIM-XPT/GASWORKS plot file only: * Open the label file and read the species names if(ipar_typ.lt.3) then open (unit = inp3, file = labelfil, status = 'old') icount = 1 do read (inp3, fmt1,end=1000) num, specnam(icount) icount = icount + 1 end do 1000 close (unit = inp3) end if * SOLVEQ file: specnam read earlier by call to ReadPlotS *-------------------------------------------------------------------- * Initialize the selection array to zero do k = 1, max5 iplt(k) = 0 label(k) = ' ' end do *-------------------------------------------------------------------- * If the last entry key read is not equal -1 then adjust * the count of species by adding 2. if (key .ne. -1) then icount = key + 2 endif *-------------------------------------------------------------------- * This part of the subroutine displays the species on * the screen and reads in the selections from the user. If (key .ne. -1) then icount = icount - 2 else icount = icount - 1 endif irows = icount / 6 inx = 1 istart = 1 do k = 1, irows if (mod(k, 21) .eq. 0) then write (*,"(//5x,'Enter list of species indexes:> '$)") read(*,"(a150)") lineinp read (lineinp,*,end=111) (iplt(in), in = istart, max5) 111 continue index = istart do if (iplt(index) .eq. 0) exit index = index + 1 end do istart = index endif write (*, fmt2) (j, specnam(j), j = inx, inx + 5) inx = inx + 6 end do icols = icount - irows * 6 if (icols .gt. 0) then write (*, fmt2) (j,specnam(j), j = inx, inx + icols - 1) endif if (mod(k, 21) .ne. 0) then write (*,"( //5x,'Enter list of species indexes:> '$)") read(*,"(a150)") lineinp read (lineinp,*,end=222) (iplt(in), in = istart, max5) 222 continue index = istart do if (iplt(index) .eq. 0) exit index = index + 1 end do index=index-1 !total number of indexes selected endif *-------------------------------------------------------------------- do j = 1, index label(j) = specnam(iplt(j)) end do return end subroutine GenLabels(ipar_typ) ************************************************************* * Subroutine GenLabels is used to generate the label * file xxxxxPLT.LAB which stores the species' name and an * index corresponding to each species. It is used in * subroutine GetLabels to display a list of species * for the user to select from. * This subroutine is called from the main program MINPLOT * and is called whenever xxxxxPLT.LAB doesn't exist. The * label file is also used when graphing the legends. * Used for CHIM-XPT and Gasworks include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations integer indx(max2) character*8 intxt character*12 species, inspec(500) character*12 fil logical exists *-------------------------------------------------------------------- * Querry the user which output file was generated by CHIM-XPT * to generate the xxxxxx.sum file. if (ipar_typ .eq. 1 ) then write (*, 30) 30 format (//4x, ' There is no label file for this data set ' 1 /4x, ' Enter output filename generated by CHIM-XPT', 2 /4x, ' (default = chillout.dat) :> ',$) fil='chillout.dat' call open_old(iunit1,fil) else if ( ipar_typ .eq. 2) then write (*, 31) 31 format ( 4x, ' There is no label file for this data set ' 1 /4x, ' Enter output filename generated by GASWORKS', 2 /4x, ' (default = workout.dat) :> ',$) fil='workout.dat' call open_old(iunit1,fil) endif *-------------------------------------------------------------------- * Open and read this output file and read the portion of * the file where the list of species is located. Create * and write this information in the label file. write(*,"(/5x,'Enter the name of the label file to create' & /5x,' (default = ',a12,'):> '$)") labelfil call open_new(inp3,labelfil) 100 continue read (iunit1, 40) intxt 40 format ( 2x, a8) if (intxt.ne.'Species '.and.intxt.ne.'SPECIES ') go to 100 read (iunit1, 40) intxt 200 continue If (ipar_typ.eq.2) then read (iunit1, 50) species, index 50 format ( 2x, a8, 2x, i3) else read (iunit1, 52) species, index 52 format ( 2x, a12, i4) endif if (species .ne. blank12) then write (inp3, 60) index, species 60 format ( i4, 1x, a12) go to 200 endif 300 continue read (iunit1, 70) intxt 70 format ( 4x, a8) if (intxt.ne.'Species '.and.intxt.ne.'SPECIES ') go to 300 read (iunit1, 70) intxt k = 0 400 k = k + 1 If (ipar_typ.eq.2) then read (iunit1, 888) species, index, inspec(k), indx(k) 888 format ( 4x, a8, i7, 24x, a8, i7) else read (iunit1, 80) species, index, inspec(k), indx(k) 80 format ( 4x, a12, i6, 23x, a12, i6) endif if (species .ne. blank12) then write (inp3, 60) index, species go to 400 endif k = k - 1 if (inspec(k) .eq. blank12) then k = k - 1 endif do 500 j = 1, k write (inp3, 60) indx(j), inspec(j) 500 continue *-------------------------------------------------------------------- * At the very end of our species list add pH with index 0 * if using CHIM-XPT data; add f(O2) if using Gasworks data. * if ( gtype .eq. 'c' ) then * write (inp3, 90) *90 format ( ' 0 pH') * else if ( gtype .eq. 'g' ) then * write (inp3, 91) *91 format ( ' 0 f(O2)' ) * endif *-------------------------------------------------------------------- close (unit = inp3) close (unit = iunit1) return end subroutine ReadPlotC(iend,ipar_typ) ******************************************************************************** * This routines reads one set of the CHIM-XPT or Gasworks plot file * (without first line and title lines, just the plot data * with blank line at the end). include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations character*20 name, dummy * Defaults for mineral parameters for this set * (in case some minerals not present at all steps) do n=1,nmin_save rmol(n) = 0. gra(n) = 0. solf(n) = 1.0 per(n) = 0. end do * Read in first record of parameters read (inp1, '( 8e10.4)',end=1000) (par(i), i=1,max1) * write (*, '( 8(1x,g10.5))') (par(i), i=1,max1) * pause * The preceding statement reads the chillplt file general parameters * assigned as defined for parlab0(i,ipar_typ) in the BLOCK DATA * (for example,temp, pfluid, act(H+) etc.. * Converts aH+ to pH (needs a 'pH' label to work) do i=1,max1 if(parlab0(i,ipar_typ)(1:2).eq.'pH') then par(i)=-alog10(par(i)) exit end if end do * Read in data for aq. species read (inp1, '(8e10.4)') (rmtry(j), j = 1, nst) read (inp1, '(8e10.4)') (rmtot(j), j = 1, ibstot) read (inp1, '(8e10.4)') (aqm(j), j = 1, ibstot) * Read in the gas fugacities (logarithm) do n = 1,ngas If (ipar_typ.eq.2) then read (inp1,'( i3,2x,a8,2x,f10.4)') j,gasnam(j),fuga(j) else read (inp1,'( i3,2x,a20,2x,f10.4)') j,gasnam(j),fuga(j) endif end do * Read in the mineral data nmin = 0 do n = 0 If (ipar_typ.eq.2) then read (inp1, '( i3,2x,a8,2x,6e10.4)') & mm, name, rmil, rgr, rpr, rsol else read (inp1, '( i3,2x,a20,2x,6e10.4)') & mm, name, rmil, rgr, rpr, rsol endif if (mm .ne. 0) then do n = n + 1 if(n.ge.nmin) nmin=n if((name.eq.rname(n)).or.(rname(n).eq.blank20))exit end do if(nmin.gt.max7) then write(*,"(/5x,'Number of minerals (max7) exceeded.'/ & ,5x,'Maximum allowed is: ',i3)") max7 stop endif rname(n) = name rmol(n) = rmil gra(n) = rgr solf(n) = rsol per(n) = rpr else exit endif end do * Read in 2nd blank line at end of mineral list (marks end of dataset) read (inp1, '( a8)') name if (name .ne. blank20) then write (*,"( //10x,'Error in reading CHIM-XPT data file //')") close (unit = inp1) stop endif return 1000 iend=1 !end-of-file at top of set return end subroutine GenParam *********************************** * This routine generates a table parameter file (.i.e. an ASCII file * with commands to automatically tabulate selected data. * include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations write(*,"(/5x,'Enter name of parameter file to create' & /5x,'(default = ',a12,') :> '$)") paramfil call open_new(inp2,paramfil) !inp2 is output in this case * write(inp2,"('ipar_typ ',i1)") ipar_typ * General parameters do i=1,npar write(inp2,"(a12)") parlab(i) end do * Gas fugacities do i=1,nfuga write(inp2,"(a24)") fugalab(i) !fugacity end do * Selected aq. species do i=1, nmtry write (inp2,"(a16)") mtrylab(i) !molality end do do i=1, nmtot write (inp2,"(a16)") mtotlab(i) !total moles end do do i=1, naqm write (inp2,"(a16)") aqmlab(i) !total aq. molality end do * Selected minerals do i=1,nmol write (inp2,"(a24,f5.0)") mollab(i), divmin(imol(i)) !moles end do do i=1,ngra write (inp2,"(a24)") gralab(i) !grams end do do i=1,nsol write (inp2,"(a24)") solflab(i) !solid sol. fraction end do do i=1,nper write (inp2,"(a24)") perlab(i) !weight percent end do do i=1,nqk write (inp2,"(a24)") qklab(i) !log q/k end do write(inp2,"('endlist')") !end record return end subroutine ReadParam(ipar_typ) ************************************* * * This routine reads the parameter file (.i.e. reads the file * containing info on the type of data to tabulate) * This file is run-independent and contains character strings * that tell mintab what data to tabulate as follows (in any order * except for ipar_typ * * Plot type (must be first entry or defaults to 1!) * ipar_typ type of plot file to read (1=CHIM-XPT,2=gasworks, and * 3=solveq) * * * General parameters. These must be identical to the labels * defined for array parlab0, currently: * Temp temperature * Pfluid fluid pressure * a(H+) hydrogen ion activity * Mix.frac mixing fraction * Tempc temperature of mixing solution * Enthalpy total enthalpy * Porosity porosity * Gas_wt.% gas weight percent of water+gas system * * Aq. species, gases and mineral to select must have the * exact same name as in the soltherm data file, and must be * followed directly with the following strings: * _aqm total aq. molality e.g. NaCl_aqm * _tot total moles e.g. SO4--_tot * (blank) molality of individual species e.g. Cl- * _gra gram mineral e.g quartz_gra * _mol moles mineral, followed by division factor (optional) * (div. factor 1 by default, separated by at least one blank) * _per wt. percent minerals * _frc mole frcation mineral solid solution * * Final record (mandatory): * endlist (or blank, or any record) * include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations character*24 string character*80 record1 integer ltrim * Oens the parameter file write(*,"(/5x,'Enter name of parameter file to use' & /5x,'(default = ',a12,') :> '$)") paramfil call open_old(inp2,paramfil) * CHIM-XPT/GASWORKS files only if(ipar_typ.lt.3) then * Opens and reads the label file (aqueous species names) * (this file is run-dependent and created from main) open (unit = inp3, file = labelfil, status = 'old') icount = 1 do read (inp3,"(1x,i3,1x,a12)",end=1000) & indx, specnam(icount) icount = icount + 1 end do 1000 close (unit = inp3) end if * SOLVEQ files: specnam was read by call to ReadPlotS earlier * Initialize all table data index parameters npar=0 !number of general parameters to select nfuga=0 !number of gas species fugacity values to select nmtry=0 !number of species molality values to select nmtot=0 !number of components total moles values to select naqm=0 !number of species aq. molality values to select nmol=0 !number of mineral moles values to select ngra=0 !number of mineral grams values nsol=0 !number of mineral mole fraction values to select nper=0 !number of mineral weight percent to select nqk=0 !number of mineral log q/k * ipar_typ=1 !default * At this stage, nothing is selected * Reads the parameter file and assigns the data to be tabulated * according to the content of the parameter file do !main loop (reads each record) iloop=0 !flag =1 to loop back to read next entry read(inp2,"(a80)",end=2000) record1 string = record1(1:24) * checks plot type first, if none defaults to 1 * if (string(1:8).eq.'ipar_typ') then * read(record1(10:80),*) ipar_typ * if(ipar_typ.lt.1.or.ipar_typ.gt.3) ipar_typ=1 * iloop=1 * end if * if(iloop.eq.1) cycle * Checks for general parameters l=len_trim(string) !length without trailing blanks do i=1,max1 if(string(1:l).eq.parlab0(i,ipar_typ)(1:l)) then npar=npar+1 ipar(npar)=i parlab(npar)=string iloop=1 exit end if end do if(iloop.eq.1) cycle * Checks for aqueous species do i = 1, nst * l=len_trim(specnam(i)) ltrim=len_trim(specnam(i)) l=len_trim(string) if(l.gt.4.and.string(l-3:l-3).eq.'_') l=l-4 if(string(1:l).eq.specnam(i)(1:ltrim)) then * Found species match, check for type of data select case (string(l+1:l+4)) case('_aqm') !total aqueous molalities naqm=naqm+1 iaqm(naqm)=i aqmlab(naqm)=string iloop=1 exit case('_tot') !total molalities nmtot=nmtot+1 imtot(nmtot)=i mtotlab(nmtot)=string iloop=1 exit case default !molalities of actual species nmtry=nmtry+1 imtry(nmtry)=i mtrylab(nmtry)=string iloop=1 exit end select end if end do if(iloop.eq.1) cycle * Checks for gases do i = 1, ngas l=len_trim(gasnam(i)) if(string(1:l).eq.gasnam(i)(1:l) & .and.string(l+1:l+1).ne.'_') then * Found gas match, check for type of data nfuga=nfuga+1 ifuga(nfuga)=i fugalab(nfuga)=string iloop=1 exit end if end do if(iloop.eq.1) cycle * Checks for minerals *ptx do i = 1, nmin nmin_save = nmin do i = 1, nmin_save l=len_trim(rname(i)) if(string(1:l).eq.rname(i)(1:l)) then * Found mineral match, check for type of data select case (string(l+1:l+4)) case('_mol') !moles of mineral nmol=nmol+1 imol(nmol)=i mollab(nmol)=string iloop=1 read(record1(l+6:80),*,end=100) divmin(i) 100 if(divmin(i).eq.0.) divmin(i)=1.0 exit case('_gra') !grams of mineral ngra=ngra+1 igra(ngra)=i gralab(ngra)=string iloop=1 exit case('_per') !wt. percent of mineral nper=nper+1 iper(nper)=i perlab(nper)=string iloop=1 exit case('_frc') !wt. percent of mineral nsol=nsol+1 isol(nsol)=i solflab(nsol)=string iloop=1 exit case('_lqk') !wt. percent of mineral nqk=nqk+1 iqk(nqk)=i qklab(nqk)=string iloop=1 exit end select end if end do if(iloop.eq.1) cycle end do 2000 close(inp2) return end subroutine ReadPlotS(iend) *********************************************************** * This routines reads one set of the SOLVEQ output file * The top of each set is assumed to start with the * string ' Temperature:' include 'mintabdfx.inc' !array dimensions, io numbers etc include 'mintabcmx.inc' !common blocks, declarations character*20 name, dummy character*80 record1, dum80 character*25 string * Defaults for mineral parameters for this set * (in case some minerals not present at all steps) do n=1,nmin qklog(n)=-500. end do * Looks for the top of the data set (temperature record) do string='Temperature field' read (inp1, "(a80)",end=1000) record1 if(record1(1:14).eq.' Temperature:') exit end do * Reads temperature read(record1(15:23),*) par(1) read(record1(47:54),*) par(2) * Looks for water do string='Kg of water' read (inp1, "(a80)",end=1001) record1 if(record1(1:8).eq.' Water:') exit end do * Reads kg H2O read(record1(9:19),*) rmtry(2) * Looks for alkalinity do string='Alkalinity field' read (inp1, "(a80)",end=1001) record1 if(record1(49:58).eq.'Alkalinity') exit end do * Reads alkalinity read(record1(61:73),*) par(3) * Looks for a(H+) do string='First H+ field' read (inp1, "(a80)",end=1001) record1 if(record1(1:8).eq.' 1 H+ ') exit end do * Reads a(H+) read(record1(43:54),*) par(4) par(4)=-alog10(par(4)) !get pH * Reads mtrys and species names read(record1(1:30),"(i4,1x,a12,1x,e12.4)") + indx,specnam(1),rmtry(1) read(inp1,"(i4,1x,a12)") indx, specnam(2) do i=3,max2 read(inp1,"(i4,1x,a12,1x,e12.4)") indx,specnam(i),rmtry(i) if(indx.eq.0) exit end do nst=i-1 * Looks for list of comp. species, starting with H+ do string='Second H+ field' read (inp1, "(a80)",end=1001) record1 if(record1(1:8).eq.' 1 H+ ') exit end do * Reads a(H+) read(record1(1:72),"(i4,1x,a12,14X,e12.4,17x,e12.4)") + indx,dummy,rmtot(1),aqm(1) * Reads total molality and moles do i=2,max4 read(inp1,"(i4,1x,a12,14X,e12.4,17x,e12.4)") + indx,dummy,rmtot(i),aqm(i) if(indx.eq.0) exit end do ibstot=i-1 * Reads the mineral data * Looks for +++++++ line do string='++++++ after species data' read (inp1, "(a80)",end=1001) record1 if(record1(1:12).eq.' ++++++++++') exit end do * Now we are either at the top of the mineral list or at the top of the * next set if there are no mineral list. Check for mineral header: read(inp1,"(///a80)") record1 if(record1(1:9).ne.' Gas ') return read(inp1,"(a1)") dummy !skips blank line nmin = 0 !initialize number of minerals do k=1,3 !three passes for each mineral set (gases+supersat+undersat) do !loop reading through minerals n = 0 string='Gas names' read (inp1,"(a80)",end=1002) record1 if (record1(1:8).eq.blank8) exit if (record1(1:4).eq.' ') exit !case if other text read(record1(1:55),"(2x,a20,23x,f10.4)") name, qkdat do n = n + 1 if(n.ge.nmin) nmin=n * ptx nmin_save added to get mineral Q/Ks dropped at high P, T because of absent data if(nmin.ge.nmin_save) nmin_save = nmin if((name.eq.rname(n)).or.(rname(n).eq.blank20)) exit end do if(nmin.gt.max7) then write(*,"(/5x,'Number of minerals (max7) exceeded.'/ & ,5x,'Maximum allowed is: ',i3)") max7 stop endif rname(n) = name qklog(n) = qkdat end do * At this point, we reached a blank line. Check to see if there is a * list of minerals (i.e. SUPERSATURATED or UNDERSATURATED or both.) If (k.eq.1) then read(inp1,"(//a80)", end=1002) record1 else read(inp1,"(a80)", end=1002) record1 endif if (record1(20:33).ne.'UNDERSATURATED'.and. + record1(20:33).ne.'SUPERSATURATED') exit IF (record1(20:33).eq.'SUPERSATURATED') + string='Supersaturated minerals' if (record1(20:33).eq.'UNDERSATURATED') + string='Undersaturated minerals' * Skip to top of mineral block, super- or undersaturated read(inp1,"(//a80)",end=1001) record1 end do * Finished reading one set 1002 return 1000 iend=1 !end-of-file when reading top of set return 1001 write(*,"(//5x,'End-of-file encountered when looking for:', + a25/5x,'in SOLVEQ output file'/)") string stop end