Skip to content

Commit

Permalink
deleted not needed background gas vars
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicholaswogan committed Feb 3, 2025
1 parent 8bbc520 commit a9f6ce2
Show file tree
Hide file tree
Showing 9 changed files with 29 additions and 205 deletions.
3 changes: 1 addition & 2 deletions src/evoatmosphere/photochem_evoatmosphere_init.f90
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,7 @@ module function create_EvoAtmosphere(mechanism_file, settings_file, flux_file, a
allocate(self%wrk)

self%var%data_dir = data_dir
call setup(mechanism_file, s, flux_file, atmosphere_txt, .false., &
self%dat, self%var, err)
call setup(mechanism_file, s, flux_file, atmosphere_txt, self%dat, self%var, err)
if (allocated(err)) return

call self%wrk%init(self%dat%nsp, self%dat%np, self%dat%nq, &
Expand Down
4 changes: 2 additions & 2 deletions src/evoatmosphere/photochem_evoatmosphere_utils.f90
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@ module subroutine set_rate_fcn(self, species, fcn, err)
ind = findloc(self%dat%species_names(1:self%dat%nq), species, 1)
if (ind == 0) then
err = 'Species "'//species//'" is not in the list of species, '// &
'or is a background or short-lived species.'
'or is a short-lived species.'
return
endif

Expand Down Expand Up @@ -607,7 +607,7 @@ subroutine properties_for_new_TOA(self, usol, top_atmos_new, &
particle_radius_new, pressure_new, err)
use photochem_enum, only: DensityBC, PressureBC
use futils, only: interp
use photochem_eqns, only: vertical_grid, molar_weight, press_and_den, gravity
use photochem_eqns, only: vertical_grid, press_and_den, gravity
use photochem_const, only: small_real, k_boltz
class(EvoAtmosphere), target, intent(inout) :: self
real(dp), intent(in) :: usol(:,:)
Expand Down
11 changes: 3 additions & 8 deletions src/input/photochem_input.f90
Original file line number Diff line number Diff line change
Expand Up @@ -24,13 +24,11 @@ module subroutine after_read_setup(dat, var, err)
character(:), allocatable, intent(out) :: err
end subroutine

module subroutine read_all_files(mechanism_file, s, flux_file, atmosphere_txt, back_gas, &
dat, var, err)
module subroutine read_all_files(mechanism_file, s, flux_file, atmosphere_txt, dat, var, err)
character(len=*), intent(in) :: mechanism_file
type(PhotoSettings), intent(in) :: s
character(len=*), intent(in) :: flux_file
character(len=*), intent(in) :: atmosphere_txt
logical, intent(in) :: back_gas
type(PhotochemData), intent(inout) :: dat
type(PhotochemVars), intent(inout) :: var
character(:), allocatable, intent(out) :: err
Expand Down Expand Up @@ -59,20 +57,17 @@ module subroutine interp2particlexsdata(dat, var, err)

contains

subroutine setup(mechanism_file, s, flux_file, atmosphere_txt, back_gas, &
dat, var, err)
subroutine setup(mechanism_file, s, flux_file, atmosphere_txt, dat, var, err)

character(len=*), intent(in) :: mechanism_file
type(PhotoSettings), intent(in) :: s
character(len=*), intent(in) :: flux_file
character(len=*), intent(in) :: atmosphere_txt
logical, intent(in) :: back_gas
type(PhotochemData), intent(inout) :: dat
type(PhotochemVars), intent(inout) :: var
character(:), allocatable, intent(out) :: err

call read_all_files(mechanism_file, s, flux_file, atmosphere_txt, back_gas, &
dat, var, err)
call read_all_files(mechanism_file, s, flux_file, atmosphere_txt, dat, var, err)
if (allocated(err)) return

call after_read_setup(dat, var, err)
Expand Down
25 changes: 5 additions & 20 deletions src/input/photochem_input_after_read.f90
Original file line number Diff line number Diff line change
Expand Up @@ -74,27 +74,12 @@ subroutine interp2atmosfile(dat, var, err)
endif
var%edd = 10.0_dp**var%edd

if (dat%back_gas) then
do i = 1,dat%nq
call interp(var%nz, dat%nzf, var%z, dat%z_file,&
log10(abs(dat%mix_file(i,:))), var%usol_init(i,:), ierr)
if (ierr /= 0) then
err = 'Subroutine interp returned an error.'
return
endif
enddo
var%usol_init = 10.0_dp**var%usol_init

if (dat%conserving_init) then
call interp2atmosfile_mixconserving(dat, var, err)
if (allocated(err)) return
else

if (dat%conserving_init) then
call interp2atmosfile_mixconserving(dat, var, err)
if (allocated(err)) return
else
call interp2atmosfile_mix(dat, var, err)
if (allocated(err)) return
endif

call interp2atmosfile_mix(dat, var, err)
if (allocated(err)) return
endif

if (dat%there_are_particles) then
Expand Down
118 changes: 13 additions & 105 deletions src/input/photochem_input_read.f90
Original file line number Diff line number Diff line change
Expand Up @@ -3,27 +3,16 @@

contains

module subroutine read_all_files(mechanism_file, s, flux_file, atmosphere_txt, back_gas, &
dat, var, err)
module subroutine read_all_files(mechanism_file, s, flux_file, atmosphere_txt, dat, var, err)
character(len=*), intent(in) :: mechanism_file
type(PhotoSettings), intent(in) :: s
character(len=*), intent(in) :: flux_file
character(len=*), intent(in) :: atmosphere_txt
logical, intent(in) :: back_gas
type(PhotochemData), intent(inout) :: dat
type(PhotochemVars), intent(inout) :: var
character(:), allocatable, intent(out) :: err

! stuff dat needs before entering get_photomech
dat%back_gas = back_gas
if (dat%back_gas) then
if (allocated(s%back_gas_name)) then
dat%back_gas_name = s%back_gas_name
else
err = 'A background gas is required but not specified in '//trim(s%filename)
return
endif
endif
dat%nsl = s%nsl
dat%SL_names = s%SL_names

Expand Down Expand Up @@ -97,7 +86,7 @@ subroutine get_rxmechanism(mapping, infile, dat, var, err)
! temporary work variables
character(len=str_len) :: tmpchar
character(len=str_len) :: tmp
character(len=:), allocatable :: rxstring, back_gas_name_tmp
character(len=:), allocatable :: rxstring
integer :: i, ii, j, k, kk, l, ind(1)
logical :: reverse
! all_species causes a small memory leak. Not sure how to free the memory properly
Expand Down Expand Up @@ -266,13 +255,7 @@ subroutine get_rxmechanism(mapping, infile, dat, var, err)
dat%ng = dat%ng + 1
enddo

if (dat%back_gas) then
dat%nll = dat%ng - dat%nsl - 1 ! minus 1 for background
back_gas_name_tmp = dat%back_gas_name
else
dat%nll = dat%ng - dat%nsl
back_gas_name_tmp = "Not a => gas!"
endif
dat%nll = dat%ng - dat%nsl

dat%ng_1 = dat%npq + 1 ! the long lived gas index
! dat%nq is the last ll gas index
Expand Down Expand Up @@ -326,8 +309,6 @@ subroutine get_rxmechanism(mapping, infile, dat, var, err)
if (ind(1) /= 0) then ! short lived species
j = dat%nq + l
l = l + 1
elseif (tmpchar == back_gas_name_tmp) then ! background gas
j = dat%nsp
else ! long lived species
j = kk
kk = kk + 1
Expand Down Expand Up @@ -372,13 +353,6 @@ subroutine get_rxmechanism(mapping, infile, dat, var, err)
err = 'IOError: One of the short lived species is not in the file '//trim(infile)
return
endif
if (dat%back_gas) then
ind = findloc(dat%species_names,dat%back_gas_name)
if (ind(1) == 0) then
err = 'IOError: The specified background gas is not in '//trim(infile)
return
endif
endif
i = check_for_duplicates(dat%species_names)
if (i /= 0) then
err = 'Species "'//trim(dat%species_names(i))//'" is a duplicate in '//trim(infile)
Expand All @@ -405,7 +379,7 @@ subroutine get_rxmechanism(mapping, infile, dat, var, err)
err = "Particle "//trim(dat%particle_names(i))// &
" can not be made from "//trim(dat%particle_gas_phase(i))// &
" because "//trim(dat%particle_gas_phase(i))//" is"// &
" is a particle, short-lived species, or background gas."
" is a particle or is short-lived species."
return
else
dat%particle_gas_phase_ind(i) = ind(1)
Expand Down Expand Up @@ -586,7 +560,7 @@ subroutine get_rxmechanism(mapping, infile, dat, var, err)

subroutine unpack_settings(infile, s, dat, var, err)
use photochem_enum, only: CondensingParticle
use photochem_enum, only: VelocityBC, DensityBC, MixingRatioBC, PressureBC
use photochem_enum, only: VelocityBC, DensityBC, PressureBC
use photochem_enum, only: DiffusionLimHydrogenEscape, ZahnleHydrogenEscape, NoHydrogenEscape
use photochem_types, only: PhotoSettings
character(len=*), intent(in) :: infile
Expand Down Expand Up @@ -622,16 +596,6 @@ subroutine unpack_settings(infile, s, dat, var, err)
!!!!!!!!!!!!!!
!!! planet !!!
!!!!!!!!!!!!!!
! dat%back_gas
! dat%back_gas_name
! already set earlier
if (dat%back_gas) then
if (.not. allocated(s%P_surf)) then
err = 'The settings file does not contain a "surface-pressure"'
return
endif
var%surface_pressure = s%P_surf
endif
dat%planet_mass = s%planet_mass
dat%planet_radius = s%planet_radius
var%surface_albedo = s%surface_albedo
Expand Down Expand Up @@ -672,13 +636,6 @@ subroutine unpack_settings(infile, s, dat, var, err)
allocate(dat%H_escape_coeff)
dat%H_escape_coeff = zahnle_Hescape_coeff(s%H_escape_S1)

if (dat%back_gas) then
if (dat%back_gas_name == "H2") then
err = "zahnle hydrogen escape does not work when the background gas is H2."
return
endif
endif

endblock; endif

endif
Expand Down Expand Up @@ -770,23 +727,12 @@ subroutine unpack_settings(infile, s, dat, var, err)
!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! boundary-conditions !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!

if (dat%back_gas) then
ind = findloc(dat%species_names,trim(dat%back_gas_name))
dat%back_gas_ind = ind(1)
dat%back_gas_mu = dat%species_mass(ind(1))
endif

allocate(var%lowerboundcond(dat%nq))
allocate(var%lower_vdep(dat%nq))
allocate(var%lower_flux(dat%nq))
allocate(var%lower_dist_height(dat%nq))
if (dat%back_gas) then
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%lower_fix_den(dat%nq))
allocate(var%lower_fix_press(dat%nq))
allocate(var%upperboundcond(dat%nq))
allocate(var%upper_veff(dat%nq))
allocate(var%upper_flux(dat%nq))
Expand Down Expand Up @@ -815,26 +761,15 @@ subroutine unpack_settings(infile, s, dat, var, err)
' in settings file is not in the reaction mechanism file.'
return
endif
if (dat%back_gas) then
if (ind(1) == dat%back_gas_ind) then ! can't be background gas
err = "IOError: Species "//trim(s%sp_names(j))// &
' in settings file is the background gas, and can not have boundary conditions.'
return
endif
endif

if (s%sp_types(j) == 'long lived') then

var%lowerboundcond(ind(1)) = s%lbcs(j)%bc_type
var%lower_vdep(ind(1)) = s%lbcs(j)%vel
var%lower_flux(ind(1)) = s%lbcs(j)%flux
var%lower_dist_height(ind(1)) = s%lbcs(j)%height
if (dat%back_gas) then
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%lower_fix_den(ind(1)) = s%lbcs(j)%den
var%lower_fix_press(ind(1)) = s%lbcs(j)%press

var%upperboundcond(ind(1)) = s%ubcs(j)%bc_type
var%upper_veff(ind(1)) = s%ubcs(j)%vel
Expand All @@ -849,20 +784,10 @@ subroutine unpack_settings(infile, s, dat, var, err)
! Make sure that upper boundary condition for H and H2 are
! effusion velocities, if diffusion limited escape
if (dat%H_escape_type == DiffusionLimHydrogenEscape) then
if (dat%back_gas) then
if (dat%back_gas_name /= "H2") then
if (var%upperboundcond(dat%LH2) /= VelocityBC) then
err = "IOError: H2 must have a have a effusion velocity upper boundary"// &
" for diffusion limited hydrogen escape"
return
endif
endif
else
if (var%upperboundcond(dat%LH2) /= VelocityBC) then
err = "IOError: H2 must have a have a effusion velocity upper boundary"// &
" for diffusion limited hydrogen escape"
return
endif
if (var%upperboundcond(dat%LH2) /= VelocityBC) then
err = "IOError: H2 must have a have a effusion velocity upper boundary"// &
" for diffusion limited hydrogen escape"
return
endif
if (var%upperboundcond(dat%LH) /= VelocityBC) then
err = "IOError: H must have a have a effusion velocity upper boundary"// &
Expand All @@ -883,23 +808,6 @@ subroutine unpack_settings(infile, s, dat, var, err)
enddo
endif

! make sure bc work for the model
if (dat%back_gas) then
if (any(var%lowerboundcond == DensityBC)) then
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 = 'Fixed mixing ratio boundary conditions are not allowed for class "EvoAtmosphere".'
return
endif
endif

! check for SL nonlinearities
call check_sl(dat, err)
if (allocated(err)) return
Expand Down
1 change: 0 additions & 1 deletion src/photochem_enum.f90
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@ module photochem_enum
enumerator :: &
MosesBC = -1, &
VelocityBC = 0, &
MixingRatioBC = 1, &
FluxBC = 2, &
VelocityDistributedFluxBC = 3, &
DensityBC = 4, &
Expand Down
20 changes: 0 additions & 20 deletions src/photochem_eqns.f90
Original file line number Diff line number Diff line change
Expand Up @@ -221,26 +221,6 @@ pure subroutine press_and_den(nz, T, grav, Psurf, dz, &

end subroutine

pure subroutine molar_weight(nll, usol_layer, sum_usol_layer, masses, background_mu, mubar_layer)
implicit none
integer, intent(in) :: nll
real(dp), intent(in) :: usol_layer(nll)
real(dp), intent(in) :: sum_usol_layer
real(dp), intent(in) :: masses(nll)
real(dp), intent(in) :: background_mu
real(dp), intent(out) :: mubar_layer
integer :: j
real(dp) :: f_background

mubar_layer = 0.0_dp
do j = 1, nll
mubar_layer = mubar_layer + usol_layer(j) * masses(j)
enddo
f_background = 1.0_dp - sum_usol_layer
mubar_layer = mubar_layer + f_background * background_mu

end subroutine

pure function dynamic_viscosity_air(T) result(eta)
real(dp), intent(in) :: T
real(dp) :: eta ! dynamic viscosity [dynes s/cm^2]
Expand Down
14 changes: 4 additions & 10 deletions src/photochem_types.f90
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,6 @@ module photochem_types ! make a giant IO object
integer :: nz

! planet
character(:), allocatable :: back_gas_name
real(dp), allocatable :: P_surf
real(dp) :: planet_mass
real(dp) :: planet_radius
real(dp) :: surface_albedo
Expand Down Expand Up @@ -266,12 +264,12 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)

! species
! Organization is as follows
! [ nsp ]
! [ nq + nsl + 1]
! [np + ng + nsl + 1]
! [ nsp ]
! [ nq + nsl]
! [np + ng + nsl]
! |_______|
! |
! Only np + ng = nq evolve through time. nsl are assumed to be in equilibrium. +1 is background gas.
! Only np + ng = nq evolve through time. nsl are assumed to be in equilibrium.
integer :: nq !! number of gases + particles which evolve over time from integration
integer :: ng_1 !! index of first gas
integer :: nll !! number of long-lived gas molecules
Expand Down Expand Up @@ -346,10 +344,6 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
real(dp), allocatable :: particle_radius_file(:,:) !! (np,nzf) cm

! settings
logical :: back_gas !! True if background gas is used
character(:), allocatable :: back_gas_name !! Normally N2, but can be most any gas.
real(dp), allocatable :: back_gas_mu !! g/mol
integer, allocatable :: back_gas_ind !! index of the background gas
real(dp) :: planet_mass !! grams
real(dp) :: planet_radius !! cm
logical :: fix_water_in_trop !! True if fixing water in troposphere
Expand Down
Loading

0 comments on commit a9f6ce2

Please sign in to comment.