optional upwind scheme for molecular advection for stability

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)
+
 
 

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)
+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)

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
 

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)

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