From da2091ccc29807ea29dc1e404b53c843d1d99f6a Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Fri, 7 Apr 2023 16:12:31 +0200
Subject: [PATCH 1/9] Implement dust calculations in
 analysis_common_envelope.F90

---
 src/main/dust_formation.F90            |  10 +-
 src/utils/analysis_common_envelope.F90 | 219 ++++++++++++++++++++-----
 2 files changed, 180 insertions(+), 49 deletions(-)

diff --git a/src/main/dust_formation.F90 b/src/main/dust_formation.F90
index cb55a2d4c..5e53dc246 100644
--- a/src/main/dust_formation.F90
+++ b/src/main/dust_formation.F90
@@ -32,7 +32,7 @@ module dust_formation
  integer, public :: idust_opacity = 0
 
  public :: set_abundances,evolve_dust,evolve_chem,calc_kappa_dust,&
-      calc_kappa_bowen,chemical_equilibrium_light,&
+      calc_kappa_bowen,chemical_equilibrium_light,psat_C,calc_nucleation,&
       read_options_dust_formation,write_options_dust_formation,&
       calc_Eddington_factor,calc_muGamma,init_muGamma,init_nucleation,&
       write_headeropts_dust_formation,read_headeropts_dust_formation
@@ -41,6 +41,7 @@ module dust_formation
 !
 
  real, public :: kappa_gas   = 2.d-4
+ real, public, parameter :: Scrit = 2. ! Critical saturation ratio
 
  private
 
@@ -82,7 +83,6 @@ module dust_formation
       -4.38897d+05, -1.58111d+05, 2.49224d+01, 1.08714d-03, -5.62504d-08, & !TiO
       -3.32351d+05, -3.04694d+05, 5.86984d+01, 1.17096d-03, -5.06729d-08, & !TiO2
        2.26786d+05, -1.43775d+05, 2.92429d+01, 1.69434d-04, -1.79867d-08], shape(coefs)) !C2
- real, parameter :: Scrit = 2. ! Critical saturation ratio
  real, parameter :: vfactor = sqrt(kboltz/(2.*pi*atomic_mass_unit*12.01))
  !real, parameter :: vfactor = sqrt(kboltz/(8.*pi*atomic_mass_unit*12.01))
 
@@ -179,7 +179,7 @@ subroutine evolve_chem(dt, T, rho_cgs, JKmuS)
     S = pC/psat_C(T)
     if (S > Scrit) then
        !call nucleation(T, pC, pC2, 0., pC2H, pC2H2, S, JstarS, taustar, taugr)
-       call nucleation(T, pC, 0., 0., 0., pC2H2, S, JstarS, taustar, taugr)
+       call calc_nucleation(T, pC, 0., 0., 0., pC2H2, S, JstarS, taustar, taugr)
        JstarS = JstarS/ nH_tot
        call evol_K(JKmuS(idJstar), JKmuS(idK0:idK3), JstarS, taustar, taugr, dt, Jstar_new, K_new)
     else
@@ -272,7 +272,7 @@ end function calc_Eddington_factor
 !  Compute nucleation rate
 !
 !----------------------------
-subroutine nucleation(T, pC, pC2, pC3, pC2H, pC2H2, S, JstarS, taustar, taugr)
+subroutine calc_nucleation(T, pC, pC2, pC3, pC2H, pC2H2, S, JstarS, taustar, taugr)
 ! all quantities are in cgs
  real, intent(in)  :: T, pC, pC2, pC3, pC2H, pC2H2, S
  real, intent(out) :: JstarS, taustar, taugr
@@ -320,7 +320,7 @@ subroutine nucleation(T, pC, pC2, pC3, pC2H, pC2H2, S, JstarS, taustar, taugr)
     taustar = 1.d-30
  endif
  taugr = kboltz*T/(A0*v1*(alpha1*pC*(1.-1./S) + 2.*alpha2/sqrt(2.)*(pC2+pC2H+pC2H2)*(1.-1./S**2)))
-end subroutine nucleation
+end subroutine calc_nucleation
 
 !------------------------------------
 !
diff --git a/src/utils/analysis_common_envelope.F90 b/src/utils/analysis_common_envelope.F90
index 89f6fd72c..d76271588 100644
--- a/src/utils/analysis_common_envelope.F90
+++ b/src/utils/analysis_common_envelope.F90
@@ -83,7 +83,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
 
  !chose analysis type
  if (dump_number==0) then
-    print "(40(a,/))", &
+    print "(41(a,/))", &
             ' 1) Sink separation', &
             ' 2) Bound and unbound quantities', &
             ' 3) Energies', &
@@ -122,9 +122,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
             '37) Planet profile',&
             '38) Velocity profile',&
             '39) Angular momentum profile',&
-            '40) Keplerian velocity profile'
+            '40) Keplerian velocity profile',&
+            '41) Total dust mass'
     analysis_to_perform = 1
-    call prompt('Choose analysis type ',analysis_to_perform,1,40)
+    call prompt('Choose analysis type ',analysis_to_perform,1,41)
  endif
 
  call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
@@ -132,7 +133,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
                                    xyzmh_ptmass,vxyz_ptmass,omega_corotate,dump_number)
 
  ! List of analysis options that require specifying EOS options
- requires_eos_opts = any((/2,3,4,6,8,9,11,13,14,15,20,21,22,23,24,25,26,29,30,31,32,33,35/) == analysis_to_perform)
+ requires_eos_opts = any((/2,3,4,6,8,9,11,13,14,15,20,21,22,23,24,25,26,29,30,31,32,33,35,41/) == analysis_to_perform)
  if (dump_number == 0 .and. requires_eos_opts) call set_eos_options(analysis_to_perform)
 
  select case(analysis_to_perform)
@@ -200,6 +201,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
     call angular_momentum_profile(time,num,npart,particlemass,xyzh,vxyzu)
  case(40) ! Keplerian velocity profile
     call vkep_profile(time,num,npart,particlemass,xyzh,vxyzu)
+  case(41) !Total dust mass
+    call total_dust_mass(time,npart,particlemass,xyzh)
  case(12) !sink properties
     call sink_properties(time,npart,particlemass,xyzh,vxyzu)
  case(13) !MESA EoS compute total entropy and other average thermodynamical quantities
@@ -302,38 +305,97 @@ end subroutine do_analysis
 
 
 subroutine total_dust_mass(time,npart,particlemass,xyzh)
- use part,           only:nucleation,idK3
+ use part,           only:nucleation,idK3,idK0,idK1, idJstar
  use dust_formation, only:set_abundances, mass_per_H
  use physcon, only:atomic_mass_unit
  real, intent(in)               :: time,particlemass,xyzh(:,:)
  integer, intent(in)            :: npart
- integer                        :: i,ncols
- real, dimension(1)            :: dust_mass
+ integer                        :: i,ncols,j
+ real, dimension(2)             :: dust_mass
  character(len=17), allocatable :: columns(:)
+ real, allocatable              :: temp(:) !npart
+ real                           :: median,mass_factor,grain_size
+ real, parameter :: a0 = 1.28e-4 !radius of a carbon atom in micron
 
- call set_abundances !without the calling, the parameter mass_per_H is zero
- dust_mass(1) = 0
- ncols = 1
+ call set_abundances !initialize mass_per_H
+ dust_mass = 0.
+ ncols = 2
  print *,'size(nucleation,1) = ',size(nucleation,1)
  print *,'size(nucleation,2) = ',size(nucleation,2)
- allocate(columns(ncols))
- columns = (/'total dust mass'/)
+ allocate(columns(ncols),temp(npart))
+ columns = (/'Dust mass [Msun]', &
+             'median size [um]'/)
+ j=0
+ mass_factor = 12.*atomic_mass_unit*particlemass/mass_per_H
  do i = 1,npart
-    if (.not. isdead_or_accreted(xyzh(4,i))) then
-       dust_mass(1) = dust_mass(1) + nucleation(idK3,i) &
-                         * 12*atomic_mass_unit*particlemass*2.0E+33/mass_per_H
-       !the factor 2.0E+33 convert particlemass from solar units to cgs
-       !12*atomic_mass_unit is the mass of a Carbon atom
-    endif
+   if (.not. isdead_or_accreted(xyzh(4,i))) then
+     dust_mass(1) = dust_mass(1) + nucleation(idK3,i) *mass_factor
+     grain_size = a0*nucleation(idK1,i)/(nucleation(idK0,i)+1.0E-99) !in micron
+     if (grain_size > a0) then
+        j = j+1
+        temp(j) = grain_size
+     endif
+   endif
  enddo
 
+ call sort(temp,j)
+ if (mod(j,2)==0) then !npart
+   median = (temp(j/2)+temp(j/2+1))/2.0 !(temp(npart/2)+temp(npart/2+1))/2.0
+ else
+   median = (temp(j/2)+temp(j/2+1))/2.0 !temp(npart/2+1)
+ endif
+
+ dust_mass(2) = median
+
  call write_time_file('total_dust_mass_vs_time', columns, time, dust_mass, ncols, dump_number)
  !after execution of the analysis routine, a file named "total_dust_mass_vs_time.ev" appears
- deallocate(columns)
+ deallocate(columns,temp)
 
 end subroutine total_dust_mass
 
-
+! --------------------------------------------------------------------
+! integer function  FindMinimum():
+!    This function returns the location of the minimum in the section
+! between Start and End.
+! --------------------------------------------------------------------
+
+integer function  FindMinimum(x, Start, Fin)
+ implicit  none
+ integer, intent(in)                   :: start, fin
+ real, dimension(Fin), intent(in) :: x
+ real                            :: minimum
+ integer                            :: location
+ integer                            :: i
+
+ minimum  = x(start)          ! assume the first is the min
+ location = start             ! record its position
+ do i = start+1, fin          ! start with next elements
+    if (x(i) < minimum) then  !   if x(i) less than the min?
+       minimum  = x(i)        !      yes, a new minimum found
+       location = i                !      record its position
+    end if
+ end do
+ findminimum = location            ! return the position
+end function  findminimum
+
+! --------------------------------------------------------------------
+! subroutine  Sort():
+!    This subroutine receives an array x() and sorts it into ascending
+! order.
+! --------------------------------------------------------------------
+
+subroutine  Sort(x, longitud)
+ implicit  none
+ integer, intent(in)                   :: longitud
+ real, dimension(longitud), intent(inout) :: x
+ integer                               :: i
+ integer                               :: location
+
+ do i = 1, longitud-1             ! except for the last
+    location = findminimum(x, i, longitud)  ! find min from this to last
+    call swap(x(i), x(location))  ! swap this and the minimum
+ end do
+end subroutine  Sort
 
 
 !----------------------------------------------------------------
@@ -465,7 +527,7 @@ subroutine planet_rvm(time,particlemass,xyzh,vxyzu)
                 mass(1), mass(2), mass(3), mass(4), mass(5), rhoprev, smin /)
  call write_time_file('planet_rvm', columns, time, data_cols, ncols, dump_number)
 
- deallocate(data_cols,columns,vthreshold,mass)
+ deallocate(data_cols,columns,mass,vthreshold)
 
 end subroutine planet_rvm
 
@@ -934,8 +996,7 @@ subroutine create_profile(time, num, npart, particlemass, xyzh, vxyzu)
 
  call write_file(name_in, 'profile', columns, profile, size(profile(1,:)), ncols, num)
 
- deallocate(profile)
- deallocate(columns)
+ deallocate(profile,columns)
 end subroutine create_profile
 
 
@@ -1143,7 +1204,7 @@ subroutine roche_lobe_values(time,npart,particlemass,xyzh,vxyzu)
              ' Fallback Jz'/)
 
  call write_time_file('roche_lobes', columns, time, MRL, ncols, dump_number)
- deallocate(columns)
+ deallocate(columns,iorder)
 
 end subroutine roche_lobe_values
 
@@ -1266,7 +1327,7 @@ subroutine star_stabilisation_suite(time,npart,particlemass,xyzh,vxyzu)
 
  star_stability(imassfracout) = star_stability(imassout) / total_mass
  call write_time_file('star_stability', columns, time, star_stability, ncols, dump_number)
- deallocate(columns)
+ deallocate(columns,star_stability,iorder,iorder_a)
 
 end subroutine star_stabilisation_suite
 
@@ -1305,12 +1366,14 @@ end subroutine print_simulation_parameters
 !+
 !----------------------------------------------------------------
 subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
- use part,              only:eos_vars,itemp
+ use part,              only:eos_vars,itemp,nucleation,idK0,idK1,idK2,idK3,idJstar,idmu,idgamma
  use eos,               only:entropy
  use eos_mesa,          only:get_eos_kappa_mesa
  use mesa_microphysics, only:getvalue_mesa
  use sortutils,         only:set_r2func_origin,r2func_origin,indexxfunc
  use ionization_mod,    only:calc_thermal_energy,ionisation_fraction
+ use dust_formation,    only:psat_C,eps,set_abundances,mass_per_H, chemical_equilibrium_light, calc_nucleation!, Scrit
+ !use dim,     only:nElements
  integer, intent(in)          :: npart
  character(len=*), intent(in) :: dumpfile
  real, intent(in)             :: time,particlemass
@@ -1324,11 +1387,15 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
  real, allocatable, save      :: init_entropy(:)
  real, allocatable            :: quant(:,:)
  real, dimension(3)           :: com_xyz,com_vxyz,xyz_a,vxyz_a
+ real                         :: pC, pC2, pC2H, pC2H2, nH_tot, epsC, S
+ real                         :: taustar, taugr, JstarS
+ real, parameter :: Scrit = 2. ! Critical saturation ratio
+ logical :: verbose = .false.
 
  allocate(quant(4,npart))
- Nquantities = 12
+ Nquantities = 13
  if (dump_number == 0) then
-    print "(12(a,/))",&
+    print "(13(a,/))",&
            '1) Total energy (kin + pot + therm)', &
            '2) Mach number', &
            '3) Opacity from MESA tables', &
@@ -1340,7 +1407,8 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
            '9) Total energy (kin + pot)', &
            '10) Mass coordinate', &
            '11) Gas omega w.r.t. CoM', &
-           '12) Gas omega w.r.t. sink 1'
+           '12) Gas omega w.r.t. sink 1',&
+           '13) JstarS' !option to calculate JstarS
 
     quantities_to_calculate = (/1,2,4,5/)
     call prompt('Choose first quantity to compute ',quantities_to_calculate(1),1,Nquantities)
@@ -1356,13 +1424,13 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
  com_vxyz = 0.
  do k=1,4
     select case (quantities_to_calculate(k))
-    case(1,2,3,6,8,9) ! Nothing to do
+    case(1,2,3,6,8,9,13) ! Nothing to do
     case(4,5,11,12) ! Fractional difference between gas and orbital omega
        if (quantities_to_calculate(k) == 4 .or. quantities_to_calculate(k) == 5) then
-          com_xyz  = (xyzmh_ptmass(1:3,1)*xyzmh_ptmass(4,1) + xyzmh_ptmass(1:3,2)*xyzmh_ptmass(4,2)) &
-                     / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2))
-          com_vxyz = (vxyz_ptmass(1:3,1)*xyzmh_ptmass(4,1)  + vxyz_ptmass(1:3,2)*xyzmh_ptmass(4,2))  &
-                     / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2))
+       com_xyz  = (xyzmh_ptmass(1:3,1)*xyzmh_ptmass(4,1) + xyzmh_ptmass(1:3,2)*xyzmh_ptmass(4,2)) &
+                  / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2))
+       com_vxyz = (vxyz_ptmass(1:3,1)*xyzmh_ptmass(4,1)  + vxyz_ptmass(1:3,2)*xyzmh_ptmass(4,2))  &
+                  / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2))
        elseif (quantities_to_calculate(k) == 11 .or. quantities_to_calculate(k) == 12) then
           com_xyz = xyzmh_ptmass(1:3,1)
           com_vxyz = vxyz_ptmass(1:3,1)
@@ -1378,16 +1446,60 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
        call set_r2func_origin(0.,0.,0.)
        allocate(iorder(npart))
        call indexxfunc(npart,r2func_origin,xyzh,iorder)
+       deallocate(iorder)
     case default
        print*,"Error: Requested quantity is invalid."
        stop
     end select
  enddo
 
+ !set initial abundances to get mass_per_H
+ call set_abundances
  ! Calculations performed in loop over particles
  do i=1,npart
     do k=1,4
        select case (quantities_to_calculate(k))
+       case(13) !to calculate JstarS
+          rhopart = rhoh(xyzh(4,i), particlemass)
+          !call equationofstate to obtain temperature and store it in tempi
+          call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i))
+          JstarS = 0.
+          !nH_tot is needed to normalize JstarS
+          nH_tot = rhopart/mass_per_H
+          epsC   = eps(3) - nucleation(idK3,i)
+          if (epsC < 0.) then
+             print *,'eps(C) =',eps(3),', K3=',nucleation(idK3,i),', epsC=',epsC,', T=',tempi,' rho=',rhopart
+             print *,'JKmuS=',nucleation(:,i)
+             stop '[S-dust_formation] epsC < 0!'
+          endif
+          if (tempi > 450.) then
+             !call chemical_equilibrium_light to obtain pC, and pC2H2
+             call chemical_equilibrium_light(rhopart, tempi, epsC, pC, pC2, pC2H, pC2H2, nucleation(idmu,i), nucleation(idgamma,i))
+             S = pC/psat_C(tempi)
+             if (S > Scrit) then
+                !call nucleation_function to obtain JstarS
+                call calc_nucleation(tempi, pC, 0., 0., 0., pC2H2, S, JstarS, taustar, taugr)
+                JstarS = JstarS/ nH_tot
+             endif
+          endif
+          !Check if the variables have meaningful values close to condensation temperatures
+          if (tempi.ge.1400. .and. tempi.le.1500. .and. verbose ) then
+             print *,'size(nucleation,1) = ',size(nucleation,1)
+             print *,'size(nucleation,2) = ',size(nucleation,2)
+             print *,'nucleation(idK3,i) = ',nucleation(idK3,i)
+             print *,'epsC = ',epsC
+             print *,'tempi = ',tempi
+             print *,'S = ',S
+             print *,'pC =',pC
+             print *,'psat_C(tempi) = ',psat_C(tempi)
+             print *,'nucleation(idmu,i) = ',nucleation(idmu,i)
+             print *,'nucleation(idgamma,i) = ',nucleation(idgamma,i)
+             print *,'taustar = ',taustar
+             print *,'eps = ',eps
+	     print *,'JstarS = ',JstarS
+          endif
+          quant(k,i) = JstarS
+
        case(1,9) ! Total energy (kin + pot + therm)
           rhopart = rhoh(xyzh(4,i), particlemass)
           call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i))
@@ -1467,6 +1579,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
     write(iu) (quant(k,i),i=1,npart)
  enddo
  close(iu)
+ deallocate(quant)
 
 end subroutine output_divv_files
 
@@ -1676,6 +1789,8 @@ subroutine tau_profile(time,num,npart,particlemass,xyzh)
  open(unit=unitnum, file=trim(adjustl(filename)), position='append')
  write(unitnum,data_formatter) time,tau_r
  close(unit=unitnum)
+ deallocate(rad_part,kappa_part,rho_part)
+ deallocate(rho_hist,kappa_hist,sepbins,tau_r)
 
 end subroutine tau_profile
 
@@ -1837,6 +1952,7 @@ subroutine recombination_tau(time,npart,particlemass,xyzh,vxyzu)
        call write_time_file("recombination_tau",(/'          tau'/),-1.,tau_recombined(i),1,i-1) ! Set num = i-1 so that header will be written for particle 1 and particle 1 only
     enddo
  endif
+ deallocate(recombined_pid,rad_part,kappa_part,rho_part)
 
 end subroutine recombination_tau
 
@@ -1899,6 +2015,7 @@ subroutine energy_hist(time,npart,particlemass,xyzh,vxyzu)
     write(unitnum,data_formatter) time,hist
     close(unit=unitnum)
  enddo
+ deallocate(filename,coord,hist,Emin,Emax,quant)
 
 end subroutine energy_hist
 
@@ -2061,6 +2178,8 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu)
     write(unitnum,data_formatter) time,hist
     close(unit=unitnum)
  enddo
+ deallocate(iorder,coord,headerline,filename,quant,hist)
+
 end subroutine energy_profile
 
 
@@ -2136,6 +2255,7 @@ subroutine rotation_profile(time,num,npart,xyzh,vxyzu)
     write(unitnum,data_formatter) time,hist_var(:)
     close(unit=unitnum)
  enddo
+ deallocate(hist_var,grid_file,dist_part,rad_part)
 
 end subroutine rotation_profile
 
@@ -2254,6 +2374,7 @@ subroutine velocity_profile(time,num,npart,particlemass,xyzh,vxyzu)
  open(newunit=iu, file=trim(adjustl(file_name)), position='append')
  write(iu,data_formatter) time,hist
  close(unit=iu)
+ deallocate(hist,dist_part,rad_part)
 
 end subroutine velocity_profile
 
@@ -2363,6 +2484,7 @@ subroutine vkep_profile(time,num,npart,particlemass,xyzh,vxyzu)
  open(newunit=iu, file=trim(adjustl(file_name)), position='append')
  write(iu,data_formatter) time,hist
  close(unit=iu)
+ deallocate(hist,dist_part,rad_part)
 
 end subroutine vkep_profile
 
@@ -2416,6 +2538,7 @@ subroutine planet_profile(num,dumpfile,particlemass,xyzh,vxyzu)
  enddo
 
  close(unit=iu)
+ deallocate(R,z,rho)
 
 end subroutine planet_profile
 
@@ -2534,6 +2657,8 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu)
 
     close(unit=unitnum)
  enddo
+ deallocate(hist_var)
+
 end subroutine unbound_profiles
 
 
@@ -3276,7 +3401,7 @@ subroutine gravitational_drag(time,npart,particlemass,xyzh,vxyzu)
     write (filename, "(A16,I0)") "sink_drag_", i
     call write_time_file(trim(adjustl(filename)), columns, time, drag_force(:,i), ncols, dump_number)
  enddo
- deallocate(columns)
+ deallocate(columns,drag_force,force_cut_vec,Rcut)
 
 end subroutine gravitational_drag
 
@@ -3314,7 +3439,7 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu)
  data(1,:) = data(1,:) / particlemass ! specific energy
 
  call write_file('JEplane','JEplane',columns,data,size(data(1,:)),ncols,num)
- deallocate(columns)
+ deallocate(columns,data)
 
 end subroutine J_E_plane
 
@@ -3332,7 +3457,7 @@ subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu)
  character(len=18)                :: filename
  real, allocatable                :: planetDestruction(:)
  integer                          :: ncols,i,j
- real, allocatable, save          :: time_old
+ real, save                       :: time_old
  real, allocatable, save          :: particleRho(:)
  character(len=50)                :: planetRadiusPromptString
  real, allocatable, save          :: planetRadii(:) !In units of Rsun
@@ -3344,8 +3469,7 @@ subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu)
  real, allocatable, save          :: currentKhAblatedMass(:)
 
  ncols=5
- allocate(columns(ncols))
- allocate(planetDestruction(ncols))
+ allocate(columns(ncols),planetDestruction(ncols))
  columns=(/"      rhoGas", &
            "  kh_rhoCrit", &
            "     kh_lmax", &
@@ -3372,7 +3496,6 @@ subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu)
           call prompt(planetRadiusPromptString,planetRadii(i),0.0,1.0)
        enddo
 
-       allocate(time_old)
        allocate(particleRho(npart))
        allocate(currentKhAblatedMass(nptmass))
 
@@ -3409,6 +3532,8 @@ subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu)
  enddo
 
  time_old=time
+
+ deallocate(columns,planetDestruction)
 end subroutine planet_destruction
 
 !-----------------------------------------------------------------------------------------
@@ -3461,8 +3586,9 @@ subroutine create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu)
     profile(3,i)=previousBindingEnergyU+currentParticleGPE-(vxyzu(4,j)*particlemass)
  enddo
 
-
  call write_file('bEnergyProfile','bEnergyProfiles',columns,profile,npart,ncols,num)
+ deallocate(columns,iorder,profile)
+
 end subroutine create_bindingEnergy_profile
 
 
@@ -3880,6 +4006,7 @@ subroutine stellar_profile(time,ncols,particlemass,npart,xyzh,vxyzu,profile,simp
     if (i > 1) profile(2,i) = profile(2,i-1) + particlemass
  enddo
 
+ deallocate(profile)
  print*, "Profile completed"
 
 end subroutine stellar_profile
@@ -3923,6 +4050,7 @@ subroutine get_interior_mass(xyzh,vxyzu,donor_xyzm,companion_xyzm,particlemass,n
  interior_mass = npart_int * particlemass
 
  call get_centreofmass(com_xyz,com_vxyz,npart_int,xyz_int,vxyz_int,nptmass,xyzmh_ptmass,vxyz_ptmass)
+ deallocate(iorder)
 
 end subroutine get_interior_mass
 
@@ -3974,6 +4102,7 @@ subroutine orbit_com(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass,com_xyz,c
  enddo
  npart_a = k - 1
  call get_centreofmass(com_xyz,com_vxyz,npart_a,xyz_a,vxyz_a,nptmass,xyzmh_ptmass,vxyz_ptmass)
+ deallocate(iorder,xyz_a,vxyz_a)
 
 end subroutine orbit_com
 
@@ -4282,9 +4411,9 @@ function sphInterpolation(npart,particlemass,particleRho,particleXyzh,interpolat
 
  integer              :: i,j
  integer, allocatable :: iorder(:)
- real                 :: currentR,currentQ,currentQ2
- real                 :: nearestSphH
- real                 :: currentParticleRho,currentSphSummandFactor
+ real                :: currentR,currentQ,currentQ2
+ real                :: nearestSphH
+ real                :: currentParticleRho,currentSphSummandFactor
 
  interpolatedData=0.0
  allocate(iorder(npart))
@@ -4309,6 +4438,8 @@ function sphInterpolation(npart,particlemass,particleRho,particleXyzh,interpolat
     currentSphSummandFactor=(particlemass/currentParticleRho)*((1.0/((nearestSphH**3.0)*pi))*wkern(currentQ2,currentQ))
     interpolatedData=interpolatedData+(currentSphSummandFactor*toInterpolate(:,j))
  enddo
+ deallocate(iorder)
+
 end function sphInterpolation
 
 !Sorting routines

From c8c7d75adc64b8e8f821f9bd79f9110e914d814b Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Mon, 26 Jun 2023 15:26:22 +0200
Subject: [PATCH 2/9] analysis_CE : bug fixes

---
 src/utils/analysis_common_envelope.f90 | 11 ++++++-----
 src/utils/moddump_sink.f90             |  2 +-
 src/utils/phantomanalysis.f90          |  3 ++-
 3 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90
index d76271588..a8822cfbb 100644
--- a/src/utils/analysis_common_envelope.f90
+++ b/src/utils/analysis_common_envelope.f90
@@ -1381,8 +1381,8 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
  integer                      :: i,k,Nquantities,ierr,iu
  integer, save                :: quantities_to_calculate(4)
  integer, allocatable         :: iorder(:)
- real                         :: ekini,einti,epoti,ethi,phii,rhopart,ponrhoi,spsoundi,tempi,&
-                                 omega_orb,kappai,kappat,kappar,pgas,mu,entropyi,&
+ real                         :: ekini,einti,epoti,ethi,phii,rhopart,rho_cgs,ponrhoi,spsoundi,tempi,&
+                                 omega_orb,kappai,kappat,kappar,pgas,mu,entropyi,rhopart,&
                                  dum1,dum2,dum3,dum4,dum5
  real, allocatable, save      :: init_entropy(:)
  real, allocatable            :: quant(:,:)
@@ -1461,20 +1461,21 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
        select case (quantities_to_calculate(k))
        case(13) !to calculate JstarS
           rhopart = rhoh(xyzh(4,i), particlemass)
+          rho_cgs = rhopart*unit_density
           !call equationofstate to obtain temperature and store it in tempi
           call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i))
           JstarS = 0.
           !nH_tot is needed to normalize JstarS
-          nH_tot = rhopart/mass_per_H
+          nH_tot = rho_cgs/mass_per_H
           epsC   = eps(3) - nucleation(idK3,i)
           if (epsC < 0.) then
-             print *,'eps(C) =',eps(3),', K3=',nucleation(idK3,i),', epsC=',epsC,', T=',tempi,' rho=',rhopart
+             print *,'eps(C) =',eps(3),', K3=',nucleation(idK3,i),', epsC=',epsC,', T=',tempi,' rho=',rho_cgs
              print *,'JKmuS=',nucleation(:,i)
              stop '[S-dust_formation] epsC < 0!'
           endif
           if (tempi > 450.) then
              !call chemical_equilibrium_light to obtain pC, and pC2H2
-             call chemical_equilibrium_light(rhopart, tempi, epsC, pC, pC2, pC2H, pC2H2, nucleation(idmu,i), nucleation(idgamma,i))
+             call chemical_equilibrium_light(rho_cgs, tempi, epsC, pC, pC2, pC2H, pC2H2, nucleation(idmu,i), nucleation(idgamma,i))
              S = pC/psat_C(tempi)
              if (S > Scrit) then
                 !call nucleation_function to obtain JstarS
diff --git a/src/utils/moddump_sink.f90 b/src/utils/moddump_sink.f90
index 880f51081..4cffdabea 100644
--- a/src/utils/moddump_sink.f90
+++ b/src/utils/moddump_sink.f90
@@ -77,7 +77,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
     xyzmh_ptmass(1,isinkpart) = newx
     print*,'x-coordinate changed to ',xyzmh_ptmass(1,isinkpart)
 
-    Lnuc = xyzmh_ptmass(1,ilum)
+    Lnuc = xyzmh_ptmass(ilum,isinkpart)
     Lnuc_cgs = Lnuc * unit_energ / utime
     call prompt('Enter new sink heating luminosity in erg/s:',Lnuc_cgs,0.)
     xyzmh_ptmass(ilum,isinkpart) = Lnuc_cgs / unit_energ * utime
diff --git a/src/utils/phantomanalysis.f90 b/src/utils/phantomanalysis.f90
index b9ac76e5e..c758f103e 100644
--- a/src/utils/phantomanalysis.f90
+++ b/src/utils/phantomanalysis.f90
@@ -58,7 +58,8 @@ program phantomanalysis
 !
     if (iarg==1) then
 
-       iloc = index(dumpfile,'_0')
+       !iloc = index(dumpfile,'_0')
+       iloc = index(dumpfile,'_',.true.) !to load dump > 9999
 
        if (iloc > 1) then
           fileprefix = trim(dumpfile(1:iloc-1))

From 8e79433331db97cae1c05ddf6773b8a46d64368b Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Mon, 26 Jun 2023 18:53:53 +0200
Subject: [PATCH 3/9] Update CE.rst

---
 docs/CE.rst | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/docs/CE.rst b/docs/CE.rst
index a3be7f23d..ce86e2e71 100644
--- a/docs/CE.rst
+++ b/docs/CE.rst
@@ -124,3 +124,12 @@ if you come from 2.10, then use as initial model (hereafter initial_nnnnn) one o
         softening length for the primary core = 1., softening length for companion = 0.1)
   2.16 vim binary.in (optional, tmax=200.00, dtmax=0.100)
   2.17 ./phantom binary.in
+
+
+**D. Setup sink properties (luminosity)**
+
+::
+
+  2.18 ./phantommoddump binary_00000.tmp dusty_binary_00000.tmp 0.0
+       option 9, 12 lum  
+  2.19 ./phantom dusty_binary.in

From cb708f8ac4bac54f481d2f8b2975da3b48899908 Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Wed, 28 Jun 2023 15:11:57 +0200
Subject: [PATCH 4/9] (set_star) fix critical typo

---
 src/setup/set_star.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90
index 7a617d3c3..ce9460408 100644
--- a/src/setup/set_star.f90
+++ b/src/setup/set_star.f90
@@ -697,7 +697,7 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
     else
        star%isinkcore = .true.
        call read_inopt(star%input_profile,'input_profile'//trim(c),db,errcount=nerr)
-       call read_inopt(star%outputfilename,'outputfilename//trim(c)',db,errcount=nerr)
+       call read_inopt(star%outputfilename,'outputfilename'//trim(c),db,errcount=nerr)
        if (star%isoftcore==1) call read_inopt(star%isofteningopt,'isofteningopt'//trim(c),&
                                               db,errcount=nerr,min=0)
        if ((star%isofteningopt==1) .or. (star%isofteningopt==3)) then

From f75105e270bfdeace6ecb67d6c0ad4c1c920bee5 Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Wed, 28 Jun 2023 16:06:30 +0200
Subject: [PATCH 5/9] (checksetup.F90) if isink_radiation > 1 check that the
 sink particle has a defined luminosity

---
 src/main/checksetup.F90 | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/main/checksetup.F90 b/src/main/checksetup.F90
index fb8b187a7..2e753ebb5 100644
--- a/src/main/checksetup.F90
+++ b/src/main/checksetup.F90
@@ -512,7 +512,8 @@ end function in_range
 !------------------------------------------------------------------
 subroutine check_setup_ptmass(nerror,nwarn,hmin)
  use dim,  only:maxptmass
- use part, only:nptmass,xyzmh_ptmass,ihacc,ihsoft,gr,iTeff,sinks_have_luminosity
+ use part, only:nptmass,xyzmh_ptmass,ihacc,ihsoft,gr,iTeff,sinks_have_luminosity,ilum
+ use ptmass_radiation, only:isink_radiation
  integer, intent(inout) :: nerror,nwarn
  real,    intent(in)    :: hmin
  integer :: i,j,n
@@ -589,6 +590,11 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
  !
  !  check that radiation properties are sensible
  !
+ if (isink_radiation > 1 .and. xyzmh_ptmass(ilum,1) < 1e-10) then
+    nerror = nerror + 1
+    print*,'ERROR: isink_radiation > 1 and sink particle has no luminosity'
+    return
+ endif
  if (sinks_have_luminosity(nptmass,xyzmh_ptmass)) then
     if (any(xyzmh_ptmass(iTeff,1:nptmass) < 100.)) then
        print*,'WARNING: sink particle temperature less than 100K'

From b603ba36d62e29e9b2ee1c3931a6bbc5eb908485 Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Wed, 28 Jun 2023 16:15:31 +0200
Subject: [PATCH 6/9] Update CE.rst

---
 docs/CE.rst | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/CE.rst b/docs/CE.rst
index ce86e2e71..382a7b017 100644
--- a/docs/CE.rst
+++ b/docs/CE.rst
@@ -73,7 +73,7 @@ Use SETUP=star or SETUP=dustystar and if not specified, the default options.
   2.2 make setup
   2.3 ./phantomsetup star (option 5 MESA star, input profile = Jan_Star_Phantom_Profile.data, desired EOS = 10,
       use constant entropy profile, Relax star automatically = yes). The core radius is the softening radius (2-3Ro)
-      the core mass is the same as the one you have measured from MESA (0.46Mo in Jan_Star_Phantom_Profile).
+      the core mass is the same as the one you have measured from MESA (0.46Mo in Jan_Star_Phantom_Profile.data).
       This produces a file called star.setup - this file has all the options so you can edit it.
   2.4 vim  star.setup, (write_rho_to_file = T)
 Relaxation
@@ -132,4 +132,5 @@ if you come from 2.10, then use as initial model (hereafter initial_nnnnn) one o
 
   2.18 ./phantommoddump binary_00000.tmp dusty_binary_00000.tmp 0.0
        option 9, 12 lum  
-  2.19 ./phantom dusty_binary.in
+  2.19 vim dusty_binary.in (adapt isink_radiation, idust_opacity)
+  2.20 ./phantom dusty_binary.in

From f35fdbb8ed07df04970382c3add2e1c6f7ac54b7 Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Wed, 28 Jun 2023 18:49:28 +0200
Subject: [PATCH 7/9] ask for sink luminosity in set_star + add lcore to the
 star_t type

---
 src/main/inject_wind.F90            |  1 -
 src/main/readwrite_dumps_common.F90 | 12 ++++-----
 src/main/units.f90                  | 14 ++++++++++-
 src/setup/set_star.f90              | 39 +++++++++++++++++++----------
 src/setup/set_star_utils.f90        | 13 +++++++---
 5 files changed, 55 insertions(+), 24 deletions(-)

diff --git a/src/main/inject_wind.F90 b/src/main/inject_wind.F90
index 5fbad9ed1..db149ff2d 100644
--- a/src/main/inject_wind.F90
+++ b/src/main/inject_wind.F90
@@ -28,7 +28,6 @@ module inject
 !   infile_utils, injectutils, io, options, part, partinject, physcon,
 !   ptmass_radiation, setbinary, timestep, units, wind, wind_equations
 !
- use physcon,           only:solarl
  use dim,               only:isothermal
 
  implicit none
diff --git a/src/main/readwrite_dumps_common.F90 b/src/main/readwrite_dumps_common.F90
index 63bed518e..38a8ca86c 100644
--- a/src/main/readwrite_dumps_common.F90
+++ b/src/main/readwrite_dumps_common.F90
@@ -125,7 +125,7 @@ subroutine check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkpropert
  use dim,  only:maxp,maxvxyzu,maxalpha,maxBevol,mhd,h2chemistry,use_dustgrowth,gr,&
                 do_radiation,store_dust_temperature,do_nucleation
  use eos,  only:ieos,polyk,gamma,eos_is_non_ideal
- use part, only:maxphase,isetphase,set_particle_type,igas,ihacc,ihsoft,imacc,&
+ use part, only:maxphase,isetphase,set_particle_type,igas,ihacc,ihsoft,imacc,ilum,&
                 xyzmh_ptmass_label,vxyz_ptmass_label,get_pmass,rhoh,dustfrac,ndusttypes,norig
  use io,   only:warning,id,master
  use options,        only:alpha,use_dustfrac,use_var_comp
@@ -322,13 +322,13 @@ subroutine check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkpropert
     endif
     if (id==master .and. i1==1) then
        print "(2(a,i4),a)",' got ',nsinkproperties,' sink properties from ',nptmass,' sink particles'
-       if (nptmass > 0) print "(1x,47('-'),/,1x,a,'|',4(a9,1x,'|'),/,1x,47('-'))",&
-                              'ID',' Mass    ',' Racc    ',' Macc    ',' hsoft   '
+       if (nptmass > 0) print "(1x,58('-'),/,1x,a,'|',5(a9,1x,'|'),/,1x,58('-'))",&
+                              'ID',' Mass    ',' Racc    ',' Macc    ',' hsoft   ',' Lsink   '
        do i=1,min(nptmass,999)
-          if (xyzmh_ptmass(4,i) > 0.) print "(i3,'|',4(1pg9.2,1x,'|'))",i,xyzmh_ptmass(4,i), &
-                                            xyzmh_ptmass(ihacc,i),xyzmh_ptmass(imacc,i),xyzmh_ptmass(ihsoft,i)
+          if (xyzmh_ptmass(4,i) > 0.) print "(i3,'|',5(1pg9.2,1x,'|'))",i,xyzmh_ptmass(4,i),xyzmh_ptmass(ihacc,i),&
+                                            xyzmh_ptmass(imacc,i),xyzmh_ptmass(ihsoft,i),xyzmh_ptmass(ilum,i)
        enddo
-       if (nptmass > 0) print "(1x,47('-'))"
+       if (nptmass > 0) print "(1x,58('-'))"
     endif
  endif
  !
diff --git a/src/main/units.f90 b/src/main/units.f90
index 14a8cd605..455ee79af 100644
--- a/src/main/units.f90
+++ b/src/main/units.f90
@@ -35,7 +35,7 @@ module units
  public :: get_G_code, get_c_code, get_radconst_code, get_kbmh_code
  public :: c_is_unity, G_is_unity, in_geometric_units
  public :: is_time_unit, is_length_unit
- public :: in_solarr, in_solarm
+ public :: in_solarr, in_solarm, in_solarl
 
 contains
 
@@ -464,5 +464,17 @@ real(kind=8) function in_solarr(val) result(rval)
  rval = val*(udist/solarr)
 
 end function in_solarr
+!---------------------------------------------------------------------------
+!+
+!  function to convert a luminosity value from code units to solar luminosity
+!+
+!---------------------------------------------------------------------------
+real(kind=8) function in_solarl(val) result(rval)
+ use physcon, only:solarl
+ real, intent(in) :: val
+
+ rval = val*(unit_luminosity/solarl)
+
+end function in_solarl
 
 end module units
diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90
index ce9460408..26d5ef5b6 100644
--- a/src/setup/set_star.f90
+++ b/src/setup/set_star.f90
@@ -41,6 +41,7 @@ module setstar
     real :: initialtemp
     real :: rcore
     real :: mcore
+    real :: lcore
     real :: hsoft
     real :: hacc   ! accretion radius if star is a sink particle
     character(len=120) :: input_profile,dens_profile
@@ -80,6 +81,7 @@ subroutine set_defaults_star(star)
  star%hacc           = 1.
  star%rcore          = 0.
  star%mcore          = 0.
+ star%lcore          = 0.
  star%isofteningopt  = 1 ! By default, specify rcore
  star%np             = 1000
  star%input_profile  = 'P12_Phantom_Profile.data'
@@ -108,7 +110,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,&
                               write_kepler_comp
  use radiation_utils,    only:set_radiation_and_gas_temperature_equal
  use relaxstar,          only:relax_star
- use part,               only:ihsoft,igas,imu,set_particle_type
+ use part,               only:ihsoft,igas,imu,set_particle_type,ilum
  use extern_densprofile, only:write_rhotab
  use unifdis,            only:mask_prototype
  use physcon,            only:pi
@@ -185,10 +187,11 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,&
  !
  ! add sink particle stellar core
  !
- if (star%isinkcore) call set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,&
-                                           ihsoft,star%mcore,star%hsoft,ierr)
+ if (star%isinkcore) call set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,&
+                                           star%mcore,star%hsoft,ilum,star%lcore,ierr)
  if (ierr==1) call fatal('set_stellar_core','mcore <= 0')
  if (ierr==2) call fatal('set_stellar_core','hsoft <= 0')
+ if (ierr==3) call fatal('set_stellar_core','lcore < 0')
  !
  ! Write the desired profile to file (do this before relaxation)
  !
@@ -440,8 +443,8 @@ end subroutine set_defaults_given_profile
 subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
  use prompting,     only:prompt
  use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar
- use units,         only:in_solarm,in_solarr,udist,umass
- use physcon,       only:solarr,solarm
+ use units,         only:in_solarm,in_solarr,in_solarl,udist,umass,unit_luminosity
+ use physcon,       only:solarr,solarm,solarl
  type(star_t), intent(out)   :: star
  integer,      intent(in)    :: id,master
  logical,      intent(out)   :: use_var_comp
@@ -449,13 +452,14 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
  integer,      intent(inout) :: ieos
  real,         intent(inout) :: polyk
  integer :: i
- real :: mstar_msun,rstar_rsun,rcore_rsun,mcore_msun,hsoft_rsun
+ real :: mstar_msun,rstar_rsun,rcore_rsun,mcore_msun,lcore_lsun,hsoft_rsun
 
  ! set defaults
  call set_defaults_star(star)
  mstar_msun = real(in_solarm(star%mstar))
  rstar_rsun = real(in_solarr(star%rstar))
  mcore_msun = real(in_solarm(star%mcore))
+ lcore_lsun = real(in_solarl(star%lcore))
  rcore_rsun = real(in_solarr(star%rcore))
  hsoft_rsun = real(in_solarr(star%hsoft))
 
@@ -506,8 +510,10 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
        if (star%isinkcore) then
           call prompt('Enter mass of the created sink particle core [Msun]',mcore_msun,0.)
           call prompt('Enter softening length of the sink particle core [Rsun]',hsoft_rsun,0.)
+          call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.)
           star%mcore = mcore_msun*real(solarm/umass)
           star%hsoft = hsoft_rsun*real(solarr/udist)
+          star%lcore = lcore_lsun*real(solarl/unit_luminosity)
        endif
     case(1)
        star%isinkcore = .true. ! Create sink particle core automatically
@@ -532,19 +538,20 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
           star%mcore = mcore_msun*real(solarm/umass)
           star%rcore = rcore_rsun*real(solarr/udist)
        end select
-
-       call prompt('Enter output file name of cored stellar profile:',star%outputfilename)
+       call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.)
+       star%lcore = lcore_lsun*real(solarl/unit_luminosity)
 
     case(2)
        star%isinkcore = .true. ! Create sink particle core automatically
        print*,'Specify core radius and initial guess for mass of sink particle core'
        call prompt('Enter core radius in Rsun : ',rcore_rsun,0.)
        call prompt('Enter guess for core mass in Msun : ',mcore_msun,0.)
+       call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.)
        call prompt('Enter output file name of cored stellar profile:',star%outputfilename)
        star%mcore = mcore_msun*real(solarm/umass)
        star%rcore = rcore_rsun*real(solarr/udist)
+       star%lcore = lcore_lsun*real(solarl/unit_luminosity)
     end select
-
  case(ievrard)
     call prompt('Enter the specific internal energy (units of GM/R) ',star%ui_coef,0.)
  case(:0)
@@ -561,7 +568,7 @@ end subroutine set_star_interactive
 subroutine write_options_star(star,iunit,label)
  use infile_utils,  only:write_inopt,get_optstring
  use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar
- use units,         only:in_solarm,in_solarr
+ use units,         only:in_solarm,in_solarr,in_solarl
  type(star_t),     intent(in) :: star
  integer,          intent(in) :: iunit
  character(len=*), intent(in), optional :: label
@@ -614,6 +621,8 @@ subroutine write_options_star(star,iunit,label)
           call write_inopt(in_solarm(star%mcore),'mcore'//trim(c),&
                'Initial guess for mass of sink particle stellar core [Msun]',iunit)
        endif
+       call write_inopt(in_solarl(star%lcore),'lcore'//trim(c),&
+                              'Luminosity of point mass stellar core [Lsun]',iunit)
     else
        call write_inopt(star%isinkcore,'isinkcore'//trim(c),&
                'Add a sink particle stellar core',iunit)
@@ -623,6 +632,8 @@ subroutine write_options_star(star,iunit,label)
           call write_inopt(in_solarr(star%hsoft),'hsoft'//trim(c),&
                'Softening length of sink particle stellar core [Rsun]',iunit)
        endif
+       call write_inopt(in_solarl(star%lcore),'lcore'//trim(c),&
+               'Luminosity of sink core particle [Lsun]',iunit)
     endif
  case (ievrard)
     call write_inopt(star%ui_coef,'ui_coef'//trim(c),&
@@ -646,8 +657,8 @@ end subroutine write_options_star
 subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
  use infile_utils,  only:inopts,read_inopt
  use setstar_utils, only:need_inputprofile,need_rstar,nprofile_opts
- use units,         only:umass,udist
- use physcon,       only:solarm,solarr
+ use units,         only:umass,udist,unit_luminosity
+ use physcon,       only:solarm,solarr,solarl
  type(star_t),              intent(out)   :: star
  type(inopts), allocatable, intent(inout) :: db(:)
  integer,                   intent(out)   :: need_iso
@@ -656,7 +667,7 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
  integer,                   intent(inout) :: nerr
  character(len=*),          intent(in), optional :: label
  character(len=10) :: c
- real :: mcore_msun,rcore_rsun,mstar_msun,rstar_rsun,hsoft_rsun
+ real :: mcore_msun,rcore_rsun,lcore_lsun,mstar_msun,rstar_rsun,hsoft_rsun
 
  ! set defaults
  call set_defaults_star(star)
@@ -709,6 +720,8 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
           call read_inopt(mcore_msun,'mcore'//trim(c),db,errcount=nerr,min=0.)
           star%mcore = mcore_msun*real(solarm/umass)
        endif
+       call read_inopt(lcore_lsun,'lcore'//trim(c),db,errcount=nerr,min=0.)
+       star%lcore = lcore_lsun*real(solarl/unit_luminosity)
     endif
  case(ievrard)
     call read_inopt(star%ui_coef,'ui_coef'//trim(c),db,errcount=nerr,min=0.)
diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90
index a07d903c0..402cea0d7 100644
--- a/src/setup/set_star_utils.f90
+++ b/src/setup/set_star_utils.f90
@@ -310,11 +310,13 @@ end subroutine set_star_density
 !  Add a sink particle as a stellar core
 !+
 !-----------------------------------------------------------------------
-subroutine set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,mcore,hsoft,ierr)
+subroutine set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,mcore,&
+                            hsoft,ilum,lcore,ierr)
  integer, intent(out) :: nptmass,ierr
  real, intent(out)    :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
- real, intent(in)     :: mcore,hsoft
- integer              :: n,ihsoft
+ real, intent(in)     :: mcore,hsoft,lcore
+ integer, intent(in)  :: ihsoft,ilum
+ integer              :: n
 
  ierr = 0
  ! Check for sensible values
@@ -326,12 +328,17 @@ subroutine set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,mcore,hsoft,
     ierr = 2
     return
  endif
+ if (lcore < 0.) then
+    ierr = 3
+    return
+ endif
 
  nptmass                = 1
  n                      = nptmass
  xyzmh_ptmass(:,n)      = 0. ! zero all quantities by default
  xyzmh_ptmass(4,n)      = mcore
  xyzmh_ptmass(ihsoft,n) = hsoft
+ xyzmh_ptmass(ilum,n)   = lcore
  vxyz_ptmass(:,n)       = 0.
 
 end subroutine set_stellar_core

From c96029e584b75132c3a987f0064d405469ae4019 Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Thu, 29 Jun 2023 10:09:27 +0200
Subject: [PATCH 8/9] (analysis_common_envelope) fix multiple assignments of
 the same variable

---
 src/utils/analysis_common_envelope.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90
index a8822cfbb..37e762f03 100644
--- a/src/utils/analysis_common_envelope.f90
+++ b/src/utils/analysis_common_envelope.f90
@@ -1381,7 +1381,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
  integer                      :: i,k,Nquantities,ierr,iu
  integer, save                :: quantities_to_calculate(4)
  integer, allocatable         :: iorder(:)
- real                         :: ekini,einti,epoti,ethi,phii,rhopart,rho_cgs,ponrhoi,spsoundi,tempi,&
+ real                         :: ekini,einti,epoti,ethi,phii,rho_cgs,ponrhoi,spsoundi,tempi,&
                                  omega_orb,kappai,kappat,kappar,pgas,mu,entropyi,rhopart,&
                                  dum1,dum2,dum3,dum4,dum5
  real, allocatable, save      :: init_entropy(:)

From 315650b92915d7958b87d7f2588a87ae9fa3af1f Mon Sep 17 00:00:00 2001
From: Lionel Siess <lionel.siess@astro.ulb.ac.be>
Date: Thu, 29 Jun 2023 15:23:24 +0200
Subject: [PATCH 9/9] (analysis_CE) fix tabulation

---
 src/utils/analysis_common_envelope.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90
index 37e762f03..979d7b09e 100644
--- a/src/utils/analysis_common_envelope.f90
+++ b/src/utils/analysis_common_envelope.f90
@@ -1497,7 +1497,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
              print *,'nucleation(idgamma,i) = ',nucleation(idgamma,i)
              print *,'taustar = ',taustar
              print *,'eps = ',eps
-	     print *,'JstarS = ',JstarS
+             print *,'JstarS = ',JstarS
           endif
           quant(k,i) = JstarS