diff --git a/build/Makefile_setups b/build/Makefile_setups
index 613f126e7..edd7170cd 100644
--- a/build/Makefile_setups
+++ b/build/Makefile_setups
@@ -177,6 +177,18 @@ ifeq ($(SETUP), grtde)
     ANALYSIS= analysis_gws.f90
 endif
 
+ifeq ($(SETUP), radiotde)
+#   radio tidal disruption event in general relativity
+    GR=yes
+    METRIC=minkowski
+    KNOWN_SETUP=yes
+    GRAVITY=no
+    IND_TIMESTEPS=no
+    ANALYSIS=analysis_radiotde.f90
+    MODFILE=moddump_radiotde.f90
+    SYSTEM=gfortran
+endif
+
 ifeq ($(SETUP), srpolytrope)
 #   polytrope in special relativity
     FPPFLAGS= -DPRIM2CONS_FIRST
diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90
index 5e4ca24b0..c54602fd5 100644
--- a/src/main/step_leapfrog.F90
+++ b/src/main/step_leapfrog.F90
@@ -823,7 +823,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
     !
     ! predictor step for external forces, also recompute external forces
     !
-    !$omp parallel do default(none) &
+    !$omp parallel do default(none) schedule(runtime)  &
     !$omp shared(npart,xyzh,vxyzu,fext,iphase,ntypes,massoftype) &
     !$omp shared(maxphase,maxp,eos_vars) &
     !$omp shared(dt,hdt,xtol,ptol) &
@@ -957,7 +957,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
     !$omp reduction(min:dtextforce_min) &
     !$omp reduction(+:accretedmass,naccreted,nlive) &
     !$omp shared(idamp,damp_fac)
-    !$omp do
+    !$omp do schedule(runtime)
     accreteloop: do i=1,npart
        if (.not.isdead_or_accreted(xyzh(4,i))) then
           if (ntypes > 1 .and. maxphase==maxp) then
diff --git a/src/utils/analysis_radiotde.f90 b/src/utils/analysis_radiotde.f90
new file mode 100644
index 000000000..4d009de92
--- /dev/null
+++ b/src/utils/analysis_radiotde.f90
@@ -0,0 +1,289 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2023 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.bitbucket.io/                                          !
+!--------------------------------------------------------------------------!
+module analysis
+!
+! Computes the outflow profile in a TDE simulation
+!
+! :References: None
+!
+! :Owner: Fitz Hu
+!
+! :Runtime parameters:
+!   - rad_cap   : *capture shell radius*
+!   - drad_cap  : *capture shell thickness*
+!   - v_max     : *max velocity*
+!   - v_min     : *min velocity*
+!   - theta_max : *max azimuthal angle*
+!   - theta_min : *min azimuthal angle*
+!   - phi_max   : *max altitude angle*
+!   - phi_min   : *min altitude angle*
+!
+! :Dependencies: infile_utils, io, physcon, readwrite_dumps, units
+!
+ implicit none
+ character(len=8), parameter, public :: analysistype = 'radiotde'
+ public :: do_analysis
+
+ private
+
+ real, dimension(:), allocatable    :: theta,plot_theta,phi,vr,vtheta,vphi
+ logical, dimension(:), allocatable :: cap
+
+ !---- These can be changed in the params file
+ real    :: rad_cap = 1.e16 ! radius where the outflow in captured (in cm)
+ real    :: drad_cap = 4.7267e14 ! thickness of the shell to capture outflow (in cm)
+ real    :: v_min = 0.
+ real    :: v_max = 1.
+ real    :: theta_min = -180.
+ real    :: theta_max = 180.
+ real    :: phi_min = -90.
+ real    :: phi_max = 90.
+ real    :: m_accum, m_cap, vr_accum_mean, vr_cap_mean, v_accum_mean, v_cap_mean, e_accum, e_cap
+ integer :: n_accum, n_cap
+
+contains
+
+subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit)
+ use readwrite_dumps, only: opened_full_dump
+ use units,           only: utime,udist,unit_energ,umass
+ use physcon,         only: solarm,days
+ character(len=*),   intent(in) :: dumpfile
+ integer,            intent(in) :: numfile,npart,iunit
+ real,               intent(in) :: xyzh(:,:),vxyzu(:,:)
+ real,               intent(in) :: pmass,time
+ character(len=120) :: output
+ character(len=30)  :: filename
+ integer            :: i,ierr
+ logical            :: iexist
+ real               :: toMsun,todays
+
+ m_accum = 0.
+ n_accum = 0
+ m_cap = 0.
+ n_cap = 0
+ e_accum = 0.
+ e_cap = 0.
+
+ toMsun = umass/solarm
+ todays = utime/days
+
+ if (.not.opened_full_dump) then
+    write(*,'("SKIPPING FILE -- (Not a full dump)")')
+    return
+ endif
+
+! Print the analysis being done
+ write(*,'("Performing analysis type ",A)') analysistype
+ write(*,'("Input file name is ",A)') dumpfile
+
+ write(output,"(a8,i5.5)") 'outflow_',numfile
+ write(*,'("Output file name is ",A)') output
+
+ ! Read black hole mass from params file
+ filename = 'analysis_'//trim(analysistype)//'.params'
+ inquire(file=filename,exist=iexist)
+ if (iexist) call read_tdeparams(filename,ierr)
+ if (.not.iexist.or.ierr/=0) then
+    call write_tdeparams(filename)
+    print*,' Edit '//trim(filename)//' and rerun phantomanalysis'
+    stop
+ endif
+
+ rad_cap = rad_cap/udist
+ if (drad_cap < 0.) then
+    drad_cap = huge(0.)
+ else
+    drad_cap = drad_cap/udist
+ endif
+ print*, 'Capture particles from', rad_cap, 'to', rad_cap+drad_cap
+
+ allocate(theta(npart),plot_theta(npart),phi(npart),vr(npart),vtheta(npart),vphi(npart),cap(npart))
+ cap = .false.
+
+ call tde_analysis(npart,pmass,xyzh,vxyzu)
+
+ if (n_cap > 0) then
+ open(iunit,file=output)
+ write(iunit,'("# ",es20.12,"   # TIME")') time
+ write(iunit,"('#',6(1x,'[',i2.2,1x,a11,']',2x))") &
+       1,'theta',      &
+       2,'thetap',  &
+       3,'phi', &
+       4,'vr',     &
+       5,'vtheta',   &
+       6,'vphi'
+
+ do i = 1,npart
+    if (cap(i)) then
+       write(iunit,'(6(es18.10,1X))') &
+          theta(i),   &
+          plot_theta(i), &
+          phi(i),    &
+          vr(i),   &
+          vtheta(i),    &
+          vphi(i)
+    endif
+ enddo
+ close(iunit)
+ endif
+
+ deallocate(theta,plot_theta,phi,vr,vtheta,vphi,cap)
+
+ inquire(file='outflows',exist=iexist)
+ if (iexist) then
+    open(iunit,file='outflows',status='old',position='append')
+ else
+    open(iunit,file='outflows',status='new')
+    write(iunit,'(9(A15,1X))') '# time', 'm_cap[msun]', 'm_accum[msun]', 'vr_accum_mean[c]', 'vr_cap_mean[c]', &
+                               'v_accum_mean[c]', 'v_cap_mean[c]', 'e_accum[erg]', 'e_cap[erg]'
+ endif
+ write(iunit,'(9(es18.10,1x))') &
+    time*todays, &
+    m_cap*toMsun, &
+    m_accum*toMsun, &
+    vr_accum_mean, &
+    vr_cap_mean, &
+    v_accum_mean, &
+    v_cap_mean, &
+    e_accum*unit_energ, &
+    e_cap*unit_energ
+ close(iunit) 
+
+ write(*,'(I8,1X,A2,1X,I8,1X,A34)') n_cap, 'of', npart, 'particles are in the capture shell'
+ write(*,'(I8,1X,A2,1X,I8,1X,A40)') n_accum, 'of', npart, 'particles are outside the capture radius'
+
+end subroutine do_analysis
+
+!--------------------------------------------------------------------------------------------------------------------
+!
+!-- Actual subroutine where the analysis is done!
+!
+!--------------------------------------------------------------------------------------------------------------------
+subroutine tde_analysis(npart,pmass,xyzh,vxyzu)
+ integer, intent(in) :: npart
+ real, intent(in)    :: pmass,xyzh(:,:),vxyzu(:,:)
+ integer :: i
+ real    :: r,v,x,y,z,xyz(1:3),vx,vy,vz,vxyz(1:3)
+ real    :: thetai,phii,vri
+ real    :: vr_accum_add,vr_cap_add,v_accum_add,v_cap_add
+
+ vr_accum_add = 0.
+ vr_cap_add = 0.
+ v_accum_add = 0.
+ v_cap_add = 0.
+ 
+ do i = 1,npart
+    x = xyzh(1,i)
+    y = xyzh(2,i)
+    z = xyzh(3,i)
+    xyz = (/x,y,z/)
+    vx = vxyzu(1,i)
+    vy = vxyzu(2,i)
+    vz = vxyzu(3,i)
+    vxyz = (/vx,vy,vz/)
+    r = sqrt(dot_product(xyz,xyz))
+    v = sqrt(dot_product(vxyz,vxyz))
+    if (r > rad_cap) then
+       m_accum = m_accum + pmass 
+       n_accum = n_accum + 1
+       e_accum = e_accum + 0.5*pmass*v**2
+       vri = dot_product(vxyz,xyz)/r
+       vr_accum_add = vr_accum_add + vri
+       v_accum_add = v_accum_add + v
+       if (r-rad_cap < drad_cap .and. (v .ge. v_min .and. v .le. v_max)) then
+          thetai = atan2d(y,x)
+          phii = atan2d(z,sqrt(x**2+y**2))
+          if ((thetai .ge. theta_min .and. thetai .le. theta_max) .and. (phii .ge. phi_min .and. phii .le. phi_max)) then
+             m_cap = m_cap + pmass
+             n_cap = n_cap + 1
+             cap(i) = .true.
+             theta(i) = thetai
+             phi(i) = phii
+             plot_theta(i) = theta(i) * sqrt(cosd(phi(i)))
+             vr(i) = vri
+             vtheta(i) = -sind(phii)*vx + cosd(phii)*vy
+             vphi(i) = cosd(thetai)*cosd(phii)*vx + cosd(thetai)*sind(phii)*vy - sind(thetai)*vz
+             e_cap = e_cap + 0.5*pmass*v**2
+             vr_cap_add = vr_cap_add + vri
+             v_cap_add = v_cap_add + v
+          endif
+       endif
+    endif
+ enddo
+ vr_accum_mean = vr_accum_add/n_accum
+ v_accum_mean = v_accum_add/n_accum
+ vr_cap_mean = vr_cap_add/n_cap
+ v_cap_mean = v_cap_add/n_cap
+
+end subroutine tde_analysis
+
+!----------------------------------------------------------------
+!+
+!  Read/write tde information from/to params file
+!+
+!----------------------------------------------------------------
+subroutine write_tdeparams(filename)
+ use infile_utils, only:write_inopt
+ character(len=*), intent(in) :: filename
+ integer, parameter :: iunit = 20
+
+ print "(a)",' writing analysis options file '//trim(filename)
+ open(unit=iunit,file=filename,status='replace',form='formatted')
+ write(iunit,"(a,/)") '# options when performing radio TDE analysis'
+
+ call write_inopt(rad_cap,'rad_cap','capture inner radius (in cm)',iunit)
+ call write_inopt(drad_cap,'drad_cap','capture thickness (in cm) (-ve for all particles at outer radius)',iunit)
+
+ call write_inopt(v_min,'v_min','min velocity (in c)',iunit)
+ call write_inopt(v_max,'v_max','max velocity (in c)',iunit)
+
+ call write_inopt(theta_min,'theta_min','min theta (in deg)',iunit)
+ call write_inopt(theta_max,'theta_max','max theta (in deg)',iunit)
+
+ call write_inopt(phi_min,'phi_min','min phi (in deg)',iunit)
+ call write_inopt(phi_max,'phi_max','max phi (in deg)',iunit)
+
+ close(iunit)
+
+end subroutine write_tdeparams
+
+subroutine read_tdeparams(filename,ierr)
+ use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db
+ use io,           only:error
+ character(len=*), intent(in)  :: filename
+ integer,          intent(out) :: ierr
+ integer, parameter        :: iunit = 21
+ integer                   :: nerr
+ type(inopts), allocatable :: db(:)
+
+ print "(a)",'reading analysis options from '//trim(filename)
+ nerr = 0
+ ierr = 0
+ call open_db_from_file(db,filename,iunit,ierr)
+ 
+ call read_inopt(rad_cap,'rad_cap',db,min=0.,errcount=nerr)
+ call read_inopt(drad_cap,'drad_cap',db,errcount=nerr)
+
+ call read_inopt(v_min,'v_min',db,min=0.,max=1.,errcount=nerr)
+ call read_inopt(v_max,'v_max',db,min=0.,max=1.,errcount=nerr)
+
+ call read_inopt(theta_min,'theta_min',db,min=-180.,max=180.,errcount=nerr)
+ call read_inopt(theta_max,'theta_max',db,min=-180.,max=180.,errcount=nerr)
+
+ call read_inopt(phi_min,'phi_min',db,min=-90.,max=90.,errcount=nerr)
+ call read_inopt(phi_max,'phi_max',db,min=-90.,max=90.,errcount=nerr)
+
+ call close_db(db)
+ if (nerr > 0) then
+    print "(1x,i2,a)",nerr,' error(s) during read of params file: re-writing...'
+    ierr = nerr
+ endif
+
+end subroutine read_tdeparams
+
+end module analysis
+
diff --git a/src/utils/moddump_radiotde.f90 b/src/utils/moddump_radiotde.f90
new file mode 100644
index 000000000..d854923a5
--- /dev/null
+++ b/src/utils/moddump_radiotde.f90
@@ -0,0 +1,412 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2023 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.bitbucket.io/                                          !
+!--------------------------------------------------------------------------!
+module moddump
+!
+! Setup a circumnuclear gas cloud around outflowing TDE
+!
+! :References: None
+!
+! :Owner: Fitz Hu
+!
+! :Runtime parameters:
+!   - temperature : *Temperature*
+!   - mu          : *mean molecular mass*
+!   - ieos_in     : *equation of state*
+!   - use_func    : *use broken power law or profile date points*
+!
+! :Dependencies: datafiles, eos, io, stretchmap, kernel,
+!   mpidomain, part, physcon, setup_params,
+!   spherical, timestep, units, infile_utils
+!
+ implicit none
+ public :: modify_dump
+ private :: rho,rho_tab,get_temp_r,uerg,calc_rhobreak,write_setupfile,read_setupfile
+
+ private
+ integer           :: ieos_in,nprof,nbreak,nbreak_old
+ real              :: temperature,mu,ignore_radius,rad_max,rad_min
+ character(len=50) :: profile_filename
+ character(len=3)  :: interpolation
+ real, allocatable :: rhof_n(:),rhof_rbreak(:),rhof_rhobreak(:)
+ real, allocatable :: rhof_n_in(:),rhof_rbreak_in(:)
+ real, allocatable :: rad_prof(:),dens_prof(:)
+ real              :: rhof_rho0
+ logical           :: use_func,use_func_old,remove_overlap
+
+contains
+
+!----------------------------------------------------------------
+!
+!  Sets up a circumnuclear gas cloud
+!
+!----------------------------------------------------------------
+subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
+ use physcon,      only:solarm,years,mass_proton_cgs
+ use setup_params, only:npart_total
+ use part,         only:igas,set_particle_type,delete_particles_inside_radius,delete_particles_outside_sphere
+ use io,           only:fatal,master,id
+ use units,        only:umass,udist,utime,set_units,unit_density
+ use timestep,     only:dtmax,tmax
+ use eos,          only:ieos,gmw
+ use kernel,       only:hfact_default
+ use stretchmap,   only:get_mass_r
+ use spherical,    only:set_sphere
+ use mpidomain,    only:i_belong
+ integer,           intent(inout)   :: npart
+ integer,           intent(inout)   :: npartoftype(:)
+ real,              intent(inout)   :: xyzh(:,:)
+ real,              intent(inout)   :: vxyzu(:,:)
+ real,              intent(inout)   :: massoftype(:)
+ integer                       :: i,ierr,iunit=12,iprof
+ integer                       :: np_sphere,npart_old
+ real                          :: totmass,delta,r
+ character(len=120)            :: fileset,fileprefix='radio'
+ logical                       :: read_temp,setexists
+ real, allocatable             :: masstab(:),temp_prof(:)
+ character(len=15), parameter  :: default_name = 'default_profile'
+ real, dimension(7), parameter :: dens_prof_default = (/8.9e-21, 5.1e-21, 3.3e-21, 2.6e-21, &
+                                                        6.6e-25, 3.4e-25, 8.1e-26/), &
+                                  rad_prof_default = (/8.7e16, 1.2e17, 1.4e17, 2.0e17, &
+                                                       4.0e17, 4.8e17, 7.1e17/) ! profile from Cendes+2021
+ procedure(rho), pointer       :: rhof
+
+ !--Check for existence of the .params files
+ fileset=trim(fileprefix)//'.params'
+ inquire(file=fileset,exist=setexists)
+
+ !--Set default values
+ temperature       = 10.           ! Temperature in Kelvin
+ mu                = 2.            ! mean molecular weight
+ ieos_in           = 2
+ ignore_radius     = 1.e14          ! in cm 
+ use_func          = .true.
+ use_func_old      = use_func
+ remove_overlap    = .true.
+  !--Power law default setups
+ rad_max           = 7.1e16        ! in cm
+ rad_min           = 8.7e15        ! in cm
+ nbreak            = 1
+ nbreak_old        = nbreak
+ rhof_rho0         = 1.e4*mu*mass_proton_cgs
+ allocate(rhof_n(nbreak),rhof_rbreak(nbreak))
+ rhof_n            = -1.7
+ rhof_rbreak       = rad_min
+  !--Profile default setups
+ read_temp         = .false.
+ profile_filename  = default_name
+ nprof             = 7
+ interpolation     = 'log'
+
+ !--Read values from .setup
+ if (setexists) call read_setupfile(fileset,ierr)
+ if (.not. setexists .or. ierr /= 0) then
+    !--Prompt to get inputs and write to file
+    call write_setupfile(fileset)
+    stop
+ elseif (nbreak /= nbreak_old) then
+    !--Rewrite setup file
+    write(*,'(a)') ' [nbreak] changed. Rewriting setup file ...'
+    deallocate(rhof_n,rhof_rbreak)
+    allocate(rhof_n(nbreak),rhof_rbreak(nbreak))
+    rhof_n = -1.7
+    rhof_rbreak = rad_min
+    call write_setupfile(fileset)
+    stop
+ elseif (use_func .neqv. use_func_old) then
+    !--Rewrite setup fi.e
+    write(*,'(a)') ' [use_func] changed. Rewriting setup file ...'
+    call write_setupfile(fileset)
+    stop
+ endif
+
+ !--allocate memory
+ if (use_func) then
+    rhof => rho
+    deallocate(rhof_n,rhof_rbreak)
+    allocate(rhof_n(nbreak),rhof_rbreak(nbreak),rhof_rhobreak(nbreak))
+    rhof_n(:) = rhof_n_in(1:nbreak)
+    rhof_rbreak(:) = rhof_rbreak_in(1:nbreak)
+    call calc_rhobreak()
+ else
+    if (temperature .le. 0) read_temp = .true.
+    rhof => rho_tab    
+
+    deallocate(rhof_n,rhof_rbreak)
+    allocate(dens_prof(nprof),rad_prof(nprof),masstab(nprof))
+    if (read_temp) allocate(temp_prof(nprof))
+
+    !--Read profile from data
+    if (profile_filename == default_name) then
+       rad_prof = rad_prof_default
+       dens_prof = dens_prof_default
+    else
+       open(iunit,file=profile_filename)
+       if (.not. read_temp) then
+          do iprof = 1,nprof
+             read(iunit,*) rad_prof(iprof), dens_prof(iprof)
+          enddo
+       else
+          do iprof = 1,nprof
+             read(iunit,*) rad_prof(iprof), dens_prof(iprof), temp_prof(iprof)
+          enddo
+       endif
+    endif
+ endif
+ ieos = ieos_in
+ gmw = mu
+
+ !--Everything to code unit
+ ignore_radius = ignore_radius/udist
+ if (use_func) then
+    rad_min = rad_min/udist
+    rad_max = rad_max/udist
+    rhof_rbreak = rhof_rbreak/udist
+    rhof_rhobreak = rhof_rhobreak/unit_density
+ else
+    rad_prof = rad_prof/udist
+    dens_prof = dens_prof/unit_density
+    rad_min = rad_prof(1)
+    rad_max = rad_prof(nprof)
+ endif
+
+ !--remove unwanted particles
+ npart_old = npart
+ call delete_particles_inside_radius((/0.,0.,0./),ignore_radius,npart,npartoftype)
+ write(*,'(I10,1X,A23,1X,E8.2,1X,A14)') npart_old - npart, 'particles inside radius', ignore_radius*udist, 'cm are deleted'
+ npart_old = npart
+ if (remove_overlap) then
+    call delete_particles_outside_sphere((/0.,0.,0./),rad_min,npart)
+    write(*,'(I10,1X,A24,1X,E8.2,1X,A14)') npart_old - npart, 'particles outside radius', rad_min*udist, 'cm are deleted'
+    npart_old = npart
+ endif
+
+ !--setup cloud
+ totmass = get_mass_r(rhof,rad_max,rad_min)
+ write(*,'(A42,1X,F5.2,1X,A10)') ' Total mass of the circumnuclear gas cloud:', totmass*umass/solarm, 'solar mass'
+ np_sphere = nint(totmass/massoftype(igas))
+ call set_sphere('random',id,master,rad_min,rad_max,delta,hfact_default,npart,xyzh, &
+                 rhofunc=rhof,nptot=npart_total,exactN=.true.,np_requested=np_sphere,mask=i_belong)
+ if (ierr /= 0) call fatal('moddump','error setting up the circumnuclear gas cloud')
+ 
+ npartoftype(igas) = npart
+ !--Set particle properties
+ do i = npart_old+1,npart
+    call set_particle_type(i,igas)
+    r = dot_product(xyzh(1:3,i),xyzh(1:3,i))
+    if (read_temp) temperature = get_temp_r(r,rad_prof,temp_prof)       
+    vxyzu(4,i) = uerg(rhof(r),temperature)
+    vxyzu(1:3,i) = 0. ! stationary for now
+ enddo
+ 
+ !--Set timesteps
+ tmax = 10.*years/utime
+ dtmax = tmax/1000.
+
+end subroutine modify_dump
+
+!--Functions
+
+real function rho(r)
+ real, intent(in) :: r
+ integer          :: i
+ logical          :: found_rad
+
+ found_rad = .false.
+ do i = 1,nbreak-1
+    if (r > rhof_rbreak(i) .and. r < rhof_rbreak(i+1)) then
+       rho = rhof_rhobreak(i)*(r/rhof_rbreak(i))**rhof_n(i)
+       found_rad = .true.
+    endif
+ enddo
+ if (.not. found_rad) rho = rhof_rhobreak(nbreak)*(r/rhof_rbreak(nbreak))**rhof_n(nbreak)
+
+end function rho
+
+real function rho_tab(r)
+ real, intent(in) :: r
+ integer          :: i
+ real             :: logr1,logr2,logr
+ real             :: logrho1,logrho2,logrho_tab
+ real             :: gradient
+
+ rho_tab = 0.
+ do i = 1,nprof-1
+    if (r > rad_prof(i) .and. r < rad_prof(i+1)) then
+       select case (interpolation)
+       case ('log')
+          logr1 = log10(rad_prof(i))
+          logr2 = log10(rad_prof(i+1))
+          logrho1 = log10(dens_prof(i))
+          logrho2 = log10(dens_prof(i+1))
+          logr = log10(r)
+          gradient = (logrho2-logrho1)/(logr2-logr1)
+          logrho_tab = logrho1 + gradient*(logr-logr1)
+          rho_tab = 10**logrho_tab
+       case ('lin')
+          gradient = (dens_prof(i+1)-dens_prof(i))/(rad_prof(i+1)-rad_prof(i))
+          rho_tab = dens_prof(i) + gradient*(r-rad_prof(i))
+       case default
+          write(*,'(a29,1x,a)') 'Unknown interpolation option:', trim(interpolation)
+          write(*,'(a53)') "Support only 'lin'ear/'log'arithmic interpolation now"
+       end select
+    endif
+ enddo
+end function rho_tab
+
+real function get_temp_r(r,rad_prof,temp_prof)
+ real, intent(in) :: r,rad_prof(nprof),temp_prof(nprof)
+ integer :: i
+ real    :: t1,r1
+
+ get_temp_r = temperature
+ do i = 1,nprof
+    if (r > rad_prof(i) .and. r < rad_prof(i+1)) then
+       t1 = temp_prof(i)
+       r1 = rad_prof(i)
+       get_temp_r = (temp_prof(i+1)-t1)/(rad_prof(i+1)-r1)*(r-r1) + t1
+       exit
+    endif
+ enddo
+
+end function get_temp_r
+
+real function uerg(rho,T)
+ use physcon, only:kb_on_mh,radconst
+ use units,   only:unit_density,unit_ergg
+ real, intent(in) :: rho,T
+ real :: ucgs_gas,ucgs_rad,rhocgs
+
+ rhocgs = rho*unit_density
+ ucgs_gas = 1.5*kb_on_mh*T/mu
+ ucgs_rad = 0. !radconst*T**4/rhocgs
+ uerg = (ucgs_gas+ucgs_rad)/unit_ergg
+
+end function uerg
+
+subroutine calc_rhobreak()
+ integer :: i
+
+ rhof_rhobreak(1) = rhof_rho0
+ if (nbreak > 1) then
+    do i = 2,nbreak
+       rhof_rhobreak(i) = rhof_rhobreak(i-1)*(rhof_rbreak(i)/rhof_rbreak(i-1))**rhof_n(i-1)
+    enddo
+ endif
+
+end subroutine calc_rhobreak
+
+!----------------------------------------------------------------
+!+
+!  write parameters to setup file
+!+
+!----------------------------------------------------------------
+subroutine write_setupfile(filename)
+ use infile_utils, only: write_inopt
+ character(len=*), intent(in) :: filename
+ integer, parameter :: iunit = 20
+ integer            :: i
+ character(len=20)  :: rstr,nstr
+ 
+ write(*,"(a)") ' writing setup options file '//trim(filename)
+ open(unit=iunit,file=filename,status='replace',form='formatted')
+ write(iunit,"(a)") '# input file for setting up a circumnuclear gas cloud'
+
+ write(iunit,"(/,a)") '# geometry'
+ call write_inopt(ignore_radius,'ignore_radius','tde particle inside this radius will be ignored',iunit)
+ call write_inopt(remove_overlap,'remove_overlap','remove outflow particles overlap with circum particles',iunit)
+ call write_inopt(use_func,'use_func','if use broken power law for density profile',iunit)
+ if (use_func) then
+    call write_inopt(rad_min,'rad_min','inner radius of the circumnuclear gas cloud',iunit)
+    call write_inopt(rad_max,'rad_max','outer radius of the circumnuclear gas cloud',iunit)
+    write(iunit,"(/,a)") '# density broken power law'
+    call write_inopt(rhof_rho0,'rhof_rho0','density at rad_min (in g/cm^3)',iunit)
+    call write_inopt(nbreak,'nbreak','number of broken power laws',iunit)
+    write(iunit,"(/,a)") '#    section 1 (from rad_min)'
+    call write_inopt(rhof_n(1),'rhof_n_1','power law index of the section',iunit)
+    if (nbreak > 1) then
+       do i=2,nbreak
+          write(iunit,"(a,1x,i1)") '#    section',i
+          write(rstr,'(a12,i1)') 'rhof_rbreak_',i
+          write(nstr,'(a7,i1)') 'rhof_n_',i
+          call write_inopt(rhof_rbreak(i),trim(rstr),'inner radius of the section',iunit)
+          call write_inopt(rhof_n(i),trim(nstr),'power law index of the section',iunit)
+       enddo
+    endif
+ else
+    call write_inopt(profile_filename,'profile_filename','filename for the cloud profile',iunit)
+    call write_inopt(nprof,'nprof','number of data points in the cloud profile',iunit)
+    call write_inopt(interpolation,'interpolation',"use 'lin'ear/'log'arithmic interpolation between data points",iunit)
+ endif
+
+ write(iunit,"(/,a)") '# eos'
+ call write_inopt(ieos_in,'ieos','equation of state used',iunit)
+ call write_inopt(temperature,'temperature','temperature of the gas cloud (-ve = read from file)',iunit)
+ call write_inopt(mu,'mu','mean molecular density of the cloud',iunit)
+
+ close(iunit)
+
+end subroutine write_setupfile
+
+!----------------------------------------------------------------
+!+
+!  Read parameters from setup file
+!+
+!----------------------------------------------------------------
+subroutine read_setupfile(filename,ierr)
+ use infile_utils, only: open_db_from_file,inopts,read_inopt,close_db
+ use io,           only: fatal
+ character(len=*), intent(in)  :: filename
+ integer,          intent(out) :: ierr
+ integer, parameter            :: iunit=21,in_num=50
+ integer                       :: i
+ type(inopts), allocatable     :: db(:)
+ character(len=20)             :: rstr,nstr
+ real                          :: use_func_test
+
+ write(*,"(a)")'  reading setup options from '//trim(filename)
+ call open_db_from_file(db,filename,iunit,ierr)
+
+ call read_inopt(ignore_radius,'ignore_radius',db,min=0.,err=ierr)
+ call read_inopt(remove_overlap,'remove_overlap',db,err=ierr)
+ call read_inopt(use_func,'use_func',db,err=ierr)
+ call read_inopt(use_func_test,'nbreak',db,err=ierr)
+ if (ierr == -1) use_func_old = .false.
+ if (use_func) then
+    call read_inopt(rad_min,'rad_min',db,min=ignore_radius,err=ierr)
+    call read_inopt(rad_max,'rad_max',db,min=rad_min,err=ierr)
+    call read_inopt(rhof_rho0,'rhof_rho0',db,min=0.,err=ierr)
+    call read_inopt(nbreak,'nbreak',db,min=1,err=ierr)
+    allocate(rhof_rbreak_in(in_num),rhof_n_in(in_num))
+    call read_inopt(rhof_n_in(1),'rhof_n_1',db,err=ierr)
+    rhof_rbreak_in(1) = rad_min
+    do i=2,nbreak+1
+       write(rstr,'(a12,i1)') 'rhof_rbreak_',i
+       write(nstr,'(a7,i1)') 'rhof_n_',i
+       call read_inopt(rhof_rbreak_in(i),trim(rstr),db,min=rhof_rbreak_in(i-1),max=rad_max,err=ierr)
+       call read_inopt(rhof_n_in(i),trim(nstr),db,err=ierr)
+       if (ierr == 0) nbreak_old = i
+    enddo
+ else
+    call read_inopt(profile_filename,'profile_filename',db,err=ierr)
+    call read_inopt(nprof,'nprof',db,min=1,err=ierr)
+    call read_inopt(interpolation,'interpolation',db,err=ierr)
+ endif
+
+ call read_inopt(ieos_in,'ieos',db,err=ierr)
+ call read_inopt(temperature,'temperature',db,err=ierr)
+ call read_inopt(mu,'mu',db,err=ierr)
+
+ call close_db(db)
+
+ if (ierr /= 0) then
+    call fatal('moddump_radiotde','Error in reading setup file')
+ endif
+
+end subroutine read_setupfile
+!----------------------------------------------------------------
+end module moddump
+