      Program gcecontam2
c
c***************************************************************************************
c
c      This program is used to create a formatted INCON file for T2VOC.  It 
c      uses the ROCKS, CHEMP, RPCAP, and ELEME (with element coordinates) blocks, along with
c      a short list of input parameters to generate the initial conditions for
c      the case of gravity-capillary equilibrium, with or without VOC contamination.  The
c      van Genucten capillary pressure curve for air/water (ICP=8) must be used, and Leverett
c      scaling of capillary pressure is assumed for heterogeneous grids using permeability
c      modifiers in the ELEME block.
c
c      If VOC contamination is desired, cylindrical volumes may be contaminated 
c      anywhere in the grid, at any contamination level (using units of total 
c      mass fraction).  The program does a flash calculation to determine the phase
c      conditions (NAPL or not), and assigns the correct primary variables in the
c      generated INCON file.  This file is named "INCON.gce", and it can be renamed
c      "INCON", or it may be pasted into a T2VOC file.
c 
c      In order to use this program, a few lines of input are added to the top of a
c      working T2VOC input file (or at least the ROCKS, CHEMP, RPCAP and ELEME blocks).  The default
c      output from this program to unit 6 is not needed for any purpose, but it may be
c      used to check the assignment of variables, etc. 
c
c      Ron Falta, 9/4/01
c--------use subroutine COWAT to get liquid water density
c--------falta 10/8/01
c
c****************************************************************************************
c
c-------Some historical comments
c
c      Now modified to allow for correlation of contaminant level
c      with permeability using the permeability modifyer
c      correlation may be postitive or negative, with any exponent
c       falta 3/15/01
c     NEW VERSION OF GRAVCAP THAT GENERATES AN INITIAL CONDITION
C     CONTAINING CONTAMINANT DISTRIBUTIONS.  THE CONTAMINANT MAY 
C     EXIST IN A SERIES OF CONSTANT SOIL MASS FRACTION CYLINDERS
C     OF ANY SIZE AT ANY LOCATION WITH ANY CONCENTRATION.
C     FALTA 8/18/00
c
c     gce
c-----Modified 8/14/00 to have gas static pressure dist.  Falta
c
c
c-----this program computes distributions of water saturation
c-----for 3-d heterogeneous tough2 grids (t2voc) assuming
c-----vertical equilibrium with no infiltration. 0/3/97
c
c********************************************************************************************
c     INPUT VARIABLES
c********************************************************************************************
c
c*****The input file contains an initial line with the variables:
c
c     gselev, wtelev, temp, patm and iunits.  
c     format (4e10.4,i1)
c
c     Gselev is the elevation of the
c-----ground surface in grid elevation units, Wtelev is the elevation of the water 
c-----table in grid elevation units, temp is the temperature,
c-----patm is the atmospheric pressure in Pascals,
c-----and iunits is a flag to indicate
c-----whether the units are meters (iunit=0) or feet (iunit=1)
c
c*****The next input line reads the number of contaminated cylinders to be read
c
c     ncyl
c     format (i2)
c
c*****Then, for n=1,ncyl, read in the geometry, location, and soil mass fraction (mg/kg)
c     for each contaminated cylinder.  They may overlap, in which case the contamination
c     is additive.  If no contamination is desired, use a dummy cylinder with dummy coordinates
c     outside of the grid.
c
c     xcyl(n),ycyl(n),zcyl(n),radius(n),thick(n),cmassf(n),expcorr(n)
c     format (7e10.4)
c
c     where xcyl, ycyl, and zcyl are the coordinates of the center of the cylinder,
c     radius is the radius of the cylinder, thick is the thickness of the cylinder, centered 
c     around zcyl, and cmassf is the contaminant soil mass fraction in mg/kg, on a dry weight
c     basis. 
c
c     For an r-z grid, xcyl is used for r, and ycyl is zero.  This approach will probably also work for
c     1-D and 2-D cartesian grids, although it has not been tested. 
c
c     The variable expcorr(n) correlates the mass fraction to a permeability modifier (either the
c     absolute value of a permeability modifier, or the ratio of element perm to average perm) in ELEME
c     by X=Xave*(permod)**expcorr(n) where Xave is the mass fraction read in for
c     each cylinder.  If expcorr is zero, there is no correlation, if it is 1 there is
c     a linear positive correlation with k, if it is -1, there is an inverse linear correlation
c     and so forth.  Note that each cylinder may have a different correlation.
c
c     A value of expcorr=.5 would correlate contaminant inversley with Pcap and -.5 would be
c     a linear correlation with Pcap.
c
c
c
c     Note that the 5 geometry variables must have the same units as the grid.
c
c**************************************************************************************************
***************************************************************************************************
c 
c-----It should work to just add these lines to the top of a 
c-----working T2VOC file
c
c-----Otherwise, the ROCKS block is pasted in, followed by the CHEMP and RPCAP blocks,
c-----finishing with the ELEME block.  They must be in this order, and
c-----both the ROCKS and ELEME blocks must end with blank lines
c-----Input file should end in ENDCY
c----
c-----The program returns a formatted INCON file for a air/water/chemical/heat
c-----system
c
c------Ron Falta 9/3/97
c
      IMPLICIT REAL*8(A-H,O-Z)
      dimension DM(27),POR(27),PER(3,27),CWET(27),SH(27)
      dimension MAT(27)
      dimension COM(27),EXPAN(27),CDRY(27),TORT(27),GK(27)
      dimension FOC(27)
      dimension IRP(27),RP(7,27),ICP(27),CP(7,27),RPD(7),CPD(7)
	dimension xcyl(50),ycyl(50),zloc(50),radius(50),thick(50)
	dimension cmassf(50),expcorr(50)
      COMMON/CRITP/TCRIT,PCRIT,ZCRIT,OMEGA,DIPOLM                       
      COMMON/VPOIL/TBOIL,VPA,VPB,VPC,VPD                                
      COMMON/HCAPL/AMO,CPA,CPB,CPC,CPDD                                 
      COMMON/DENOIL/RHOREF,TDENRF,DIFV0,TDIFRF,TEXPO                    
      COMMON/VISOIL/VLOA,VLOB,VLOC,VLOD,VOLCRT                          
      COMMON/HCOIL/SOLA,SOLB,SOLC,SOLD                                  
      DIMENSION VER(5),XIN(4,27),NAM(27)
c
      CHARACTER*5 ELEM1,ELEM2,MAT,ELST,VER,WORD,MA12,NAM
      CHARACTER*3 EL,MA1,EL1,EL2,SL
      CHARACTER MA2*2,TYPE*4,ITAB*1
C
      DATA VER /'CHEMP','ELEME','ROCKS','RPCAP','ENDCY'/
C
C
c-----read in ground surface elevation,water table elevation, 
c-----temperature, atmospheric pressure, and units flag
c
      read 8002, gselev,wtelev, temp, patm,iunits
 8002 format (4e10.4,i1)
c   read in the contaminated cylinders
      read 8018, ncyl
 8018 format (i2)
c
      do 8020 n=1,ncyl
	read 8017, xcyl(n),ycyl(n),zloc(n),radius(n),
     &thick(n),cmassf(n),expcorr(n)
 8017 format (7e10.4)
 8020 continue	      

c---- calculate gas density at atmospheric pressure
c
      rhog=((patm)*29.0d0)/(8314.4d0*(temp+273.15d0))
c
c---- water density
      call cowat(temp,patm,rhow,enthalpy)
c
c---- calculate gas pressure at water table
c
      if(iunits.eq.0) pwt=patm+(gselev-wtelev)*9.806d0*rhog
	if(iunits.eq.1) pwt=patm+(gselev-wtelev)*0.3048d0*9.806*rhog
c
 5019 READ 5020,WORD
 5020 FORMAT(A5)
C
      DO900 K=1,5
  900 IF(WORD.EQ.VER(K)) GOTO920
c     PRINT 901,WORD
c 901 FORMAT(' HAVE READ UNKNOWN BLOCK LABEL "',A5,
c    X'" --- IGNORE THIS, AND CONTINUE READING INPUT DATA')
      GOTO 5019
C
  920 GOTO(2600,1100,1700,2200,1500),K
C
C*****READ ROCK PROPERTIES.*********************************************
C                                                                      *
 1700 NM=1
    2 READ 1,MAT(NM),NAD,DM(NM),POR(NM),(PER(I,NM),I=1,3),
     ACWET(NM),SH(NM)
    1 FORMAT(A5,I5,7E10.4)
C
      IF(MAT(NM).EQ.'     ') GOTO 3
C
      COM(NM)=0.D0
      EXPAN(NM)=0.D0
      CDRY(NM)=0.D0
      TORT(NM)=0.D0
      GK(NM)=0.D0
      IRP(NM)=0
      ICP(NM)=0
      IF(NAD.GE.1) READ 7,COM(NM),EXPAN(NM),CDRY(NM),TORT(NM),GK(NM)
     X,FOC(NM)
    7 FORMAT(8E10.4)
      IF(CDRY(NM).EQ.0.D0) CDRY(NM)=CWET(NM)
      IF(NAD.LE.1) GOTO 1701
C-----READ PARAMETERS FOR RELATIVE PERMEABILITY AND CAPILLARY PRESSURE.
      READ 1702,IRP(NM),(RP(I,NM),I=1,7)
      READ 1702,ICP(NM),(CP(I,NM),I=1,7)
 1702 FORMAT(I5,5X,7E10.4)
 1701 CONTINUE
C
      NAD=0
C
      PRINT 5,NM,MAT(NM)
    5 FORMAT(/11H DOMAIN NO.,I3,22H     MATERIAL NAME -- ,A5)
      NM=NM+1
      GOTO2
    3 CONTINUE
      NM=NM-1
      GOTO5019
C
C-----END OF ROCK PROPERTIES.-------------------------------------------
C
 2200 CONTINUE
C-----COME HERE TO READ DEFAULT ASSIGNMENTS FOR RELATIVE PERMEABILITY
C     AND CAPILLARY PRESSURE PARAMETERS.
      READ 1702,IRPD,(RPD(I),I=1,7)
      READ 1702,ICPD,(CPD(I),I=1,7)
      GOTO 5019
C
C
 2600 CONTINUE
C-----ASSIGN VOC CHEMICAL PARAMETERS ('CHEMP')
      READ 2002,TCRIT,PCRIT,ZCRIT,OMEGA,DIPOLM
      READ 2002,TBOIL,VPA,VPB,VPC,VPD
      READ 2002,AMO,CPA,CPB,CPC,CPDD
      READ 2002,RHOREF,TDENRF,DIFV0,TDIFRF,TEXPO
      READ 2002,VLOA,VLOB,VLOC,VLOD,VOLCRT
      READ 2002,SOLA,SOLB,SOLC,SOLD
      READ 2002,OCK,FOX,ALAM
 2002 FORMAT(8E10.4)
c
c     calculate aqueous solubility (mole fr.)
      call hco(temp,solub)
c     calculate NAPL vapor pressure
      call sato(temp,pvap)
c     calculate NAPL density
      call cowato(temp,dloil)
c-----calculate max gas concen
      print *,'tcrit pcrit zcrit omega dipolm'
	print 2002,tcrit,pcrit,zcrit,omega,dipolm
	print *,'tboil,vpa,vpb,vpc,vpd'
	print 2002, tboil,vpa,vpb,vpc,vpd
	print *, 'amo=', amo
      print *,' temp = ',temp
      print *, 'pvap =',pvap
      cmaxg=pvap*amo/(8314.4*(temp+273.15))
c     calculate max aqueous concen
      xmassc=(solub*amo)/(solub*amo+(1.-solub)*18.0)
	denw=1./((1-xmassc)/999.0+xmassc/dloil)
	cmaxw=xmassc*denw
	print *, 'cmaxg=',cmaxg,'cmaxw=',cmaxw
c
      GOTO5019
C
C*****READ ELEMENT DATA, and compute primary 
c-----variables for  INCON block*****
C
 1100 CONTINUE
      open(unit=4,file='INCON.gce',status='unknown')
      PRINT 1103
 1103 FORMAT(' WRITE FILE *INCON.gce* assuming gravity/capillary',
     &' equilibrium')
      WRITE(4,1101)
 1101 FORMAT(5HINCON)
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c   the main element loop starts here
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      nele=0
 1102 READ 10,EL,NE,NSEQ,NADD,MA1,MA2,VOLX,AHTX,permod,X,Y,Z
      nele=nele+1
c
   10 FORMAT(A3,I2,2I5,A3,A2,2E10.4,4E10.4)
      IF(EL.EQ.'   '.AND.NE.EQ.0) GOTO40
C
      NSEQ1=NSEQ+1
C
c------now for correlation of contamination with permod----falta 3/15/01
c------ratio is the ratio of perm to average perm
c
       if(permod.lt.0.) ratio=-permod
	 if(permod.gt.0.) ratio=permod/per(1,m)
	 if(permod.eq.0.) ratio=1.0
c
c **** now loop over the contaminated cylinders, to see which ones
c      (if any) this gridblock resides in.  Add up the total mass fraction
c      of contaminant in this gridblock
c     
      tmassf=0.0
	icflag=0
      do 8060 i=1,ncyl
       rad=sqrt((x-xcyl(i))**2+(y-ycyl(i))**2)
       if(rad.gt.radius(i)) go to 8060
	 vert=abs(z-zloc(i))
	 if(vert.gt.thick(i)/2.) go to 8060
	 cormassf=cmassf(i)*ratio**expcorr(i)
	 tmassf=tmassf+cormassf
	 icflag=1
 8060 continue
      print *, 'element number', nele,'tmassf=',tmassf
c
c
      IF(MA1.EQ.'   ') GOTO15
C
C-----FIND MATERIAL INDEX.
      MA12=MA1//MA2
      DO11 M=1,NM
      IF(MA12.EQ.MAT(M)) GOTO 12
   11 CONTINUE
      PRINT 9,MA1,MA2,EL,NE
    9 FORMAT(31H REFERENCE TO UNKNOWN MATERIAL ,A3,A2,12H AT ELEMENT ,A2
     A,I3,24H --- WILL IGNORE ELEMENT)
      GOTO1102
C
   12 CONTINUE
C-----COME HERE FOR ELEMENTS WITH MATERIAL NAME.
C
c-----assume constant z for nseq not equal to zero
      DO13 I=1,NSEQ1
      N=NE+(I-1)*NADD
c-----check for position of water table
      if(wtelev.ge.z) then
c--------we are below the water table
	 if(iunits.eq.0) p=(wtelev-z)*9.806*rhow+pwt
	 if(iunits.ne.0) p=(wtelev-z)*9.806*rhow*0.3048+pwt
	 xairw=1.e-8
	 write(4,8010)el,n,por(m)
 8010    format(a3,i2,10x,e15.8)
c   no contamination, below wt.
	 if(icflag.eq.0) then
	   write(4,8011)p,xairw+50.0,0.0,temp
 8011    format(4e20.13)
         go to 1102
       endif
c   we do have contamination
c   calculate ctmax and check for napl
       bulkdens=dm(m)*(1.0-por(m))
	 if(foc(m).eq.0.0) foc(m)=fox
	 cmaxads=bulkdens*ock*foc(m)*cmaxw
	 ctmax=por(m)*cmaxw+cmaxads
	 ctot=bulkdens*tmassf/1.e6
	print *, 'ctmax=',ctmax,'ctot=',ctot
c  we have napl
	 if(ctot.gt.ctmax) then
	  sn=(ctot-ctmax)/(por(m)*dloil)
	print *, ' sn = ',sn
	  write (4,8011) p,1.0-sn,xairw+50.0,temp
       else
c  no napl
	  xmolow=solub*ctot/ctmax
	  write (4,8011) p,xairw+50.0,xmolow,temp
       endif
c go to next element
c
      else 
c
c-----   we are above the water table
	 if(iunits.eq.0) pcap=(z-wtelev)*9.806*rhow
	 if(iunits.eq.0) p=(wtelev-z)*rhog*9.806+pwt
	 if(iunits.ne.0) pcap=(z-wtelev)*9.806*rhow*0.3048
	 if(iunits.ne.0) p=(wtelev-z)*rhog*9.806*.3048+pwt
	 print 9002, el,n,pcap
 9002    format(' element ',a3,i2,' pcap = ',f15.5)
c--------get capillary parameters -check for default values
	 if(icp(m).eq.0) then
	      icp(m)=icpd
	      cp(1,m)=cpd(1)
	      cp(2,m)=cpd(2)
	      cp(3,m)=cpd(3)
	      cp(4,m)=cpd(4)
	 endif
	 if(icp(m).lt.8) then
	    print 8050
 8050       format (' you must use ICP>=8 option here ')
	    stop
	 endif
	  if(icp(m).eq.9) then
		   sw=.1
	       goto8900
	    endif
	  sm=cp(1,m)
	 an=cp(2,m)
	 alphgn=cp(3,m)
	 alphnw=cp(4,m)
	 alphgw=1./(1/alphgn+1/alphnw)
	 if(permod.lt.0.) scale=dsqrt(-1./permod)
	 if(permod.gt.0.) scale=dsqrt(per(1,m)/permod)
	 if(permod.eq.0.) scale=1.0
	
	 const=scale*9796./alphgw
	 am=1.-1./an
	 sw=(1.-sm)*((pcap/const)**an+1)**(-am)+sm
c	  print 9003, el,n,permod,scale,sm,an,alphgn,alphnw,sw  
c9003     format(1x,a3,i2,'pmod=',f10.4,'scale=',f10.4,5f10.4)         
 8900    write(4,8010)el,n,por(m)
c  no contamination, above water table
       if(icflag.eq.0) then
	  write(4,8011)p,sw,0.0d0,temp
	  go to 1102
	 endif
c  we do have contamination
c   calculate ctmax and check for napl
       bulkdens=dm(m)*(1.0-por(m))
	 if(foc(m).eq.0.0) foc(m)=fox
	 cmaxads=bulkdens*ock*foc(m)*cmaxw
	 ctmax=por(m)*sw*cmaxw+por(m)*(1.-sw)*cmaxg+cmaxads
	 ctot=bulkdens*tmassf/1.e6
	print *, 'ctmax=',ctmax,'ctot=',ctot
c  we have napl
	 if(ctot.gt.ctmax) then
	  sn=(ctot-ctmax)/(por(m)*dloil)
	print *, ' sn = ',sn
c  displace the gas and water in proportion to relative distance
c  from the water residual sat and sw=1
        swfrac=(sw-sm)/(1.0-sm)
	  sgfrac=1.0-swfrac
	  swnew=sw-swfrac*sn
	  sgnew=(1.0-sw)-sgfrac*sn
	  write (4,8011) p,swnew,sgnew+10.0d0,temp
       else
c  no napl
	  xmolog=pvap/p*ctot/ctmax
	  write (4,8011) p,sw,xmolog,temp
       endif
c
c  go to next element
      endif
c
   13 continue 
c-----read the next element
      GOTO1102
C
c
c
C-----COME HERE FOR ELEMENTS WITH MATERIAL NUMBER.
C
   15 CONTINUE
c------convert ma2 to integer
      read (ma2,'(i2)') m
      DO113 I=1,NSEQ1
      N=NE+(I-1)*NADD
c-----check for position of water table and write to incon
      if(wtelev.ge.z) then
c--------we are below the water table
	 if(iunits.eq.0) p=(wtelev-z)*9.806*rhow+pwt
	 if(iunits.ne.0) p=(wtelev-z)*9.806*rhow*0.3048+pwt
	 xairw=1.e-8
	 write(4,8010)el,n,por(m)
c no contamination
	 if(icflag.eq.0) then
	  write(4,8011)p,xairw+50.0,0.0d0,temp
	  go to 1102
       endif
c we have contamination
c   calculate ctmax and check for napl
       bulkdens=dm(m)*(1.0-por(m))
	 if(foc(m).eq.0.0) foc(m)=fox
	 cmaxads=bulkdens*ock*foc(m)*cmaxw
	 ctmax=por(m)*cmaxw+cmaxads
	 ctot=bulkdens*tmassf/1.e6
	print *, 'ctmax=',ctmax,'ctot=',ctot
c  we have napl
	 if(ctot.gt.ctmax) then
	  sn=(ctot-ctmax)/(por(m)*dloil)
	print *, ' sn = ',sn
	  write (4,8011) p,1.0-sn,xairw+50.0,temp
       else
c  no napl
	  xmolow=solub*ctot/ctmax
	  write (4,8011) p,xairw+50.0,xmolow,temp
       endif 
c go to next element	     
      else
c
c-----   we are above the water table
	 if(iunits.eq.0) pcap=(z-wtelev)*9.806*rhow
	 if(iunits.eq.0) p=(wtelev-z)*rhog*9.806+pwt
	 if(iunits.ne.0) pcap=(z-wtelev)*9.806*rhow*0.3048
	 if(iunits.ne.0) p=(wtelev-z)*rhog*9.806*0.3048+pwt
c--------get capillary parameters -check for default values
	 if(icp(m).eq.0) then
	      icp(m)=icpd
	      cp(1,m)=cpd(1)
	      cp(2,m)=cpd(2)
	      cp(3,m)=cpd(3)
	      cp(4,m)=cpd(4)
	 endif
	 if(icp(m).lt.8) then
	    print 8050
	    stop
	 endif
	 if(icp(m).eq.9) then
		   sw=.1
	       goto8901
	    endif
	 sm=cp(1,m)
	 an=cp(2,m)
	 alphgn=cp(3,m)
	 alphnw=cp(4,m)
	 alphgw=1./(1./alphgn+1./alphnw)
	 if(permod.lt.0.) scale=dsqrt(-1./permod)
	 if(permod.gt.0.) scale=dsqrt(per(1,m)/permod)
	 if(permod.eq.0.) scale=1.0
	 const=scale*9796./alphgw
	 am=1.-1./an
	 sw=(1.-sm)*((pcap/const)**an+1)**(-am)+sm
 8901    write(4,8010)el,n,por(m)
c  no contamination
	 if(icflag.eq.0) then
	  write(4,8011)p,sw,0.0d0,temp
	  go to 1102
	 endif
c  we have contamination
c
c   calculate ctmax and check for napl
       bulkdens=dm(m)*(1.0-por(m))
	 if(foc(m).eq.0.0) foc(m)=fox
	 cmaxads=bulkdens*ock*foc(m)*cmaxw
	 ctmax=por(m)*sw*cmaxw+por(m)*(1.-sw)*cmaxg+cmaxads
	 ctot=bulkdens*tmassf/1.e6
	print *, 'ctmax=',ctmax,'ctot=',ctot
c  we have napl
	 if(ctot.gt.ctmax) then
	  sn=(ctot-ctmax)/(por(m)*dloil)
	print *, ' sn = ',sn
c  displace the gas and water in proportion to relative distance
c  from the water residual sat and sw=1
        swfrac=(sw-sm)/(1.0-sm)
	  sgfrac=1.0-swfrac
	  swnew=sw-swfrac*sn
	  sgnew=(1.0-sw)-sgfrac*sn
	  write (4,8011) p,swnew,sgnew+10.0d0,temp
       else
c  no napl
	  xmolog=pvap/p*ctot/ctmax
	  write (4,8011) p,sw,xmolog,temp
       endif
c
c
      endif
c
  113 continue 
c
      GOTO1102
C
C-----END OF ELEMENT DATA.----------------------------------------------
C
   40 WRITE(4,41)
   41 FORMAT(5H     )
      GOTO5019
 1500 print 8051, nele
 8051 format (' done computing initial conditions for ',i5,'elements')
      stop
      end
c
c
      SUBROUTINE SATO(T,PSATO)
      IMPLICIT REAL*8(A-H,O-Z)
C
C-----CALCULATE HYDROCARBON VAPOR PRESSURE FROM WAGNER EQN USING
C-----CONSTANTS FOUND IN REID ET AL., 1987----------C
C
      COMMON/CRITP/TCRIT,PCRIT,ZCRIT,OMEGA,DIPOLM
      COMMON/VPOIL/TBOIL,VPA,VPB,VPC,VPD
      TR=(T+273.15D0)/TCRIT
      X=1.D0-TR
      IF(X.LE.0.D0) GO TO 10
      IF(T.LE.0.D0) GO TO 10
      X2=VPA*X+VPB*X**1.5D0+VPC*X**3+VPD*X**6
      X3=X2/TR
      PSATO=1.D5*PCRIT*dEXP(X3)
      RETURN
  10  IGOOD=2
      PRINT 11,T
   11 FORMAT(' TEMPERATURE = ',E12.5,'  OUT OF RANGE IN SATO')
      RETURN
      END
C
      SUBROUTINE COWATO(T,DLO)
      IMPLICIT REAL*8(A-H,O-Z)
C
C-----DENSITY OF LIQUID HYDROCARBON
C----REID ET AL., 1987
C
      COMMON/CRITP/TCRIT,PCRIT,ZCRIT,OMEGA,DIPOLM
      COMMON/DENOIL/RHOREF,TDENRF,DIFV0,TDIFRF,TEXPO
      COMMON/HCAPL/AMO,CPA,CPB,CPC,CPDD
      TKEL=T+273.15D0
      TR=TKEL/TCRIT
      TRREF=TDENRF/TCRIT
      ZRA=0.29056D0-0.08775D0*OMEGA
      PHI=(1.D0-TRREF)**0.2857D0-(1.D0-TR)**0.2857D0
      DLO=RHOREF*ZRA**PHI
      RETURN
      END
C
c
      SUBROUTINE HCO(T,sol)
      IMPLICIT REAL*8(A-H,O-Z)
C
C-----CALCULATE solubility
      COMMON/HCOIL/SOLA,SOLB,SOLC,SOLD
      SOL=SOLA+SOLB*T+SOLC*T**2+SOLD*T**3
      RETURN
      END
c
      SUBROUTINE COWAT(TF,PP,D,U)
C--------- Fast COWAT M.J.O'Sullivan - 17 SEPT 1990 ---------
C     20 September 1990.  VAX needs double precision, CRAY does not.
      IMPLICIT REAL*8 (A-H,O-Z)
c      IMPLICIT REAL (A-H,O-Z)
c      COMMON/KONIT/KON,DELT,IGOOD
      save icall
      DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12 /
     16.824687741E3,-5.422063673E2,-2.096666205E4,3.941286787E4,
     2-13.466555478E4,29.707143084E4,-4.375647096E5,42.954208335E4,
     3-27.067012452E4,9.926972482E4,-16.138168904E3,7.982692717E0/
      DATA A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23 /
     4-2.616571843E-2,1.522411790E-3,2.284279054E-2,2.421647003E2,
     51.269716088E-10,2.074838328E-7,2.174020350E-8,1.105710498E-9,
     61.293441934E1,1.308119072E-5,6.047626338E-14/
      DATA SA1,SA2,SA3,SA4,SA5,SA6,SA7,SA8,SA9,SA10,SA11,SA12 /
     18.438375405E-1,5.362162162E-4,1.720000000E0,7.342278489E-2,
     24.975858870E-2,6.537154300E-1,1.150E-6,1.51080E-5,
     31.41880E-1,7.002753165E0,2.995284926E-4,2.040E-1   /
C
      DATA ICALL/0/
      ICALL=ICALL+1
c      IF(ICALL.EQ.1) WRITE(11,899)
c  899 FORMAT(6X,'COWAT    1.0 S    17 SEPTEMBER 1990',6X,
c     X'LIQUID WATER DENSITY AND INT.',
c     X' ENERGY VERSUS TEMPERATURE AND PRESSURE (M. OS.)')
C
      TKR=(TF+273.15)/647.3
      TKR2=TKR*TKR
      TKR3=TKR*TKR2
      TKR4=TKR2*TKR2
      TKR5=TKR2*TKR3
      TKR6=TKR4*TKR2
      TKR7=TKR4*TKR3
      TKR8=TKR4*TKR4
      TKR9=TKR4*TKR5
      TKR10=TKR4*TKR6
      TKR11=TKR*TKR10
      TKR19=TKR8*TKR11
      TKR18=TKR8*TKR10
      TKR20=TKR10*TKR10
      PNMR=PP/2.212E7
      PNMR2=PNMR*PNMR
      PNMR3=PNMR*PNMR2
      PNMR4=PNMR*PNMR3
      Y=1.-SA1*TKR2-SA2/TKR6
      ZP=SA3*Y*Y-2.*SA4*TKR+2.*SA5*PNMR
      IF(ZP.LT.0.) GOTO 1
C     19 September 1990.  double on VAX, single on CRAY
      Z=Y+ DSQRT(ZP)
c      Z=Y+SQRT(ZP)
      CZ=Z**(5./17.)
      PAR1=A12*SA5/CZ
      CC1=SA6-TKR
      CC2=CC1*CC1
      CC4=CC2*CC2
      CC8=CC4*CC4
      CC10=CC2*CC8
      AA1=SA7+TKR19
      PAR2=A13+A14*TKR+A15*TKR2+A16*CC10+A17/AA1
      PAR3=(A18+2.*A19*PNMR+3.*A20*PNMR2)/(SA8+TKR11)
      DD1=SA10+PNMR
      DD2=DD1*DD1
      DD4=DD2*DD2
      PAR4=A21*TKR18*(SA9+TKR2)*(-3./DD4+SA11)
      PAR5=3.*A22*(SA12-TKR)*PNMR2+4.*A23/TKR20*PNMR3
      VMKR=PAR1+PAR2-PAR3-PAR4+PAR5
      V=VMKR*3.17E-3
      D=1./V
      YD=-2.*SA1*TKR+6.*SA2/TKR7
      SNUM= A10+A11*TKR
      SNUM=SNUM*TKR + A9
      SNUM=SNUM*TKR + A8
      SNUM=SNUM*TKR + A7
      SNUM=SNUM*TKR + A6
      SNUM=SNUM*TKR + A5
      SNUM=SNUM*TKR + A4
      SNUM=SNUM*TKR2 - A2
      PRT1=A12*(Z*(17.*(Z/29.-Y/12.)+5.*TKR*YD/12.)+SA4*TKR-
     1(SA3-1.)*TKR*Y*YD)/CZ
      PRT2=PNMR*(A13-A15*TKR2+A16*(9.*TKR+SA6)*CC8*CC1
     2+A17*(19.*TKR19+AA1)/(AA1*AA1))
      BB1=SA8+TKR11
      BB2=BB1*BB1
      PRT3=(11.*TKR11+BB1)/BB2*(A18*PNMR+A19*
     3PNMR2+A20*PNMR3)
      EE1=SA10+PNMR
      EE3=EE1*EE1*EE1
      PRT4=A21*TKR18*(17.*SA9+19.*TKR2)*(1./EE3+SA11*PNMR)
      PRT5=A22*SA12*PNMR3+21.*A23/TKR20*PNMR4
      ENTR= A1*TKR - SNUM +PRT1+PRT2-PRT3+PRT4+PRT5
      H=ENTR*70120.4
      U=H-PP*V
      RETURN
    1 IGOOD=2
c      WRITE(6,2)TF
c      WRITE(7,2)TF
c    2 FORMAT(' TEMPERATURE = ',E12.6,'  OUT OF RANGE IN COWAT')
c  100 FORMAT(1H ,5X,A6,2X,E20.10)
c  102 FORMAT(1H ,5X,A6,5X,I2,2X,E20.10)
      RETURN
      END
