From 3f0d2358512e3b5aca28eaea0a47d5ed894e9e64 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Mon, 22 Jul 2024 07:56:13 -0700 Subject: [PATCH] optional upwind scheme for molecular advection for stability --- photochem/cython/PhotochemVars.pyx | 12 ++++ photochem/cython/PhotochemVars_pxd.pxd | 5 +- photochem/fortran/PhotochemVars_wrapper.f90 | 16 ++++++ .../photochem_evoatmosphere_rhs.f90 | 55 ++++++++++++++++++- src/photochem_types.f90 | 4 ++ 5 files changed, 89 insertions(+), 3 deletions(-) diff --git a/photochem/cython/PhotochemVars.pyx b/photochem/cython/PhotochemVars.pyx index 8c222ab..ac953dd 100644 --- a/photochem/cython/PhotochemVars.pyx +++ b/photochem/cython/PhotochemVars.pyx @@ -345,5 +345,17 @@ cdef class PhotochemVars: def __set__(self, double val): var_pxd.photochemvars_fast_arbitrary_rate_set(self._ptr, &val) + property upwind_molec_diff: + """bool. If True, then the code uses a 1st order upwind method for the advective molecular + diffusion terms instead of a centered scheme. This permits stability (at the cost + of accuracy) for atmospheres with strong molcular advection in the upper atmosphere. + """ + def __get__(self): + cdef bool val + var_pxd.photochemvars_upwind_molec_diff_get(self._ptr, &val) + return val + def __set__(self, bool val): + var_pxd.photochemvars_upwind_molec_diff_set(self._ptr, &val) + \ No newline at end of file diff --git a/photochem/cython/PhotochemVars_pxd.pxd b/photochem/cython/PhotochemVars_pxd.pxd index f362cd6..1fd176e 100644 --- a/photochem/cython/PhotochemVars_pxd.pxd +++ b/photochem/cython/PhotochemVars_pxd.pxd @@ -98,4 +98,7 @@ cdef extern void photochemvars_verbose_get(PhotochemVars *ptr, int *val) cdef extern void photochemvars_verbose_set(PhotochemVars *ptr, int *val) cdef extern void photochemvars_fast_arbitrary_rate_get(PhotochemVars *ptr, double *val) -cdef extern void photochemvars_fast_arbitrary_rate_set(PhotochemVars *ptr, double *val) \ No newline at end of file +cdef extern void photochemvars_fast_arbitrary_rate_set(PhotochemVars *ptr, double *val) + +cdef extern void photochemvars_upwind_molec_diff_get(PhotochemVars *ptr, bool *val) +cdef extern void photochemvars_upwind_molec_diff_set(PhotochemVars *ptr, bool *val) \ No newline at end of file diff --git a/photochem/fortran/PhotochemVars_wrapper.f90 b/photochem/fortran/PhotochemVars_wrapper.f90 index 8bf83a5..5105c85 100644 --- a/photochem/fortran/PhotochemVars_wrapper.f90 +++ b/photochem/fortran/PhotochemVars_wrapper.f90 @@ -487,4 +487,20 @@ subroutine photochemvars_fast_arbitrary_rate_set(ptr, val) bind(c) call c_f_pointer(ptr, var) var%fast_arbitrary_rate = val end subroutine + + subroutine photochemvars_upwind_molec_diff_get(ptr, val) bind(c) + type(c_ptr), value, intent(in) :: ptr + logical(c_bool), intent(out) :: val + type(PhotochemVars), pointer :: var + call c_f_pointer(ptr, var) + val = var%upwind_molec_diff + end subroutine + + subroutine photochemvars_upwind_molec_diff_set(ptr, val) bind(c) + type(c_ptr), value, intent(in) :: ptr + logical(c_bool), intent(in) :: val + type(PhotochemVars), pointer :: var + call c_f_pointer(ptr, var) + var%upwind_molec_diff = val + end subroutine \ No newline at end of file diff --git a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 index 79f5e36..81e1fe1 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 @@ -242,6 +242,7 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & real(dp) :: scale_height_av(var%nz-1),scale_height_i_av real(dp) :: b1x2av(dat%nll,var%nz-1) real(dp) :: gamma_i_gas_av(dat%nll,var%nz-1), gamma_i_part_av(dat%np,var%nz-1) + real(dp) :: gamma_i_gas_D(dat%nll,var%nz) real(dp) :: grav_av, mubar_av real(dp) :: bx1x2 @@ -262,6 +263,33 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & binary_diffusion_param => default_binary_diffusion_param endif + if (var%upwind_molec_diff) then; block + real(dp) :: dTdz_tmp, scale_height_i_tmp, b1x2_tmp, grav_tmp, mubar_tmp + ! If using upwind scheme for molecular diffusion, then compute the needed + ! advective velocities. + do j = 1,var%nz-1 + dTdz_tmp = (var%temperature(j+1) - var%temperature(j))/var%dz(j) + grav_tmp = var%grav(j) + mubar_tmp = mubar(j) + do i = dat%ng_1,dat%nq + k = i - dat%np + scale_height_i_tmp = (N_avo*k_boltz*var%temperature(j))/(grav_tmp*dat%species_mass(i)) + b1x2_tmp = binary_diffusion_param(dat%species_mass(i), mubar_tmp, var%temperature(j)) + gamma_i_gas_D(k,j) = (b1x2_tmp/den(j))*(1.0_dp/scale_height_i_tmp + (1.0_dp/var%temperature(j))*dTdz_tmp) + enddo + enddo + j = var%nz + dTdz_tmp = (var%temperature(var%nz) - var%temperature(var%nz-1))/var%dz(var%nz) + grav_tmp = var%grav(j) + mubar_tmp = mubar(j) + do i = dat%ng_1,dat%nq + k = i - dat%np + scale_height_i_tmp = (N_avo*k_boltz*var%temperature(j))/(grav_tmp*dat%species_mass(i)) + b1x2_tmp = binary_diffusion_param(dat%species_mass(i), mubar_tmp, var%temperature(j)) + gamma_i_gas_D(k,j) = (b1x2_tmp/den(j))*(1.0_dp/scale_height_i_tmp + (1.0_dp/var%temperature(j))*dTdz_tmp) + enddo + endblock; 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)) @@ -279,8 +307,15 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & scale_height_i_av = (N_avo*k_boltz*tav(j))/(grav_av*dat%species_mass(i)) b1x2av(k,j) = binary_diffusion_param(dat%species_mass(i), mubar_av, tav(j)) - gamma_i_gas_av(k,j) = eddav(j)*(1.0_dp/scale_height_av(j) + (1.0_dp/tav(j))*dTdz(j)) + & - (b1x2av(k,j)/denav(j))*(1.0_dp/scale_height_i_av + (1.0_dp/tav(j))*dTdz(j)) + gamma_i_gas_av(k,j) = eddav(j)*(1.0_dp/scale_height_av(j) + (1.0_dp/tav(j))*dTdz(j)) + + ! If we don't use upwind scheme for molecular diffusion, then we add the centered + ! molecular diffusion terms here. + if (.not.var%upwind_molec_diff) then + gamma_i_gas_av(k,j) = gamma_i_gas_av(k,j) + & + (b1x2av(k,j)/denav(j))*(1.0_dp/scale_height_i_av + (1.0_dp/tav(j))*dTdz(j)) + endif + enddo enddo @@ -298,6 +333,12 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & ADU(i,j) = gamma_i_gas_av(k,j)/(2.0_dp*var%dz(j)) ADL(i,j) = - gamma_i_gas_av(k,j-1)/(2.0_dp*var%dz(j)) ADD(i,j) = ADU(i,j) + ADL(i,j) + + if (var%upwind_molec_diff) then + ADU(i,j) = ADU(i,j) + gamma_i_gas_D(k,j+1)/var%dz(j) + ADL(i,j) = ADL(i,j) + 0.0_dp + ADD(i,j) = ADD(i,j) + (- gamma_i_gas_D(k,j)/var%dz(j)) + endif enddo enddo ! lower boundary @@ -309,6 +350,11 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & ADU(i,j) = gamma_i_gas_av(k,j)/(2.0_dp*var%dz(j)) ADD(i,j) = ADU(i,j) + + if (var%upwind_molec_diff) then + ADU(i,j) = ADU(i,j) + gamma_i_gas_D(k,j+1)/var%dz(j) + ADD(i,j) = ADD(i,j) + 0.0_dp + endif enddo ! upper boundary j = var%nz @@ -319,6 +365,11 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & ADL(i,j) = - gamma_i_gas_av(k,j-1)/(2.0_dp*var%dz(j)) ADD(i,j) = ADL(i,j) + + if (var%upwind_molec_diff) then + ADL(i,j) = ADL(i,j) + 0.0_dp + ADD(i,j) = ADD(i,j) + (- gamma_i_gas_D(k,j)/var%dz(j)) + endif enddo ! ! particles (eddy diffusion) diff --git a/src/photochem_types.f90 b/src/photochem_types.f90 index 35b9336..7549187 100644 --- a/src/photochem_types.f90 +++ b/src/photochem_types.f90 @@ -497,6 +497,10 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) integer :: verbose = 1 !! 0 == no printing. 1 == some printing. 2 == bunch of printing. !> Arbitrary rate that is fast (1/s). Used for keeping H2O at saturation in troposphere real(dp) :: fast_arbitrary_rate = 1.0e-2_dp + !> If True, then the code uses a 1st order upwind method for the advective molecular + !> diffusion terms instead of a centered scheme. This permits stability (at the cost + !> of accuracy) for atmospheres with strong molcular advection in the upper atmosphere. + logical :: upwind_molec_diff = .false. end type type :: SundialsData