MODULE enxen_module !### !***************************************************************************** ! * ! Module Name: enxen_module * ! * ! File Name: enxen_module.f90 * ! * ! Purpose: Contains interface for subroutine enxen * ! * ! Author: Robert L. ***** * ! * ! Date: March 03, 2004 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Cosmpatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 03/04/04 0 Initial Release * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! USE enxen_module * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! None * ! * ! Global Variables: * ! ----------------- * ! * ! None * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! None * ! * ! Modules: * ! -------- * ! None * ! * ! Interfaces: * ! ----------- * ! * ! enxen - M R BOLDUC subroutine to solve an n x n system of simultaneous* ! linear equations by the process of diagonalization. * ! * !***************************************************************************** !### IMPLICIT none INTERFACE enxen_sub SUBROUTINE enxen( & lv, & stop_flag) INTEGER,INTENT(inout) :: lv LOGICAL,INTENT(inout) :: stop_flag END SUBROUTINE enxen END INTERFACE enxen_sub END MODULE enxen_module MODULE read_error_module !### !***************************************************************************** ! * ! Module Name: read_error_module * ! * ! File Name: read_error_module.f90 * ! * ! Purpose: Contains interface for subroutine read_error * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 2003 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Cosmpatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/03 0 Initial Release * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! USE read_error_module * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! None * ! * ! Global Variables: * ! ----------------- * ! * ! None * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! None * ! * ! Modules: * ! -------- * ! None * ! * ! Interfaces: * ! ----------- * ! * ! read_error - Process the status returned from a read * ! * !***************************************************************************** !### IMPLICIT none INTERFACE error_read SUBROUTINE read_error(unit, & message, & status) IMPLICIT none INTEGER,INTENT(in) :: message INTEGER,INTENT(inout) :: status INTEGER,INTENT(in) :: unit END SUBROUTINE read_error END INTERFACE error_read END MODULE read_error_module MODULE lease_module !### !***************************************************************************** ! * ! Module Name: lease_module * ! * ! File Name: lease_module.f90 * ! * ! Purpose: Contains constants and globals for lease * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 2003 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Cosmpatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/03 0 Initial Release * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! USE lease_module * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! None * ! * ! Global Variables: * ! ----------------- * ! * ! a - Real array containing elements of the matrix * ! ac - Real array containing coefficients of an n degree polynominal * ! ans - Real array containing coefficients derived from matrix * ! solution * ! av - Logical variable used to determine if points are to be * ! constrained to fit at 1000 degrees Kelvin * ! a_store- Real array used to store sum_all for transfer to subroutine * ! enxen * ! b - Real array containing stoichiometric coefficient of elements * ! in current specie * ! c - Real array containing answer region of the matrix solution * ! cp - Real array of computed heat capacity * ! cpr - Real array containing cp/r calculated from the coefficients * ! cpsum - Real array containing summation of heat capacity for nd * ! points * ! date1 - Character*3 variable containing reference code of source of * ! data used in the calculations * ! date2 - Character*3 variable containing reference date of source of * ! data used in the calculations * ! hcalc - Real array containing computed heat content * ! hhrt - Real array containing heat content of species divided by rt * ! hsum - Real array containing summation of computed heat content for * ! nd points * ! hth - Real array containing heat content of specie at T * ! htho - Real array containing heat content converted to * ! calories * ! hthoe - Real array containing percentage of error between computed * ! heat content and heat content read from input records * ! infile - Character*80 variable containing input file name * ! ispeci - Character*4 variable containing alphanumeric name of specie * ! being considered * ! kk - Integer array used as position of transition T * ! lp - Integer variable used as upper limit of values being stored * ! m - Integer value used to determine which npt value is being * ! considered * ! mp - Integer variable used as location of minimum value of points * ! being considered * ! mt - Character*2 array containing alphanumeric name of element in * ! current specie * ! nd - Integer value used as number of readings being considered in * ! one pass * ! ndeg - Integer value used as number of degrees in polynominal * ! npts - Integer array containing total number of points being * ! considered * ! outfile- Character*80 variable containing output file name * ! phaz - Character*1 variable containing phase of species being * ! considered * ! re_read- Character*80 buffer used to read previous record * ! sr - Real array used as temporary storage of st/r * ! stop_it- Logical variable use to test if stop condition met * ! sum_all- Real array containing summation of pi * ! st - Real array containing entropy read from input records * ! t1 - Real variable containing lower limit of T range of specie * ! t2 - Real variable containing upper limit of T range of specie * ! tcp - Real array used as temporary storage of calculated heat * ! capacity * ! thcalc - Real array used as temporary storage of calculated heat * ! content * ! tt - Real array used as multiplier factor in final computation of * ! heat capacity * ! vol - Real variable containing molar volume of the specie * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! av_kount - Integer value of page line count after ave listing * ! av_value - Real value used to compute averager value * ! coeff_set - Integer value of number of coefficients in a set * ! date_error - Integer value of error message for date error * ! high_coeff - Integer index to high set of coefficients * ! init_kount - Integer initial value of page line count * ! inp - Integer value of logical input unit = 5 * ! kount_inc - Integer value of kount increment and new init value * ! last - Character*4 parameter set to 'LAST' * ! ldeg_start - Integer value used as start of ldeg index in eqfind * ! ls_inc - Integer value to add to ls in routine eqfind * ! ls_start - Integer value used as start of ls index in eqfind * ! lv_start - Integer value used as start of lv index in eqfind * ! low_coeff - Integer index to low set of coefficients * ! matrix_size - Integer value of number of rows and columns in matrix * ! max_kount - Integer value of maximum lines per page * ! message - Chracter array containing error messages to output * ! message_offset - Integer value for offset of message number to use * ! ndeg_start - Integer value used as start of ndeg index in eqfind * ! newpage - Character containing new page control character * ! npts_error - Integer value of error message for npts error * ! out - Integer value of logical output unit = 6 * ! percent - Real parameter containing percentage conversion value * ! r - Real parameter containing universal gas constant * ! save_size - Integer value of size of saved arrays in subroutine * ! ave * ! specie_error - Integer value of error message for specie error * ! temp_break - Real value of temperature break over point for sets of* ! coefficients * ! temp_conv - Real parameter containing temperature conversion * ! divisor * ! temp_error - Integer value of error message for temp error * ! * ! Modules: * ! -------- * ! None * ! * ! Interfaces: * ! ----------- * ! * ! None * ! * !***************************************************************************** !### IMPLICIT none LOGICAL :: av LOGICAL :: stop_it REAL,DIMENSION(10,11),SAVE:: a REAL,DIMENSION(100),SAVE :: a_store REAL,DIMENSION(60),SAVE :: ac REAL,DIMENSION(2,7),SAVE :: ans REAL,DIMENSION(10),SAVE :: b REAL,DIMENSION(10),SAVE :: c REAL,DIMENSION(60),SAVE :: cp REAL,DIMENSION(60),SAVE :: cpr REAL,DIMENSION(60),SAVE :: cpsum REAL,DIMENSION(60),SAVE :: hcalc REAL,DIMENSION(60),SAVE :: hhrt REAL,DIMENSION(60),SAVE :: hsum REAL,DIMENSION(60),SAVE :: hth REAL,DIMENSION(60),SAVE :: htho REAL,DIMENSION(60),SAVE :: hthoe REAL,DIMENSION(60),SAVE :: sr REAL,DIMENSION(60),SAVE :: st REAL,DIMENSION(60),SAVE :: sum_all REAL,DIMENSION(5,60),SAVE :: tcp REAL,DIMENSION(5,60),SAVE :: thcalc REAL,DIMENSION(60),SAVE :: tt REAL,SAVE :: t1 REAL,SAVE :: t2 REAL,SAVE :: vol INTEGER,SAVE :: lp INTEGER,SAVE :: m INTEGER,SAVE :: mp INTEGER,SAVE :: nd INTEGER,SAVE :: ndeg INTEGER,DIMENSION(2),SAVE :: kk INTEGER,DIMENSION(2),SAVE :: npts REAL,PARAMETER :: av_value = 2.0 REAL,PARAMETER :: percent = 100.0 REAL,PARAMETER :: r = 1.98726 REAL,PARAMETER :: temp_break= 1000.0 REAL,PARAMETER :: temp_conv = 1000.0 CHARACTER(len=4),PARAMETER :: last = 'LAST' CHARACTER(len=1),PARAMETER :: newpage = '1' INTEGER,PARAMETER :: av_kount = 28 INTEGER,PARAMETER :: coeff_set = 7 INTEGER,PARAMETER :: date_error = 2 INTEGER,PARAMETER :: high_coeff = 1 INTEGER,PARAMETER :: inp = 5 INTEGER,PARAMETER :: init_kount = 24 INTEGER,PARAMETER :: kount_inc = 3 INTEGER,PARAMETER :: ldeg_start = 2 INTEGER,PARAMETER :: low_coeff = 2 INTEGER,PARAMETER :: ls_inc = 2 INTEGER,PARAMETER :: ls_start = 2 INTEGER,PARAMETER :: lv_start = 2 INTEGER,PARAMETER :: matrix_size = 10 INTEGER,PARAMETER :: max_kount = 60 INTEGER,PARAMETER :: message_offset= 2 INTEGER,PARAMETER :: ndeg_start = 2 INTEGER,PARAMETER :: npts_error = 3 INTEGER,PARAMETER :: out = 6 INTEGER,PARAMETER :: save_size = 60 INTEGER,PARAMETER :: specie_error = 1 INTEGER,PARAMETER :: temp_error = 4 CHARACTER(len=4),DIMENSION(8),SAVE :: ispeci CHARACTER(len=2),DIMENSION(10),SAVE :: mt CHARACTER(len = 80) :: re_read CHARACTER(len=3),SAVE :: date1 CHARACTER(len=3),SAVE :: date2 CHARACTER(len=80),SAVE:: infile CHARACTER(len=80),SAVE:: outfile CHARACTER(len=1),SAVE :: phaz CHARACTER(len = *),DIMENSION(8),PARAMETER :: message = & (/' ERROR IN ABOVE SPECIE RECORD ', & ' EOF READ WHILE READING SPECIE RECORD ', & ' ERROR IN ABOVE DATE RECORD ', & ' EOF READ WHILE READING DATE RECORD ', & ' ERROR IN ABOVE NPTS RECORD ', & ' EOF READ WHILE READING NPTS RECORD ', & ' ERROR IN ABOVE TEMPERATURE RECORD ', & ' EOF READ WHILE READING TEMPERATURE RECORD'/) END MODULE lease_module PROGRAM lease !### !*********************************************************************** ! * ! Program Name: lease * ! * ! File Name: lease.f90 * ! * ! Purpose: Program to calculate thermodynamic data required for * ! input into a program which computes chemical * ! equilibrium compositions (writeen by Bonnie J. McBride * ! and Sanford Gordon at Lewis Research Center, Cleveland,* ! Ohio). This program will either read heat capacity * ! values or compute the values. Written by Robert L. * ! ***** under DR. John Holloway, at Arizona State * ! University. * ! * ! Author: Robert L. ***** * ! * ! Date: JUNE 1, 1973 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Compatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/73 0 Initial Release * ! R.L. ***** 06/01/03 0 Converted to FORTRAN 95* ! * !----------------------------------------------------------------------* ! * ! Execution Command: * ! * ! lease Where: * ! * ! No Parameters * ! * ! Program will request name of input parameter file and output * ! listing file. * ! * !----------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! --------------- * ! * ! an - Real variable used to store the value of the 6th * ! coefficient minus enthalpy divided by r * ! bypass_code - Logical variable used to by pass code * ! ccpr - Real variable containing cp/r calculated from the * ! coefficients * ! chhrt - Real variable containing ht-h298/rt calculated * ! from the coefficients * ! cpre - Real variable containing percentage of error * ! between calculated cp/r and cp/r used in the * ! calculations * ! csr - Real variable containing s/r calculated from the * ! coefficients * ! day - Integer value of current day of month * ! dh - Real variable containing enthalpy value of specie * ! at 298 degrees Kelvin * ! exit_loop - Logical variable used to exit a loop * ! file_error - Logical variable used to denote a file open error * ! hhrte - Real variable containing percentage of error * ! between calculated ht-h298/rt and ht-h298/rt used * ! in the calculations * ! hoo - Real variable containing enthalpy at zero degree * ! in Kelvin * ! hoom - Real variable containing hoo converted to calories* ! hour - Integer value of current hour of day * ! hundred - Integer value of current hundredth second * ! i - Integer varaible used as do loop index * ! ii - Integer varaible used as do loop index * ! j - Integer varaible used as do loop index * ! jj - Integer varaible used as do loop index * ! k - Integer varaible used as temporary storage of kk * ! in calculation of matrix elements. Also do loop * ! limit in matrix reflection and solution * ! kount - Integer varaible used to count number of lines * ! printed. Used for page control * ! l - Integer varaible used as limit of the set of input* ! values being considered * ! ll - Integer varaible used as limit of the set of input* ! values being considered * ! lt - Real variable equal to log T of T being considered* ! main_loop - Logical variable used to control looping for the * ! main loop * ! minute - Integer value of current minute of day * ! month - Integer value of current month of the year * ! mm - Integer varaible used to determine which set of * ! coefficients are being stored in ans array * ! n - Integer varaible used as do loop index * ! not_open - Integer varaible used to return status of opening * ! a file * ! np - Integer varaible used as temporary storage of npts* ! pts - Real variable containing npts converted to real * ! rcp - Logical variable used to test if heat capacity * ! should be calculated from ht - h298 values * ! rt - Real variable used as temporary storage of r*t * ! second - Integer value of current second of day * ! skip_code - Logical variable used to test for skipping code * ! sre - Real variable used as percentage of error between * ! st/r calculated from the coefficiens and st/r used* ! in the calculations * ! status - Integer value of status of current read * ! t - Real variable containing T value read for each * ! specie * ! to - Real variable containing T of transition point * ! year - Integer value of current year * ! * ! * ! Global Variables: * ! ---------------- * ! * ! lease_module * ! a - Real array containing elements of the matrix * ! ans - Real array containing coefficients derived from * ! matrix solution * ! av - Logical variable used to test if values should be * ! constrained to be equal at 1000 degrees Kelvin * ! b - Real array containing stoichiometric coefficients * ! of elements of species * ! c - Real array containing answer region of the matrix * ! solution * ! cp - Real array containing heat capacity of species in * ! calories * ! cpr - Real array containing cp/r calculated from the * ! coefficients * ! date1 - Character*3 containing reference code and date of * ! data used to calculate the coefficients * ! date2 - Character*3 containing reference code and date of * ! data used to calculate the coefficients * ! hhrt - Real array containing heat content of species * ! divided by rt * ! hth - Real array containing heat content of specie in * ! Kcal * ! htho - Real array containing heat content converted to * ! calories * ! hthoe - Real array containing percentage of error between * ! heat content calculated from heat capacity and * ! heat content input * ! infile - Character*80 variable containing input file name * ! ispeci - Character*4 variable containing alphameric name of* ! specie being considered * ! kk - Integer array used as position of transition T * ! lp - Integer variable used to locate last point of * ! input values being considered * ! m - Integer value used to determine which set of input* ! values is being considered * ! mp - Integer value used to locate minimum point of * ! input values being considered * ! mt - Character*2 array containing alphameric name of * ! element of species being considered * ! nd - Integer variable containing number of points to * ! consider in averaging corresponding values in * ! calculating heat capacity * ! ndeg - Integer value of number of degrees for polynomial * ! fit * ! npts - Integer array containing number of points before * ! and/or after transition T * ! outfile - Character*80 variable containing output file name * ! phaz - Character*1 variable containing phase of specie * ! being considered * ! sr - Real array used as temporary storage of st/r * ! st - Real array containing entropy values of species at* ! T in calories * ! stop_it - Logical flag used to test if stop condition met * ! t1 - Real value of inital T of specie being considered * ! t2 - Real value of final T of specie being considered * ! tt - Real array containing T being considered * ! vol - Real variable containing molar volume of specie * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! ** Note - Some magic numbers left in to ease examination of formulas * ! lease_module * ! av_kount - Integer value of page line count after ave * ! listing * ! date_error - Integer value of error message for date error * ! high_coeff - Integer index to high set of coefficients * ! init_kount - Integer initial value of page line count * ! inp - Integer value of logical input unit = 5 * ! kount_inc - Integer value of kount increment and new init * ! value * ! last - Character*4 parameter set to 'LAST' * ! low_coeff - Integer index to low set of coefficients * ! matrix_size - Integer value of number of rows and columns in * ! matrix * ! max_kount - Integer value of maximum lines per page * ! newpage - Character containing new page control character * ! npts_error - Integer value of error message for npts error * ! out - Integer value of logical output unit = 6 * ! r - Real parameter containing universal gas constant* ! specie_error - Integer value of error message for specie error * ! temp_conv - Real parameter containing temperature conversion* ! divisor * ! temp_break - Real value of temperature break over point for * ! sets of coefficients * ! temp_error - Integer value of error message for temp error * ! * ! Modules: * ! -------- * ! * ! lease_module - Contains constants and globals for lease * ! read_error_module - Contains interface for read_error * ! * ! External Software: * ! ----------------- * ! * ! alog - FORTRAN 95 math library routine for log value * ! eqfind - Computes the coefficients of a second degree * ! polynomial * ! getdat - FORTRAN 95 routine to get current date * ! gettim - FORTRAN 95 routine to get current time * ! punch - Outputs species coefficient records for input into* ! a chemical equilibrium program * ! read_error - Checks status of last read and outputs error * ! error message if error occurred * ! * ! Files: * ! ----- * ! * ! Input Files: * ! * ! infile - Input parameter file * ! * ! Intermediate Files: * ! * ! None * ! * ! Output Files: * ! * ! outfile - Output listing file * ! * !*********************************************************************** !### USE lease_module USE read_error_module IMPLICIT none LOGICAL :: bypass_code LOGICAL :: exit_loop LOGICAL :: file_error LOGICAL :: main_loop LOGICAL :: rcp LOGICAL :: skip_code INTEGER :: i INTEGER :: ii INTEGER :: j INTEGER :: jj INTEGER :: k INTEGER :: kount INTEGER :: l INTEGER :: ll INTEGER :: mm INTEGER :: n INTEGER :: np INTEGER :: not_open INTEGER :: status INTEGER :: day INTEGER :: hour INTEGER :: hundred INTEGER :: minute INTEGER :: month INTEGER :: second INTEGER :: year REAL :: an REAL :: ccpr REAL :: chhrt REAL :: cpre REAL :: csr REAL :: dh REAL :: hhrte REAL :: hoom REAL :: hoo REAL :: lt REAL :: pts REAL :: rt REAL :: sre REAL :: t REAL :: to a = 0.0 stop_it = .false. WRITE(*,fmt = "(a)",advance = 'NO')' Input name of Input file: ' READ(*,*) infile WRITE(*,fmt = "(a,a)")' INPUT FILE : ',infile WRITE(*,fmt = "(a)",advance = 'NO')' Input name of Output file: ' READ(*,*) outfile WRITE(*,fmt = "(a,a)") ' OUTPUT FILE : ',outfile file_error = .false. main_loop = .true. OPEN(unit = out,file = outfile,status = 'unknown',iostat = not_open) lease_01 : IF (not_open /= 0) THEN WRITE(*,fmt = "(a)") ' UNABLE TO OPEN LISTING FILE' WRITE(*,fmt = "(a)") ' PROGRAM TERMINATED' main_loop = .false. file_error = .true. ELSE lease_01 REWIND out OPEN(unit = inp,file = infile,status = 'old',iostat = not_open) lease_02 : IF (not_open /= 0) THEN WRITE(*,fmt = "(a)") ' UNABLE TO OPEN INPUT FILE' WRITE(*,fmt = "(a)") ' PROGRAM TERMINATED' main_loop = .false. file_error = .true. END IF lease_02 END IF lease_01 lease_03 : IF(.not.file_error)THEN REWIND inp END IF lease_03 CALL getdat( & year, & month,day) CALL gettim( & hour, & minute, & second, & hundred) WRITE(unit = out,fmt = "(a,i2.2,a,i2.2,a,i4)")& ' DATE = ',month,'/',day,'/',year WRITE(unit = out,fmt = "(a,i2.2,a,i2.2)")& ' TIME = ',hour,':',minute WRITE(unit = out,fmt = "(a)")" " WRITE(unit = out,fmt="(a,a)")" Input File = ",infile WRITE(unit = out,fmt="(a,a)")" Output File = ",outfile WRITE(unit = out,fmt = "(a)")" " lease_04 : DO WHILE(main_loop) status = 0 tt = 0.0 rcp = .true. av = .false. m = 0 kount = init_kount lp = 0 READ(unit = inp,fmt = "(8a4)",iostat = status)(ispeci(i),i = 1,8) CALL read_error( & unit = inp, & message = specie_error, & status = status) lease_05 : IF(ispeci(1) == last.or.status /= 0)THEN EXIT lease_04 END IF lease_05 BACKSPACE(unit = inp) READ(unit = inp,fmt = "(8a4,l1,2x,l1,4f8.0)",iostat = status) & (ispeci(i),i = 1,8),av,rcp,hoom,dh,t1,t2 CALL read_error( & unit = inp, & message = specie_error, & status = status) lease_06 : IF(status /=0)THEN EXIT lease_04 END IF lease_06 WRITE(unit = out,fmt = "(a,2x,8a4,/,/)") & 'THE SPECIE BEING CONSIDERED IS',(ispeci(i),i = 1,8) READ(unit = inp,fmt = "(2a3,2x,a1,10(a2,f3.0),e15.0)",iostat = status) & date1,date2,phaz,(mt(i),b(i),i = 1,10),vol CALL read_error( & unit = inp, & message = date_error, & status = status) lease_07 : IF(status /=0)THEN EXIT lease_04 END IF lease_07 lease_08 : DO bypass_code = .false. skip_code = .false. m = m+1 READ(unit = inp,fmt = "(i4,f8.0,i4,i4,i4)",iostat = status) & npts(m),to,kk(m),nd,ndeg CALL read_error( & unit = inp, & message = npts_error, & status = status) lease_09 : IF(status /= 0)THEN EXIT lease_08 END IF lease_09 WRITE(unit = out,fmt = "(a,6x,a,8x,a,5x,a,7x,a,2x,a,2x,a,3x,a,3x,a,/)") & ' NPTS','TO','K','HOO','RCP','DEG','ND','DELTA H','AV' WRITE(unit = out,fmt = "(i4,f12.3,i6,f12.3,4x,l1,i6,i5,f10.3,4x,l1,/)") & npts(m),to,kk(m),hoom,rcp,ndeg,nd,dh,av mp = lp+1 lp = lp+npts(m) ! ! If rcp =.true., cp is computed. If rcp=.false., cp is read in. ! lease_10 : IF(rcp)THEN bypass_code = .true. ELSE lease_10 READ(unit = inp,fmt = "(f6.0,f9.0,f9.0,f9.0)",iostat = status) & (tt(n),cp(n),st(n),hth(n),n = mp,lp) CALL read_error( & unit = inp, & message = temp_error, & status = status) lease_11 : IF(status /= 0)THEN EXIT lease_08 END IF lease_11 lease_12 : IF(mp > 1) THEN EXIT lease_08 END IF lease_12 END IF lease_10 lease_13 : IF(.not.bypass_code) THEN lease_14 : IF(av)THEN skip_code = .true. ELSE lease_14 EXIT lease_08 END IF lease_14 END IF lease_13 lease_15 : IF(.not.skip_code)THEN np = npts(m) READ(unit = inp,fmt = "(f6.0,f8.0,f8.0)",iostat = status) & (tt(n),hth(n),st(n),n = 1,np) CALL read_error( & unit = inp, & message = temp_error, & status = status) lease_16 : IF(status /= 0)THEN EXIT lease_08 END IF lease_16 CALL eqfind ! ! Exit main_loop i stop condition met in enxen ! lease_17 : IF(stop_it)THEN EXIT lease_04 END IF lease_17 lease_18 : IF(mp > 1)THEN EXIT lease_08 END IF lease_18 END IF lease_15 END DO lease_08 lease_19 : IF(status /= 0) THEN EXIT lease_04 END IF lease_19 hoo = hoom*temp_conv dh = dh*temp_conv lease_20 : DO i = 1,lp htho(i) = (hth(i)*temp_conv)-hoo END DO lease_20 lease_21 : IF(av) THEN kount = av_kount END IF lease_21 ! ! Clear the matrix region ! lease_22 : IF(m > 1) THEN m = m-1 END IF lease_22 lease_23 : DO a = 0.0 ! ! Set up matrix elements for diagonal and above for first 10 rows ! pts = npts(m) skip_code = .false. lease_24 : IF(av)THEN skip_code = .true. l = mp ll = lp lease_25 : IF(m == 1) THEN ll = mp-1 END IF lease_25 lease_26 : IF(m == 1) THEN l = mp-npts(1) END IF lease_26 k = kk(m)+(l-1) END IF lease_24 lease_27 : IF(.not.skip_code)THEN l = 1 k = kk(m) ll = npts(1) END IF lease_27 lease_28 : DO n = l,ll cpr(n) = cp(n)/r t = tt(n) rt = r * t hhrt(n) = htho(n)/rt sr(n) = st(n)/r lt = alog(t) a(1,1) = a(1,1) + 2. + lt*lt a(1,2) = a(1,2) + (3./2.+lt)*t a(2,2) = a(2,2)+t*t a(1,3) = a(1,3)+(4./3.+(1./2.)*lt)*t*t a(2,3) = a(2,3)+t**3 a(3,3) = a(3,3)+t**4 a(1,4) = a(1,4)+(5./4.+(1./3.)*lt)*t**3 a(2,4) = a(2,4)+t**4 a(3,4) = a(3,4)+t**5 a(4,4) = a(4,4)+t**6 a(1,5) = a(1,5)+(6./5.+(1./4.)*lt)*t**4 a(2,5) = a(2,5)+t**5 a(3,5) = a(3,5)+t**6 a(4,5) = a(4,5)+t**7 a(5,5) = a(5,5)+t**8 a(1,6) = a(1,6)+ 1./t a(3,6) = a(3,6)+t a(4,6) = a(4,6)+t*t a(5,6) = a(5,6)+t**3 a(6,6) = a(6,6) +1./(t*t) a(1,7) = a(1,7) + lt a(2,7) = a(2,7) + t a(3,7) = a(3,7) + t*t a(4,7) = a(4,7) + t**3 a(5,7) = a(5,7) + t**4 a(1,11)= a(1,11)+ cpr(n)+hhrt(n)+(sr(n)*lt) a(2,11)= a(2,11)+ t * (cpr(n)+(1./2.)*hhrt(n) +sr(n)) a(3,11)= a(3,11)+ t*t*(cpr(n)+(1./3.)*hhrt(n)+(1./2.)*sr(n)) a(4,11)= a(4,11)+ t**3 * (cpr(n)+ (1./4.)*hhrt(n)+ (1./3.)*sr(n)) a(5,11)= a(5,11)+ t**4 * (cpr(n)+ (1./5.) * hhrt(n)+(1./4.)*sr(n)) a(6,11)= a(6,11)+(1./t)*hhrt(n) a(7,11)= a(7,11)+sr(n) END DO lease_28 a(2,2) = a(2,2)*(9./4.) a(2,3) = a(2,3)*(5./3.) a(3,3) = a(3,3)*(49./36.) a(2,4) = a(2,4)*(35./24.) a(3,4) = a(3,4)*(5./4.) a(4,4) = a(4,4)*(169./144.) a(2,5) = a(2,5)*(27./20.) a(3,5) = a(3,5)*(143./120.) a(4,5) = a(4,5)*(17./15.) a(5,5) = a(5,5)*(441./400.) a(2,6) = pts/2. a(3,6) = a(3,6)* (1./3.) a(4,6) = a(4,6)* (1./4.) a(5,6) = a(5,6)* (1./5.) a(3,7) = a(3,7)* (1./2.) a(4,7) = a(4,7)* (1./3.) a(5,7) = a(5,7)* (1./4.) a(7,7) = pts a(1,8) = 1. a(2,8) = to a(3,8) = to * to a(4,8) = to**3 a(5,8) = to**4 a(1,9) = 1. a(2,9) = to/2. a(3,9) = (to*to)/3. a(4,9) = (to**3)/4. a(5,9) = (to**4)/5. a(6,9) = 1./to a(1,10) = alog(to) a(2,10) = to a(3,10) = (to*to)/2. a(4,10) = to**3/3. a(5,10) = to**4/4. a(7,10) = 1. a(8,11) = cp(k)/r a(9,11) = htho(k)/(r*tt(k)) a(10,11)= st(k)/r ! ! Compute the matrix by reflecting symetrical elements above the ! diagonal ! k = 2 lease_29 : DO i = 1,(matrix_size-1) lease_30 : DO j = k,matrix_size a(j,i) = a(i,j) END DO lease_30 k = k+1 END DO lease_29 ! ! Solve the matrix ! n = matrix_size c = 0.0 lease_31 : DO i = 1,matrix_size lease_32 : DO j = i,matrix_size a(i,j+1) = a(i,j+1)/a(i,i) exit_loop = .false. lease_33 : IF ((i-n) == 0)THEN exit_loop = .true. EXIT lease_32 END IF lease_33 END DO lease_32 lease_34 : IF(exit_loop)THEN EXIT lease_31 END IF lease_34 k = i+1 lease_35 : DO ii = k,matrix_size lease_35a : DO jj = i,matrix_size a(ii,jj+1) = -a(ii,i)*a(i,jj+1)+a(ii,jj+1) END DO lease_35a END DO lease_35 END DO lease_31 c(n) = a(i,j+1) j = (matrix_size - 1) lease_36 : DO i = 1,(matrix_size -1) k = j+1 lease_37 : DO mm = 1,i c(j) = c(j)+c(k)*a(j,k) k = k+1 END DO lease_37 c(j) = a(j,k)-c(j) j = j-1 END DO lease_36 ! ! Store the coefficients in ans(i,j) array. ! c(6) = c(6)+(dh/r) mm = high_coeff lease_38 : IF(tt(l) < temp_break) THEN mm = low_coeff END IF lease_38 lease_39 : DO i = 1,7 ans(mm,i) = c(i) END DO lease_39 lease_40 : IF(l > 1) THEN EXIT lease_23 END IF lease_40 m = m+ 1 lease_41 : IF(.not.av)THEN EXIT lease_23 END IF lease_41 END DO lease_23 lease_42 : IF(status /= 0)THEN EXIT lease_04 END IF lease_42 ! ! The coefficients for the species have been computed, call ! subroutine punch to output records for the species being considered. ! CALL punch skip_code = .false. lease_43 : IF(rcp)THEN skip_code = .true. WRITE(unit = out,fmt = "(/,4x,a,5x,a,11x,a,8x,a,5x,a,8x,a,/)") & 'TEMP.','C','S','H-H298','H-H298E','H-HO' lease_44 : DO n = 1,lp kount = kount+1 WRITE(unit = out,fmt = "(f9.2,f8.3,4f12.3)") tt(n),cp(n),st(n), & hth(n),hthoe(n),htho(n) lease_45 : IF(kount >= max_kount)THEN WRITE(unit = out,fmt = "(a)") newpage kount = 0 END IF lease_45 END DO lease_44 END IF lease_43 lease_46 : IF(.not.skip_code)THEN WRITE(unit = out,fmt = "(/,4x,a,5x,a,11x,a,8x,a,7x,a,/)") & 'TEMP.','C','S','H-H298','H-HO' lease_47 : DO n = 1,lp kount = kount+1 WRITE(unit = out,fmt = "(f9.2,f8.3,3f12.3)") tt(n),cp(n),st(n),& hth(n),htho(n) lease_48 : IF(kount >= max_kount)THEN WRITE(unit = out,fmt = "(a)")newpage kount = 0 END IF lease_48 END DO lease_47 END IF lease_46 kount = kount+ kount_inc lease_49 : IF(kount >= max_kount) THEN WRITE(unit = out,fmt = "(a)")newpage kount = kount_inc END IF lease_49 WRITE(unit = out,fmt = "(/,4x,a,7x,a,7x,a,6x,a,6x,a,5x,a,5x,a,6x,a,7x,a,7x,a,/)") & 'TEMP.','CPR','CCPR','CPRE','HHRT','CHHRT','HHRTE','SR','CSR','SRE' ! ! Calculate from the least square coefficients values of cp/r,H-HO/ ! rt,sr, and the errors in these at each temperature. ! lease_50 : DO j = 1,lp skip_code = .false. kount = kount + 1 n = low_coeff lease_51 : IF(j > npts(1))THEN n = high_coeff END IF lease_51 lease_52 : IF(tt(1) >= temp_break)THEN n = high_coeff END IF lease_52 ccpr = ans(n,1) + ans(n,2)* tt(j)+ans(n,3)* tt(j)**2+ans(n,4)*tt(j)**3 & + ans(n,5)*tt(j)**4 an = ans(n,6)-(dh/r) chhrt = ans(n,1)+ans(n,2)*(tt(j)/2.)+ ans(n,3)*(tt(j)**2/3.) & +ans(n,4)*(tt(j)**3/4.) +ans(n,5)*(tt(j)**4/5.)+an/tt(j) csr = ans(n,1)*alog(tt(j))+ans(n,2)*tt(j)+ans(n,3)*(tt(j)**2/2.) & + ans(n,4)*(tt(j)**3/3.)+ans(n,5)*(tt(j)**4/4.)+ans(n,7) cpre = ((cpr(j)-ccpr)/cpr(j))*100. sre = ((sr(j)-csr)/sr(j))*100. lease_53 : IF(hhrt(j) /= 0.0)THEN hhrte = ((hhrt(j)-chhrt)/hhrt(j))*100. skip_code = .true. END IF lease_53 lease_54 : IF(.not.skip_code)THEN hhrte = 0. END IF lease_54 WRITE(unit = out,fmt = "(f9.2,1x,9f10.4)")tt(j),cpr(j),ccpr,cpre, & hhrt(j),chhrt,hhrte,sr(j),csr,sre lease_55 : IF(kount >= max_kount)THEN WRITE(unit = out,fmt = "(a)")newpage kount = 0 END IF lease_55 END DO lease_50 END DO lease_04 END PROGRAM lease SUBROUTINE ave !### !***************************************************************************** ! * ! Subroutine Name: ave * ! * ! File Name: ave.f90 * ! * ! Purpose: A subroutine to find the average heat capacity of the values * ! generated in a program to find the coefficients of a second * ! degree plynomial * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 1973 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Compatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/73 0 Initial Release * ! R.L. ***** 06/01/03 A Converted to FORTRAN 95 * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! CALL ave where: * ! * ! No arguments * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! aa - Real variable containing number of cp and ht-h298 values for * ! each T * ! acp - Real variable containing average heat capacity value at 1000 * ! degrees Kelvin * ! cpp - Real array used as temporary storage for heat capacity * ! hc - Real array used as temporary storage for calculated heat * ! content * ! he - Real array used as temporary storage for heat content * ! computation error * ! ht - Real array used as temporary storae for input heat content * ! i - Integer variable used as do loop index * ! ii - Integer variable used as counter to locate heat capacity and * ! heat content values being averaged * ! j - Integer variable used as do loop index * ! jc - Integer variable used as counter to locate heat capacity and * ! heat content values being averaged * ! jj - Integer variable used as lower column number of values being * ! averaged * ! k - Integer variable used as upper column numbeer of values being * ! averaged * ! kkk - Integer variable used as position of values being stored * ! mdeg - Integer variable used as cross over point where number of rows* ! begins to decrease * ! np - Integer variable used as upper limit of do loop (=npts) * ! s - Real array used as temporary storage for entropy * ! t - Real array used as temporary storage for temperatures * ! * ! Global Variables: * ! ----------------- * ! lease_module * ! av - Logical variable used to determine if points are to be * ! constrained to fit at 1000 degrees Kelvin * ! cp - Real array of computed heat capacity * ! cpsum - Real array containing summation of heat capacity for nd * ! points * ! hcalc - Real array containing computed heat content * ! hsum - Real array containing summation of computed heat content for* ! nd points * ! hth - Real array containing heat content from input records * ! hthoe - Real array containing percentage of error between computed * ! heat content and heat content read from input records * ! lp - Integer variable used as upper limit of values being stored * ! m - Integer variable used as position of npts value being * ! considered * ! mp - Integer variable used as location of minimum value of points* ! being considered * ! nd - Integer variable used to store current npts value * ! npts - Integer variable used as number of values below and above * ! 1000 degrees Kelvin * ! st - Real array containing entropy read from input records * ! tcp - Real array used as temporary storage for heat capacity * ! thcalc - Real array used as temporary storage for calculated heat * ! content * ! tt - Real array used to store t from input records * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! lease_module * ! av_value - Real value used to compute averager value * ! percent - Real parameter containing percentage conversion value * ! save_size - Integer value of size of saved arrays in subroutine ave * ! temp_conv - Real parameter containing temperature conversion divisior * ! * ! Modules: * ! -------- * ! * ! lease_module - Module containing globals and constants for lease * ! * ! External Software: * ! ----------------- * ! * ! None * ! * ! Files: * ! ----- * ! * ! Input Files: * ! * ! None * ! * ! Intermediate Files: * ! * ! None * ! * ! Output Files: * ! * ! None * ! * !***************************************************************************** !### USE lease_module IMPLICIT none INTEGER :: i INTEGER :: ii INTEGER :: j INTEGER :: jc INTEGER :: jj INTEGER :: k INTEGER :: kkk INTEGER :: mdeg INTEGER :: np REAL :: aa REAL :: acp REAL,SAVE,DIMENSION(save_size) :: cpp REAL,SAVE,DIMENSION(save_size) :: hc REAL,SAVE,DIMENSION(save_size) :: he REAL,SAVE,DIMENSION(save_size) :: ht REAL,SAVE,DIMENSION(save_size) :: s REAL,SAVE,DIMENSION(save_size) :: t mdeg = npts(m) - (nd-1) np = npts(m) hsum = 0.0 cpsum = 0.0 ! ! Sum the rows and compute the average ! jj = 1 k = 1 ave_01 : DO i = 1,np ii = 1 jc = i ave_02 : IF(i > mdeg) THEN jc = mdeg ii = i - (mdeg-1) END IF ave_02 ave_03 : DO j = jj,k cpsum(i) = cpsum(i) + tcp(ii,jc) hsum(i) = hsum(i) + thcalc(ii,jc) ii = ii + 1 jc = jc - 1 END DO ave_03 aa = k - (jj-1) ave_04 : IF(aa > 0.0)THEN cp(i) = cpsum(i)/aa hcalc(i) = hsum(i)/aa END IF ave_04 ave_05 : IF(i >= nd) THEN jj = jj + 1 END IF ave_05 ave_06 : IF(mdeg >= i)THEN ave_07 : IF(k < mdeg)THEN k = k + 1 END IF ave_07 END IF ave_06 END DO ave_01 ! ! Compute the relative error of ht-ho ! ave_08 : DO i = 1,np ave_09 : IF(hth(i) /= 0.0) THEN hthoe(i) = ((hth(i)-hcalc(i))/hth(i))* percent ELSE ave_09 hthoe(i) = 0.0 END IF ave_09 hth(i) = hth(i) / temp_conv END DO ave_08 ! ! Store the cp, st, delta h and delta h error values. ! ave_10 : IF(av)THEN kkk = mp - 1 ave_11 : DO i = 1,np kkk = kkk + 1 t(kkk) = tt(i) ht(kkk) = hth(i) cpp(kkk)= cp(i) hc(kkk) = hcalc(i) he(kkk) = hthoe(i) s(kkk) = st(i) j = j + 1 END DO ave_11 ave_12 : IF(mp > 1) THEN ! ! Remove the cp, st, delta h and delta h error values from storage ! and put into their arrays. ! ave_13 : DO i = 1,lp tt(i) = t(i) hth(i) = ht(i) cp(i) = cpp(i) hcalc(i) = hc(i) hthoe(i) = he(i) st(i) = s(i) END DO ave_13 ! ! The second interval has been read in. Compute the average gp at ! 1000 degrees Kelvin. ! acp = (cp(mp-1) + cp(mp))/av_value cp(mp) = acp cp(mp-1) = acp END IF ave_12 END IF ave_10 END SUBROUTINE ave SUBROUTINE enxen( & lv, & stop_flag) !### !***************************************************************************** ! * ! Subroutine Name: enxen * ! * ! File Name: enxen.f90 * ! * ! Purpose: M R BOLDUC subroutine to solve an n x n system of simultaneous* ! linear equations by the process of diagonalization. * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 1973 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Compatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/73 0 Initial Release * ! R.L. ***** 06/01/03 A Converted to FORTRAN 95 * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! CALL enxen(lv,stop_flag) where: * ! * ! lv - Integer variable used to transfer value of ndeg from eqfind * ! stop_flag - Logical variable used to determine if stop condition met * ! * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! i - Integer variable used as do loop index * ! ibidig - Integer variable used as control index to locate value of a * ! being considered * ! idiag - Integer variable used as control index to locate value of a * ! being considered * ! ijk - Integer variable used as control index to locate value of a * ! being considered * ! isub - Integer variable used as control index to locate value of a * ! being considered * ! j - Integer variable used as do loop index * ! k - Integer variable used as do loop index * ! l - Integer variable used to store value of ndeg transfered from * ! n - Integer variable used as do loop upper limit * ! neils - Integer variable used as counter to locate value of a being * ! considered * ! * ! Global Variables: * ! ----------------- * ! lease_module * ! a_store - Real array used to transfer value of SUM from eqfind * ! ac - Real array containing coefficients of an n degree polynomial* ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! None * ! * ! Modules: * ! -------- * ! * ! lease_module - Contains constants and global variables for lease * ! * ! External Software: * ! ----------------- * ! * ! None * ! * ! Files: * ! ----- * ! * ! Input Files: * ! * ! None * ! * ! Intermediate Files: * ! * ! None * ! * ! Output Files: * ! * ! Outputs to the screen * ! * !***************************************************************************** !### USE lease_module IMPLICIT none INTEGER :: i INTEGER :: ibidig INTEGER :: idiag INTEGER :: ijk INTEGER :: isub INTEGER :: j INTEGER :: k INTEGER :: l INTEGER :: n INTEGER :: neils INTEGER,INTENT(inout) :: lv LOGICAL,INTENT(inout) :: stop_flag l = lv n = lv-1 idiag = 1 enxen_01 : DO i = 1,n enxen_02 : IF(a_store(idiag) == 0)THEN WRITE(*,fmt = "(a)")' DIAGONAL ELEMENT EQUAL TO ZERO' stop_flag = .true. EXIT enxen_01 END IF enxen_02 enxen_03 : DO j = i,n ibidig = j*l+i a_store(ibidig) = a_store(ibidig)/a_store(idiag) END DO enxen_03 ac(i) = ac(i)/a_store(idiag) enxen_04 : DO k = 1,l enxen_05 : IF ((k-i) /=0)THEN neils = k+(i-1) * l enxen_06 : DO j = i,n isub = j*l+k ijk = j*l+i a_store(isub) = a_store(isub)-a_store(ijk)*a_store(neils) END DO enxen_06 ac(k) = ac(k)-ac(i)*a_store(neils) END IF enxen_05 END DO enxen_04 idiag = i + l * i + 1 END DO enxen_01 enxen_07 : IF(.not.stop_flag)THEN i = l*l ac(l) = ac(l)/a_store(i) enxen_08 : DO i = 1,n j = n*l+i ac(i) = ac(i)-a_store(j)*ac(l) END DO enxen_08 END IF enxen_07 END SUBROUTINE enxen SUBROUTINE eqfind !### !***************************************************************************** ! * ! Subroutine Name: eqfind * ! * ! File Name: eqfind.f90 * ! * ! Purpose: A subroutine compute the coefficients of a second degree * ! polynomial * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 1973 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Compatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/73 0 Initial Release * ! R.L. ***** 06/01/03 A Converted to FORTRAN 95 * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! CALL eqfind where: * ! * ! No arguments * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! i - Integer variable used as do loop index * ! ii - Integer variable used as upper position of ht-h298 value * ! being considered * ! j - Integer variable used as do loop index * ! jc - Integer variable used as counter to determine which t * ! interval is being considered * ! k - Integer variable used as do loop index * ! kkk - Integer variable used as lower position of ht-h298 and cp * ! value being considered * ! l - Integer variable used to locate position of ht-h298 and cp * ! value being stored * ! lb - Integer variable used as number of degrees plus one * ! ldeg - Integer variable used to store value of ndeg * ! ls - Integer variable used as number of degrees plus two * ! lv - Integer variable used to store value of number of degrees * ! mmno - Integer variable used as counter to determine which ht-h298 * ! and t is being stored in a * ! mno - Integer variable used to determine which npt is being * ! considered * ! np - Integer variable used to store value of current npts * ! npind - Integer variable used to calculate value of mmno * ! p - Real variable containing ti (i=1,lv) or (i=lb,ls) * ! x - Real variable used as multiplier factor in final computation * ! of heat capacity * ! * ! Global Variables: * ! ----------------- * ! lease_module * ! ac - Real array containing coefficients of an n degree * ! polynominal * ! a_store - Real array used to store sum_all for transfer to subroutine* ! enxen * ! hth - Real array containing heat content of specie at T * ! m - Integer value used to determine which npt value is being * ! considered * ! nd - Integer value used as number of readings being considered * ! in one pass * ! ndeg - Integer value used as number of degrees in polynominal * ! npts - Integer array containing total number of points being * ! considered * ! stop_it - Logical variable used to determine if stop condition met * ! sum_all - Real array containing summation of pi * ! tcp - Real array used as temporary storage of calculated heat * ! capacity * ! thcalc - Real array used as temporary storage of calculated heat * ! content * ! tt - Real array used as multiplier factor in final computation * ! of heat capacity * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! lease_module * ! ldeg_start- Integer value used as start of ldeg index in eqfind * ! ls_inc - Integer value to add to ls in routine eqfind * ! ls_start - Integer value used as start of ls index in eqfind * ! lv_start - Integer value used as start of lv index in eqfind * ! ndeg_start- Integer value used as start of ndeg index in eqfind * ! temp_conv - Real parameter containing temperature conversion divisor * ! * ! Modules: * ! -------- * ! * ! enxen_module - Contains interface for subroutine enxen * ! lease_module - Contains global variables and constants for lease * ! * ! External Software: * ! ----------------- * ! * ! ave - Find the average heat capacity of the values generated in a * ! program to find the coefficients of a second degree plynomial* ! enxen - M R BOLDUC subroutine to solve an n x n system of * ! simultaneous linear equations by the process of * ! diagonalization. * ! * ! Files: * ! ----- * ! * ! Input Files: * ! * ! None * ! * ! Intermediate Files: * ! * ! None * ! * ! Output Files: * ! * ! None * ! * !***************************************************************************** !### USE enxen_module USE lease_module IMPLICIT none INTEGER :: i INTEGER :: ii INTEGER :: j INTEGER :: jc INTEGER :: k INTEGER :: kkk INTEGER :: l INTEGER :: lb INTEGER :: ldeg INTEGER :: ls INTEGER :: lv INTEGER :: mno INTEGER :: mmno INTEGER :: np INTEGER :: npind REAL :: p REAL :: x lv = ndeg lb = lv + 1 ls = lv + ls_inc jc = 0 kkk = 0 ii = nd - 1 np = npts(m) eqfind_01 : DO i = 1,np hth(i) = hth(i) * temp_conv END DO eqfind_01 eqfind_02 : DO kkk = kkk+1 ii = ii + 1 eqfind_03 : IF(ii > np) THEN EXIT eqfind_02 END IF eqfind_03 eqfind_04 : DO j = ls_start,ls sum_all(1) = nd sum_all(j) = 0.0 END DO eqfind_04 eqfind_05 : DO j = 1,lv ac(j) = 0.0 END DO eqfind_05 npind = 0 eqfind_06 : DO i = kkk,ii p = 1.0 mmno = npind + i ac(1) = ac(1)+ hth(mmno) eqfind_07 : DO j = lv_start,lv p = tt(mmno)* p sum_all(j) = sum_all(j)+p ac(j) = ac(j)+hth(mmno)*p END DO eqfind_07 eqfind_08 : DO j = lb,ls p = tt(mmno) * p sum_all(j) = sum_all(j)+p END DO eqfind_08 END DO eqfind_06 eqfind_09 : DO k = 1,lv eqfind_10 : DO i = 1,lv j = i+k mno = (k-1)* lv + i a_store(mno) = sum_all(j-1) END DO eqfind_10 END DO eqfind_09 CALL enxen( & lv = lv, & stop_flag = stop_it) eqfind_11 : IF(stop_it)THEN EXIT eqfind_02 END IF eqfind_11 ldeg = ndeg + 1 jc = jc+1 ! Compute ht-h298 l = 0 eqfind_12 : DO i = kkk,ii l = l+1 thcalc(l,jc) = 0.0 eqfind_13 : DO j = ldeg_start,ldeg thcalc(l,jc) = thcalc(l,jc)+ac(j-1)*tt(i)**(j-ldeg_start) END DO eqfind_13 END DO eqfind_12 ! Compute heat capacity l = 0 eqfind_14 : DO i = kkk,ii l = l+1 tcp(l,jc) = 0.0 eqfind_15 : DO j = ndeg_start,ndeg x = j-1 tcp(l,jc) = tcp(l,jc)+ x * ac(j)* tt(i)**(j-ndeg_start) END DO eqfind_15 END DO eqfind_14 END DO eqfind_02 ! All the heat capacity values have been calculated. Call subroutine ! ave to compute the average value for each temperature. eqfind_16 : IF(.not.stop_it)THEN CALL ave END IF eqfind_16 END SUBROUTINE eqfind SUBROUTINE punch !### !***************************************************************************** ! * ! Subroutine Name: punch * ! * ! File Name: punch.f90 * ! * ! Purpose: A subroutine to output species coefficient records for input * ! into a chemical equilibrium program * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 1973 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Compatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/73 0 Initial Release * ! R.L. ***** 06/01/03 A Converted to FORTRAN 95 * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! CALL punch where: * ! * ! No arguments * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! i - Integer variable used as do loop index * ! j - Integer variable used as do loop index * ! k - Integer variable used as do loop index * ! mm - Integer variable used to locate set of coefficients to be set to * ! zero * ! * ! Global Variables: * ! ----------------- * ! lease_module * ! ans - Real array containing coefficients derived from matrix * ! solution * ! av - Logical variable used to determine if one or two sets of * ! coefficients are to be output * ! b - Real array containing stoichiometric coefficient of elements * ! in current specie * ! date1 - Character*3 variable containing reference code of source of * ! data used in the calculations * ! date2 - Character*3 variable containing reference date of source of * ! data used in the calculations * ! ispeci - Character*4 variable containing alphanumeric name of specie * ! being considered * ! mt - Character*2 array containing alphanumeric name of element in * ! current specie * ! phaz - Character*1 variable containing phase of species being * ! considered * ! t1 - Real variable containing lower limit of T range of specie * ! t2 - Real variable containing upper limit of T range of specie * ! tt - Real variable containing first T used in the calculations * ! vol - Real variable containing molar volume of the specie * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! lease_module * ! coeff_set - Integer value of number of coefficients in a set * ! high_coeff - Integer index to high set of coefficients * ! low_coeff - Integer index to low set of coefficients * ! out - Integer value of logical output unit = 6 * ! * ! Modules: * ! -------- * ! * ! lease_module - Contains global variables and constants for lease * ! * ! External Software: * ! ----------------- * ! * ! None * ! * ! Files: * ! ----- * ! * ! Input Files: * ! * ! None * ! * ! Intermediate Files: * ! * ! None * ! * ! Output Files: * ! * ! outfile - Output listing file * ! * !***************************************************************************** !### USE lease_module IMPLICIT none INTEGER :: i INTEGER :: j INTEGER :: k INTEGER :: mm punch_1 : IF(.not.av)THEN mm = high_coeff punch_2 : IF(tt(1) >= 1000.)THEN mm = low_coeff END IF punch_2 ! ! av = .false. Zero the unused coefficients. ! punch_3 : DO i = 1,coeff_set ans(mm,i) = 0.0 END DO punch_3 END IF punch_1 WRITE(unit = out,fmt = "(/,/,a,/)")' COEFFICIENT OUTPUT' WRITE(unit = out,fmt = "(8a4,6x,2a3,2x,a1,2f10.3)")(ispeci(i),i = 1,8),& date1,date2,phaz,t1,t2 WRITE(unit = out,fmt = "(10(a2,f3.0),e15.8)") (mt(j),b(j),j = 1,10),vol WRITE(unit = out,fmt = "(5e15.8)")(ans(1,j),j = 1,5) WRITE(unit = out,fmt = "(5e15.8)")(ans(1,j),j = 6,7),(ans(2,k),k = 1,3) WRITE(unit = out,fmt = "(4e15.8)")(ans(2,k),k = 4,7) END SUBROUTINE punch SUBROUTINE getdat(iyr, & imon, & iday) ! Temporary fix for Salford Compiler which does not have getdat IMPLICIT none INTEGER :: iday INTEGER :: imon INTEGER :: iyr INTEGER,PARAMETER :: date_all_size = 8 INTEGER,PARAMETER :: day_index = 3 INTEGER,PARAMETER :: mon_index = 2 INTEGER,PARAMETER :: time_date_size = 8 INTEGER,PARAMETER :: time_size = 10 INTEGER,PARAMETER :: year_index = 1 INTEGER,PARAMETER :: zone_size = 5 INTEGER,DIMENSION(time_date_size) :: time_date CHARACTER(len=date_all_size) :: date_all CHARACTER(len=time_size) :: time CHARACTER(len=zone_size) :: zone CHARACTER(len=time_size) :: dummy ! !-----Get current data and time and output title line. ! CALL date_and_time(date_all,time,zone,time_date) imon = time_date(mon_index) iday = time_date(day_index) iyr = time_date(year_index) ! !-----Satisfy Lahey Compiler warning messages ! dummy(1:date_all_size) = date_all date_all = dummy(1:date_all_size) dummy(1:time_size) = time time = dummy(1:time_size) dummy(1:zone_size) = zone zone = dummy(1:zone_size) END SUBROUTINE getdat SUBROUTINE gettim(ihr, & imin, & isec, & ihun) ! Temporary fix for Salford compiler which does not have gettim IMPLICIT none INTEGER :: ihr INTEGER :: ihun INTEGER :: imin INTEGER :: isec INTEGER,PARAMETER :: date_all_size = 8 INTEGER,PARAMETER :: hour_index = 5 INTEGER,PARAMETER :: min_index = 6 INTEGER,PARAMETER :: sec_index = 7 INTEGER,PARAMETER :: tenth_index = 8 INTEGER,PARAMETER :: time_date_size = 8 INTEGER,PARAMETER :: time_size = 10 INTEGER,PARAMETER :: zone_size = 5 INTEGER,DIMENSION(time_date_size) :: time_date CHARACTER(len=date_all_size) :: date_all CHARACTER(len=time_size) :: time CHARACTER(len=zone_size) :: zone CHARACTER(len=time_size) :: dummy ! !-----Get current data and time and output title line. ! CALL date_and_time(date_all,time,zone,time_date) ihr = time_date(hour_index) imin = time_date(min_index) isec = time_date(sec_index) ihun = time_date(tenth_index) ! !-----Satisfy Lahey Compiler warning messages ! dummy(1:date_all_size) = date_all date_all = dummy(1:date_all_size) dummy(1:time_size) = time time = dummy(1:time_size) dummy(1:zone_size) = zone zone = dummy(1:zone_size) END SUBROUTINE gettim SUBROUTINE read_error(unit, & message_number, & status) !### !***************************************************************************** ! * ! Subroutine Name: read_error * ! * ! File Name: read_error.f90 * ! * ! Purpose: Checks the status of a read and if an error occurred, * ! outputs the appropiate error message. * ! * ! Author: Robert L. ***** * ! * ! Date: June 01, 2003 * ! * ! Language: Salford FORTRAN 95, Version 0.9.0 * ! * ! Hardware: IBM Compatible PC * ! * ! Operating System: Microsoft Windows 98 * ! * ! Update History: * ! * ! Name Date Revision Changes * ! ---- ---- --------- ------- * ! * ! R.L. ***** 06/01/03 0 Initial Release * ! * !----------------------------------------------------------------------------* ! * ! Usage: * ! * ! CALL read_error(unit,message_number,status) where * ! * ! unit - integer value of logical output unit * ! message_number - integer value of message to output * ! status - integer value of status returned from read * ! <0 - end of file read * ! 0 - no error * ! >0 - read error * ! * ! * !----------------------------------------------------------------------------* ! * ! Definitions: * ! * ! Local Variables: * ! ---------------- * ! * ! message_index - Integer variable used to find which error message to * ! output * ! * ! Global Variables: * ! ----------------- * ! * ! lease_module * ! re_read - Character*80 buffer used to read previous record * ! * ! Local Common Blocks: * ! ------------------- * ! * ! None * ! * ! Constants: * ! --------- * ! * ! lease_module * ! message - Chracter array containing error messages to output * ! message_offset - Integer value for offset of message number to use * ! * ! Modules: * ! -------- * ! * ! lease_module - Contains constants and globals for lease routines * ! * ! * ! External Software: * ! ----------------- * ! * ! None * ! * ! Files: * ! ----- * ! * ! Input Files: * ! * ! unit - Input logical unit where read error occurred * ! * ! Intermediate Files: * ! * ! None * ! * ! Output Files: * ! * ! Outputs to computer screen * ! * !***************************************************************************** !### USE lease_module IMPLICIT none INTEGER :: message_index INTEGER,INTENT(in) :: message_number INTEGER,INTENT(inout) :: status INTEGER,INTENT(in) :: unit read_err : IF(status >0)THEN BACKSPACE(unit) READ(unit,fmt="(a80)")re_read WRITE(*,fmt="(1x,a80)")re_read message_index = (message_number*message_offset) - 1 WRITE(*,fmt="(a)")message(message_index) ELSE IF(status < 0)THEN read_err message_index = message_number*message_offset WRITE(*,fmt="(a)")message(message_index) END IF read_err END SUBROUTINE read_error