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