Skip to content

Commit

Permalink
added ability to set cumstom binary diffusion parameter function
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicholaswogan committed Nov 10, 2023
1 parent c44e8c2 commit 8bb2926
Show file tree
Hide file tree
Showing 6 changed files with 76 additions and 7 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
cmake_minimum_required(VERSION "3.14")
cmake_policy(SET CMP0148 OLD)
project(Photochem LANGUAGES Fortran C VERSION "0.4.4")
project(Photochem LANGUAGES Fortran C VERSION "0.4.5")

include(FortranCInterface)
FortranCInterface_VERIFY()
Expand Down
12 changes: 11 additions & 1 deletion src/atmosphere/photochem_atmosphere_rhs.f90
Original file line number Diff line number Diff line change
Expand Up @@ -169,9 +169,10 @@ subroutine dochem(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol, density
subroutine diffusion_coefficients(dat, var, den, mubar, &
DU, DD, DL, ADU, ADL, ADD, wfall, VH2_esc, VH_esc)
use photochem_eqns, only: dynamic_viscosity_air, fall_velocity, slip_correction_factor, &
binary_diffusion_param
default_binary_diffusion_param
use photochem_const, only: k_boltz, N_avo
use photochem_enum, only: DiffusionLimHydrogenEscape
use photochem_types, only: binary_diffusion_fcn

type(PhotochemData), intent(in) :: dat
type(PhotochemVars), intent(in) :: var
Expand All @@ -195,7 +196,16 @@ subroutine diffusion_coefficients(dat, var, den, mubar, &
real(dp) :: viscosity_pp, viscosity
real(dp) :: FF2, FF1

! molecular diffusion function
procedure(binary_diffusion_fcn), pointer :: binary_diffusion_param

integer :: i, j

if (associated(var%custom_binary_diffusion_fcn)) then
binary_diffusion_param => var%custom_binary_diffusion_fcn
else
binary_diffusion_param => default_binary_diffusion_param
endif

! eddy diffusion. particles and gases
! middle
Expand Down
11 changes: 10 additions & 1 deletion src/evoatmosphere/photochem_evoatmosphere_rhs.f90
Original file line number Diff line number Diff line change
Expand Up @@ -171,9 +171,10 @@ subroutine dochem(self, usol, rx_rates, &
subroutine diffusion_coefficients_evo(dat, var, den, mubar, &
DU, DD, DL, ADU, ADL, ADD, wfall, VH2_esc, VH_esc)
use photochem_eqns, only: dynamic_viscosity_air, fall_velocity, slip_correction_factor, &
binary_diffusion_param
default_binary_diffusion_param
use photochem_const, only: k_boltz, N_avo
use photochem_enum, only: DiffusionLimHydrogenEscape
use photochem_types, only: binary_diffusion_fcn

type(PhotochemData), intent(in) :: dat
type(PhotochemVars), intent(in) :: var
Expand All @@ -199,9 +200,17 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, &
real(dp) :: viscosity_pp, viscosity
real(dp) :: FF2, FF1

! molecular diffusion function
procedure(binary_diffusion_fcn), pointer :: binary_diffusion_param

integer :: j, i, k

if (associated(var%custom_binary_diffusion_fcn)) then
binary_diffusion_param => var%custom_binary_diffusion_fcn
else
binary_diffusion_param => default_binary_diffusion_param
endif

! compute relevant parameters at the edges of the grid cells
do j = 1,var%nz-1
eddav(j) = sqrt(var%edd(j)*var%edd(j+1))
Expand Down
7 changes: 4 additions & 3 deletions src/photochem_eqns.f90
Original file line number Diff line number Diff line change
Expand Up @@ -267,9 +267,10 @@ pure function slip_correction_factor(partical_radius, density) result(correct_fa
(1.257e0_dp + 0.4e0_dp*exp((-1.1e0_dp*partical_radius)/(lambda)))
end function

pure function binary_diffusion_param(mu_i, mubar, T) result(b)
real(dp), intent(in) :: mu_i, mubar, T
real(dp) :: b
function default_binary_diffusion_param(mu_i, mubar, T) result(b)
use iso_c_binding, only: c_double
real(c_double), value, intent(in) :: mu_i, mubar, T
real(c_double) :: b
! Banks and Kockarts 1973, Eq 15.29
! also Catling and Kasting 2017, Eq B.4 (although Catling has a typo,
! and is missing a power of 0.5)
Expand Down
14 changes: 13 additions & 1 deletion src/photochem_types.f90
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ module photochem_types ! make a giant IO object
public :: Reaction, Efficiencies, BaseRate, PhotolysisRate
public :: ElementaryRate, ThreeBodyRate, FalloffRate, ProdLoss
public :: SundialsDataFinalizer
public :: time_dependent_flux_fcn, time_dependent_rate_fcn
public :: time_dependent_flux_fcn, time_dependent_rate_fcn, binary_diffusion_fcn

!!!!!!!!!!!!!!!!
!!! Settings !!!
Expand Down Expand Up @@ -107,6 +107,16 @@ module function create_PhotoSettings(filename, err) result(s)
!!!!!!!!!!!!!!!!!

abstract interface
!> Custom binary diffusion function
function binary_diffusion_fcn(mu_i, mubar, T) result(b)
use iso_c_binding, only: c_double
real(c_double), value, intent(in) :: mu_i !! molar weight of species i (g/mol)
real(c_double), value, intent(in) :: mubar !! molar weight of background gas (g/mol)
real(c_double), value, intent(in) :: T !! Temperature (K)
real(c_double) :: b !! binary diffusion parameter of species i
!! with respect to the background gas (molecules cm^-1 s^1)
end function

!> Sets the time-dependent photon flux
subroutine time_dependent_flux_fcn(tn, nw, photon_flux)
use iso_c_binding, only: c_double, c_int
Expand Down Expand Up @@ -416,6 +426,8 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
real(dp), allocatable :: z(:) !! (nz) cm
real(dp), allocatable :: dz(:) !! (nz) cm
real(dp), allocatable :: edd(:) !! (nz) cm2/s
!> A function for specifying a custom binary diffusion parameter (b_ij)
procedure(binary_diffusion_fcn), nopass, pointer :: custom_binary_diffusion_fcn => null()
real(dp), allocatable :: grav(:) !! (nz) cm/s2
real(dp), allocatable :: usol_init(:,:) !! (nq,nz) mixing ratio (Atmosphere) or densities (EvoAtmosphere).
real(dp), allocatable :: particle_radius(:,:) !! (np,nz) cm
Expand Down
37 changes: 37 additions & 0 deletions tests/memtest.f90
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ program memtest
call test_atom_conservation(pcs)
call test_evolve(pcs)
call test_set_temperature(pcs)
call test_custom_binary_diffusion(pcs)

contains

Expand Down Expand Up @@ -245,5 +246,41 @@ subroutine test_set_temperature(pc)
endif

end subroutine

subroutine test_custom_binary_diffusion(pc)
type(Atmosphere), intent(inout) :: pc
character(:), allocatable :: err
real(dp) :: tn

pc%var%custom_binary_diffusion_fcn => custom_binary_diffusion_param

call pc%initialize_stepper(pc%var%usol_init, err)
if (allocated(err)) then
print*,trim(err)
stop 1
endif

tn = pc%step(err)
if (allocated(err)) then
print*,trim(err)
stop 1
endif

call pc%destroy_stepper(err)
if (allocated(err)) then
print*,trim(err)
stop 1
endif

pc%var%custom_binary_diffusion_fcn => null()

end subroutine

function custom_binary_diffusion_param(mu_i, mubar, T) result(b)
use iso_c_binding, only: c_double
real(c_double), value, intent(in) :: mu_i, mubar, T
real(c_double) :: b
b = 1.52e18_dp*((1.0_dp/mu_i+1.0_dp/mubar)**0.5e0_dp)*(T**0.5e0_dp)
end function

end program

0 comments on commit 8bb2926

Please sign in to comment.