Added supercritical saturation to gases to try to create more physica…

…l behavior

src/adiabat/clima_adiabat_general.f90

@@ -352,9 +352,7 @@ subroutine make_profile(T_surf, P_i_surf, &
       ! compute saturation vapor pressures
       P_sat = huge(1.0_dp)
       if (allocated(d%sp%g(i)%sat)) then
-        if (T_surf < d%sp%g(i)%sat%T_critical) then
-          P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T_surf)
-        endif
+        P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T_surf)
       endif
 
       ! determine if species are condensing, or not
@@ -640,9 +638,7 @@ subroutine root_fcn(d, P, T, gout)
       ! compute saturation vapor pressures
       P_sat = huge(1.0_dp)
       if (allocated(d%sp%g(i)%sat)) then
-        if (T < d%sp%g(i)%sat%T_critical) then
-          P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)
-        endif
+        P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)
       endif
 
       if (d%sp_type(i) == CondensingSpeciesType) then
@@ -703,13 +699,13 @@ subroutine mixing_ratios(d, P, T, f_i_layer, f_dry)
     f_moist = 0.0_dp
     do i = 1,d%sp%ng
       if (d%sp_type(i) == CondensingSpeciesType) then
-        f_i_layer(i) = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)/P
+        f_i_layer(i) = min(d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)/P,1.0_dp)
         f_moist = f_moist + f_i_layer(i)
       endif
     enddo
 
     ! fraction of the atmosphere that is dry
-    f_dry = 1.0_dp - f_moist
+    f_dry = max(1.0_dp - f_moist, 1.0e-40_dp)
     ! mixing ratios of dry species
     do i = 1,d%sp%ng
       if (d%sp_type(i) == DrySpeciesType) then

src/clima_saturationdata.f90

@@ -20,11 +20,15 @@ module clima_saturationdata
     real(dp) :: a_v, b_v
     !> Latent heat fit parameters for sublimation (no units)
     real(dp) :: a_s, b_s
+    !> Non-physical latent heat fit parameters for super-critical gas (no units)
+    real(dp) :: a_c, b_c
 
   contains
+    procedure :: latent_heat_crit
     procedure :: latent_heat_vap
     procedure :: latent_heat_sub
     procedure :: latent_heat
+    procedure :: sat_pressure_crit
     procedure :: sat_pressure_vap
     procedure :: sat_pressure_sub
     procedure :: sat_pressure
@@ -35,6 +39,16 @@ module clima_saturationdata
 
 contains
 
+  !> Latent heat of condensation above the critical point.
+  !> This is non-physical. Its to approximate transitions in atmospheres
+  !> between super-critical and sub-critical gases.
+  function latent_heat_crit(self, T) result(L)
+    class(SaturationData), intent(inout) :: self
+    real(dp), intent(in) :: T !! K
+    real(dp) :: L !! erg/g
+    L = self%a_c + self%b_c*T
+  end function  
+
   !> Latent heat of vaporization
   function latent_heat_vap(self, T) result(L)
     class(SaturationData), intent(inout) :: self
@@ -56,13 +70,29 @@ function latent_heat(self, T) result(L)
     class(SaturationData), intent(inout) :: self
     real(dp), intent(in) :: T !! K
     real(dp) :: L !! erg/g
-    if (T > self%T_triple) then
+    if (T >= self%T_critical) then
+      L = self%latent_heat_crit(T) ! non-physical approximation
+    elseif (T > self%T_triple .and. T < self%T_critical) then
       L = self%latent_heat_vap(T)
     else ! (T <= T_triple) then
       L = self%latent_heat_sub(T)
     endif
   end function
 
+  !> Saturation pressure of a super-critical gas.
+  !> This is non-physical. Its to approximate transitions in atmospheres
+  !> between super-critical and sub-critical gases.
+  function sat_pressure_crit(self, T) result(p_sat)
+    use clima_const, only: Rgas
+    class(SaturationData), intent(inout) :: self
+    real(dp), intent(in) :: T !! K
+    real(dp) :: p_sat !! dynes/cm2
+    real(dp) :: tmp
+    tmp = (integral_fcn(self%a_v, self%b_v, self%T_critical) - integral_fcn(self%a_v, self%b_v, self%T_ref)) + &
+          (integral_fcn(self%a_c, self%b_c, T) - integral_fcn(self%a_c, self%b_c, self%T_critical))
+    p_sat = self%P_ref*exp((self%mu/Rgas)*(tmp))
+  end function
+
   !> Saturation pressure over liquid
   function sat_pressure_vap(self, T) result(p_sat)
     use clima_const, only: Rgas
@@ -91,7 +121,9 @@ function sat_pressure(self, T) result(p_sat)
     class(SaturationData), intent(inout) :: self
     real(dp), intent(in) :: T !! K
     real(dp) :: p_sat !! dynes/cm2
-    if (T > self%T_triple) then
+    if (T >= self%T_critical) then
+      p_sat = self%sat_pressure_crit(T) ! non-physical
+    elseif (T > self%T_triple .and. T < self%T_critical) then
       p_sat = self%sat_pressure_vap(T)
     else ! (T <= T_triple) then
       p_sat = self%sat_pressure_sub(T)
@@ -189,6 +221,17 @@ function create_SaturationData(s, name, filename, err) result(sat)
     sat%b_s = tmpdict%get_real('b',error = io_err)
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
 
+    ! non-physical parameters to approximate atmospheric transitions from super-critical
+    ! to sub-critical states.
+    tmpdict => s%get_dictionary("super-critical",.true.,error = io_err)
+    if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
+
+    sat%a_c = tmpdict%get_real('a',error = io_err)
+    if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
+
+    sat%b_c = tmpdict%get_real('b',error = io_err)
+    if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
+
     ! test evaluations
     block
       real(dp) :: L(3), P(3)

templates/EarlyMars/species.yaml

@@ -22,6 +22,7 @@ species:
     T-critical: 647.0
     vaporization: {a: 2.841421e+10, b: -1.399732e+07}
     sublimation: {a: 2.746884e+10, b: 4.181527e+06}
+    super-critical: {a: 1.793161e+12, b: 0.0}
   note: From the NIST database
 - name: CO2
   composition: {C: 1, O: 2}
@@ -40,6 +41,7 @@ species:
     T-critical: 304.13
     vaporization: {a: 4.656475e+09, b: -3.393595e+06}
     sublimation: {a: 6.564668e+09, b: -3.892217e+06}
+    super-critical: {a: 1.635908e+11, b: 0.0}
   note: From the NIST database
 - name: N2
   composition: {N: 2}

templates/runaway_greenhouse/species.yaml

@@ -26,6 +26,7 @@ species:
     T-critical: 647.0
     vaporization: {a: 2.841421e+10, b: -1.399732e+07}
     sublimation: {a: 2.746884e+10, b: 4.181527e+06}
+    super-critical: {a: 1.793161e+12, b: 0.0}
   note: From the NIST database
 - name: CO2
   composition: {C: 1, O: 2}
@@ -44,6 +45,7 @@ species:
     T-critical: 304.13
     vaporization: {a: 4.656475e+09, b: -3.393595e+06}
     sublimation: {a: 6.564668e+09, b: -3.892217e+06}
+    super-critical: {a: 1.635908e+11, b: 0.0}
   note: From the NIST database
 - name: N2
   composition: {N: 2}