set_temperature EvoAtmosphere

photochem/cython/EvoAtmosphere.pyx

@@ -263,6 +263,31 @@ cdef class EvoAtmosphere:
     ea_pxd.evoatmosphere_set_rate_fcn_wrapper(&self._ptr, species_c, fcn_c, err)
     if len(err.strip()) > 0:
        raise PhotoException(err.decode("utf-8").strip())
+
+  def set_temperature(self, ndarray[double, ndim=1] temperature, trop_alt = None):
+    """Changes the temperature profile.
+
+    Parameters
+    ----------
+    temperature : ndarray[double,ndim=1]
+        new temperature at each atomspheric layer
+    trop_alt : float, optional
+        Tropopause altitude (cm). Only necessary if rainout == True, 
+        or fix_water_in_trop == True.
+    """
+    cdef char err[ERR_LEN+1]
+    cdef int nz = temperature.size
+
+    cdef double trop_alt_ = 0.0
+    cdef bool trop_alt_present = False
+    if trop_alt != None:
+      trop_alt_present = True
+      trop_alt_ = trop_alt
+
+    ea_pxd.evoatmosphere_set_temperature_wrapper(&self._ptr, &nz, <double *>temperature.data, 
+                                       &trop_alt_, &trop_alt_present, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
 
   def regrid_prep_atmosphere(self, ndarray[double, ndim=2] usol, double top_atmos):
     """This subroutine calculates re-grids the model so that the top of the model domain 

photochem/cython/EvoAtmosphere_pxd.pxd

@@ -23,6 +23,8 @@ cdef extern void evoatmosphere_set_lower_bc_wrapper(void *ptr, char *species, ch
 cdef extern void evoatmosphere_set_upper_bc_wrapper(void *ptr, char *species, 
                                                     char *bc_type, double *veff, double *flux, bool *missing, char *err)
 cdef extern void evoatmosphere_set_rate_fcn_wrapper(void *ptr, char *species_c, time_dependent_rate_fcn fcn, char *err)
+cdef extern void evoatmosphere_set_temperature_wrapper(void *ptr, int *nz, double *temperature, 
+                                                    double *trop_alt, bool *trop_alt_present, char *err)
 cdef extern void evoatmosphere_regrid_prep_atmosphere_wrapper(void *ptr, int *nq, int *nz, double *usol, double *top_atmos, char *err)
 
 cdef extern void evoatmosphere_evolve_wrapper(void *ptr, char *filename, 

photochem/fortran/EvoAtmosphere_wrapper.f90

@@ -219,6 +219,32 @@ subroutine evoatmosphere_set_rate_fcn_wrapper(ptr, species_c, fcn_c, err) bind(c
     endif 
   end subroutine
 
+  subroutine evoatmosphere_set_temperature_wrapper(ptr, nz, temperature, &
+                                                trop_alt, trop_alt_present, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: nz
+    real(c_double), intent(in) :: temperature(nz)
+    real(c_double), intent(in) :: trop_alt
+    logical(c_bool), intent(in) :: trop_alt_present
+    character(kind=c_char), intent(out) :: err(err_len+1)
+
+    character(:), allocatable :: err_f
+    type(EvoAtmosphere), pointer :: pc
+
+    call c_f_pointer(ptr, pc)
+
+    if (trop_alt_present) then
+      call pc%set_temperature(temperature, trop_alt, err_f)
+    else
+      call pc%set_temperature(temperature, err=err_f)
+    endif
+    err(1) = c_null_char
+    if (allocated(err_f)) then
+      call copy_string_ftoc(err_f, err)
+    endif
+
+  end subroutine
+
   subroutine evoatmosphere_regrid_prep_atmosphere_wrapper(ptr, nq, nz, usol, top_atmos, err) bind(c)
     type(c_ptr), intent(in) :: ptr
     integer(c_int), intent(in) :: nq, nz

src/evoatmosphere/photochem_evoatmosphere.f90

@@ -60,6 +60,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     procedure :: set_lower_bc
     procedure :: set_upper_bc
     procedure :: set_rate_fcn
+    procedure :: set_temperature
     procedure :: rebin_update_vertical_grid
     procedure :: regrid_prep_atmosphere
 
@@ -259,14 +260,26 @@ module subroutine set_upper_bc(self, species, bc_type, veff, flux, err)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> Sets a function describing a custom rate for a species.
+    !> This could be useful for modeling external processes not in the
+    !> model.
     module subroutine set_rate_fcn(self, species, fcn, err)
       use photochem_types, only: time_dependent_rate_fcn
       class(EvoAtmosphere), target, intent(inout) :: self
-      character(*), intent(in) :: species
+      character(*), intent(in) :: species !! Species name
       procedure(time_dependent_rate_fcn), pointer :: fcn
       character(:), allocatable, intent(inout) :: err
     end subroutine
 
+    !> Changes the temperature profile.
+    module subroutine set_temperature(self, temperature, trop_alt, err)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      real(dp), intent(in) :: temperature(:) !! new temperature at each atomspheric layer
+      real(dp), optional, intent(in) :: trop_alt !! Tropopause altitude (cm). Only necessary if
+                                                 !! rainout == True, or fix_water_in_trop == True.
+      character(:), allocatable, intent(out) :: err
+    end subroutine
+
     module subroutine rebin_update_vertical_grid(self, usol_old, top_atmos, usol_new, err)
       class(EvoAtmosphere), target, intent(inout) :: self
       real(dp), intent(in) :: usol_old(:,:)

src/evoatmosphere/photochem_evoatmosphere_utils.f90

@@ -256,6 +256,80 @@ module subroutine set_rate_fcn(self, species, fcn, err)
 
   end subroutine
 
+  module subroutine set_temperature(self, temperature, trop_alt, err)
+    use photochem_input, only: interp2xsdata, compute_gibbs_energy
+
+    class(EvoAtmosphere), target, intent(inout) :: self
+    real(dp), intent(in) :: temperature(:)
+    real(dp), optional, intent(in) :: trop_alt
+    character(:), allocatable, intent(out) :: err
+
+    type(PhotochemData), pointer :: dat
+    type(PhotochemVars), pointer :: var
+    type(PhotochemVars) :: var_save
+
+    dat => self%dat
+    var => self%var
+
+    if (size(temperature) /= var%nz) then
+      err = "temperature has the wrong input dimension"
+      return
+    endif
+    if (self%evolve_climate) then
+      err = "You can not set the temperature when evolving climate"
+      return
+    endif
+
+    ! save in case there is an issue
+    var_save = var
+
+    var%temperature = temperature
+
+    ! xsections and gibbs energy needs updating
+    call interp2xsdata(dat, var, err)
+    if (allocated(err)) then
+      var = var_save
+      return
+    endif
+    if (dat%reverse) then
+      call compute_gibbs_energy(dat, var, err)
+      if (allocated(err)) then
+        var = var_save
+        return
+      endif
+    endif
+
+    ! if water is fixed in troposhere or gas rainout, and trop_alt present
+    ! then we need to change trop_ind, reallocate some stuff
+    ! in wrk, then we will re-prep the atmosphere
+    if ((dat%fix_water_in_trop .or. dat%gas_rainout) .and. present(trop_alt)) then
+      if (trop_alt < var%bottom_atmos .or. trop_alt > var%top_atmos) then
+        var = var_save
+        err = "trop_alt is above or bellow the atmosphere!"
+        return
+      endif
+
+      var%trop_alt = trop_alt
+      var%trop_ind = max(minloc(abs(var%z - var%trop_alt), 1) - 1, 1)
+
+      if (var%trop_ind < 3) then
+        var = var_save
+        err = 'Tropopause is too low.'
+        return
+      elseif (var%trop_ind > var%nz-2) then
+        var = var_save
+        err = 'Tropopause is too high.'
+        return
+      endif
+
+    endif
+
+    ! Fill wrk with new values
+    call self%prep_atmosphere(self%wrk%usol, err)
+    if (allocated(err)) return
+
+  end subroutine
+
   module subroutine rebin_update_vertical_grid(self, usol_old, top_atmos, usol_new, err)
     class(EvoAtmosphere), target, intent(inout) :: self
     real(dp), intent(in) :: usol_old(:,:)