Fixed pressure lower bc for EvoAtmosphere

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION "3.14")
 cmake_policy(SET CMP0148 OLD)
-project(Photochem LANGUAGES Fortran C VERSION "0.4.5")
+project(Photochem LANGUAGES Fortran C VERSION "0.4.6")
 
 include(FortranCInterface)
 FortranCInterface_VERIFY()

src/evoatmosphere/photochem_evoatmosphere.f90

@@ -23,11 +23,13 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     type(PhotochemWrkEvo), allocatable :: wrk
 
     logical :: evolve_climate
+    ! Below are only relevant for evolve_climate = .true.
     real(dp) :: T_surf
     real(dp) :: T_trop = 200.0_dp
     !> Callback function to set a temperature dependent albedo (e.g. ice-albedo feedback).
     procedure(temp_dependent_albedo_fcn), nopass, pointer :: albedo_fcn => null()
     type(Radtran), allocatable :: rad
+    ! Above are only relevant for evolve_climate = .true.
 
     ! Modern Earth has a pressure of 4e-7 at 100 km
     ! so it makes sense to try to keep pressure between these values
@@ -42,9 +44,9 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     procedure :: prep_atm_evo_gas
     procedure :: prep_atmosphere => prep_all_evo_gas
     procedure :: right_hand_side_chem
+    procedure :: production_and_loss
     procedure :: right_hand_side => rhs_evo_gas
     procedure :: jacobian => jac_evo_gas
-    procedure :: production_and_loss
 
     !!! photochem_evoatmosphere_integrate.f90 !!!
     procedure :: evolve

src/evoatmosphere/photochem_evoatmosphere_init.f90

@@ -9,6 +9,7 @@ module function create_EvoAtmosphere(mechanism_file, settings_file, flux_file, a
     use iso_c_binding, only : c_associated
     use photochem_input, only: setup
     use photochem_types, only: PhotoSettings
+    use photochem_enum, only: PressureBC
     use clima_types, only: ClimaSettings
 
     character(len=*), intent(in) :: mechanism_file
@@ -45,8 +46,16 @@ module function create_EvoAtmosphere(mechanism_file, settings_file, flux_file, a
                        self%dat%nw)
 
     if (s%evolve_climate) then
-      ! allocate
       self%evolve_climate = .true.
+
+      ! Make sure there are no pressure boundary conditions.
+      if (any(self%var%lowerboundcond == PressureBC)) then
+        err = 'Fixed pressure boundary conditions are not allowed for class "EvoAtmosphere" '// &
+              'when evolve_climate is true.'
+        return
+      endif
+
+      ! allocate
       if (allocated(self%rad)) deallocate(self%rad)
       allocate(self%rad)
 

src/evoatmosphere/photochem_evoatmosphere_integrate.f90

@@ -201,7 +201,8 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, re
     use fsundials_linearsolver_mod, only: SUNLinearSolver, FSUNLinSolFree
     use fsunlinsol_band_mod, only: FSUNLinSol_Band
 
-    use photochem_enum, only: DensityBC
+    use photochem_enum, only: DensityBC, PressureBC
+    use photochem_const, only: k_boltz
     use photochem_types, only: SundialsDataFinalizer
 
     ! in/out
@@ -330,6 +331,8 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, re
     do i = 1,dat%nq
       if (var%lowerboundcond(i) == DensityBC) then
         wrk%sun%yvec(i) = var%lower_fix_den(i)
+      elseif (var%lowerboundcond(i) == PressureBC) then
+        wrk%sun%yvec(i) = var%lower_fix_press(i)/(k_boltz*var%temperature(1))
       endif
     enddo
     ! set abstol
@@ -697,7 +700,8 @@ module subroutine initialize_stepper(self, usol_start, err)
     use fsundials_linearsolver_mod, only: SUNLinearSolver, FSUNLinSolFree
     use fsunlinsol_band_mod, only: FSUNLinSol_Band
 
-    use photochem_enum, only: DensityBC
+    use photochem_enum, only: DensityBC, PressureBC
+    use photochem_const, only: k_boltz
 
     class(EvoAtmosphere), target, intent(inout) :: self
     real(dp), intent(in) :: usol_start(:,:)
@@ -757,6 +761,8 @@ module subroutine initialize_stepper(self, usol_start, err)
     do i = 1,dat%nq
       if (var%lowerboundcond(i) == DensityBC) then
         wrk%sun%yvec(i) = var%lower_fix_den(i)
+      elseif (var%lowerboundcond(i) == PressureBC) then
+        wrk%sun%yvec(i) = var%lower_fix_press(i)/(k_boltz*var%temperature(1))
       endif
     enddo
     ! set abstol

src/evoatmosphere/photochem_evoatmosphere_rhs.f90

@@ -395,7 +395,7 @@ module subroutine prep_atm_evo_gas(self, usol_in, usol, &
     use photochem_eqns, only: press_and_den
     use photochem_common, only: molec_per_particle
     use photochem_const, only: small_real, N_avo, k_boltz
-    use photochem_enum, only: DensityBC
+    use photochem_enum, only: DensityBC, PressureBC
     use photochem_input, only: compute_gibbs_energy, interp2xsdata
     class(EvoAtmosphere), target, intent(inout) :: self
     real(dp), intent(in) :: usol_in(:,:)
@@ -427,6 +427,8 @@ module subroutine prep_atm_evo_gas(self, usol_in, usol, &
     do i = 1,dat%nq
       if (var%lowerboundcond(i) == DensityBC) then
         usol(i,1) = var%lower_fix_den(i)
+      elseif (var%lowerboundcond(i) == PressureBC) then
+        usol(i,1) = var%lower_fix_press(i)/(k_boltz*var%temperature(1))
       endif
     enddo
 
@@ -611,7 +613,7 @@ module subroutine prep_all_evo_gas(self, usol_in, err)
   end subroutine
 
   module subroutine rhs_evo_gas(self, neqs, usol_flat, rhs, err)
-    use photochem_enum, only: MosesBC, VelocityBC, DensityBC, FluxBC, VelocityDistributedFluxBC
+    use photochem_enum, only: MosesBC, VelocityBC, DensityBC, PressureBC, FluxBC, VelocityDistributedFluxBC
     use photochem_enum, only: ZahnleHydrogenEscape
     use iso_c_binding, only: c_ptr, c_f_pointer
     use photochem_const, only: pi, small_real  
@@ -668,7 +670,8 @@ module subroutine rhs_evo_gas(self, neqs, usol_flat, rhs, err)
         rhs(i) = rhs(i) + wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) &
                         + wrk%DD(i,1)*wrk%usol(i,1) + wrk%ADD(i,1)*wrk%usol(i,1) &
                         - wrk%lower_vdep_copy(i)*wrk%usol(i,1)/var%dz(1)
-      elseif (var%lowerboundcond(i) == DensityBC) then
+      elseif (var%lowerboundcond(i) == DensityBC .or. &
+              var%lowerboundcond(i) == PressureBC) then
         rhs(i) = 0.0_dp
       elseif (var%lowerboundcond(i) == FluxBC) then
         rhs(i) = rhs(i) + wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) &
@@ -727,7 +730,7 @@ module subroutine rhs_evo_gas(self, neqs, usol_flat, rhs, err)
   end subroutine
 
   module subroutine jac_evo_gas(self, lda_neqs, neqs, usol_flat, jac, err)
-    use photochem_enum, only: MosesBC, VelocityBC, DensityBC, FluxBC, VelocityDistributedFluxBC
+    use photochem_enum, only: MosesBC, VelocityBC, DensityBC, PressureBC, FluxBC, VelocityDistributedFluxBC
     use photochem_enum, only: ZahnleHydrogenEscape
     use iso_c_binding, only: c_ptr, c_f_pointer
     use photochem_const, only: pi, small_real
@@ -821,7 +824,8 @@ module subroutine jac_evo_gas(self, lda_neqs, neqs, usol_flat, jac, err)
 
         djac(dat%ku,i+dat%nq) = wrk%DU(i,1) + wrk%ADU(i,1)
         djac(dat%kd,i) = djac(dat%kd,i) + wrk%DD(i,1) + wrk%ADD(i,1) - wrk%lower_vdep_copy(i)/var%dz(1)
-      elseif (var%lowerboundcond(i) == DensityBC) then
+      elseif (var%lowerboundcond(i) == DensityBC .or. &
+              var%lowerboundcond(i) == PressureBC) then
 
         do m=1,dat%nq
           mm = dat%kd + i - m

src/input/photochem_input_read.f90

@@ -623,7 +623,7 @@ subroutine H2SO4_interpolator(var, s2, err)
   subroutine unpack_settings(infile, s, dat, var, err)
     use photochem_enum, only: H2SO4Saturation
     use photochem_enum, only: CondensingParticle
-    use photochem_enum, only: VelocityBC, DensityBC, MixingRatioBC
+    use photochem_enum, only: VelocityBC, DensityBC, MixingRatioBC, PressureBC
     use photochem_enum, only: DiffusionLimHydrogenEscape, ZahnleHydrogenEscape, NoHydrogenEscape
     use photochem_types, only: PhotoSettings
     character(len=*), intent(in) :: infile
@@ -843,6 +843,7 @@ subroutine unpack_settings(infile, s, dat, var, err)
       allocate(var%lower_fix_mr(dat%nq))
     else
       allocate(var%lower_fix_den(dat%nq))
+      allocate(var%lower_fix_press(dat%nq))
     endif
     allocate(var%upperboundcond(dat%nq))
     allocate(var%upper_veff(dat%nq))
@@ -890,6 +891,7 @@ subroutine unpack_settings(infile, s, dat, var, err)
           var%lower_fix_mr(ind(1)) = s%lbcs(j)%mix
         else
           var%lower_fix_den(ind(1)) = s%lbcs(j)%den
+          var%lower_fix_press(ind(1)) = s%lbcs(j)%press
         endif
 
         var%upperboundcond(ind(1)) = s%ubcs(j)%bc_type
@@ -942,12 +944,16 @@ subroutine unpack_settings(infile, s, dat, var, err)
     ! make sure bc work for the model
     if (dat%back_gas) then
       if (any(var%lowerboundcond == DensityBC)) then
-        err = 'Fixing density boundary conditions are not allowed for class "Atmosphere".'
+        err = 'Fixed density boundary conditions are not allowed for class "Atmosphere".'
+        return
+      endif
+      if (any(var%lowerboundcond == PressureBC)) then
+        err = 'Fixed pressure boundary conditions are not allowed for class "Atmosphere".'
         return
       endif
     else
       if (any(var%lowerboundcond == MixingRatioBC)) then
-        err = 'Fixing mixing ratio boundary conditions are not allowed for class "EvoAtmosphere".'
+        err = 'Fixed mixing ratio boundary conditions are not allowed for class "EvoAtmosphere".'
         return
       endif
     endif

src/photochem_enum.f90

@@ -39,7 +39,8 @@ module photochem_enum
     MixingRatioBC = 1, &
     FluxBC = 2, &
     VelocityDistributedFluxBC = 3, &
-    DensityBC = 4
+    DensityBC = 4, &
+    PressureBC = 5
 
   enumerator :: &
     DiffusionLimHydrogenEscape, &

src/photochem_types.f90

@@ -35,6 +35,7 @@ module photochem_types ! make a giant IO object
     real(dp) :: flux
     real(dp) :: height
     real(dp) :: den
+    real(dp) :: press
   end type
 
   type :: PhotoSettings
@@ -384,6 +385,7 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
     real(dp), allocatable :: lower_dist_height(:)
     real(dp), allocatable :: lower_fix_mr(:)
     real(dp), allocatable :: lower_fix_den(:)
+    real(dp), allocatable :: lower_fix_press(:)
     integer, allocatable :: upperboundcond(:) ! 0 or 2
     real(dp), allocatable :: upper_veff(:)
     real(dp), allocatable :: upper_flux(:)

src/photochem_types_create.f90

@@ -401,7 +401,7 @@ function unpack_PhotoSettings(root, filename, err) result(s)
   subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
     use photochem_types, only: SettingsBC
     use photochem_enum, only: MosesBC, VelocityBC, MixingRatioBC, FluxBC
-    use photochem_enum, only: VelocityDistributedFluxBC, DensityBC
+    use photochem_enum, only: VelocityDistributedFluxBC, DensityBC, PressureBC
     type(type_dictionary), intent(in) :: bc
     character(*), intent(in) :: bc_kind
     character(*), intent(in) :: sp_name
@@ -417,6 +417,14 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
     elseif (bc_kind == "lower") then
       vel = 'vdep'
     endif
+
+    ! initialize all values
+    sbc%vel = -huge(1.0_dp)
+    sbc%mix = -huge(1.0_dp)
+    sbc%flux = -huge(1.0_dp)
+    sbc%height = -huge(1.0_dp)
+    sbc%den = -huge(1.0_dp)
+    sbc%press = -huge(1.0_dp)
 
     bctype = bc%get_string("type",error = io_err)
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
@@ -425,11 +433,6 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
       sbc%bc_type = VelocityBC
       sbc%vel = bc%get_real(vel,error = io_err)
       if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-
-      sbc%mix = -huge(1.0_dp)
-      sbc%flux = -huge(1.0_dp)
-      sbc%height = -huge(1.0_dp)
-      sbc%den = -huge(1.0_dp)
 
       if (sbc%vel < 0.0_dp) then
         err = 'Velocity '//trim(bc_kind)//' boundary condition for '//trim(sp_name)// &
@@ -446,11 +449,6 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
       sbc%mix = bc%get_real("mix",error = io_err)
       if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
 
-      sbc%vel = -huge(1.0_dp)
-      sbc%flux = -huge(1.0_dp)
-      sbc%height = -huge(1.0_dp)
-      sbc%den = -huge(1.0_dp)
-
       if (sbc%mix < 0.0_dp .or. sbc%mix > 1.0_dp) then
         err = 'Fixed '//trim(bc_kind)//' boundary condition for '//trim(sp_name)// &
               ' must be between 0 and 1.'
@@ -461,10 +459,6 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
       sbc%flux = bc%get_real("flux",error = io_err)
       if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
 
-      sbc%vel = -huge(1.0_dp)
-      sbc%mix = -huge(1.0_dp)
-      sbc%height = -huge(1.0_dp)
-      sbc%den = -huge(1.0_dp)
     elseif (bctype == "vdep + dist flux") then
       if (bc_kind == "upper") then
         err = 'Upper boundary conditions can not be "vdep + dist flux" for '//trim(sp_name)
@@ -480,9 +474,6 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
 
       sbc%height = bc%get_real("height",error = io_err)
       if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-
-      sbc%mix = -huge(1.0_dp)
-      sbc%den = -huge(1.0_dp)
 
       if (sbc%vel < 0.0_dp) then
         err = 'Velocity '//trim(bc_kind)//' boundary condition for '//trim(sp_name)// &
@@ -508,25 +499,29 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err)
       sbc%bc_type = DensityBC
       sbc%den = bc%get_real("den",error = io_err)
       if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-
-      sbc%vel = -huge(1.0_dp)
-      sbc%mix = -huge(1.0_dp)
-      sbc%flux = -huge(1.0_dp)
-      sbc%height = -huge(1.0_dp)
 
-      if (sbc%den < 0.0_dp) then
+      if (sbc%den <= 0.0_dp) then
         err = 'Fixed density '//trim(bc_kind)//' boundary condition for '//trim(sp_name)// &
-              ' must be greater than 1.'
+              ' must be greater than 0.'
+        return
+      endif
+    elseif (bctype == "press") then
+      if (bc_kind == "upper") then
+        err = 'Upper boundary conditions can not be "press" for '//trim(sp_name)
+        return
+      endif
+
+      sbc%bc_type = PressureBC
+      sbc%press = bc%get_real("press",error = io_err)
+      if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
+
+      if (sbc%press <= 0.0_dp) then
+        err = 'Fixed pressure '//trim(bc_kind)//' boundary condition for '//trim(sp_name)// &
+              ' must be greater than 0.'
         return
       endif
     elseif (bctype == "Moses") then
       sbc%bc_type = MosesBC
-
-      sbc%vel = -huge(1.0_dp)
-      sbc%mix = -huge(1.0_dp)
-      sbc%flux = -huge(1.0_dp)
-      sbc%height = -huge(1.0_dp)
-      sbc%den = -huge(1.0_dp)
     else
       err = 'IOError: "'//trim(bctype)//'" is not a valid lower boundary condition for '//trim(sp_name)
       return

tests/testevo_settings1.yaml

@@ -2,7 +2,7 @@
 atmosphere-grid:
   bottom: 0.0 # cm
   top: 1.0e7 # cm
-  number-of-layers: 100
+  number-of-layers: 80
 
 photolysis-grid:
   regular-grid: true
@@ -13,7 +13,7 @@ photolysis-grid:
 planet:
   planet-mass: 5.972e27 # grams
   planet-radius: 6.371e8 # cm. Radius to bottom of atmosphere-grid
-  surface-albedo: 0.2 # cm
+  surface-albedo: 0.3 # cm
   solar-zenith-angle: 60.0
   hydrogen-escape:
     type: diffusion limited
@@ -56,17 +56,15 @@ particles:
 
 optical-properties:
   ir:
-    k-method: RandomOverlapResortRebin
-    number-of-bins: 8
+    k-method: AdaptiveEquivalentExtinction
     opacities:
       k-distributions: [H2O, CO2, CH4, O3]
       CIA: [N2-N2, O2-O2, N2-O2, CH4-CH4, CO2-CO2, H2-CH4, H2-H2]
       rayleigh: true
       photolysis-xs: true
       water-continuum: MT_CKD
   solar:
-    k-method: RandomOverlapResortRebin
-    number-of-bins: 8
+    k-method: AdaptiveEquivalentExtinction
     opacities:
       k-distributions: [H2O, CO2, CH4, O3]
       CIA: [N2-N2, O2-O2, N2-O2, CH4-CH4, CO2-CO2, H2-CH4, H2-H2]

tests/testevo_settings2.yaml

@@ -15,6 +15,7 @@ planet:
   planet-radius: 6.371e8 # cm. Radius to bottom of atmosphere-grid
   surface-albedo: 0.2 # cm
   solar-zenith-angle: 60.0
+  # evolve-climate: true
   hydrogen-escape:
     type: diffusion limited
   water:
@@ -55,10 +56,10 @@ particles:
 
 boundary-conditions:
 - name: N2
-  lower-boundary: {type: den, den: 1.9091467479365386e+19}
+  lower-boundary: {type: press, press: 0.78e6}
   upper-boundary: {type: veff, veff: 1.0e-8}
 - name: O2
-  lower-boundary: {type: flux, flux: 529988353109.4235}
+  lower-boundary: {type: press, press: 0.21e6}
   upper-boundary: {type: veff, veff: 0.0}
 - name: H2
   lower-boundary: {type: vdep, vdep: 0.002101353876412517}