diff --git a/CMakeLists.txt b/CMakeLists.txt
index b3d7599..868dd98 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION "3.14")
 cmake_policy(SET CMP0148 OLD)
-project(Photochem LANGUAGES Fortran C VERSION "0.5.4")
+project(Photochem LANGUAGES Fortran C VERSION "0.5.5")
 
 include(FortranCInterface)
 FortranCInterface_VERIFY()
diff --git a/photochem/cython/Atmosphere.pyx b/photochem/cython/Atmosphere.pyx
index 2435270..420dd6a 100644
--- a/photochem/cython/Atmosphere.pyx
+++ b/photochem/cython/Atmosphere.pyx
@@ -3,7 +3,7 @@ cimport Atmosphere_pxd as a_pxd
   
 cdef class Atmosphere:
   """A photochemical model which assumes a background gas (e.g. N2). Once initialized,
-  this class can integrate an atmosphere forward in time, and can investigate
+  this class can integrate an atmosphere forward in time and can investigate
   the reactions producing and destroying each molecule.
   """
 
@@ -27,9 +27,26 @@ cdef class Atmosphere:
       raise PhotoException('The "__init__" method of Atmosphere has not been called.')
     PyObject_GenericSetAttr(self, name, value)
   
-  def __init__(self, mechanism_file = None, settings_file = None, 
-               flux_file = None, atmosphere_txt = None, data_dir = None):           
-    
+  def __init__(self, str mechanism_file, str settings_file, 
+               str flux_file, str atmosphere_txt, data_dir = None):           
+    """Initializes the photochemical model.
+
+    Parameters
+    ----------
+    mechanism_file : str
+        Path to the reaction mechanism file (yaml format).
+    settings_file : str
+        Path to the settings file (yaml format).
+    flux_file : str
+        Path to the file describing the stellar flux.
+    atmosphere_txt : str
+        Path to the file containing altitude, total number density, temperature, 
+        eddy diffusion, initial concentrations of each gas (mixing ratios), 
+        and particle radii.
+    data_dir : str, optional
+        Path to the data directory containing photolysis cross sections and other data
+        needed to run the model
+    """
     self._init_called = True
 
     if data_dir == None:
@@ -142,7 +159,7 @@ cdef class Atmosphere:
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
       
-  def out2atmosphere_txt(self,filename = None, int number_of_decimals=5, bool overwrite = False, bool clip = True):
+  def out2atmosphere_txt(self,str filename, int number_of_decimals=5, bool overwrite = False, bool clip = True):
     """Saves state of the atmosphere using the mixing ratios in self.wrk.usol.
 
     Parameters
@@ -201,7 +218,7 @@ cdef class Atmosphere:
     -------
     tuple
         First element are the surface fluxes, and the second are top-of-atmosphere
-        fluxes.
+        fluxes. Units are molecules/cm^2/s
     """
     cdef ndarray surf_fluxes = np.empty(self.dat.nq, np.double)
     cdef ndarray top_fluxes = np.empty(self.dat.nq, np.double)
@@ -220,7 +237,7 @@ cdef class Atmosphere:
       top[names[i]] = top_fluxes[i]
     return surface, top
     
-  def set_lower_bc(self, species = None, bc_type = None, vdep = None, mix = None,
+  def set_lower_bc(self, str species, str bc_type, vdep = None, mix = None,
                             flux = None, height = None):
     """Sets a lower boundary condition.
 
@@ -279,7 +296,7 @@ cdef class Atmosphere:
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
   
-  def set_upper_bc(self, species = None, bc_type = None, veff = None,flux = None):
+  def set_upper_bc(self, str species, str bc_type, veff = None, flux = None):
     """Sets upper boundary condition.
 
     Parameters
diff --git a/photochem/cython/EvoAtmosphere.pyx b/photochem/cython/EvoAtmosphere.pyx
index ef7cefb..eeb804e 100644
--- a/photochem/cython/EvoAtmosphere.pyx
+++ b/photochem/cython/EvoAtmosphere.pyx
@@ -3,9 +3,9 @@ cimport EvoAtmosphere_pxd as ea_pxd
 
 cdef class EvoAtmosphere:
   """A photochemical model which assumes no background gas. Once initialized,
-  this class can integrate an atmosphere forward in time, and can investigate
-  the reactions producing and destroying each molecule. The model can also
-  optionally self-consistently evolve climate.
+  this class can integrate an atmosphere forward in time to a steady
+  state. The model can also, optionally, self-consistently evolve climate for
+  terrestrial planets.
   """
 
   cdef ea_pxd.EvoAtmosphere *_ptr
@@ -28,9 +28,26 @@ cdef class EvoAtmosphere:
       raise PhotoException('The "__init__" method of EvoAtmosphere has not been called.')
     PyObject_GenericSetAttr(self, name, value)
 
-  def __init__(self, mechanism_file = None, settings_file = None, 
-               flux_file = None, atmosphere_txt = None, data_dir = None):           
+  def __init__(self, str mechanism_file, str settings_file, 
+               str flux_file, str atmosphere_txt, data_dir = None):           
+    """Initializes the photochemical model.
 
+    Parameters
+    ----------
+    mechanism_file : str
+        Path to the reaction mechanism file (yaml format).
+    settings_file : str
+        Path to the settings file (yaml format).
+    flux_file : str
+        Path to the file describing the stellar flux.
+    atmosphere_txt : str
+        Path to the file containing altitude, total number density, temperature, 
+        eddy diffusion, initial concentrations of each gas (mixing ratios), 
+        and particle radii.
+    data_dir : str, optional
+        Path to the data directory containing photolysis cross sections and other data
+        needed to run the model
+    """
     self._init_called = True
 
     if data_dir == None:
@@ -94,7 +111,7 @@ cdef class EvoAtmosphere:
     Parameters
     ----------
     usol : ndarray[double,ndim=2]
-        densities
+        Number densities (molecules/cm^3)
     """
     cdef char err[ERR_LEN+1]
     cdef int nq = self.dat.nq
@@ -107,8 +124,8 @@ cdef class EvoAtmosphere:
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
 
-  def out2atmosphere_txt(self,filename = None, int number_of_decimals=5, bool overwrite = False, bool clip = True):
-    """Saves state of the atmosphere using the mixing ratios in self.wrk.usol.
+  def out2atmosphere_txt(self,str filename, int number_of_decimals=5, bool overwrite = False, bool clip = True):
+    """Saves state of the atmosphere using the concentrations in self.wrk.usol.
 
     Parameters
     ----------
@@ -138,7 +155,7 @@ cdef class EvoAtmosphere:
     -------
     tuple
         First element are the surface fluxes, and the second are top-of-atmosphere
-        fluxes.
+        fluxes. Units molecules/cm^2/s.
     """
     cdef ndarray surf_fluxes = np.empty(self.dat.nq, np.double)
     cdef ndarray top_fluxes = np.empty(self.dat.nq, np.double)
@@ -157,7 +174,7 @@ cdef class EvoAtmosphere:
       top[names[i]] = top_fluxes[i]
     return surface, top
 
-  def set_lower_bc(self, species = None, bc_type = None, vdep = None, den = None, press = None,
+  def set_lower_bc(self, str species, str bc_type, vdep = None, den = None, press = None,
                             flux = None, height = None):
     """Sets a lower boundary condition.
 
@@ -224,7 +241,7 @@ cdef class EvoAtmosphere:
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
   
-  def set_upper_bc(self, species = None, bc_type = None, veff = None,flux = None):
+  def set_upper_bc(self, str species, str bc_type, veff = None,flux = None):
     """Sets upper boundary condition.
 
     Parameters
@@ -437,6 +454,8 @@ cdef class EvoAtmosphere:
         times to evaluate the solution
     overwrite : bool
         If true, then overwrites pre-existing files with `filename`
+    restart_from_file : bool
+        If true, then the integration restarts from the input file.
 
     Returns
     -------
@@ -474,7 +493,7 @@ cdef class EvoAtmosphere:
     Parameters
     ----------
     usol_start : ndarray[double,ndim=2]
-        Initial mixing ratios
+        Initial number densities (molecules/cm^3)
     """   
     cdef char err[ERR_LEN+1]
     cdef int nq = self.dat.nq
diff --git a/src/atmosphere/photochem_atmosphere.f90 b/src/atmosphere/photochem_atmosphere.f90
index 9803015..aa061bb 100644
--- a/src/atmosphere/photochem_atmosphere.f90
+++ b/src/atmosphere/photochem_atmosphere.f90
@@ -169,11 +169,11 @@ module subroutine initialize_stepper(self, usol_start, err)
     end subroutine
     
     !> Takes one internal integration step. Function `initialize_stepper`
-    !> must have been called befe this
+    !> must have been called before this
     module function step(self, err) result(tn)
       class(Atmosphere), target, intent(inout) :: self
       character(:), allocatable, intent(out) :: err
-      real(dp) :: tn
+      real(dp) :: tn !! Current time in the integration.
     end function
     
     !> Deallocates memory created during `initialize_stepper`
diff --git a/src/dependencies/CMakeLists.txt b/src/dependencies/CMakeLists.txt
index ad0e6b5..1814acd 100644
--- a/src/dependencies/CMakeLists.txt
+++ b/src/dependencies/CMakeLists.txt
@@ -1,9 +1,9 @@
 
 CPMAddPackage(
   NAME futils
-  VERSION 0.1.9
+  VERSION 0.1.11
   GITHUB_REPOSITORY "Nicholaswogan/futils"
-  GIT_TAG "v0.1.9"
+  GIT_TAG "v0.1.11"
   EXCLUDE_FROM_ALL ON
 )
 
@@ -49,7 +49,7 @@ CPMAddPackage(
 
 CPMAddPackage(
   NAME clima 
-  VERSION 0.4.5
+  VERSION 0.4.6
   OPTIONS 
     "BUILD_EXECUTABLE OFF"
     "BUILD_WITH_OPENMP ${BUILD_WITH_OPENMP}"
@@ -57,7 +57,7 @@ CPMAddPackage(
     "BUILD_PYTHON_CLIMA ${BUILD_PYTHON_PHOTOCHEM}"
     "PYTHON_CLIMA_DESTINATION photochem"
   GITHUB_REPOSITORY "Nicholaswogan/clima"
-  GIT_TAG "v0.4.5"
+  GIT_TAG "v0.4.6"
   EXCLUDE_FROM_ALL OFF
 )
 
diff --git a/src/evoatmosphere/photochem_evoatmosphere.f90 b/src/evoatmosphere/photochem_evoatmosphere.f90
index a085238..51a5832 100644
--- a/src/evoatmosphere/photochem_evoatmosphere.f90
+++ b/src/evoatmosphere/photochem_evoatmosphere.f90
@@ -39,7 +39,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
 
   contains
 
-    !!! photochem_evoatmosphere_rhs.f90 !!!
+    !~~ photochem_evoatmosphere_rhs.f90 ~~!
     procedure :: set_trop_ind
     procedure :: prep_atm_evo_gas
     procedure :: prep_atmosphere => prep_all_evo_gas
@@ -48,14 +48,14 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     procedure :: right_hand_side => rhs_evo_gas
     procedure :: jacobian => jac_evo_gas
 
-    !!! photochem_evoatmosphere_integrate.f90 !!!
+    !~~ photochem_evoatmosphere_integrate.f90 ~~!
     procedure :: evolve
     procedure :: check_for_convergence
     procedure :: initialize_stepper
     procedure :: step
     procedure :: destroy_stepper
 
-    !!! photochem_evoatmosphere_utils.f90 !!!
+    !~~ photochem_evoatmosphere_utils.f90 ~~!
     procedure :: out2atmosphere_txt
     procedure :: gas_fluxes
     procedure :: set_lower_bc
@@ -73,23 +73,28 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
   end interface
 
   interface
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    !!! photochem_evoatmosphere_init.f90 !!!
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+    !~~ photochem_evoatmosphere_init.f90 ~~!
+
+    !> Initializes the photochemical model.
     module function create_EvoAtmosphere(mechanism_file, &
                                          settings_file, flux_file, atmosphere_txt, data_dir, err) result(self)
-      character(len=*), intent(in) :: mechanism_file
-      character(len=*), intent(in) :: settings_file
-      character(len=*), intent(in) :: flux_file
+      character(len=*), intent(in) :: mechanism_file !! Path to the reaction mechanism file (yaml format).
+      character(len=*), intent(in) :: settings_file !! Path to the settings file (yaml format).
+      character(len=*), intent(in) :: flux_file !! Path to the file describing the stellar flux.
+      !> Path to the file containing altitude, total number density, temperature, 
+      !> eddy diffusion, initial concentrations of each gas (mixing ratios), 
+      !> and particle radii.
       character(len=*), intent(in) :: atmosphere_txt
+      !> Path to the data directory containing photolysis cross sections and other data
+      !> needed to run the model
       character(len=*), intent(in) :: data_dir
       character(:), allocatable, intent(out) :: err
       type(EvoAtmosphere) :: self
     end function
 
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    !!! photochem_atmosphere_rhs.f90 !!!
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    !~~ photochem_atmosphere_rhs.f90 ~~!
+
     module subroutine set_trop_ind(self, usol_in, err)
       class(EvoAtmosphere), target, intent(inout) :: self
       real(dp), intent(in) :: usol_in(:,:)
@@ -108,9 +113,13 @@ module subroutine prep_atm_evo_gas(self, usol_in, usol, &
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> Given `usol`, the densities of each species in the atmosphere,
+    !> this subroutine calculates reaction rates, photolysis rates, etc.
+    !> and puts this information into self.wrk. self.wrk contains all the
+    !> information needed for `dochem` to compute chemistry.
     module subroutine prep_all_evo_gas(self, usol_in, err)
       class(EvoAtmosphere), target, intent(inout) :: self
-      real(dp), intent(in) :: usol_in(:,:)
+      real(dp), intent(in) :: usol_in(:,:) !! Number densities (molecules/cm^3)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
@@ -121,6 +130,8 @@ module subroutine right_hand_side_chem(self, usol, rhs, err)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> Computes the right-hand-side of the ODEs describing atmospheric chemistry
+    !> and transport.
     module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err)
       class(EvoAtmosphere), target, intent(inout) :: self
       integer, intent(in) :: neqs
@@ -128,17 +139,15 @@ module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err)
       real(dp), target, intent(in) :: usol_flat(neqs)
       real(dp), intent(out) :: rhs(neqs)
       character(:), allocatable, intent(out) :: err
-      ! Computes the right-hand-side of the ODEs describing atmospheric chemistry
-      ! and transport.
     end subroutine
     
+    !> The jacobian of the rhs_background_gas.
     module subroutine jac_evo_gas(self, lda_neqs, neqs, usol_flat, jac, err)
       class(EvoAtmosphere), target, intent(inout) :: self
       integer, intent(in) :: lda_neqs, neqs
       real(dp), target, intent(in) :: usol_flat(neqs)
       real(dp), intent(out), target :: jac(lda_neqs)
       character(:), allocatable, intent(out) :: err
-      ! The jacobian of the rhs_background_gas.
     end subroutine
 
     module subroutine production_and_loss(self, species, usol, pl, err)  
@@ -150,21 +159,20 @@ module subroutine production_and_loss(self, species, usol, pl, err)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    !!! photochem_evoatmosphere_integrate.f90 !!!
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    !~~ photochem_evoatmosphere_integrate.f90 ~~!
+
+    !> Evolve atmosphere through time, and saves output in a binary Fortran file.
     module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, restart_from_file, err) result(success)
       use, intrinsic :: iso_c_binding
       class(EvoAtmosphere), target, intent(inout) :: self
-      character(len=*), intent(in) :: filename
-      real(c_double), intent(inout) :: tstart
-      real(dp), intent(inout) :: usol_start(:,:)
-      real(c_double), intent(in) :: t_eval(:)
-      logical, optional, intent(in) :: overwrite
-      logical, optional, intent(in) :: restart_from_file
-      logical :: success
+      character(len=*), intent(in) :: filename !! Filename to save results.
+      real(c_double), intent(inout) :: tstart !! start time in seconds
+      real(dp), intent(inout) :: usol_start(:,:) !! Initial number densities (molecules/cm^3)
+      real(c_double), intent(in) :: t_eval(:) !! times to evaluate the solution
+      logical, optional, intent(in) :: overwrite !! If true, then overwrites pre-existing files with `filename`
+      logical, optional, intent(in) :: restart_from_file !! If true, then the integration restarts from the input file.
+      logical :: success !! If True, then integration was successful.
       character(:), allocatable, intent(out) :: err
-      ! Evolve atmosphere through time, and saves output in a binary Fortran file.
     end function
 
     !> Determines if integration has converged to photochemical steady-state.
@@ -177,7 +185,7 @@ module function check_for_convergence(self, err) result(converged)
     !> Initializes an integration starting at `usol_start`
     module subroutine initialize_stepper(self, usol_start, err)      
       class(EvoAtmosphere), target, intent(inout) :: self
-      real(dp), intent(in) :: usol_start(:,:) !! Initial mixing ratios
+      real(dp), intent(in) :: usol_start(:,:) !! Initial number densities (molecules/cm^3)
       character(:), allocatable, intent(out) :: err
     end subroutine
     
@@ -186,7 +194,7 @@ module subroutine initialize_stepper(self, usol_start, err)
     module function step(self, err) result(tn)
       class(EvoAtmosphere), target, intent(inout) :: self
       character(:), allocatable, intent(out) :: err
-      real(dp) :: tn
+      real(dp) :: tn !! Current time in the integration.
     end function
     
     !> Deallocates memory created during `initialize_stepper`
@@ -224,16 +232,18 @@ module function JacFn_evo(tn, sunvec_y, sunvec_f, sunmat_J, user_data, &
       integer(c_int)        :: ierr
     end function
 
+    !~~ photochem_evoatmosphere_utils.f90 ~~!
 
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    !!! photochem_evoatmosphere_utils.f90 !!!
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    !> Saves state of the atmosphere using the concentrations in self%wrk%usol.
     module subroutine out2atmosphere_txt(self, filename, number_of_decimals, overwrite, clip, err)
       use photochem_common, only: out2atmosphere_txt_base
       class(EvoAtmosphere), target, intent(inout) :: self
-      character(len=*), intent(in) :: filename
-      integer, intent(in) :: number_of_decimals
-      logical, intent(in) :: overwrite, clip
+      character(len=*), intent(in) :: filename !! Output filename
+      integer, intent(in) :: number_of_decimals !! Number of decimals
+      logical, intent(in) :: overwrite !! If true, then output file can be overwritten, by default False
+      !> If true, then mixing ratios are clipped at a very small 
+      !> positive number, by default False
+      logical, intent(in) :: clip
       character(:), allocatable, intent(out) :: err
     end subroutine
 
@@ -325,10 +335,14 @@ module subroutine rebin_update_vertical_grid(self, usol_old, top_atmos, usol_new
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> This subroutine calculates re-grids the model so that the top of the model domain 
+    !> is at `top_atmos` the computes reaction rates, photolysis rates, etc.
+    !> and puts this information into self.wrk. self.wrk contains all the
+    !> information needed for `dochem` to compute chemistry.
     module subroutine regrid_prep_atmosphere(self, usol_new, top_atmos, err)
       class(EvoAtmosphere), target, intent(inout) :: self
-      real(dp), intent(in) :: usol_new(:,:)
-      real(dp), intent(in) :: top_atmos
+      real(dp), intent(in) :: usol_new(:,:) !! The number densities (molecules/cm^3)
+      real(dp), intent(in) :: top_atmos !! The top of the model domain (cm)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
diff --git a/src/evoatmosphere/photochem_evoatmosphere_utils.f90 b/src/evoatmosphere/photochem_evoatmosphere_utils.f90
index 7d253df..2a03072 100644
--- a/src/evoatmosphere/photochem_evoatmosphere_utils.f90
+++ b/src/evoatmosphere/photochem_evoatmosphere_utils.f90
@@ -9,7 +9,8 @@ module subroutine out2atmosphere_txt(self, filename, number_of_decimals, overwri
     class(EvoAtmosphere), target, intent(inout) :: self
     character(len=*), intent(in) :: filename
     integer, intent(in) :: number_of_decimals
-    logical, intent(in) :: overwrite, clip
+    logical, intent(in) :: overwrite
+    logical, intent(in) :: clip
     character(:), allocatable, intent(out) :: err
     
     real(dp) :: rhs(self%var%neqs)  
diff --git a/src/input/photochem_input_read.f90 b/src/input/photochem_input_read.f90
index a63d621..e0e0767 100644
--- a/src/input/photochem_input_read.f90
+++ b/src/input/photochem_input_read.f90
@@ -2002,7 +2002,7 @@ subroutine get_photorad(dat, var, err)
   ! also returns the radii of particles for that file.
   subroutine read_mie_data_file(filename, nw, wavl, &
                                  nrad_file, radii_file, w0_file, qext_file, g_file, err)
-    use futils, only: addpnt, inter2
+    use futils, only: addpnt, inter2, FileCloser
     
     character(len=*), intent(in) :: filename
     integer, intent(in) :: nw
@@ -2020,9 +2020,11 @@ subroutine read_mie_data_file(filename, nw, wavl, &
     integer :: nw_tmp
     real(dp) :: dum
     integer :: i, j, io, ierr
+    type(FileCloser) :: file
     
     
     open(2,file=filename,form="unformatted",status='old',iostat=io)
+    file%unit = 2
     if (io /= 0) then
       err = "Was unable to open mie data file "//trim(filename)
       return
@@ -2079,7 +2081,6 @@ subroutine read_mie_data_file(filename, nw, wavl, &
             ". Should have reached the end of the file, but did not."
       return
     endif
-    close(2)
     
     ! now lets interpolate a bunch
     allocate(w0_file(nrad_file,nw))
@@ -2504,7 +2505,7 @@ subroutine rayleigh_vardavas(A, B, Delta, lambda, sigray)
   end subroutine
   
   subroutine read_stellar_flux(star_file, nw, wavl, photon_flux, err)
-    use futils, only: inter2, addpnt
+    use futils, only: inter2, addpnt, FileCloser
     use photochem_const, only: c_light, plank
     
     character(len=*), intent(in) :: star_file
@@ -2518,8 +2519,10 @@ subroutine read_stellar_flux(star_file, nw, wavl, photon_flux, err)
     real(dp) :: dum1, dum2
     integer :: io, i, n, ierr
     real(dp), parameter :: rdelta = 1.0e-4_dp
+    type(FileCloser) :: file
     
     open(1,file=star_file,status='old',iostat=io)
+    file%unit = 1
     if (io /= 0) then
       err = "The input file "//star_file//' does not exist.'
       return
@@ -2545,7 +2548,6 @@ subroutine read_stellar_flux(star_file, nw, wavl, photon_flux, err)
         return
       endif
     enddo
-    close(1)
     
     i = n
     ! interpolate 
@@ -2573,6 +2575,7 @@ subroutine read_stellar_flux(star_file, nw, wavl, photon_flux, err)
   end subroutine
   
   subroutine read_atmosphere_file(atmosphere_txt, dat, var, err)
+    use futils, only: FileCloser
     character(len=*), intent(in) :: atmosphere_txt
     type(PhotochemData), intent(inout) :: dat
     type(PhotochemVars), intent(inout) :: var
@@ -2585,8 +2588,10 @@ subroutine read_atmosphere_file(atmosphere_txt, dat, var, err)
     integer :: ind(1)
     real(dp), allocatable :: temp(:,:)
     integer :: io, i, n, nn, ii
+    type(FileCloser) :: file
     
     open(4, file=trim(atmosphere_txt),status='old',iostat=io)
+    file%unit = 4
     if (io /= 0) then
       err = 'Can not open file '//trim(atmosphere_txt)
       return
@@ -2648,7 +2653,6 @@ subroutine read_atmosphere_file(atmosphere_txt, dat, var, err)
         return
       endif
     enddo
-    close(4)
     
     ! reads in mixing ratios
     do i=1,dat%nq
diff --git a/src/photochem_common.f90 b/src/photochem_common.f90
index 5d67f85..ecc65a3 100644
--- a/src/photochem_common.f90
+++ b/src/photochem_common.f90
@@ -620,6 +620,7 @@ subroutine molec_per_particle(dat, var, molecules_per_particle)
   subroutine out2atmosphere_txt_base(dat, var, &
                                      pressure, density, densities, molecules_per_particle, &
                                      filename, number_of_decimals, overwrite, clip, err)
+    use futils, only: FileCloser
     type(PhotochemData), pointer :: dat
     type(PhotochemVars), pointer :: var
     real(dp), intent(in) :: pressure(:), density(:), densities(:,:), molecules_per_particle(:,:)
@@ -634,6 +635,7 @@ subroutine out2atmosphere_txt_base(dat, var, &
     character(len=10) :: number_of_decimals_str, number_of_spaces_str
     character(:), allocatable :: fmt_label, fmt_number
     real(dp) :: val, clip_value
+    type(FileCloser) :: file
 
     if (clip) then
       clip_value = 1.0e-40_dp
@@ -643,12 +645,14 @@ subroutine out2atmosphere_txt_base(dat, var, &
     
     if (overwrite) then
       open(1, file=filename, form='formatted', status='replace', iostat=io)
+      file%unit = 1
       if (io /= 0) then
         err = "Unable to overwrite file "//trim(filename)
         return
       endif
     else
       open(1, file=filename, form='formatted', status='new', iostat=io)
+      file%unit = 1
       if (io /= 0) then
         err = "Unable to create file "//trim(filename)//" because it already exists"
         return
@@ -730,8 +734,6 @@ subroutine out2atmosphere_txt_base(dat, var, &
         enddo
       endif
     enddo
-    
-    close(1)
 
   end subroutine