diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py
index 6d4f86dc..16f97d7e 100644
--- a/hyperion/conf/conf_files.py
+++ b/hyperion/conf/conf_files.py
@@ -18,6 +18,7 @@ def __init__(self):
         self.output_density = 'none'
         self.output_density_diff = 'none'
         self.output_specific_energy = 'last'
+        self.output_specific_energy_nu = 'last'
         self.output_n_photons = 'none'
         self._freeze()
 
@@ -27,6 +28,7 @@ def read(cls, group):
         self.output_density = group.attrs['output_density'].decode('utf-8')
         self.output_density_diff = group.attrs['output_density_diff'].decode('utf-8')
         self.output_specific_energy = group.attrs['output_specific_energy'].decode('utf-8')
+        self.output_specific_energy_nu = group.attrs['output_specific_energy_nu'].decode('utf-8')
         self.output_n_photons = group.attrs['output_n_photons'].decode('utf-8')
         return self
 
@@ -34,6 +36,7 @@ def write(self, group):
         group.attrs['output_density'] = np.string_(self.output_density.encode('utf-8'))
         group.attrs['output_density_diff'] = np.string_(self.output_density_diff.encode('utf-8'))
         group.attrs['output_specific_energy'] = np.string_(self.output_specific_energy.encode('utf-8'))
+        group.attrs['output_specific_energy_nu'] = np.string_(self.output_specific_energy_nu.encode('utf-8'))
         group.attrs['output_n_photons'] = np.string_(self.output_n_photons.encode('utf-8'))
 
 
@@ -52,6 +55,8 @@ def __init__(self):
         self.set_max_reabsorptions(1000000)
         self.set_pda(False)
         self.set_mrw(False)
+        self.compute_isrf(False)
+
         self.set_convergence(False)
         self.set_kill_on_absorb(False)
         self.set_kill_on_scatter(False)
@@ -391,6 +396,25 @@ def _read_pda(self, group):
     def _write_pda(self, group):
         group.attrs['pda'] = bool2str(self.pda)
 
+    def _read_isrf(self,group):
+        self.isrf = str2bool(group.attrs['isrf'])
+        
+    def _write_isrf(self,group):
+        group.attrs['isrf'] = bool2str(self.isrf)
+
+    def compute_isrf(self,isrf):
+
+        '''
+
+        Set whether or not to compute and save the ISRF in each cell
+
+        If enabled, the ISRF is computed in every cell at the
+        frequencies of the dust opacity tables.
+
+        '''
+        self.isrf = isrf
+
+
     def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True):
         '''
         Set whether to use the Modified Random Walk (MRW) approximation
@@ -719,6 +743,7 @@ def read_run_conf(self, group):  # not a class method because inherited
         self._read_max_reabsorptions(group)
         self._read_pda(group)
         self._read_mrw(group)
+        self._reada_isrf(group)
         self._read_convergence(group)
         self._read_kill_on_absorb(group)
         self._read_kill_on_scatter(group)
@@ -747,6 +772,7 @@ def write_run_conf(self, group):
         self._write_max_reabsorptions(group)
         self._write_pda(group)
         self._write_mrw(group)
+        self._write_isrf(group)
         self._write_convergence(group)
         self._write_kill_on_absorb(group)
         self._write_kill_on_scatter(group)
diff --git a/hyperion/conf/tests/test_conf_io.py b/hyperion/conf/tests/test_conf_io.py
index da688d32..738ef480 100644
--- a/hyperion/conf/tests/test_conf_io.py
+++ b/hyperion/conf/tests/test_conf_io.py
@@ -249,6 +249,16 @@ def test_io_run_conf_mrw(value):
     r2.read_run_conf(v)
     assert r2.mrw == r1.mrw
 
+@pytest.mark.parametrize(('value'), [True, False])
+def test_io_run_conf_isrf(value):
+    r1 = RunConf()
+    r1.compute_isrf(value)
+    r1.set_n_photons(1, 2)
+    v = virtual_file()
+    r1.write_run_conf(v)
+    r2 = RunConf()
+    r2.read_run_conf(v)
+    assert r2.isrf == r1.isrf
 
 @pytest.mark.parametrize(('value'), [False, True])
 def test_io_run_conf_convergence(value):
diff --git a/hyperion/grid/grid_helpers.py b/hyperion/grid/grid_helpers.py
index cadb3496..55194e84 100644
--- a/hyperion/grid/grid_helpers.py
+++ b/hyperion/grid/grid_helpers.py
@@ -22,7 +22,7 @@ def single_grid_dims(data, ndim=3):
     '''
 
     if type(data) in [list, tuple]:
-
+        
         n_pop = len(data)
         shape = None
         for item in data:
@@ -34,14 +34,22 @@ def single_grid_dims(data, ndim=3):
         if shape is not None and len(shape) != ndim:
             raise ValueError("Grids should be %i-dimensional" % ndim)
 
-    elif isinstance(data, np.ndarray):
 
+
+    elif isinstance(data, np.ndarray):
         if data.ndim == ndim:
             n_pop = None
             shape = data.shape
         elif data.ndim == ndim + 1:
             n_pop = data.shape[0]
             shape = data[0].shape
+            
+        #DN CRAZY ADDITIONS
+        elif data.ndim == ndim+2:
+            n_pop = data.shape[1]
+            shape = data.shape[-1]
+            shape = tuple(map(int,str(shape).split(' ')))
+
         else:
             raise Exception("Unexpected number of dimensions: %i" % data.ndim)
 
diff --git a/hyperion/model/model.py b/hyperion/model/model.py
index 2a0ce3d4..b5e1442c 100644
--- a/hyperion/model/model.py
+++ b/hyperion/model/model.py
@@ -223,7 +223,7 @@ def use_geometry(self, filename):
         # Close the file
         f.close()
 
-    def use_quantities(self, filename, quantities=['density', 'specific_energy'],
+    def use_quantities(self, filename, quantities=['density', 'specific_energy', 'specific_energy_nu'],
                        use_minimum_specific_energy=True, use_dust=True, copy=True,
                        only_initial=False):
         '''
@@ -297,6 +297,14 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy'],
                 else:
                     quantities_path['specific_energy'] = last_iteration
 
+            if self.compute_isrf == True:
+                if 'specific_energy_nu' in quantities:
+                    if only_initial or last_iteration is None:
+                        if 'specific_energy_nu' in f['/Input/Grid/Quantities']:
+                            quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities'
+                    else:
+                        quantities_path['specific_energy_nu'] = last_iteration
+
             # Minimum specific energy
             if use_minimum_specific_energy:
                 minimum_specific_energy_path = '/Input/Grid/Quantities'
@@ -313,6 +321,11 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy'],
                 if 'specific_energy' in f['/Grid/Quantities']:
                     quantities_path['specific_energy'] = '/Grid/Quantities'
 
+            if self.compute_isrf == True:
+                if 'specific_energy_nu' in quantities:
+                    if 'specific_energy_nu' in f['/Grid/Quantities']:
+                        quantities_path['specific_energy_nu'] = '/Grid/Quantities'
+                        
             # Minimum specific energy
             if use_minimum_specific_energy:
                 minimum_specific_energy_path = '/Grid/Quantities'
@@ -808,6 +821,9 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...).
             self.grid['density'] = []
             if specific_energy is not None:
                 self.grid['specific_energy'] = []
+                
+                if self.compute_isrf == True:
+                    self.grid['specific_energy_nu'] = []
 
         # Check whether the density can be added to an existing one
         if merge_if_possible:
@@ -849,6 +865,9 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...).
         # Set specific energy if specified
         if specific_energy is not None:
             self.grid['specific_energy'].append(specific_energy)
+            
+            if self.compute_isrf == True:
+                self.grid['specific_energy_nu'].append(specific_energy)
 
     def set_cartesian_grid(self, x_wall, y_wall, z_wall):
         self.set_grid(CartesianGrid(x_wall, y_wall, z_wall))
diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90
index d050bca0..2dc2457c 100644
--- a/src/grid/grid_generic.f90
+++ b/src/grid/grid_generic.f90
@@ -3,11 +3,14 @@ module grid_generic
   use core_lib, only : sp, dp, hid_t, warn, error
   use mpi_hdf5_io, only : mp_create_group
   use mpi_core, only : main_process
-
-  use grid_io, only : write_grid_3d, write_grid_4d
+  
+  use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d
   use grid_geometry, only : geo
-  use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy, density, density_original
-  use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type
+  use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, density, density_original
+  use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf
+
+  use dust_main, only: d
+  
 
   implicit none
   save
@@ -23,6 +26,7 @@ subroutine grid_reset_energy()
     if(allocated(n_photons)) n_photons = 0
     if(allocated(last_photon_id)) last_photon_id = 0
     if(allocated(specific_energy_sum)) specific_energy_sum = 0._dp
+    if(allocated(specific_energy_sum_nu)) specific_energy_sum_nu = 0._dp
   end subroutine grid_reset_energy
 
   subroutine output_grid(group, iter, n_iter)
@@ -31,6 +35,10 @@ subroutine output_grid(group, iter, n_iter)
 
     integer(hid_t),intent(in) :: group
     integer,intent(in) :: iter, n_iter
+    integer :: n_cells, n_dust, n_isrf_lam
+    integer :: i,j,k
+    real, dimension(d(1)%n_nu) :: energy_frequency_bins
+
 
     if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")')
 
@@ -61,6 +69,58 @@ subroutine output_grid(group, iter, n_iter)
        end if
     end if
 
+
+    
+    ! WRITE THE ISRF
+    if (compute_isrf) then 
+       
+       n_cells = size(specific_energy_nu, 1)
+       n_dust = size(specific_energy_nu, 2)
+       n_isrf_lam = size(specific_energy_nu,3)
+       
+       do i=1,d(1)%n_nu
+          energy_frequency_bins(i) = d(1)%nu(i)
+       end do
+       
+       if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then
+          !if(allocated(energy_frequency_bins)) then
+          call write_grid_3d(group, 'ISRF_frequency_bins',energy_frequency_bins, geo)
+          !else
+          !call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated")
+          !end if
+       end if
+       
+       if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then
+          if(allocated(specific_energy_nu)) then
+             
+
+             if (compute_isrf) then
+                print *,'[GRID_GENERIC.F90] ENTERING BLOCK FOR WRITING SPECIFIC_ENERGY_NU'
+                
+                select case(physics_io_type)
+                   
+                case(sp)
+                   print *,'[GRID_GENERIC.F90] diving into write_grid_5d sp'
+                   call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo)
+
+                case(dp)
+                   print *,'[GRID_GENERIC.F90] diving into write_grid_5d dp'
+                   call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo)
+                   print *,'[GRID_GENERIC.F90] finished with write_grid_5d dp'
+
+                case default
+                   call error("output_grid","unexpected value of physics_io_type (should be sp or dp)")
+                end select
+             else
+                call warn("output_grid","specific_energy_nu array is not allocated")
+
+             end if
+          end if
+       end if
+       
+    endif
+
+
     ! DENSITY
 
     if(trim(output_density)=='all' .or. (trim(output_density)=='last'.and.iter==n_iter)) then
diff --git a/src/grid/grid_io.f90 b/src/grid/grid_io.f90
index 8cfc82f5..0cba3f1d 100644
--- a/src/grid/grid_io.f90
+++ b/src/grid/grid_io.f90
@@ -14,6 +14,7 @@ module grid_io
   public :: read_grid_4d
   public :: write_grid_3d
   public :: write_grid_4d
+  public :: write_grid_5d
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
@@ -43,6 +44,14 @@ module grid_io
      module procedure write_grid_4d_int8
   end interface write_grid_4d
 
+  interface write_grid_5d
+     module procedure write_grid_5d_sp
+     module procedure write_grid_5d_dp
+     module procedure write_grid_5d_int
+     module procedure write_grid_5d_int8
+  end interface write_grid_5d
+
+
 contains
 
   logical function grid_exists(group, name)
@@ -108,6 +117,28 @@ subroutine read_grid_3d_int8(group, path, array, geo)
 
   end subroutine read_grid_3d_int8
 
+
+
+
+  subroutine write_grid_5d_int8(group, path, array, geo)
+
+    implicit none
+
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    integer(idp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+
+    call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/)))
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_int8
+
+
+
+
   subroutine write_grid_4d_int8(group, path, array, geo)
 
     implicit none
@@ -192,6 +223,24 @@ subroutine read_grid_3d_int(group, path, array, geo)
 
   end subroutine read_grid_3d_int
 
+
+
+
+  subroutine write_grid_5d_int(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    integer, intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/)))
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_int
+
+
   subroutine write_grid_4d_int(group, path, array, geo)
 
     implicit none
@@ -276,6 +325,24 @@ subroutine read_grid_3d_dp(group, path, array, geo)
 
   end subroutine read_grid_3d_dp
 
+
+
+  subroutine write_grid_5d_dp(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    real(dp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+
+    call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/)))
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_dp
+
+
   subroutine write_grid_4d_dp(group, path, array, geo)
 
     implicit none
@@ -360,6 +427,23 @@ subroutine read_grid_3d_sp(group, path, array, geo)
 
   end subroutine read_grid_3d_sp
 
+
+  !CRAZY DN ADDITIONS
+  subroutine write_grid_5d_sp(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    real(sp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/)))
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_sp
+
+
   subroutine write_grid_4d_sp(group, path, array, geo)
 
     implicit none
diff --git a/src/grid/grid_io_1d.f90 b/src/grid/grid_io_1d.f90
index 5d87aaa2..7388046d 100644
--- a/src/grid/grid_io_1d.f90
+++ b/src/grid/grid_io_1d.f90
@@ -14,6 +14,8 @@ module grid_io
   public :: read_grid_4d
   public :: write_grid_3d
   public :: write_grid_4d
+  public :: write_grid_5d
+  
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
@@ -43,6 +45,14 @@ module grid_io
      module procedure write_grid_4d_int8
   end interface write_grid_4d
 
+  interface write_grid_5d
+     module procedure write_grid_5d_sp
+     module procedure write_grid_5d_dp
+     module procedure write_grid_5d_int
+     module procedure write_grid_5d_int8
+  end interface write_grid_5d
+
+
 contains
 
   logical function grid_exists(group, name)
@@ -95,6 +105,23 @@ subroutine read_grid_3d_int8(group, path, array, geo)
 
   end subroutine read_grid_3d_int8
 
+
+
+  subroutine write_grid_5d_int8(group, path, array, geo)
+    
+    implicit none
+    
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    integer(idp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, array)
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_int8
+
+
   subroutine write_grid_4d_int8(group, path, array, geo)
 
     implicit none
@@ -166,6 +193,25 @@ subroutine read_grid_3d_int(group, path, array, geo)
 
   end subroutine read_grid_3d_int
 
+
+
+
+
+  subroutine write_grid_5d_int(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    integer, intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, array)
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_int
+
+
   subroutine write_grid_4d_int(group, path, array, geo)
 
     implicit none
@@ -194,7 +240,7 @@ subroutine write_grid_3d_int(group, path, array, geo)
 
   end subroutine write_grid_3d_int
 
-
+ 
   subroutine read_grid_4d_dp(group, path, array, geo)
 
     implicit none
@@ -237,6 +283,23 @@ subroutine read_grid_3d_dp(group, path, array, geo)
 
   end subroutine read_grid_3d_dp
 
+
+
+  subroutine write_grid_5d_dp(group, path, array, geo)
+    
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    real(dp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, array)
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_dp
+
+
   subroutine write_grid_4d_dp(group, path, array, geo)
 
     implicit none
@@ -265,7 +328,7 @@ subroutine write_grid_3d_dp(group, path, array, geo)
 
   end subroutine write_grid_3d_dp
 
-
+  
   subroutine read_grid_4d_sp(group, path, array, geo)
 
     implicit none
@@ -308,6 +371,23 @@ subroutine read_grid_3d_sp(group, path, array, geo)
 
   end subroutine read_grid_3d_sp
 
+
+
+  subroutine write_grid_5d_sp(group, path, array, geo)
+    
+    implicit none
+    
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    real(sp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+    
+    call mp_write_array(group, path, array)
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+    
+  end subroutine write_grid_5d_sp
+
+
   subroutine write_grid_4d_sp(group, path, array, geo)
 
     implicit none
diff --git a/src/grid/grid_io_1d_template.f90 b/src/grid/grid_io_1d_template.f90
index 1abbe776..f36aa0b7 100644
--- a/src/grid/grid_io_1d_template.f90
+++ b/src/grid/grid_io_1d_template.f90
@@ -41,6 +41,15 @@ module grid_io
      module procedure write_grid_4d_int
      module procedure write_grid_4d_int8
   end interface write_grid_4d
+  
+  interface write_grid_5d
+     module procedure write_grid_5d_sp
+     module procedure write_grid_5d_dp
+     module procedure write_grid_5d_int
+     module procedure write_grid_5d_int8
+  end interface write_grid_5d
+
+
 
 contains
 
@@ -53,6 +62,7 @@ end function grid_exists
 
   !!@FOR real(sp):sp real(dp):dp integer:int integer(idp):int8
 
+
   subroutine read_grid_4d_<T>(group, path, array, geo)
 
     implicit none
@@ -95,6 +105,22 @@ subroutine read_grid_3d_<T>(group, path, array, geo)
 
   end subroutine read_grid_3d_<T>
 
+
+  subroutine write_grid_5d_<T>(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    @T, intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, array)
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_<T>
+
+
   subroutine write_grid_4d_<T>(group, path, array, geo)
 
     implicit none
diff --git a/src/grid/grid_io_amr.f90 b/src/grid/grid_io_amr.f90
index 832e1725..b49281c4 100644
--- a/src/grid/grid_io_amr.f90
+++ b/src/grid/grid_io_amr.f90
@@ -14,6 +14,8 @@ module grid_io
   public :: read_grid_4d
   public :: write_grid_3d
   public :: write_grid_4d
+  public :: write_grid_5d
+
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
@@ -43,6 +45,14 @@ module grid_io
      module procedure write_grid_4d_int8
   end interface write_grid_4d
 
+  interface write_grid_5d
+     module procedure write_grid_5d_sp
+     module procedure write_grid_5d_dp
+     module procedure write_grid_5d_int
+     module procedure write_grid_5d_int8
+  end interface write_grid_5d
+  
+
 contains
 
   logical function grid_exists(group, name)
@@ -146,6 +156,52 @@ end subroutine test
 
   end subroutine read_grid_3d_int8
 
+
+
+    subroutine write_grid_5d_int8(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    integer(idp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in),target :: geo
+    integer(hid_t) :: g_level, g_grid
+    character(len=100) :: name
+    integer :: ilevel, igrid
+    type(level_desc), pointer :: level
+    type(grid_desc), pointer :: grid
+
+    do ilevel=1,size(geo%levels)
+       level => geo%levels(ilevel)
+       write(name, '("level_", I5.5)') ilevel
+       if(mp_path_exists(group, name)) then
+          g_level = mp_open_group(group, name)
+       else
+          g_level = mp_create_group(group, name)
+       end if
+       do igrid=1,size(level%grids)
+          grid => level%grids(igrid)
+          write(name, '("grid_", I5.5)') igrid
+          if(mp_path_exists(g_level, name)) then
+             g_grid = mp_open_group(g_level, name)
+          else
+             g_grid = mp_create_group(g_level, name)
+          end if
+
+          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :,:), &
+               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/)))
+          call mp_close_group(g_grid)
+       end do
+       call mp_close_group(g_level)
+
+    end do
+
+  end subroutine write_grid_5d_int8
+
+
+
+
   subroutine write_grid_4d_int8(group, path, array, geo)
 
     implicit none
@@ -226,6 +282,7 @@ subroutine write_grid_3d_int8(group, path, array, geo)
   end subroutine write_grid_3d_int8
 
 
+
   subroutine read_grid_4d_int(group, path, array, geo)
 
     implicit none
@@ -304,6 +361,47 @@ end subroutine test
 
   end subroutine read_grid_3d_int
 
+
+  subroutine write_grid_5d_int(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    integer, intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in),target :: geo
+    integer(hid_t) :: g_level, g_grid
+    character(len=100) :: name
+    integer :: ilevel, igrid
+    type(level_desc), pointer :: level
+    type(grid_desc), pointer :: grid
+
+    do ilevel=1,size(geo%levels)
+       level => geo%levels(ilevel)
+       write(name, '("level_", I5.5)') ilevel
+       if(mp_path_exists(group, name)) then
+          g_level = mp_open_group(group, name)
+       else
+          g_level = mp_create_group(group, name)
+       end if
+       do igrid=1,size(level%grids)
+          grid => level%grids(igrid)
+          write(name, '("grid_", I5.5)') igrid
+          if(mp_path_exists(g_level, name)) then
+             g_grid = mp_open_group(g_level, name)
+          else
+             g_grid = mp_create_group(g_level, name)
+          end if
+          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :,:), &
+               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/)))
+          call mp_close_group(g_grid)
+       end do
+       call mp_close_group(g_level)
+    end do
+    
+  end subroutine write_grid_5d_int
+
+
   subroutine write_grid_4d_int(group, path, array, geo)
 
     implicit none
@@ -462,6 +560,48 @@ end subroutine test
 
   end subroutine read_grid_3d_dp
 
+
+
+  subroutine write_grid_5d_dp(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    real(dp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in),target :: geo
+    integer(hid_t) :: g_level, g_grid
+    character(len=100) :: name
+    integer :: ilevel, igrid
+    type(level_desc), pointer :: level
+    type(grid_desc), pointer :: grid
+
+    do ilevel=1,size(geo%levels)
+       level => geo%levels(ilevel)
+       write(name, '("level_", I5.5)') ilevel
+       if(mp_path_exists(group, name)) then
+          g_level = mp_open_group(group, name)
+       else
+          g_level = mp_create_group(group, name)
+       end if
+       do igrid=1,size(level%grids)
+          grid => level%grids(igrid)
+          write(name, '("grid_", I5.5)') igrid
+          if(mp_path_exists(g_level, name)) then
+             g_grid = mp_open_group(g_level, name)
+          else
+             g_grid = mp_create_group(g_level, name)
+          end if
+          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :), &
+               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/)))
+          call mp_close_group(g_grid)
+       end do
+       call mp_close_group(g_level)
+    end do
+
+  end subroutine write_grid_5d_dp
+
+
   subroutine write_grid_4d_dp(group, path, array, geo)
 
     implicit none
@@ -620,6 +760,49 @@ end subroutine test
 
   end subroutine read_grid_3d_sp
 
+
+
+  subroutine write_grid_5d_sp(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    real(sp), intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in),target :: geo
+    integer(hid_t) :: g_level, g_grid
+    character(len=100) :: name
+    integer :: ilevel, igrid
+    type(level_desc), pointer :: level
+    type(grid_desc), pointer :: grid
+
+    do ilevel=1,size(geo%levels)
+       level => geo%levels(ilevel)
+       write(name, '("level_", I5.5)') ilevel
+       if(mp_path_exists(group, name)) then
+          g_level = mp_open_group(group, name)
+       else
+          g_level = mp_create_group(group, name)
+       end if
+       do igrid=1,size(level%grids)
+          grid => level%grids(igrid)
+          write(name, '("grid_", I5.5)') igrid
+          if(mp_path_exists(g_level, name)) then
+             g_grid = mp_open_group(g_level, name)
+          else
+             g_grid = mp_create_group(g_level, name)
+          end if
+          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :), &
+               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2),size(array,3)/)))
+          call mp_close_group(g_grid)
+       end do
+       call mp_close_group(g_level)
+    end do
+
+  end subroutine write_grid_5d_sp
+
+
+
   subroutine write_grid_4d_sp(group, path, array, geo)
 
     implicit none
diff --git a/src/grid/grid_io_amr_template.f90 b/src/grid/grid_io_amr_template.f90
index f80eb96c..5e06a8c1 100644
--- a/src/grid/grid_io_amr_template.f90
+++ b/src/grid/grid_io_amr_template.f90
@@ -13,6 +13,8 @@ module grid_io
   public :: read_grid_4d
   public :: write_grid_3d
   public :: write_grid_4d
+  public :: write_grid_5d
+
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
@@ -42,6 +44,15 @@ module grid_io
      module procedure write_grid_4d_int8
   end interface write_grid_4d
 
+  interface write_grid_5d
+     module procedure write_grid_5d_sp
+     module procedure write_grid_5d_dp
+     module procedure write_grid_5d_int
+     module procedure write_grid_5d_int8
+  end interface write_grid_5d
+  
+  
+
 contains
 
   logical function grid_exists(group, name)
@@ -146,6 +157,48 @@ end subroutine test
 
   end subroutine read_grid_3d_<T>
 
+
+
+  subroutine write_grid_5d_<T>(group, path, array, geo)
+
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    @T, intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in),target :: geo
+    integer(hid_t) :: g_level, g_grid
+    character(len=100) :: name
+    integer :: ilevel, igrid
+    type(level_desc), pointer :: level
+    type(grid_desc), pointer :: grid
+
+    do ilevel=1,size(geo%levels)
+       level => geo%levels(ilevel)
+       write(name, '("level_", I5.5)') ilevel
+       if(mp_path_exists(group, name)) then
+          g_level = mp_open_group(group, name)
+       else
+          g_level = mp_create_group(group, name)
+       end if
+       do igrid=1,size(level%grids)
+          grid => level%grids(igrid)
+          write(name, '("grid_", I5.5)') igrid
+          if(mp_path_exists(g_level, name)) then
+             g_grid = mp_open_group(g_level, name)
+          else
+             g_grid = mp_create_group(g_level, name)
+          end if
+          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), &
+               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/)))
+          call mp_close_group(g_grid)
+       end do
+       call mp_close_group(g_level)
+
+    end do
+
+  end subroutine write_grid_5d_<T>
+
   subroutine write_grid_4d_<T>(group, path, array, geo)
 
     implicit none
diff --git a/src/grid/grid_io_template.f90 b/src/grid/grid_io_template.f90
index f42e6eab..696575d0 100644
--- a/src/grid/grid_io_template.f90
+++ b/src/grid/grid_io_template.f90
@@ -13,6 +13,8 @@ module grid_io
   public :: read_grid_4d
   public :: write_grid_3d
   public :: write_grid_4d
+  public :: write_grid_5d
+
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
@@ -42,6 +44,14 @@ module grid_io
      module procedure write_grid_4d_int8
   end interface write_grid_4d
 
+  interface write_grid_5d
+     module procedure write_grid_5d_sp
+     module procedure write_grid_5d_dp
+     module procedure write_grid_5d_int
+     module procedure write_grid_5d_int8
+  end interface write_grid_5d
+  
+
 contains
 
   logical function grid_exists(group, name)
@@ -53,6 +63,7 @@ end function grid_exists
 
   !!@FOR real(sp):sp real(dp):dp integer:int integer(idp):int8
 
+
   subroutine read_grid_4d_<T>(group, path, array, geo)
 
     implicit none
@@ -107,6 +118,22 @@ subroutine read_grid_3d_<T>(group, path, array, geo)
     array = reshape(array3d, (/n_cells/))
 
   end subroutine read_grid_3d_<T>
+  
+  subroutine write_grid_5d_<T>(group, path, array, geo)
+    
+    implicit none
+
+    integer(hid_t), intent(in) :: group
+    character(len=*), intent(in) :: path
+    @T, intent(in) :: array(:,:,:)
+    type(grid_geometry_desc),intent(in) :: geo
+
+    call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/)))
+    call mp_write_keyword(group, path, 'geometry', geo%id)
+
+  end subroutine write_grid_5d_<T>
+
+
 
   subroutine write_grid_4d_<T>(group, path, array, geo)
 
diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90
index 90674116..a9a28fdf 100644
--- a/src/grid/grid_physics_3d.f90
+++ b/src/grid/grid_physics_3d.f90
@@ -36,8 +36,12 @@ module grid_physics
   integer(idp),allocatable, public :: n_photons(:)
   integer(idp),allocatable, public :: last_photon_id(:)
   real(dp),allocatable, public :: specific_energy(:,:)
+  real(dp),allocatable, public :: specific_energy_nu(:,:,:) 
   real(dp),allocatable, public :: specific_energy_sum(:,:)
+  real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:)
+
   real(dp),allocatable, public :: specific_energy_additional(:,:)
+  real(dp),allocatable, public :: specific_energy_additional_nu(:,:,:)
   real(dp),allocatable, public :: energy_abs_tot(:)
   real(dp),allocatable, public :: minimum_specific_energy(:)
 
@@ -51,12 +55,16 @@ module grid_physics
   ! Temporary variables (used for convenience in external modules)
   real(dp), allocatable, public :: tmp_column_density(:)
 
+  !Counting indices
   integer :: id
+  integer :: idx 
 
   logical :: debug = .false.
 
   type(pdf_discrete_dp) :: absorption
 
+
+
 contains
 
   real(dp) function tau_inv_planck_to_closest_wall(p) result(tau)
@@ -89,16 +97,19 @@ integer function select_dust_specific_energy_rho(icell) result(id_select)
     id_select = sample_pdf(absorption)
   end function select_dust_specific_energy_rho
 
-  subroutine setup_grid_physics(group, use_mrw, use_pda)
+  subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf)
 
     implicit none
 
     integer(hid_t),intent(in) :: group
-    logical,intent(in) :: use_mrw, use_pda
+    logical,intent(in) :: use_mrw, use_pda, compute_isrf
+    integer :: n_isrf_wavelengths
 
+    n_isrf_wavelengths = d(1)%n_nu
     ! Density
     allocate(density(geo%n_cells, n_dust))
     allocate(specific_energy(geo%n_cells, n_dust))
+    allocate(specific_energy_nu(geo%n_cells,n_dust,n_isrf_wavelengths))
 
     if(n_dust > 0) then
 
@@ -146,6 +157,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda)
           ! Check number of dust types for specific_energy
           if(size(specific_energy, 2).ne.n_dust) call error("setup_grid","specific_energy array has wrong number of dust types")
 
+
+
           ! Reset specific energy to zero in masked cells
           if(geo%masked) then
              if(main_process()) write(*, '(" [grid_physics] applying mask to specific_energy grid")')
@@ -153,19 +166,39 @@ subroutine setup_grid_physics(group, use_mrw, use_pda)
                 where(.not.geo%mask)
                    specific_energy(:, id) = 0.
                 end where
+                
+                if (compute_isrf) then
+                   do idx=1,n_isrf_wavelengths
+                      where(.not.geo%mask)
+                         specific_energy_nu(:, id, idx) = 0.
+                      endwhere
+                   end do
+                end if
+
              end do
           end if
 
           if(trim(specific_energy_type) == 'additional') then
              allocate(specific_energy_additional(geo%n_cells, n_dust))
+             if (compute_isrf) then 
+                allocate(specific_energy_additional_nu(geo%n_cells,n_dust,n_isrf_wavelengths))
+             end if
              ! We store a copy of the initial specific energy in a separate
              ! array, and we set the specific energy to the minimum specific
              ! energy. After the first iteration, specific_energy will get
              ! re-calculated and we will then add specific_energy_additional
              specific_energy_additional = specific_energy
+             specific_energy_additional_nu = specific_energy_nu
              do id=1,n_dust
                 specific_energy(:,id) = minimum_specific_energy(id)
-             end do
+                
+                if (compute_isrf) then
+                   do idx=1,n_isrf_wavelengths
+                      specific_energy_nu(:,id,idx) = minimum_specific_energy(id)
+                   end do
+                end if
+
+            end do
           end if
 
 
@@ -178,6 +211,13 @@ subroutine setup_grid_physics(group, use_mrw, use_pda)
           ! Set all specific_energy to minimum requested
           do id=1,n_dust
              specific_energy(:,id) = minimum_specific_energy(id)
+
+             if (compute_isrf) then
+                do idx = 1,n_isrf_wavelengths
+                   specific_energy_nu(:,id,idx) = minimum_specific_energy(id)
+                end do
+             end if
+
           end do
 
        end if
@@ -191,6 +231,11 @@ subroutine setup_grid_physics(group, use_mrw, use_pda)
     allocate(specific_energy_sum(geo%n_cells, n_dust))
     specific_energy_sum = 0._dp
 
+    ! Set up basics for ISRF calculation
+    !n_isrf_wavelengths = d(1)%n_nu
+    allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths))
+    specific_energy_sum_nu = 0._dp
+    
     ! Total energy absorbed
     allocate(energy_abs_tot(n_dust))
     energy_abs_tot = 0._dp
@@ -257,6 +302,9 @@ subroutine sublimate_dust()
     implicit none
     integer :: ic, id
     integer :: reset
+    integer :: n_isrf_wavelengths
+
+    n_isrf_wavelengths = d(1)%n_nu
 
     reset = 0
 
@@ -270,7 +318,15 @@ subroutine sublimate_dust()
                 density(ic, id) = 0.
                 specific_energy(ic, id) = minimum_specific_energy(id)
                 reset = reset + 1
+
+                if (compute_isrf) then
+                   do idx=1,n_isrf_wavelengths
+                      specific_energy_nu(ic,id,idx) = minimum_specific_energy(id)
+                   end do
+                end if
+
              end if
+             
           end do
           if(reset > 0) write(*,'(" [sublimate_dust] dust removed in ",I0," cells")') reset
 
@@ -284,6 +340,14 @@ subroutine sublimate_dust()
                      & / chi_rosseland(id, d(id)%sublimation_specific_energy))**2
                 specific_energy(ic,id) = d(id)%sublimation_specific_energy
                 reset = reset + 1
+
+
+                if (compute_isrf) then
+                   do idx=1,n_isrf_wavelengths
+                      specific_energy_nu(ic,id,idx) = minimum_specific_energy(id)
+                   end do
+                end if
+                   
              end if
           end do
 
@@ -296,6 +360,14 @@ subroutine sublimate_dust()
                 specific_energy(ic, id) = d(id)%sublimation_specific_energy
                 reset = reset + 1
              end if
+
+             
+             if (compute_isrf) then
+                do idx=1,n_isrf_wavelengths
+                   specific_energy_nu(ic,id,idx) = minimum_specific_energy(id)
+                end do
+             end if
+
           end do
 
           if(reset > 0) write(*,'(" [sublimate_dust] capping dust specific_energy in ",I0," cells")') reset
@@ -316,12 +388,23 @@ subroutine update_energy_abs(scale)
 
     real(dp), intent(in) :: scale
 
-    integer :: id
+    integer :: id,idx
+    integer :: n_isrf_wavelengths
+    
+    n_isrf_wavelengths = d(1)%n_nu
 
     if(main_process()) write(*,'(" [grid_physics] updating energy_abs")')
 
     do id=1,n_dust
        specific_energy(:,id) = specific_energy_sum(:,id) * scale / geo%volume
+       
+       if (compute_isrf) then
+          do idx=1,n_isrf_wavelengths
+             specific_energy_nu(:,id,idx) = specific_energy_sum_nu(:,id,idx) * scale/geo%volume
+          end do
+       end if
+          
+          
        where(geo%volume == 0._dp)
           specific_energy(:,id) = 0._dp
        end where
@@ -331,10 +414,17 @@ subroutine update_energy_abs(scale)
        write(*,'(" [update_energy_abs] ",I0," cells have no energy")') count(specific_energy==0.and.density>0.)
     end if
 
+
     ! Add in additional source of heating
     if(trim(specific_energy_type) == 'additional') then
        if(main_process()) write(*,'(" [grid_physics] adding additional heating source")')
        specific_energy = specific_energy + specific_energy_additional
+
+       
+       if (compute_isrf) then 
+          specific_energy_nu = specific_energy_nu + specific_energy_additional_nu
+       end if
+
     end if
 
     call update_energy_abs_tot()
@@ -357,7 +447,10 @@ subroutine check_energy_abs()
           if(main_process()) call warn("update_energy_abs","specific_energy below minimum requested in some cells - resetting")
           where(specific_energy(:,id) < minimum_specific_energy(id))
              specific_energy(:,id) = minimum_specific_energy(id)
+             
+ 
           end where
+
        end if
 
        if(any(specific_energy(:,id) < d(id)%specific_energy(1))) then
@@ -365,7 +458,8 @@ subroutine check_energy_abs()
              if(main_process()) call warn("update_energy_abs","specific_energy below minimum allowed in some cells - resetting")
              where(specific_energy(:,id) < d(id)%specific_energy(1))
                 specific_energy(:,id) = d(id)%specific_energy(1)
-             end where
+                
+              end where
           else
              if(main_process()) call warn("update_energy_abs","specific_energy below minimum allowed in some cells - will pick closest emissivities")
           end if
@@ -376,6 +470,7 @@ subroutine check_energy_abs()
              if(main_process()) call warn("update_energy_abs","specific_energy above maximum allowed in some cells - resetting")
              where(specific_energy(:,id) > d(id)%specific_energy(d(id)%n_e))
                 specific_energy(:,id) = d(id)%specific_energy(d(id)%n_e)
+
              end where
           else
              if(main_process()) call warn("update_energy_abs","specific_energy above maximum allowed in some cells - will pick closest emissivities")
diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90
index 3c345161..1f10ddcb 100644
--- a/src/grid/grid_propagate_3d.f90
+++ b/src/grid/grid_propagate_3d.f90
@@ -3,12 +3,14 @@ module grid_propagate
   use core_lib
   use type_photon, only : photon
   use type_grid_cell
-  use dust_main, only : n_dust
+  use dust_main, only : n_dust,d
   use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall
-  use grid_physics, only : specific_energy_sum, density, n_photons, last_photon_id
+  use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, density, n_photons, last_photon_id
   use sources
   use counters
-  use settings, only : frac_check => propagation_check_frequency
+  use settings, only : frac_check => propagation_check_frequency, compute_isrf
+
+  use grid_geometry, only : geo
 
   implicit none
   save
@@ -57,6 +59,15 @@ subroutine grid_integrate(p,tau_required,tau_achieved)
 
     integer :: source_id
 
+    integer :: idx
+    real, dimension(d(1)%n_nu) :: energy_frequency_bins
+
+    
+    do id=1,d(1)%n_nu
+       energy_frequency_bins(id) = d(1)%nu(id)
+    end do
+   
+
     radial = (p%r .dot. p%v) > 0.
 
     if(debug) write(*,'(" [debug] start grid_integrate")')
@@ -130,14 +141,29 @@ subroutine grid_integrate(p,tau_required,tau_achieved)
           p%r = p%r + tmin * p%v
           tau_achieved = tau_achieved + tau_cell
 
-          ! NEED INDIVIDUAL ALPHA HERE
 
-          do id=1,n_dust
-             if(density(p%icell%ic, id) > 0._dp) then
-                specific_energy_sum(p%icell%ic, id) = &
-                     & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy
-             end if
-          end do
+          ! Compute the ISRF 
+          if (compute_isrf) then 
+
+             idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1)
+
+             
+             do id=1,n_dust
+                
+                if(density(p%icell%ic, id) > 0._dp) then
+                   specific_energy_sum(p%icell%ic, id) = &
+                        & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy
+                end if
+                
+                
+                if(density(p%icell%ic,id) > 0._dp) then
+                   specific_energy_sum_nu(p%icell%ic,id,idx) = &
+                        & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy
+                end if
+             end do
+          endif
+
+
 
           p%on_wall = .true.
           p%icell = next_cell(p%icell, id_min, intersection=p%r)
@@ -221,6 +247,15 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved)
 
     integer :: source_id
 
+    integer :: id
+    integer :: idx
+    real, dimension(d(1)%n_nu) :: energy_frequency_bins
+
+    do id=1,d(1)%n_nu
+       energy_frequency_bins(id) = d(1)%nu(id)
+       !print *,'[grid_propagate_3d last iteration] energy_frequency_bins(id) = ',energy_frequency_bins(id)
+    end do
+
     radial = (p%r .dot. p%v) > 0.
 
     if(debug) write(*,'(" [debug] start grid_integrate_noenergy")')
@@ -268,9 +303,25 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved)
        end if
 
        chi_rho_total = 0._dp
+
+
+       !Compute the ISRF
+
+       !Figure out what frequency bin in the ISRF calculation the current photon's frequency is closest to
+       idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1)
+
        do id=1,n_dust
           chi_rho_total = chi_rho_total + p%current_chi(id) * density(p%icell%ic, id)
+
+          if(density(p%icell%ic,id) > 0._dp) then
+             specific_energy_sum_nu(p%icell%ic,id,idx) = &
+                  & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy
+             !print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu(p%icell%ic,id,idx)
+          end if
+
        end do
+
+
        tau_cell = chi_rho_total * tmin
 
        if(tau_cell < tau_needed) then
diff --git a/src/main/settings.f90 b/src/main/settings.f90
index daf4ad23..d8038eef 100644
--- a/src/main/settings.f90
+++ b/src/main/settings.f90
@@ -20,7 +20,7 @@ module settings
   integer(idp),public :: n_last_photons_dust = 0
   integer(idp),public :: n_raytracing_photons_sources = 0
   integer(idp),public :: n_raytracing_photons_dust = 0
-  logical,public :: use_raytracing, use_mrw, use_pda
+  logical,public :: use_raytracing, use_mrw, use_pda, compute_isrf
   logical, public :: kill_on_absorb, kill_on_scatter
   real(dp),public :: mrw_gamma
   logical, public :: forced_first_interaction
diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90
index 649799fa..5be59417 100644
--- a/src/main/setup_rt.f90
+++ b/src/main/setup_rt.f90
@@ -74,6 +74,7 @@ subroutine setup_initial(input_handle)
 
     call mp_read_keyword(input_handle, '/', 'pda', use_pda)
     call mp_read_keyword(input_handle, '/', 'mrw', use_mrw)
+    call mp_read_keyword(input_handle, '/', 'isrf', compute_isrf)
 
     if(use_mrw) then
        call mp_read_keyword(input_handle, '/', 'mrw_gamma', mrw_gamma)
@@ -173,7 +174,7 @@ subroutine setup_initial(input_handle)
     call mp_close_group(g_geometry)
 
     g_physics = mp_open_group(input_handle, '/Grid/Quantities')
-    call setup_grid_physics(g_physics, use_mrw, use_pda)
+    call setup_grid_physics(g_physics, use_mrw, use_pda, compute_isrf)
     call mp_close_group(g_physics)
 
     call mp_read_keyword(input_handle, '/', 'physics_io_bytes', physics_io_bytes)
diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90
index fb0d6e7b..32e1c777 100644
--- a/src/mpi/mpi_routines.f90
+++ b/src/mpi/mpi_routines.f90
@@ -266,6 +266,8 @@ subroutine mp_collect_physical_arrays()
     real(dp) :: tmp
     real(dp) :: dummy_dp
     integer(idp) :: dummy_idp
+
+    real(dp),allocatable :: tmp_3d(:,:,:)
     real(dp),allocatable :: tmp_2d(:,:)
     integer(idp),allocatable :: tmp_int_1d(:)
 
@@ -278,6 +280,17 @@ subroutine mp_collect_physical_arrays()
        call mpi_reduce(specific_energy_sum, dummy_dp, size(specific_energy_sum), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr)
     end if
 
+    if(main_process()) then
+       allocate(tmp_3d(size(specific_energy_sum_nu,1),size(specific_energy_sum_nu,2),size(specific_energy_sum_nu,3)))
+       call mpi_reduce(specific_energy_sum_nu, tmp_3d, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr)
+       specific_energy_sum_nu = tmp_3d
+       deallocate(tmp_3d)
+    else
+       call mpi_reduce(specific_energy_sum_nu, dummy_dp, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr)
+    end if
+
+
+
     if(allocated(n_photons)) then
        if(main_process()) then
           allocate(tmp_int_1d(size(n_photons,1)))