Changed H2O condensation to new scheme. Removed unnecessary particle …

…rates in settings file.
Nicholaswogan · May 21, 2024 · 9ad07b0 · 9ad07b0
1 parent ad3e7f2
commit 9ad07b0
Show file tree

Hide file tree

Showing 14 changed files with 50 additions and 275 deletions.
diff --git a/src/atmosphere/photochem_atmosphere_rhs.f90 b/src/atmosphere/photochem_atmosphere_rhs.f90
@@ -36,7 +36,7 @@ subroutine dochem_${NAME}$(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol
     ${TYPE1}$, intent(inout) :: densities(nsp+1,nz), xp(nz), xl(nz)
     ${TYPE1}$, intent(inout) :: rhs(neqs)
 
-    real(dp) :: H2O_cold_trap, cond_rate0
+    real(dp) :: cond_rate0
     ${TYPE1}$ :: rh, df_gas_dt, df_particle_dt, cond_rate
     integer :: i, ii, j, k, kk
     type(PhotochemData), pointer :: dat
@@ -112,15 +112,21 @@ subroutine dochem_${NAME}$(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol
         if (var%temperature(j) < T_crit_H2O) then
           ! water will condense if it is below the critical point.
 
-          k = dat%LH2O + (j - 1) * dat%nq
-          H2O_cold_trap = var%H2O_condensation_rate(2)*H2O_sat_mix(j)
-          if (usol(dat%LH2O,j) >= H2O_cold_trap) then
-
-            cond_rate = damp_condensation_rate(var%H2O_condensation_rate(1), &
-                                              var%H2O_condensation_rate(2), &
-                                              var%H2O_condensation_rate(3), &
-                                              usol(dat%LH2O,j)/H2O_sat_mix(j))
-            rhs(k) = rhs(k) - cond_rate*(usol(dat%LH2O,j) - H2O_cold_trap)
+          k = dat%LH2O + (j - 1) * dat%nq ! gas phase rhs index
+
+          rh = max(usol(dat%LH2O,j)/H2O_sat_mix(j),small_real)
+
+          if (rh > var%H2O_cond_params%RHc) then
+
+            cond_rate0 = var%H2O_cond_params%k_cond*(var%edd(j)/scale_height(j)**2.0_dp)
+            cond_rate = damp_condensation_rate(cond_rate0, &
+                                               var%H2O_cond_params%RHc, &
+                                               (1.0_dp + var%H2O_cond_params%smooth_factor)*var%H2O_cond_params%RHc, &
+                                               rh)
+
+            ! Rate H2O gas is destroyed (1/s)
+            df_gas_dt = - cond_rate*usol(dat%LH2O,j)
+            rhs(k) = rhs(k) + df_gas_dt
 
           endif
         endif

diff --git a/src/atmosphere/photochem_atmosphere_utils.f90 b/src/atmosphere/photochem_atmosphere_utils.f90
@@ -100,6 +100,7 @@ module function atom_conservation(self, atom, err) result(con)
     use photochem_enum, only: VelocityDistributedFluxBC
     use photochem_eqns, only: damp_condensation_rate
     use photochem_types, only: AtomConservation
+    use photochem_const, only: small_real
     class(Atmosphere), target, intent(inout) :: self
     character(len=*), intent(in) :: atom
     character(:), allocatable, intent(out) :: err
@@ -108,7 +109,7 @@ module function atom_conservation(self, atom, err) result(con)
     real(dp) :: surf_fluxes(self%dat%nq)
     real(dp) :: top_fluxes(self%dat%nq)
     real(dp) :: integrated_rainout(self%dat%nq)
-    real(dp) :: cond_rate, con_evap_rate
+    real(dp) :: rh, df_gas_dt, cond_rate0, cond_rate, con_evap_rate
 
     integer :: ind(1), i, j, kk
 
@@ -200,14 +201,20 @@ module function atom_conservation(self, atom, err) result(con)
         i = 1
       endif
       do j = i,var%nz
-        if (wrk%usol(dat%LH2O,j) >= var%H2O_condensation_rate(2)*wrk%H2O_sat_mix(j)) then
-          cond_rate = damp_condensation_rate(var%H2O_condensation_rate(1), &
-                                             var%H2O_condensation_rate(2), &
-                                             var%H2O_condensation_rate(3), &
-                                             wrk%usol(dat%LH2O,j)/wrk%H2O_sat_mix(j))
-          con%out_other = con%out_other + cond_rate*(wrk%usol(dat%LH2O,j)  &
-                        - var%H2O_condensation_rate(2)*wrk%H2O_sat_mix(j)) &
-                        *wrk%density(j)*var%dz(j)*dat%species_composition(kk,dat%LH2O)
+
+        rh = max(wrk%usol(dat%LH2O,j)/wrk%H2O_sat_mix(j),small_real)
+
+        if (rh > var%H2O_cond_params%RHc) then
+
+          cond_rate0 = var%H2O_cond_params%k_cond*(var%edd(j)/wrk%scale_height(j)**2.0_dp)
+          cond_rate = damp_condensation_rate(cond_rate0, &
+                                             var%H2O_cond_params%RHc, &
+                                             (1.0_dp + var%H2O_cond_params%smooth_factor)*var%H2O_cond_params%RHc, &
+                                             rh)
+
+          ! Rate H2O gas is destroyed (1/s)
+          df_gas_dt = - cond_rate*wrk%usol(dat%LH2O,j)
+          con%out_other = con%out_other - df_gas_dt*wrk%density(j)*var%dz(j)*dat%species_composition(kk,dat%LH2O)
 
         endif
       enddo

diff --git a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90
@@ -48,7 +48,7 @@ subroutine dochem_${NAME}$(self, usol, rx_rates, &
     ${TYPE1}$, intent(inout) :: densities(:,:), xp(:), xl(:)
     ${TYPE1}$, intent(inout) :: rhs(:) ! neqs
 
-    real(dp) :: H2O_cold_trap, cond_rate0
+    real(dp) :: cond_rate0
     ${TYPE1}$ :: rh, dn_gas_dt, dn_particle_dt, cond_rate
     integer :: i, ii, j, k, kk
 
@@ -121,15 +121,22 @@ subroutine dochem_${NAME}$(self, usol, rx_rates, &
         if (var%temperature(j) < T_crit_H2O) then
           ! water will condense if it is below the critical point.
 
-          k = dat%LH2O + (j - 1) * dat%nq
-          H2O_cold_trap = var%H2O_condensation_rate(2)*H2O_sat_mix(j)
-          if (mix(dat%LH2O,j) >= H2O_cold_trap) then
+          k = dat%LH2O + (j - 1) * dat%nq ! gas phase rhs index
+
+          ! compute the relative humidity
+          rh = max(mix(dat%LH2O,j)/H2O_sat_mix(j),small_real)
+
+          if (rh > var%H2O_cond_params%RHc) then
+
+            cond_rate0 = var%H2O_cond_params%k_cond*(var%edd(j)/scale_height(j)**2.0_dp)
+            cond_rate = damp_condensation_rate(cond_rate0, &
+                                               var%H2O_cond_params%RHc, &
+                                               (1.0_dp + var%H2O_cond_params%smooth_factor)*var%H2O_cond_params%RHc, &
+                                               rh)
 
-            cond_rate = damp_condensation_rate(var%H2O_condensation_rate(1), &
-                                              var%H2O_condensation_rate(2), &
-                                              var%H2O_condensation_rate(3), &
-                                              mix(dat%LH2O,j)/H2O_sat_mix(j))
-            rhs(k) = rhs(k) - cond_rate*(mix(dat%LH2O,j) - H2O_cold_trap)*density(j)
+            ! Rate H2O gas is destroyed (molecules/cm^3/s)
+            dn_gas_dt = - cond_rate*usol(dat%LH2O,j)
+            rhs(k) = rhs(k) + dn_gas_dt
 
           endif
         endif

diff --git a/src/input/photochem_input_read.f90 b/src/input/photochem_input_read.f90
@@ -805,10 +805,6 @@ subroutine unpack_settings(infile, s, dat, var, err)
       endif
     endif
 
-    if (dat%water_cond) then        
-      var%H2O_condensation_rate = s%H2O_condensation_rate
-    endif
-
     !!!!!!!!!!!!!!!!!
     !!! particles !!!
     !!!!!!!!!!!!!!!!!

diff --git a/src/photochem_types.f90 b/src/photochem_types.f90
@@ -70,7 +70,6 @@ module photochem_types ! make a giant IO object
     real(dp) :: rainfall_rate
     character(s_str_len), allocatable :: rainout_species(:)
     real(dp) :: trop_alt
-    real(dp) :: H2O_condensation_rate(3)
 
     ! boundary-conditions
     type(SettingsBC), allocatable :: ubcs(:)
@@ -423,7 +422,7 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
     integer :: trop_ind !! index of troposphere (only for fix_water_in_trop == true or gas_rainout == true)
     logical :: use_manabe !! use manabe formula
     real(dp) :: relative_humidity !! relative humidity if no manabe
-    real(dp) :: H2O_condensation_rate(3) !! H2O condesation rate parameters
+    type(CondensationParameters) :: H2O_cond_params !! H2O condesation rate parameters
 
     ! radiative transfer
     real(dp), allocatable :: photon_flux(:) !! (nw) photon/cm^2/s in each wavelength bin hitting planet.

diff --git a/src/photochem_types_create.f90 b/src/photochem_types_create.f90
@@ -39,7 +39,7 @@ function unpack_PhotoSettings(root, filename, err) result(s)
 
     type(PhotoSettings) :: s
 
-    type(type_dictionary), pointer :: dict, tmp2, tmp3
+    type(type_dictionary), pointer :: dict, tmp2
     type(type_list), pointer :: list, bcs
     type(type_list_item), pointer :: item
     type(type_scalar), pointer :: scalar
@@ -258,24 +258,6 @@ function unpack_PhotoSettings(root, filename, err) result(s)
 
     endif
 
-    if (s%water_cond) then        
-
-      tmp3 => tmp2%get_dictionary('condensation-rate',.true.,error = io_err)
-      if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-
-      s%H2O_condensation_rate(1) = tmp3%get_real('A',error = io_err)
-      if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-      s%H2O_condensation_rate(2) = tmp3%get_real('rhc',error = io_err)
-      if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-      s%H2O_condensation_rate(3) = tmp3%get_real('rh0',error = io_err)
-      if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-      if (s%H2O_condensation_rate(3) <= s%H2O_condensation_rate(2)) then
-        err = 'IOError: Rate constant "rh0" for H2O condensation must be > "rhc". See '//trim(filename)
-        return
-      endif
-
-    endif
-
     !!!!!!!!!!!!!!!!!!!!!!!!!!!
     !!! boundary-conditions !!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!

diff --git a/templates/Hadean/settings_Hadean.yaml b/templates/Hadean/settings_Hadean.yaml
@@ -26,33 +26,6 @@ planet:
     rainfall-rate: 1
     tropopause-altitude: 1.1e6
     water-condensation: true
-    condensation-rate: {A: 1.0e-5, rhc: 0.01, rh0: 0.015}
-
-particles:
-- name: H2SO4aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S2aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S8aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: HCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: HCCCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: CH3CNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H2aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H6aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C4H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: N2Oaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: NH3aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}  
 
 boundary-conditions:
 - name: H2O

diff --git a/templates/Jupiter/settings_Jupiter.yaml b/templates/Jupiter/settings_Jupiter.yaml
@@ -25,33 +25,6 @@ planet:
     fix-water-in-troposphere: false
     gas-rainout: false
     water-condensation: true
-    condensation-rate: {A: 1.0e-8, rhc: 1, rh0: 1.05}
-
-particles:
-- name: H2SO4aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S2aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S8aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: HCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: HCCCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: CH3CNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H2aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H6aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C4H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: N2Oaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: NH3aer
-  condensation-rate: {A: 1.0e-7, rhc: 1.0, rh0: 1.05}
 
 boundary-conditions:
 - name: O1D

diff --git a/templates/Mars/settings_Mars.yaml b/templates/Mars/settings_Mars.yaml
@@ -26,33 +26,6 @@ planet:
     gas-rainout: false
     tropopause-altitude: 1.5e6
     water-condensation: true
-    condensation-rate: {A: 1.0e-05, rhc: 0.01, rh0: 0.0105}
-
-particles:
-- name: H2SO4aer
-  condensation-rate: {A: 1.0e-05, rhc: 1.0, rh0: 1.05}
-- name: S2aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S8aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: HCNaer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: HCCCNaer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: CH3CNaer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: C2H2aer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: C2H4aer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: C2H6aer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: C4H4aer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: N2Oaer
-  condensation-rate: {A: 1.0e-08, rhc: 1.0, rh0: 1.05}
-- name: NH3aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}  
 
 boundary-conditions:
 - name: H2O

diff --git a/templates/ModernEarth/settings_ModernEarth.yaml b/templates/ModernEarth/settings_ModernEarth.yaml
@@ -36,40 +36,7 @@ planet:
     tropopause-altitude: 1.1e6 # cm. required if gas-rainout or fix-water-in-troposphere
     # If true, then water will condense.
     water-condensation: true
-    # This is the rate of H2O condensation, and the conditions
-    # where condensation occurs. relative humidity >`rhc`, then condensation will begin to 
-    # occur. `rh0` is a value bigger than `rhc` which helps to smooth out the condensation rate,
-    # preventing discontinuities
-    condensation-rate: {A: 1.0e-5, rhc: 0.01, rh0: 0.015}
 
-# Specifies the condensation rates and RH thresholds for all particles 
-# which form from condensation.
-particles:
-- name: H2SO4aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S2aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S8aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: HCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: HCCCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: CH3CNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H2aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H6aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C4H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: N2Oaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: NH3aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-
 # Specifies boundary conditions. If a species is not specified, then
 # the model assumes zero-flux boundary conditions at the top and
 # bottom of the atmosphere

diff --git a/templates/Saturn/settings_Saturn.yaml b/templates/Saturn/settings_Saturn.yaml
@@ -25,34 +25,7 @@ planet:
     fix-water-in-troposphere: false
     gas-rainout: false
     water-condensation: true
-    condensation-rate: {A: 1.0e-8, rhc: 1, rh0: 1.05}
 
-particles:
-- name: H2SO4aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S2aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S8aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: HCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: HCCCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: CH3CNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H2aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H6aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C4H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: N2Oaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: NH3aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}  
-
 boundary-conditions:
 - name: O1D
   type: short lived

diff --git a/templates/Titan/settings_Titan.yaml b/templates/Titan/settings_Titan.yaml
@@ -24,33 +24,6 @@ planet:
     fix-water-in-troposphere: false
     gas-rainout: false
     water-condensation: true
-    condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-
-particles:
-- name: H2SO4aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S2aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: S8aer
-  condensation-rate: {A: 1.0e-5, rhc: 1.0, rh0: 1.05}
-- name: HCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: HCCCNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: CH3CNaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H2aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C2H6aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: C4H4aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: N2Oaer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}
-- name: NH3aer
-  condensation-rate: {A: 1.0e-8, rhc: 1.0, rh0: 1.05}  
 
 boundary-conditions:
 - name: CH4