Edits to enable blob setup to build

alisonkyoung1 · Oct 21, 2024 · 55eff86 · 55eff86
1 parent 9e14e04
commit 55eff86
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 27 deletions.
diff --git a/src/main/cooling_radapprox.f90 b/src/main/cooling_radapprox.f90
@@ -77,7 +77,7 @@ subroutine radcool_update_energ(i,xi,yi,zi,rhoi,ui,Tfloor,dt,dudti_cool)
  real,intent(out)::dudti_cool
  real            :: coldensi,kappaBari,kappaParti,ri2
  real            :: gmwi,Tmini4,Ti,dudti_rad,Teqi,Hstam,HLom,du_tot
- real            :: cs2,Om2,Hmod2
+ real            :: cs2,Om2,Hmod2,rhoi_cgs,ui_cgs
  real            :: opaci,ueqi,umini,tthermi,poti,presi,Hcomb,du_FLDi
 
  coldensi = huge(coldensi)
@@ -99,7 +99,9 @@ subroutine radcool_update_energ(i,xi,yi,zi,rhoi,ui,Tfloor,dt,dudti_cool)
  endif
 
  ! get opacities & Ti for ui
- call getopac_opdep(ui*unit_ergg,rhoi*unit_density,kappaBari,kappaParti,&
+ ui_cgs = ui*unit_ergg
+ rhoi_cgs = rhoi*unit_density
+ call getopac_opdep(ui_cgs,rhoi_cgs,kappaBari,kappaParti,&
            Ti,gmwi)
  presi = kb_on_mh*rhoi*unit_density*Ti/gmwi ! cgs
  presi = presi/unit_pressure !code units
@@ -146,7 +148,7 @@ subroutine radcool_update_energ(i,xi,yi,zi,rhoi,ui,Tfloor,dt,dudti_cool)
     Tmini4 = Tfloor**4d0
  endif
 
- call getintenerg_opdep(Tmini4**(1.0/4.0),rhoi*unit_density,umini)
+ call getintenerg_opdep(Tmini4**(1.0/4.0),rhoi_cgs,umini)
  umini = umini/unit_ergg
 
  opaci = (coldensi**2d0)*kappaBari + (1.d0/kappaParti) ! physical units
@@ -195,7 +197,8 @@ subroutine radcool_update_energ(i,xi,yi,zi,rhoi,ui,Tfloor,dt,dudti_cool)
             "Ti=", Ti, "poti=",poti, "rhoi=", rhoi
  endif
 
- call getintenerg_opdep(Teqi,rhoi*unit_density,ueqi)
+ rhoi_cgs = rhoi*unit_density
+ call getintenerg_opdep(Teqi,rhoi_cgs,ueqi)
  ueqi = ueqi/unit_ergg
 
  ! calculate thermalization timescale
@@ -248,7 +251,7 @@ subroutine radcool_update_energ_loop(dtsph,npart,xyzh,energ,dudt_sph,Tfloor)
  real,intent(inout) :: energ(:),dudt_sph(:)
  real            :: ui,rhoi,coldensi,kappaBari,kappaParti,ri2,dti
  real            :: gmwi,Tmini4,Ti,dudti_rad,Teqi,Hstam,HLom,du_tot
- real            :: cs2,Om2,Hmod2
+ real            :: cs2,Om2,Hmod2,ui_cgs,rhoi_cgs
  real            :: opaci,ueqi,umini,tthermi,poti,presi,Hcomb,du_FLDi
  integer         :: i,ratefile,n_uevo
 
@@ -262,8 +265,8 @@ subroutine radcool_update_energ_loop(dtsph,npart,xyzh,energ,dudt_sph,Tfloor)
  !$omp shared(isink_star,doFLD,ttherm_store,teqi_store,od_method,unit_pressure,ratefile) &
  !$omp shared(opac_store,Tfloor,dtsph,dudt_sph,utime,udist,umass,unit_ergg,gradP_cool,Lstar) &
  !$omp private(i,poti,du_FLDi,ui,rhoi,ri2,coldensi,kappaBari,Ti,iphase) &
- !$omp private(kappaParti,gmwi,Tmini4,dudti_rad,Teqi,Hstam,HLom,du_tot) &
- !$omp private(cs2,Om2,Hmod2,opaci,ueqi,umini,tthermi,presi,Hcomb,dti) &
+ !$omp private(kappaParti,gmwi,Tmini4,dudti_rad,Teqi,Hstam,HLom,du_tot,ui_cgs) &
+ !$omp private(cs2,Om2,Hmod2,opaci,ueqi,umini,tthermi,presi,Hcomb,dti,rhoi_cgs) &
  !$omp shared(maxp,maxphase,ibin) reduction(+:n_uevo)
 
  overpart: do i=1,npart
@@ -292,7 +295,9 @@ subroutine radcool_update_energ_loop(dtsph,npart,xyzh,energ,dudt_sph,Tfloor)
     endif
 
     ! get opacities & Ti for ui
-    call getopac_opdep(ui*unit_ergg,rhoi*unit_density,kappaBari,kappaParti,&
+    ui_cgs = ui*unit_ergg
+    rhoi_cgs = rhoi*unit_density
+    call getopac_opdep(ui_cgs,rhoi_cgs,kappaBari,kappaParti,&
            Ti,gmwi)
     presi = kb_on_mh*rhoi*unit_density*Ti/gmwi ! cgs
     presi = presi/unit_pressure !code units
@@ -377,10 +382,11 @@ subroutine radcool_update_energ_loop(dtsph,npart,xyzh,energ,dudt_sph,Tfloor)
             "Ti=", Ti, "poti=",poti, "rhoi=", rhoi
     endif
 
-    call getintenerg_opdep(Teqi,rhoi*unit_density,ueqi)
+	rhoi_cgs = rhoi*unit_density
+    call getintenerg_opdep(Teqi,rhoi_cgs,ueqi)
     ueqi = ueqi/unit_ergg
 
-    call getintenerg_opdep(Tmini4**(1.0/4.0),rhoi*unit_density,umini)
+    call getintenerg_opdep(Tmini4**(1.0/4.0),rhoi_cgs,umini)
     umini = umini/unit_ergg
 
     ! calculate thermalization timescale

diff --git a/src/main/dens.F90 b/src/main/dens.F90
@@ -302,6 +302,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
 !$omp reduction(max:rhomax) &
 !$omp private(i)
 
+
  call init_cell_exchange(xrecvbuf,irequestrecv,thread_complete,ncomplete_mpi,mpitype)
 
  !$omp master
@@ -1689,14 +1690,14 @@ subroutine calc_lambda_cell(cell,listneigh,nneigh,xyzh,xyzcache,vxyzu,iphase,gra
  integer      :: icell,i,iamtypei,iamtypej,j,n
  logical      :: iactivei,iamgasi,iamdusti,ignoreself
  logical      :: iactivej,iamgasj,iamdustj
- real(kind=8) :: hi,hi1,hi21,hi31,hi41
- real         :: rhoi,rho1i,dhdrhoi,pmassi,kappabari,kappaparti,Ti,gmwi
+ real(kind=8) :: hi1,hi21,hi31,hi41,hj1,hj21
+ real         :: hi,hj,rhoi,rho1i,dhdrhoi,pmassi,kappabari,kappaparti,Ti,gmwi
  real         :: xj,yj,zj,dx,dy,dz
- real         :: rij2,rij,q2i,qi,hj1,hj,hj21,q2j
+ real         :: rij2,rij,q2i,qi,q2j
  real         :: wabi,grkerni,gradhi,wkerni,dwkerni
  real         :: pmassj,rhoj,rho1j,dhdrhoj,kappabarj,kappaPartj,Tj,gmwj
  real         :: uradi,dradi,dradxi,dradyi,dradzi,runix,runiy,runiz
- real         :: dT4,R_rad
+ real         :: dT4,R_rad,u_cgs,rho_cgs
  integer      :: ngradh_err
 
  ngradh_err = 0
@@ -1735,8 +1736,9 @@ subroutine calc_lambda_cell(cell,listneigh,nneigh,xyzh,xyzcache,vxyzu,iphase,gra
        print *, "u=0 in FLD calc", vxyzu(4,i), i,rhoi*unit_density,Ti,&
             cell%xpartvec(ixi,icell),cell%xpartvec(iyi,icell)
     endif
-    call getopac_opdep(vxyzu(4,i)*unit_ergg,rhoi*unit_density,kappabari, &
-     kappaparti,Ti,gmwi)
+    u_cgs = vxyzu(4,i)*unit_ergg
+    rho_cgs = rhoi*unit_density
+    call getopac_opdep(u_cgs,rho_cgs,kappabari,kappaparti,Ti,gmwi)
 
     loop_over_neighbours: do n=1,nneigh
        j = abs(listneigh(n))
@@ -1791,7 +1793,9 @@ subroutine calc_lambda_cell(cell,listneigh,nneigh,xyzh,xyzcache,vxyzu,iphase,gra
           dwkerni = grkerni*cnormk*hi21*hi21*gradh(1,i)
           pmassj = massoftype(iamtypej)
           call rhoanddhdrho(hj,hj1,rhoj,rho1j,dhdrhoj,pmassj)
-          call getopac_opdep(vxyzu(4,j)*unit_ergg,rhoj*unit_density,&
+          u_cgs = vxyzu(4,j)*unit_ergg
+          rho_cgs = rhoj*unit_density
+          call getopac_opdep(u_cgs,rho_cgs,&
                  kappaBarj,kappaPartj,Tj,gmwj)
           uradi = uradi + get_radconst_code()*(Tj**4.0d0)*wkerni*pmassj/rhoj
 

diff --git a/src/main/eos_stamatellos.f90 b/src/main/eos_stamatellos.f90
@@ -283,13 +283,15 @@ subroutine get_k_fld(rhoi,eni,i,ki,Ti)
  use units,    only:unit_density,unit_ergg,unit_opacity,get_radconst_code
  real,intent(in)    :: rhoi,eni
  integer,intent(in) :: i
- real               :: kappaBar,gmwi,kappaPart
+ real               :: kappaBar,gmwi,kappaPart,eni_ergg,rhoi_g
  real,intent(out)   :: ki,Ti
 
  if (lambda_FLD(i) == 0d0) then
     ki = 0.
  else
-    call getopac_opdep(eni*unit_ergg,rhoi*unit_density,kappaBar,kappaPart,Ti,gmwi)
+ 	eni_ergg = eni*unit_ergg
+ 	rhoi_g = rhoi*unit_density
+    call getopac_opdep(eni_ergg,rhoi_g,kappaBar,kappaPart,Ti,gmwi)
     kappaPart = kappaPart/unit_opacity
     ! steboltz constant = 4pi/c * arad
     ki = 16d0*(fourpi/c)*get_radconst_code()*lambda_FLD(i)*Ti**3 /rhoi/kappaPart

diff --git a/src/main/radiation_utils.f90 b/src/main/radiation_utils.f90
@@ -419,7 +419,7 @@ subroutine get_opacity(opacity_type,density,temperature,kappa,u)
  real, intent(in), optional :: u
  real, intent(out) :: kappa
  integer, intent(in) :: opacity_type
- real :: kapt,kapr,rho_cgs,Ti,gmwi,kapBar,kappaPart
+ real :: kapt,kapr,rho_cgs,Ti,gmwi,kapBar,kappaPart,u_cgs
 
  select case(opacity_type)
  case(1)
@@ -440,7 +440,9 @@ subroutine get_opacity(opacity_type,density,temperature,kappa,u)
     !
     ! opacity for Stamatellos/Lombardi EOS
     !
-    call getopac_opdep(u*unit_ergg,density*unit_density,kapBar,kappaPart,Ti,gmwi)
+    rho_cgs = density*unit_density
+    u_cgs = u*unit_ergg
+    call getopac_opdep(u_cgs,rho_cgs,kapBar,kappaPart,Ti,gmwi)
     kappa = kappaPart/unit_opacity
  case default
     !

diff --git a/src/utils/analysis_disc_stresses.f90 b/src/utils/analysis_disc_stresses.f90
@@ -35,7 +35,7 @@ module analysis
  real,    allocatable,dimension(:)   :: H, toomre_q,epicyc,part_scaleheight
  real,    allocatable,dimension(:)   :: alpha_reyn,alpha_grav,alpha_mag,alpha_art
  real,    allocatable,dimension(:)   :: rpart,phipart,vrpart,vphipart, gr,gphi,Br,Bphi
- real,    allocatable,dimension(:,:) :: gravxyz
+ real,    allocatable,dimension(:,:) :: gravxyz,zsetgas
 
  logical :: write_neighbour_list = .true.  ! Write the neighbour list to file, if true
 
@@ -356,7 +356,7 @@ end subroutine transform_to_cylindrical
 !+
 !---------------------------------------------------------------
 
-subroutine radial_binning(npart,xyzh,vxyzu,pmass)
+subroutine radial_binning(npart,xyzh,vxyzu,pmass,eos_vars)
  use physcon, only:pi
  use eos,     only:get_spsound,ieos
  use part,    only:rhoh,isdead_or_accreted
@@ -365,7 +365,7 @@ subroutine radial_binning(npart,xyzh,vxyzu,pmass)
  real,intent(in) :: pmass
  real,intent(in) :: xyzh(:,:),vxyzu(:,:),eos_vars(:,:)
 
- integer :: ibin,ipart,nbinned
+ integer :: ibin,ipart,nbinned,iallocerr
  real :: area,csi
 
  print '(a,I4)', 'Carrying out radial binning, number of bins: ',nbins
@@ -464,8 +464,8 @@ subroutine calc_stresses(npart,xyzh,vxyzu,pmass)
  use physcon, only: pi,gg,kb_on_mh
  use units,   only: print_units, umass,udist,utime,unit_velocity,unit_density,unit_Bfield
  use dim,     only: gravity
- use part,    only: mhd,rhoh,alphaind,eos_vars,imu,itemp
- use eos,     only: gamma,ieos
+ use part,    only: mhd,rhoh,alphaind,imu,itemp
+ use eos,     only: ieos
 
  implicit none
 
@@ -500,7 +500,7 @@ subroutine calc_stresses(npart,xyzh,vxyzu,pmass)
 
  sigma(:) = sigma(:)*umass/(udist*udist)
  if (ieos /= 21) then
-	 csbin(:) = csbin(:)*unit_velocity
+    csbin(:) = csbin(:)*unit_velocity
  endif
  omega(:) = omega(:)/utime