PROGRAM calculate entropy potential T specific heat and density * Program to calculate the potential temperature, entropy and density C ##################################################################### C # # C # 12/12 2007 # C # FORTRAN program to compute potential temperature # C # density and entropy of geothermal fluid # C # # c # For help on the program, contact # C # Dr. Hongbing Sun # C # Rider University # C # New Jersey # C # # C # Tel: 609-896-5185 # C # Email:hsun@rider.edu # C # # C ##################################################################### C c S: unitless or, psu, t:degC, p: Mpa, pr:Mpa, pr is the reference pressure c Cp: specific heat, kj/kgK C======================================================= C IMPLICIT REAL*8 (a-g,p-t) doubleprecision s,t,p,pr,pottemp,E_freshwater, vapor_s_p character answer*1 C open (22, file='potout.dat',status='unknown') * write (*,*) write (*,*)'#####################################################' write (*,*)'# This program will compute entropy, specific heat,#' write (*,*)'# potential temperature, density of thermal saline #' write (*,*)'# fluid for LIQUID from 0.1-100 MPA, T:0-374oC and #' write (*,*)'# Salinity 0-40. #' write (*,*)'# #' write (*,*)'#####################################################' 2 write (*,*) Write (*,*)' Enter the salinity value' read (*,*) S write(*,*)' Enter the temperature value in degree C' read (*,*) t write (*,*) 'Enter the in-situ pressure value in Mpa' read (*,*) p write (*,*) 'Enter the reference pressure value in Mpa' read (*,*) pr * Check to see whether the input data is within the validity range of the equation if (t.gt.0.0.and.t.le.374.and.p.gt.0.0.and.p.le.100.and.s.ge.0. + and.s.le.40) then if (t.le.100.) goto 3 vapor_s_p =0.00000103*t**3 - 0.00030145*t**2 + 0.02805110*t + -0.45762969 if (p.gt.vapor_s_p) goto 3 endif write (*,*) '==========================================' write (*,*) 'Input data is out of the validity range.' write (*,*) 'Re-etner the data' write (*,*) '==========================================' goto 2 *Reference Pressure cannot be larger than the in-situ pressure 3 if (pr.gt.p) then write (*,*) '==========================================' write (*,*) 'Reference pressure can not be larger than the' write (*,*) 'insitu pressure. Reenter the data' write (*,*) '==========================================' goto 2 endif 4 e1=-entropy(s,t,p) theta=t C* Newton Iternation 5 dt=-(e1+entropy(S,theta,pr))/gtt(s,theta, pr) theta=theta+dt *Loop for potential temperature theta kcount=kcount+1 *End loop when it runs 200 times if (kcount.gt.200) then pottemp=theta kcount=1 goto 100 endif * if (abs(dt).lt.0.00001) then pottemp=theta goto 100 endif goto 5 100 continue C Specifc heat Cp=(t+273.15)*gtt(s,t,p) write (22,*) 'Your Input Data:' write (22,*)'s-psu ','p-Mpa ','pr-Mpa ','t-C ' write (22,*) 'S ','P ',' Pr ', 'T-c ' write (22,30)S,p,pr,T 30 format(F10.5 , 3f12.4) write (22,*) write (22,*) 'Output Data:' write (22,*) 'Potential Temperature',pottemp,' degree C' write (22,*) 'entropy:',-e1, ' Kj/kg/K' e2=entropy(s,pottemp,pr) den=density(s,t,p) den2=density(s,pottemp,pr) write(22,*)'Density (kg/m^3)' write(22,*) den, & ' density with potential temperature of',pottemp,' =', den2 * Print out the data on Screen write (*,*) write (*,*)'====================================================' write (*,*)' Your input data:' write (*,*)' s-psu ',' p-Mpa ',' pr-Mpa ',' t-C ' write (*,50)S,p,pr,T write (*,*) write (*,*) write (*,*)' Output Data:' write (*,*) write (*,*)'====================================================' write (*,*)' Potential Temperature', pottemp,' degree C ' write (*,*) write (*,*)' Entropy=',-e1, ' Kj/kg/K' 50 format(F10.5,3f12.4) write(*,*) write(*,*) ' Density (kg/m^3) ', den write(*,*) write (*,*)' Density with potential temperature of',pottemp & ,' =', den2, ' (Kg/m^3)' write (*,*) write (*,*)' Specific heat Cp= ',cp,' Kj/kg/K' write (*,*) write (*,*)'=====================================================' write (*,*) write (*,*)' Run another set of data? Select Y/N' read (*,'(a1)')answer if (answer.eq.'y'.or.answer.eq.'Y') goto 2 Write (*,*) ' The output data is also written in file potout & .dat' stop END C DOUBLEPRECISION FUNCTION entropy(S,T,P) C======================================================= C Computes specific entropy C C S unitless or in psu, T in C, and P in MPA, Etropy in kJ/kgK C======================================================= C IMPLICIT REAL*8 (A-H,O-Z) doubleprecision E_freshwater,E_freshsaltDiff C Array of polynomial coefficients C error flag DATA ErrorReturn/-9.999999999999D+99/ C C preset: flag for error return Entropy = ErrorReturn C C exclude crudely erratic calls: IF (S.LT. 0 .OR. T .LT. 0 .OR. P .LT. 0) RETURN C** Coefficients fitted from freshwater of Wagner and Pruß(2002) E_freshwater=0.007711828828327622+t*0.015011735629817753 & +t**2*-0.000023741029280375374+t**3*3.754448555306611e-8 & +t**4*-1.5222796931997174e-11+t**5*6.071845579228779e-15 & +p*(-0.00014393952866372566+t*-5.019345250631407e-6 & +t**2*2.5403141497396395e-9 & +t**3*-8.249480473915166e-11+t**4*-5.01693186251528E-15) & +p**2*(-1.219250663011731e-6+t*1.0713947224433905e-8 & +t**2*7.972446351589141e-11 +t**3*1.3348922731422378e-13) & +p**3*(4.384235177142296e-9+t*-4.7087670076904986e-11 & +t**2*-1.684701640927438e-13+t**3*-6.160796224467091e-17) & +p**4*(-6.28067180916623e-12+t*8.181511316370727e-14 & +t**2*8.890263615291897e-17) & +p**5*(2.940936732231923e-15-4.20129614025281e-17*t) C* Salt coefficients derived from the cp fit of Bromely,1974 and Fiestel 2003 E_freshsaltDiff=s*-4.679909745020024e-4 & + s**2*2.645857888308969e-5 & +s**3*-3.5050395334951934e-7+s**4*1.3550718491092872e-9 & +s*t*1.8388961252264425e-5 & +s*t**2*-0.8138308574476893e-7 & +s*t**3*2.5472417662293707e-10 & +s**2*t*-3.648801938210335e-8+s**3*t*2.496295470052084e-10 & +s*p*t*4.3458579758899757e-9 * Total Entropy =Freshwater - Freshsaltwater difference entropy=E_freshwater-E_freshsaltDiff RETURN END DOUBLEPRECISION FUNCTION GTT(S,T,P) C======================================================= C Computes the second order parital differential of Gibb's eneryg. C with respect to temperature. C S in unitless or psu, T in C, and P in MPA, Entropy in kJ/kgK C======================================================= C IMPLICIT REAL*8 (a-g,p-t) C Array of polynomial coefficients C error flag DATA ErrorReturn/-9.999999999999D+99/ C C preset: flag for error return GTT = ErrorReturn C C exclude crudely erratic calls: IF (S.LT.0 .OR. T .LT. 0 .OR. P .LT. 0) RETURN C** Coefficients derived from the freshwater portion of entropy de/dt dE_freshwater=(4.192843058828495 & +t*-0.00022732541214930105+t**2*2.3686269437287123e-6 & +t**4*1.6700924814293252e-10 & +p*(-0.0039782283405605644 & +t*0.00003229142315404513 +t**3*-1.0725210683124268e-9) & +p**2*(0.000019129676509112072 & +t*-4.175829266705576e-7+t**2*2.306273964197821e-9))/(273.15+t) C** Coefficients derived from the salt portion of entropy de/dt dE_freshsaltDiff= s*1.8388961252264425e-5 & +2*s*t*-0.8138308574476893e-7 & +3*s*t**2*2.5472417662293707e-10 & +s**2*-3.648801938210335e-8+s**3*2.496295470052084e-10 & +s*p*4.3458579758899757e-9 gtt=dE_freshwater-dE_freshsaltDiff RETURN END DOUBLEPRECISION FUNCTION density(S,t,P) c===================================================== C Computes density of the thermal fluid. C S in psu, T in C, and P in MPA, density in kg/m^3. C======================================================= C IMPLICIT REAL*8 (a-g,p-t) C Array of polynomial coefficients C error flag DATA ErrorReturn/-9.999999999999D+99/ C C preset: flag for error return density = ErrorReturn C C exclude crudely erratic calls: IF (S.LT.0 .OR. T .LT. 0 .OR. P .LT. 0) RETURN C** dens_freshwater polynomial fitted from Wagner and Pruß(2002) dens_freshwater=9.9920571E+02 & +t*9.5390097E-02+t**2*-7.6186636E-03 & +t**3*3.1305828E-05+t**4*-6.1737704E-08 & +p*(4.3368858E-01+t**2*2.5495667E-05 & +t**3*-2.8988021E-07+t**4*9.5784313E-10) & +p**2*(1.7627497E-03+t*-1.2312703E-04+t**2*1.3659381E-06 & +t**3*-4.0454583E-09) & +p**3*(-1.4673241E-05+t*8.8391585E-07+t**2*-1.1021321E-08 & +t**3*4.2472611E-11+t**4*-3.9591772E-14) C** dens_freshsaltDiff polynomial fitted from Fiestel (2003) and Grunberg,1970 dens_freshsaltDiff=s*-0.7999922297101534 & +s*t*0.002409364998382315+s*t**2*-0.000025805277543352034 & +s*t**3*6.856084054148796e-8 & +p*s*0.0006297611058924313 +p**2*s*-9.362637132680633e-7 density=dens_freshwater-dens_freshsaltDiff RETURN END