From a9f6ce21da69926a6bc32f54830b41ba8204831f Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Mon, 3 Feb 2025 10:22:47 -0800 Subject: [PATCH] deleted not needed background gas vars --- .../photochem_evoatmosphere_init.f90 | 3 +- .../photochem_evoatmosphere_utils.f90 | 4 +- src/input/photochem_input.f90 | 11 +- src/input/photochem_input_after_read.f90 | 25 +--- src/input/photochem_input_read.f90 | 118 ++---------------- src/photochem_enum.f90 | 1 - src/photochem_eqns.f90 | 20 --- src/photochem_types.f90 | 14 +-- src/photochem_types_create.f90 | 38 +----- 9 files changed, 29 insertions(+), 205 deletions(-) diff --git a/src/evoatmosphere/photochem_evoatmosphere_init.f90 b/src/evoatmosphere/photochem_evoatmosphere_init.f90 index e6a97c61..3e5f3e76 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_init.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_init.f90 @@ -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, & diff --git a/src/evoatmosphere/photochem_evoatmosphere_utils.f90 b/src/evoatmosphere/photochem_evoatmosphere_utils.f90 index 2a03072b..e1482049 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_utils.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_utils.f90 @@ -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 @@ -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(:,:) diff --git a/src/input/photochem_input.f90 b/src/input/photochem_input.f90 index e2232129..7bbd013f 100644 --- a/src/input/photochem_input.f90 +++ b/src/input/photochem_input.f90 @@ -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 @@ -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) diff --git a/src/input/photochem_input_after_read.f90 b/src/input/photochem_input_after_read.f90 index 5ba12363..7255f1ae 100644 --- a/src/input/photochem_input_after_read.f90 +++ b/src/input/photochem_input_after_read.f90 @@ -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 diff --git a/src/input/photochem_input_read.f90 b/src/input/photochem_input_read.f90 index 82eb7ac7..29364a4a 100644 --- a/src/input/photochem_input_read.f90 +++ b/src/input/photochem_input_read.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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) @@ -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 @@ -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 @@ -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 @@ -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)) @@ -815,13 +761,6 @@ 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 @@ -829,12 +768,8 @@ subroutine unpack_settings(infile, s, dat, var, err) 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 @@ -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"// & @@ -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 diff --git a/src/photochem_enum.f90 b/src/photochem_enum.f90 index 4315cff4..bc9c3dc8 100644 --- a/src/photochem_enum.f90 +++ b/src/photochem_enum.f90 @@ -32,7 +32,6 @@ module photochem_enum enumerator :: & MosesBC = -1, & VelocityBC = 0, & - MixingRatioBC = 1, & FluxBC = 2, & VelocityDistributedFluxBC = 3, & DensityBC = 4, & diff --git a/src/photochem_eqns.f90 b/src/photochem_eqns.f90 index 9815e534..06db42f7 100644 --- a/src/photochem_eqns.f90 +++ b/src/photochem_eqns.f90 @@ -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] diff --git a/src/photochem_types.f90 b/src/photochem_types.f90 index 6f14af56..e1a91ddf 100644 --- a/src/photochem_types.f90 +++ b/src/photochem_types.f90 @@ -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 @@ -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 @@ -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 diff --git a/src/photochem_types_create.f90 b/src/photochem_types_create.f90 index cddb2cb3..cd2a6a6c 100644 --- a/src/photochem_types_create.f90 +++ b/src/photochem_types_create.f90 @@ -42,10 +42,8 @@ function unpack_PhotoSettings(root, filename, err) result(s) type(type_dictionary), pointer :: dict, tmp2 type(type_list), pointer :: list, bcs type(type_list_item), pointer :: item - type(type_scalar), pointer :: scalar type(type_error), allocatable :: io_err character(:), allocatable :: temp_char - logical :: success integer :: i, j real(dp) :: tmp_real @@ -79,25 +77,6 @@ function unpack_PhotoSettings(root, filename, err) result(s) !!!!!!!!!!!!!! dict => root%get_dictionary('planet',.true.,error = io_err) if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - - scalar => dict%get_scalar('background-gas',required=.false.,error = io_err) - if (associated(scalar)) then - s%back_gas_name = trim(scalar%string) - endif - - scalar => dict%get_scalar('surface-pressure',required=.false.,error = io_err) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - if (associated(scalar)) then - s%P_surf = scalar%to_real(-1.0_dp, success=success) - if (.not. success) then - err = trim(filename)//trim(scalar%path)//' can not be converted to a real number.' - return - endif - if (s%P_surf <= 0.0_dp) then - err = 'IOError: Planet surface pressure must be greater than zero.' - return - endif - endif s%planet_mass = dict%get_real('planet-mass',error = io_err) if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif @@ -363,7 +342,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: MosesBC, VelocityBC, FluxBC use photochem_enum, only: VelocityDistributedFluxBC, DensityBC, PressureBC type(type_dictionary), intent(in) :: bc character(*), intent(in) :: bc_kind @@ -402,21 +381,6 @@ subroutine unpack_SettingsBC(bc, bc_kind, sp_name, filename, sbc, err) ' must be positive.' return endif - elseif (bctype == "mix") then - if (bc_kind == "upper") then - err = 'Upper boundary conditions can not be "mix" for '//trim(sp_name) - return - endif - - sbc%bc_type = MixingRatioBC - sbc%mix = bc%get_real("mix",error = io_err) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - - 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.' - return - endif elseif (bctype == "flux") then sbc%bc_type = FluxBC sbc%flux = bc%get_real("flux",error = io_err)