Merge pull request danieljprice#350 from sahl95/master

Adding two addition input options with MCFOST=yes

src/main/options.f90

@@ -50,6 +50,8 @@ module options
 ! mcfost
  logical, public :: use_mcfost, use_Voronoi_limits_file, use_mcfost_stellar_parameters, mcfost_computes_Lacc
  logical, public :: mcfost_uses_PdV
+ integer, public :: ISM
+ real(kind=4), public :: mcfost_keep_part
  character(len=80), public :: Voronoi_limits_file
 
  ! radiation
@@ -147,6 +149,8 @@ subroutine set_default_options
  use_mcfost_stellar_parameters = .false.
  mcfost_computes_Lacc = .false.
  mcfost_uses_PdV = .true.
+ mcfost_keep_part = 0.999
+ ISM = 0
 
  ! radiation
  if (do_radiation) then

src/main/readwrite_infile.F90

@@ -76,7 +76,7 @@ module readwrite_infile
                      icooling,psidecayfac,overcleanfac,hdivbbmax_max,alphamax,calc_erot,rhofinal_cgs, &
                      use_mcfost,use_Voronoi_limits_file,Voronoi_limits_file,use_mcfost_stellar_parameters,&
                      exchange_radiation_energy,limit_radiation_flux,iopacity_type,mcfost_computes_Lacc,&
-                     mcfost_uses_PdV,implicit_radiation
+                     mcfost_uses_PdV,implicit_radiation,mcfost_keep_part,ISM
  use timestep,  only:dtwallmax,tolv,xtol,ptol
  use viscosity, only:irealvisc,shearparam,bulkvisc
  use part,      only:hfact,ien_type
@@ -245,6 +245,10 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint)
       'Should mcfost compute the accretion luminosity',iwritein)
  call write_inopt(mcfost_uses_PdV,'mcfost_uses_PdV',&
       'Should mcfost use the PdV work and shock heating?',iwritein)
+ call write_inopt(mcfost_keep_part,'mcfost_keep_part',&
+      'Fraction of particles to keep for MCFOST',iwritein)
+      call write_inopt(ISM,'ISM',&
+      'ISM heating : 0 -> no ISM radiation field, 1 -> ProDiMo, 2 -> Bate & Keto',iwritein)
 #endif
 
  ! only write sink options if they are used, or if self-gravity is on
@@ -519,6 +523,10 @@ subroutine read_infile(infile,logfile,evfile,dumpfile)
        read(valstring,*,iostat=ierr) mcfost_computes_Lacc
     case('mcfost_uses_PdV')
        read(valstring,*,iostat=ierr) mcfost_uses_PdV
+    case('mcfost_keep_part')
+       read(valstring,*,iostat=ierr) mcfost_keep_part
+    case('ISM')
+         read(valstring,*,iostat=ierr) ISM
 #endif
     case('implicit_radiation')
        read(valstring,*,iostat=ierr) implicit_radiation

src/utils/analysis_mcfost.f90

@@ -43,7 +43,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
  use dim,            only:use_dust,lightcurve,maxdusttypes,use_dustgrowth,do_radiation
  use eos,            only:temperature_coef,gmw,gamma
  use options,        only:use_dustfrac,use_mcfost,use_Voronoi_limits_file,Voronoi_limits_file, &
-                             use_mcfost_stellar_parameters, mcfost_computes_Lacc, mcfost_uses_PdV
+                             use_mcfost_stellar_parameters, mcfost_computes_Lacc, mcfost_uses_PdV,&
+                             mcfost_keep_part, ISM
  use physcon,        only:cm,gram,c,steboltz
 
  character(len=*), intent(in)    :: dumpfile
@@ -58,11 +59,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
  integer         :: ierr,ntypes,dustfluidtype,ilen,nlum,i,nerr
  integer(kind=1) :: itype(maxp)
  logical         :: compute_Frad
+ logical         :: ISM_heating = .false.
  real(kind=8), dimension(6), save            :: SPH_limits
  real,         dimension(:),     allocatable :: dudt
  real,    parameter :: Tdefault = 1.
  logical, parameter :: write_T_files = .false. ! ask mcfost to write fits files with temperature structure
- integer, parameter :: ISM = 2 ! ISM heating : 0 -> no ISM radiation field, 1 -> ProDiMo, 2 -> Bate & Keto
  character(len=len(dumpfile) + 20) :: mcfost_para_filename
  real :: a_code,rhoi,pmassi,Tmin,Tmax,default_kappa,kappa_diffusion
 
@@ -73,12 +74,17 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
     call growth_to_fake_multi(npart)
  endif
 
+ if (ISM > 0) then
+   ISM_heating = .true.
+ endif
+
  if (.not.init_mcfost) then
     ilen = index(dumpfile,'_',back=.true.) ! last position of the '_' character
     mcfost_para_filename = dumpfile(1:ilen-1)//'.para'
     call init_mcfost_phantom(mcfost_para_filename,ndusttypes,use_Voronoi_limits_file,&
          Voronoi_limits_file,SPH_limits,ierr, fix_star = use_mcfost_stellar_parameters, &
-         turn_on_Lacc = mcfost_computes_Lacc)
+         turn_on_Lacc = mcfost_computes_Lacc, keep_particles = mcfost_keep_part, &
+         use_ISM_heating = ISM_heating)
     if (ierr /= 0) call fatal('mcfost-phantom','error in init_mcfost_phantom')
     init_mcfost = .true.
  endif
@@ -126,7 +132,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
 
  Tmin = minval(Tdust, mask=(Tdust > 1.))
  Tmax = maxval(Tdust)
-
  write(*,*) ''
  write(*,*) 'Minimum temperature = ', Tmin
  write(*,*) 'Maximum temperature = ', Tmax