diff --git a/src/main/cooling_stamatellos.f90 b/src/main/cooling_stamatellos.f90 index ec7d6bf14..a9aaedf31 100644 --- a/src/main/cooling_stamatellos.f90 +++ b/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,7 +66,7 @@ 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 @@ -74,7 +74,7 @@ subroutine cooling_S07(rhoi,ui,dudti_cool,xi,yi,zi,Tfloor,dudti_sph,dt,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 diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 28362db81..b8f1bff67 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/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) diff --git a/src/utils/moddump_sphNG2phantom_disc.f90 b/src/utils/moddump_sphNG2phantom_disc.f90 index 318bac487..b2db777da 100644 --- a/src/utils/moddump_sphNG2phantom_disc.f90 +++ b/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