update equilibrate and nasa7

Nicholaswogan · Oct 6, 2024 · 293382f · 293382f
1 parent 97a50c1
commit 293382f
Show file tree

Hide file tree

Showing 4 changed files with 54 additions and 10 deletions.
diff --git a/src/dependencies/CMakeLists.txt b/src/dependencies/CMakeLists.txt
@@ -73,13 +73,13 @@ CPMAddPackage(
 
 CPMAddPackage(
   NAME equilibrate
-  VERSION 0.1.6
+  VERSION 0.1.8
   OPTIONS 
     "BUILD_EXECUTABLE OFF"
     "SKBUILD ${SKBUILD}"
     "BUILD_PYTHON_EQUILIBRATE ${BUILD_PYTHON_PHOTOCHEM}"
     "PYTHON_EQUILIBRATE_DESTINATION photochem"
   GITHUB_REPOSITORY "Nicholaswogan/Equilibrate"
-  GIT_TAG "v0.1.6"
+  GIT_TAG "v0.1.8"
   EXCLUDE_FROM_ALL OFF
 )
diff --git a/src/input/photochem_input_read.f90 b/src/input/photochem_input_read.f90
@@ -1245,7 +1245,7 @@ subroutine compare_rxtype_string(tmp, eqr, eqp, reverse, rxtype_int, err)
 
   subroutine get_thermodata(molecule, molecule_name, infile, &
                             thermo, err)
-    use photochem_enum, only: ShomatePolynomial, Nasa9Polynomial
+    use photochem_enum, only: ShomatePolynomial, Nasa9Polynomial, Nasa7Polynomial
     use photochem_types, only: ThermodynamicData
     class(type_dictionary), intent(in) :: molecule
     character(len=*), intent(in) :: molecule_name
@@ -1274,6 +1274,8 @@ subroutine get_thermodata(molecule, molecule_name, infile, &
       thermo%dtype = ShomatePolynomial
     elseif (model == "NASA9") then
       thermo%dtype = Nasa9Polynomial
+    elseif (model == "NASA7") then
+      thermo%dtype = Nasa7Polynomial
     else
       err = "IOError: Thermodynamic data must be in Shomate format for "//trim(molecule_name)
       return
@@ -1317,12 +1319,12 @@ subroutine get_thermodata(molecule, molecule_name, infile, &
       return
     endif
 
-    if (thermo%dtype == 1) then
-      ! Shomate
+    if (thermo%dtype == ShomatePolynomial) then
       allocate(thermo%data(7,thermo%ntemps))
-    elseif (thermo%dtype == 2) then
-      ! NASA9
+    elseif (thermo%dtype == Nasa9Polynomial) then
       allocate(thermo%data(9,thermo%ntemps))
+    elseif (thermo%dtype == Nasa7Polynomial) then
+      allocate(thermo%data(7,thermo%ntemps))
     endif
 
     k = 1

diff --git a/src/photochem_enum.f90 b/src/photochem_enum.f90
@@ -7,7 +7,8 @@ module photochem_enum
   ! dtype
   enumerator :: &
     ShomatePolynomial = 1, &
-    Nasa9Polynomial = 2
+    Nasa9Polynomial = 2, &
+    Nasa7Polynomial = 3
 
   ! falloff_type
   enumerator :: &

diff --git a/src/photochem_eqns.f90 b/src/photochem_eqns.f90
@@ -79,9 +79,48 @@ pure function gibbs_energy_nasa9(coeffs, T) result(gibbs)
 
     gibbs = enthalpy - T*entropy
   end function
+
+  pure function gibbs_energy_nasa7(coeffs, T) result(gibbs)
+    use photochem_const, only: Rgas
+    real(dp), intent(in) :: coeffs(:)
+    real(dp), intent(in) :: T
+    real(dp) :: gibbs
+
+    real(dp) :: enthalpy, entropy
+    real(dp) :: a0, a1, a2, a3, a4, a5, a6
+
+    a0 = coeffs(1)
+    a1 = coeffs(2)
+    a2 = coeffs(3)
+    a3 = coeffs(4)
+    a4 = coeffs(5)
+    a5 = coeffs(6)
+    a6 = coeffs(7)
+
+    enthalpy = &
+      a0 + &
+      a1*T/2.0_dp + &
+      a2*T**2.0_dp/3.0_dp + &
+      a3*T**3.0_dp/4.0_dp + &
+      a4*T**4.0_dp/5.0_dp + &
+      a5/T
+    enthalpy = enthalpy*Rgas*T
+
+    entropy = &
+      a0*log(T) + &
+      a1*T + &
+      a2*T**2.0_dp/2.0_dp + &
+      a3*T**3.0_dp/3.0_dp + &
+      a4*T**4.0_dp/4.0_dp + &
+      a6
+    entropy = entropy*Rgas
+
+    gibbs = enthalpy - T*entropy
+
+  end function
 
   pure subroutine gibbs_energy_eval(thermo, T, found, gibbs_energy)
-    use photochem_enum, only: ShomatePolynomial, Nasa9Polynomial
+    use photochem_enum, only: ShomatePolynomial, Nasa9Polynomial, Nasa7Polynomial
     use photochem_types, only: ThermodynamicData
 
     type(ThermodynamicData), intent(in) :: thermo
@@ -100,7 +139,9 @@ pure subroutine gibbs_energy_eval(thermo, T, found, gibbs_energy)
         if (thermo%dtype == ShomatePolynomial) then
           gibbs_energy = gibbs_energy_shomate(thermo%data(1:7,k), T)
         elseif (thermo%dtype == Nasa9Polynomial) then
-          gibbs_energy = gibbs_energy_nasa9(thermo%data(1:9,k), T)          
+          gibbs_energy = gibbs_energy_nasa9(thermo%data(1:9,k), T)
+        elseif (thermo%dtype == Nasa7Polynomial) then
+          gibbs_energy = gibbs_energy_nasa7(thermo%data(1:7,k), T)          
         endif
 
         exit