Adjust stellar heating estimate in cooling_stamatellos.f90 and write …

…extra quantities to dump files

src/main/cooling_stamatellos.f90

@@ -22,7 +22,7 @@ module cooling_stamatellos
  implicit none
  real, public :: Lstar ! in units of L_sun
  integer :: isink_star ! index of sink to use as illuminating star
- integer :: od_method = 1 ! default = Stamatellos+ 2007 method
+ integer :: od_method = 4 ! default = Stamatellos+ 2007 method
  integer :: fld_opt = 1 ! by default FLD is switched on
  public :: cooling_S07,write_options_cooling_stamatellos,read_options_cooling_stamatellos
  public :: init_star
@@ -66,15 +66,15 @@ subroutine cooling_S07(rhoi,ui,dudti_cool,xi,yi,zi,Tfloor,dudti_sph,dt,i)
  use physcon,  only:steboltz,pi,solarl,Rg,kb_on_mh,piontwo,rpiontwo
  use units,    only:umass,udist,unit_density,unit_ergg,utime,unit_pressure
  use eos_stamatellos, only:getopac_opdep,getintenerg_opdep,gradP_cool,Gpot_cool,&
-          duFLD,doFLD
+          duFLD,doFLD,ttherm_store,teqi_store,floor_energy
  use part,       only:xyzmh_ptmass
 
  real,intent(in) :: rhoi,ui,dudti_sph,xi,yi,zi,Tfloor,dt
  integer,intent(in) :: i
  real,intent(out) :: dudti_cool
  real            :: coldensi,kappaBari,kappaParti,ri2
  real            :: gmwi,Tmini4,Ti,dudt_rad,Teqi,Hstam,HLom,du_tot
- real            :: cs2,Om2,Q3D,Hmod2
+ real            :: cs2,Om2,Hmod2
  real            :: tcool,ueqi,umini,tthermi,poti,presi,Hcomb,du_FLDi
 
  poti = Gpot_cool(i)
@@ -86,15 +86,6 @@ subroutine cooling_S07(rhoi,ui,dudti_cool,xi,yi,zi,Tfloor,dudti_sph,dt,i)
          + (zi-xyzmh_ptmass(3,isink_star))**2d0  
  endif
 
-!    Tfloor is from input parameters and is background heating
-!    Stellar heating
- if (isink_star > 0 .and. Lstar > 0.d0) then
-! Tfloor + stellar heating
-    Tmini4 = Tfloor**4d0 + (Lstar*solarl/(16d0*pi*steboltz*ri2*udist*udist))
- else
-    Tmini4 = Tfloor**4d0
- endif
-
 ! get opacities & Ti for ui
  call getopac_opdep(ui*unit_ergg,rhoi*unit_density,kappaBari,kappaParti,&
            Ti,gmwi)
@@ -140,6 +131,15 @@ subroutine cooling_S07(rhoi,ui,dudti_cool,xi,yi,zi,Tfloor,dudti_sph,dt,i)
     coldensi = 1.014d0 * Hcomb *rhoi*umass/udist/udist ! physical units
  end select
 
+!    Tfloor is from input parameters and is background heating
+!    Stellar heating
+ if (isink_star > 0 .and. Lstar > 0.d0) then
+! Tfloor + stellar heating
+    Tmini4 = Tfloor**4d0 + exp(-coldensi*kappaBari)*(Lstar*solarl/(16d0*pi*steboltz*ri2*udist*udist))
+ else
+    Tmini4 = Tfloor**4d0
+ endif
+
  tcool = (coldensi**2d0)*kappaBari + (1.d0/kappaParti) ! physical units
  dudt_rad = 4.d0*steboltz*(Tmini4 - Ti**4.d0)/tcool/unit_ergg*utime! code units
 
@@ -154,12 +154,12 @@ subroutine cooling_S07(rhoi,ui,dudti_cool,xi,yi,zi,Tfloor,dudti_sph,dt,i)
 
  Teqi = Teqi/4.d0/steboltz
  Teqi = Teqi + Tmini4
- if (Teqi < Tmini4) then
-    Teqi = Tmini4**(1.0/4.0)
- else
+! if (Teqi < Tmini4) then
+!    Teqi = Tmini4**(1.0/4.0)
+! else
     Teqi = Teqi**(1.0/4.0)
- endif
-
+! endif
+ teqi_store(i) = Teqi
  call getintenerg_opdep(Teqi,rhoi*unit_density,ueqi)
  ueqi = ueqi/unit_ergg
 
@@ -174,6 +174,8 @@ subroutine cooling_S07(rhoi,ui,dudti_cool,xi,yi,zi,Tfloor,dudti_sph,dt,i)
     tthermi = abs((ueqi - ui)/(du_tot))
  endif
 
+ ttherm_store(i) = tthermi
+
  if (tthermi == 0d0) then
     dudti_cool = 0.d0 ! condition if denominator above is zero
  else

src/main/readwrite_dumps_fortran.F90

@@ -244,7 +244,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG)
  use krome_user, only:krome_nmols
  use part,       only:gamma_chem,mu_chem,T_gas_cool
 #endif
- use eos_stamatellos, only:gradP_cool,Gpot_cool,doFLD,urad_FLD
+ use eos_stamatellos, only:gradP_cool,doFLD,urad_FLD,ttherm_store,teqi_store
  real,             intent(in) :: t
  character(len=*), intent(in) :: dumpfile
  integer,          intent(in), optional :: iorder(:)
@@ -409,9 +409,11 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG)
           endif
        endif
        ! write urad to file (stamatellos + FLD)
-!       if (icooling == 8 .and. doFLD) then
+       if (icooling == 8) then ! .and. doFLD) then
 !         call write_array(1,urad_FLD,'urad',npart,k,ipass,idump,nums,ierrs(13))
-!      endif
+         call write_array(1,teqi_store,'teqi',npart,k,ipass,idump,nums,ierrs(13))
+         call write_array(1,ttherm_store,'ttherm',npart,k,ipass,idump,nums,ierrs(13))
+      endif
 
        ! smoothing length written as real*4 to save disk space
        call write_array(1,xyzh,xyzh_label,1,npart,k,ipass,idump,nums,ierrs(14),use_kind=4,index=4)

src/utils/moddump_sphNG2phantom_disc.f90

@@ -32,7 +32,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  use prompting, only:prompt
  use physcon,  only:au,gg
  use readwrite_dumps_fortran, only:dt_read_in_fortran
- use timestep, only:time,dt,dtmax_max,dtmax_min,dtmax0
+ use timestep, only:time,dt,dtmax_max,dtmax_min
  use centreofmass, only: reset_centreofmass
  integer, intent(inout) :: npart
  integer, intent(inout) :: npartoftype(:)
@@ -97,7 +97,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  enddo
  close(iunit)
 
- print *, 'dtmax0, dtmax_max,dtmax_min',dtmax0,dtmax_max,dtmax_min
+ print *, 'dtmax_max,dtmax_min',dtmax_max,dtmax_min
  newutime = sqrt(au**3/(gg*umass))
  print *, "newutime/old", newutime/utime
  time = time * utime / newutime
@@ -177,7 +177,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
       'nptmass:', nptmass
  print *, 'gamma=', gamma
  print *, 'Timestep info:'
- print *, 'dtmax0, dtmax_max,dtmax_min', dtmax0,dtmax_max,dtmax_min
+ print *, 'dtmax_max,dtmax_min', dtmax_max,dtmax_min
  print *, 'utime=', utime
 
  return