From 486427078a3a8139b80552366d697dd941853b2c Mon Sep 17 00:00:00 2001
From: Jeff Curtis <jcurtis2@illinois.edu>
Date: Wed, 23 Oct 2024 13:31:52 -0500
Subject: [PATCH] add function to convert rh to water vapor mixing ratio

---
 src/env_state.F90       | 22 ++++++++++++++++++++++
 src/gas_state.F90       | 21 ++++++---------------
 src/tchem_interface.F90 |  9 +--------
 3 files changed, 29 insertions(+), 23 deletions(-)

diff --git a/src/env_state.F90 b/src/env_state.F90
index ef11f0f58..ce7dd9e28 100644
--- a/src/env_state.F90
+++ b/src/env_state.F90
@@ -197,6 +197,28 @@ real(kind=dp) function env_state_air_molar_den(env_state)
 
   end function env_state_air_molar_den
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Converts relative humidity (1) to water vapor mixing ratio (ppb).
+  !! 
+  !! Uses equation (1.10) of Seinfeld and Pandis Atmospheric Chemistry and
+  !! Physics From Air Pollution to Climate Change Second Edition, 2006.
+  real(kind=dp) function env_state_rel_humid_to_mix_rat(env_state)
+
+    !> Environment state.
+    type(env_state_t), intent(in) :: env_state
+
+    real(kind=dp), parameter :: t_steam = 373.15 ! steam temperature (K)
+    real(kind=dp) :: a, water_vp
+
+    a = 1.0 - t_steam / env_state%temp
+    a = (((-0.1299 * a - 0.6445) * a - 1.976) * a + 13.3185) * a
+    water_vp = 101325.0 * exp(a)  ! (Pa)
+    env_state_rel_humid_to_mix_rat = env_state%rel_humid * water_vp * 1.0e9 &
+         / env_state%pressure ! (ppb)
+
+  end function env_state_rel_humid_to_mix_rat
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   !> Condensation \f$A\f$ parameter.
diff --git a/src/gas_state.F90 b/src/gas_state.F90
index 3ce7eb465..9e3a69cf8 100644
--- a/src/gas_state.F90
+++ b/src/gas_state.F90
@@ -201,24 +201,15 @@ subroutine gas_state_set_camp_conc(gas_state, camp_state, gas_data)
     !> Gas data
     type(gas_data_t), intent(in) :: gas_data
 
-    real(kind=dp), parameter :: t_steam = 373.15 ! steam temperature (K)
-    real(kind=dp) :: a, water_vp
-
-    camp_state%state_var(1:size(gas_state%mix_rat)) = gas_state%mix_rat(:) &
-         / 1000.0d0
-
-    ! Convert relative humidity (1) to [H2O] (ppm)
-    ! From MOSAIC code - reference to Seinfeld & Pandis page 181
-    ! TODO Figure out how to have consistent RH<->ppm conversions
-    ! (There is only one environmental state for PartMC runs
     call assert(590005048, associated(camp_state%env_states(1)%val))
-    a = 1.0 - t_steam / camp_state%env_states(1)%val%temp
-    a = (((-0.1299 * a - 0.6445) * a - 1.976) * a + 13.3185) * a
-    water_vp = 101325.0 * exp(a)  ! (Pa)
 
+    ! Convert relative humidity (1) to [H2O] (ppb)
     camp_state%state_var(gas_data%i_camp_water) = &
-         camp_state%env_states(1)%val%rel_humid * water_vp * 1.0e6 &
-         / camp_state%env_states(1)%val%pressure ! (ppm)
+         env_state_rel_humid_to_mix_rat(camp_state%env_states(1))
+
+    ! Convert from ppb to ppm
+    camp_state%state_var(1:size(gas_state%mix_rat)) = gas_state%mix_rat(:) &
+         / 1000.0d0
 
   end subroutine gas_state_set_camp_conc
 
diff --git a/src/tchem_interface.F90 b/src/tchem_interface.F90
index 48ebc89ea..136f75856 100644
--- a/src/tchem_interface.F90
+++ b/src/tchem_interface.F90
@@ -223,8 +223,6 @@ subroutine tchem_from_partmc(aero_data, aero_state, gas_data, gas_state, &
 
     real(kind=dp), allocatable :: stateVector(:)
     integer :: stateVecDim
-    real(kind=dp), parameter :: t_steam = 373.15 ! steam temperature (K)
-    real(kind=dp) :: a, water_vp
     integer :: i_part
     integer :: i_water
 
@@ -239,12 +237,7 @@ subroutine tchem_from_partmc(aero_data, aero_state, gas_data, gas_state, &
     ! PartMC uses relative humidity and not H2O mixing ratio.
     ! Equation 1.10 from Seinfeld and Pandis - Second Edition.
     i_water = gas_data_spec_by_name(gas_data, "H2O")
-    a = 1.0 - t_steam / env_state%temp
-    a = (((-0.1299 * a - 0.6445) * a - 1.976) * a + 13.3185) * a
-    water_vp = 101325.0 * exp(a)  ! (Pa)
-    gas_state%mix_rat(i_water) = env_state%rel_humid * water_vp * 1.0e9 &
-         / env_state%pressure ! (ppb)
-
+    gas_state%mix_rat(i_water) = env_state_rel_humid_to_mix_rat(env_state)
     ! Add gas species to state vector. Convert from ppb to ppm.
     stateVector(4:gas_data_n_spec(gas_data)+3) = gas_state%mix_rat / 1000.d0