noc.F
``````#define UO   9
#define CVEL 3.063912d5

#define X1C 0.0
#define X2C 0.0
#define X3C 0.0

C #define RANDOM
#define RAND_DEN     0.01
#define RAND_EN      0.01
#define RAND_VEL_KMS 2.00

c=======================================================================
c
c    \\\\\\\\\\      B E G I N   S U B R O U T I N E      //////////
c    //////////               C L U S T E R               \\\\\\\\\\
c
c=======================================================================

subroutine cluster_co09

C     Total mass density is SIS
C     Gas density given by deprojected beta profile
C     Solve for temperature in HSE
C     Temp is isothermal in the outer parts with a peak in the center

use real_prec
use config
use param
use grid
use field
use root
use bndry
use gravmod
#ifdef MPI_USED
use mpiyes
#else
use mpino
#endif
use mpipar
C
use star

implicit none
REAL(rl) ::  d1d(in), e1d(in), gp1d(in)
REAL(rl) ::  sig,rc,d0,beta,rmin,ff_rot
REAL(rl) ::  mstar, reff, rstar
REAL(rl) ::  xl, xh, vol, frac, mshell
INTEGER  ::  i, j, k

#ifdef RANDOM
call srand(9)
#endif /* RANDOM */

.          (sig,rc,d0,beta,rmin,ff_rot,mstar,reff)

cc=====================================================
#if defined(ACCRETION) || defined(AGNFB)
C  BLACK-HOLE MASS IS NEEDED TO DETERMINE THE POINT MASS
C  GRAVITATIONAL POTENIAL. SO IT HAS TO READ "AGN PROFILE"
C  (IF ANY) BEFORE INITIALIZING THE GALACTIC PROFILE
call agn_init
#endif
cc=====================================================

C     ONE-DIMENSIONAL MODEL:
call cluster_1d(d1d, e1d, gp1d, sig, rc, d0, beta)

do k=ks-2,ke+2
do j=js-2,je+2
do i=is-2,ie+2

d (i,j,k) = d1d(i)
v1(i,j,k) = zro
v2(i,j,k) = zro
v3(i,j,k) = zro
e (i,j,k) = e1d(i)

C            Jaffe profile.  reff is half light radius. Compute
C            cell densities by getting the analytic mass in the
C            cell and then dividing by the volume.
rstar = reff/0.7447
xl = x1a(i  )/rstar
xh = x1a(i+1)/rstar
mshell = mstar*(1/(1+xl)-1/(1+xh))
frac = dvl2b(j)*dvl3b(k)/(4*pi)
vol  = dvl1b(i)*dvl2b(j)*dvl3b(k)

dstar(i,j,k) = frac*mshell/vol

#ifdef BACKGROUND_GRAVITY
C            The initial backgroud gravitational field due to
C            dark matter, stellar bulge and the tenius gas.
C            I put it into "force" manually.
gp_bg(i,j,k) = gp1d(i)
#endif /* BKGRND_GRAV */

#ifdef    RANDOM
d(i,j,k) = d(i,j,k) + RAND_DEN*(rand()-0.5)*d(i,j,k)
e(i,j,k) = e(i,j,k) + RAND_EN *(rand()-0.5)*e(i,j,k)
v1(i,j,k) = RAND_VEL_KMS*(rand()-0.5)
v2(i,j,k) = RAND_VEL_KMS*(rand()-0.5)
#endif /* RANDOM */

enddo
enddo
enddo

C Boundary values
do k=ks-2,ke+2
do j=js-2,je+2
do i=1,2
doib(i,j,k)  = d1d(i)
eoib(i,j,k)  = e1d(i)
v1oib(i,j,k) = zro
v2oib(i,j,k) = zro
v3oib(i,j,k) = zro
enddo
enddo
enddo

ccccccccccccccccccccccccccccccccccccccccccccccccc
C INITIALIZE THE "NEW PHYSICS" MODULES

#if defined(ACCRETION) || defined(AGNFB)
C     call agn_init    ! moved foreward
if(rmin.gt.TINY) then
call cluster_rotation(sig,rmin,ff_rot)
endif
#endif /* AGN */

#if defined(STAREV) || defined(STARFORM)
#endif /* STAREV || STARFORM */

#ifdef    STARFORM
call sf_init
#endif /* STARFORM */

#ifdef    ENSURE_GOOD
call ensure_init
#endif /* ENSURE_GOOD */

ccccccccccccccccccccccccccccccccccccccccccccccccc

do k=ks-2,ke+2
do j=js-2,je+2
do i=is-2,ie+2
enddo
enddo
enddo

return
end

c=======================================================================
c    //////////        E N D    OF   C L U S T E R        \\\\\\\\\\
c=======================================================================

C     TO get the "galactic profile" parameters

use real_prec
use param
#ifdef MPI_USED
use mpiyes
#else
use mpino
#endif
use mpipar
C
use star

implicit none
INTEGER num
parameter (num=8)
REAL(rl) ::  buffer(num)
REAL(rl) ::  sig, rc, d0, beta, rmin, ff, mstar, reff
namelist /cluster/  sig,rc,d0,beta,rmin,ff,mstar,reff

if (myid .eq. 0) then
rewind(1)
rewind(1)
write (2,cluster)

buffer(1) = sig
buffer(2) = rc
buffer(3) = d0
buffer(4) = beta
buffer(5) = rmin
buffer(6) = ff
buffer(7) = mstar
buffer(num) = reff
endif

call MPI_BCAST(buffer,num,MPI_FLOAT,0,comm3d,ierr)

sig   = buffer(1)
rc    = buffer(2)
d0    = buffer(3)
beta  = buffer(4)
rmin  = buffer(5)
ff    = buffer(6)
mstar = buffer(7)
reff  = buffer(num)

C cross dependency here: mass loss rate needs to know
C about velocity dispersion and total stellar mass
#if defined(STARFORM) || defined(STAREV)
sigma_sf = sig
mstar_sf = mstar
reff_sf  = reff
#endif /*    STARFORM || STAREV    */

return

c=======================================================================
c    \\\\\\\\\\      B E G I N   S U B R O U T I N E      //////////
c    //////////             C L U S T E R _ 1 D           \\\\\\\\\\
c=======================================================================

subroutine cluster_1d(d,e,gp,sig, rc, d0, beta)
C     1D solution for pressure given mass profile

use real_prec
use config
use param
use grid
use root
use cons
use gravmod
#ifdef MPI_USED
use mpiyes
#else
use mpino
#endif
use mpipar

implicit none
REAL(rl) ::  d(in), e(in), gp(in)
REAL(rl) ::  sig, rc, d0, beta
REAL(rl) ::  p(in), ma, den
REAL(rl) ::  dp, rr, pmin, pmin_local
REAL(rl) ::  plow, phigh
INTEGER  ::  i

C     specify gas densities and potential
do i=is-2,ie+2
rr = x1b(i)
gp(i) = -2*sig**2*log(rr/rc)
d(i)  = d0*(1+rr**2/rc**2)**(-3*beta)
enddo

C     pressure integral
p(is-2) = 0
do i=is-1,ie+2
den = 0.5*(d(i) + d(i+1))
ma = 2*sig**2*x1a(i)/guniv
C        Include BH mass in mass enclosed
if(xptmass) ma = ma + ptmass
dp = - guniv*ma*den*dx1a(i) /x1a(i)**2
p(i) = p(i-1) + dp
enddo

C     Connect integrals in different processors
dp = 0
if (ntiles(1).gt.1) then
if (coords(1).eq.0) then
phigh = p(ie)
call MPI_SEND(phigh, 1, MPI_FLOAT, n1p
.           ,9000, comm3d, ierr)
else
call MPI_RECV(plow , 1, MPI_FLOAT, n1m
.           ,9000, comm3d, stat, ierr)

dp = plow - p(is-1)
phigh = p(ie) + dp

if (coords(1).lt.ntiles(1)-1)
.      call MPI_SEND(phigh, 1, MPI_FLOAT, n1p
.           ,9000, comm3d, ierr)
endif
endif

C     Adjust based on lower processors
do i=is-2,ie+2
p(i) = p(i) + dp
enddo

C     Find local min pressure
pmin_local = huge
do i=is,ie
if (p(i).lt.pmin_local) pmin_local = p(i)
enddo

C     Find the global min pressure
call MPI_REDUCE(pmin_local, pmin, 1, MPI_FLOAT,
.     MPI_MIN, 0, comm3d, ierr)
call MPI_BCAST(pmin,1,MPI_FLOAT,0,comm3d,ierr)

C     Ensure that the pressure is nonnegative
dp = - pmin
do i=is-2,ie+2
p(i) = p(i) + dp
enddo

C     Ensure that last two bins have the same sound speed
if (ntiles(1).gt.1) then
if (coords(1).eq.ntiles(1)-1) then
dp = (d(ie-1)*p(ie) - d(ie)*p(ie-1))
.        / (d(ie) - d(ie-1))
call MPI_SEND(dp, 1, MPI_FLOAT, n1m
.          ,9001, comm3d, ierr)
else
call MPI_RECV(dp, 1, MPI_FLOAT, n1p
.          ,9001, comm3d, stat, ierr)
if (coords(1).gt.0)
.     call MPI_SEND(dp, 1, MPI_FLOAT, n1m
.          ,9001, comm3d, ierr)
endif
endif

do i=is-2,ie+2
p(i) = p(i) + dp
enddo

C     Set energies
do i=is-2,ie+2
e(i) = p(i)/gamm1
enddo

return
end

c=======================================================================
c    //////////      E N D    OF   C L U S T E R _ 1 D      \\\\\\\\\\
c=======================================================================

subroutine mass_enclosed(ma,mb)
C     Calculate mass enclosed within some radius.
C     Fill the values into the passed arrays

use real_prec
use config
use param
use grid
use field
use cons
use gravmod

implicit none
REAL(rl) ::  ma(in), mb(in)
INTEGER  ::  i

C     NOTE -- assumes static, spherically symmetric potential, neglects gas
C     But it does include the central point mass
do i=is,ie+1
#ifdef BACKGROUND_GRAVITY
ma(i) = -(gp_bg(i,js,ks)-gp_bg(i-1,js,ks))
.             /dx1b(i) *  x1a(i)**2/guniv
if(XPTMASS) ma(i) = ma(i) + ptmass
#else
ma(i) = zro
#endif /* BKGRND_GRAV */
enddo

do i=is,ie
C     NOTE -- could come up with a more elaborate scheme here
mb(i) = 0.5*(ma(i+1)+ma(i))
enddo

c      do i=is,ie
c#ifdef BACKGROUND_GRAVITY
c         ma(i) = 2*sig**2*x1a(i)/guniv
c         mb(i) = 2*sig**2*x1b(i)/guniv
c         if(XPTMASS) then
c            ma(i) = ma(i) + ptmass
c            mb(i) = mb(i) + ptmass
c         endif
c#else
c          ma(i) = zro
c          mb(i) = zro
c#endif /* BKGRND_GRAV */
c      enddo

return
end ! of subroutine mass_enclosed
cc
cc
cc
subroutine cluster_rotation(sig, rmin, ff)

use real_prec
c      use config
use param
use grid
use field
use root
use cons
use gravmod
use agn

implicit none
C     This gives solid body rotation in the inner part and constant
C     specific angular momentum in the outer part.
C     The specific angular momentum set by the BH accretion disk.
C     The switchover between the two is when the const angular momentum
C     implies a velocity FC_ROTATION the sound speed of the cluster.
C     The relevant radius is the radius perp to the z axis, not the r coordinate
C     The velocities go into the v_phi coordinate
C     ll = angular momentum per unit mass
REAL(rl) ::  sig,rmin, ff
REAL(rl) ::  cs, x, r2d
REAL(rl) ::  ex_rmin, rs, conc, delta, lambda, th, zf, rvir
INTEGER i,j,k

rs = 2*guniv*mbh/CVEL**2
cs = sqrt(gamma)*sig

write(2,*) "Used rmin = ", rmin
C     JPO suggestion that set specific angular momentum according
C     to the angular-momentum profile of accretion disk
ex_rmin = 0.5*(CVEL/cs)*sqrt(2*nscho)*rs
write(2,*) "JPO rmin = ", ex_rmin

C     Set rmin appropriate to halo:
C     NOTE -- expression for virial radius depends on units
th = 13.6
zf = 1.0
delta = 200
lambda = 0.05
conc = 10
rvir = sig*th*sqrt(2.0/delta)/(1+zf)**1.5

C     Assuming relevant r is virial radius
ex_rmin = lambda * rvir
write(2,*) "rvir rmin = ", ex_rmin

ex_rmin = lambda * rvir / conc
write(2,*) "rscale rmin = ", ex_rmin

do k=ks,ke
do j=js,je
do i=is,ie
r2d = x1a(i)*sin(x2a(j))
x = r2d/rmin
v3(i,j,k) = cs/(x + 1/(ff*x))
enddo
enddo
enddo

return
end ! of subroutine cluster_rotation

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

use real_prec
use config
use param
use grid

implicit none
REAL(rl) ::  pzpt, pzpq, prpt, prpq
INTEGER  ::  i, j, k
CHARACTER    flag
if(flag.eq."a") then      ! a-mesh

if(LGEOM .eq. 1) then
r2   = ( x1a(i) - X1C )**2
.          + ( x2a(j) - X2C )**2
.          + ( x3a(k) - X3C )**2
endif

if(LGEOM .eq. 3) then
if(X1C .eq. zro)then
r2 = x1a(i)**2
else
prpq  = x1a(i) * sin(x2a(j))
pzpq  = x1a(i) * cos(x2a(j))
prpt  = X1C  * sin(X2C)
pzpt  = X1C  * cos(X2C)

r2   = ( pzpq - pzpt )**2
.           + ( prpq - prpt )**2
.           +   2.0 * prpq * prpt
.           * ( 1.0 - cos(x3a(k) - X3C) )
endif
endif

elseif(flag.eq."b") then  ! b-mesh

if(LGEOM .eq. 1) then
r2   = ( x1b(i) - X1C )**2
.          + ( x2b(j) - X2C )**2
.          + ( x3b(k) - X3C )**2
endif

if(LGEOM .eq. 3) then
if(X1C .eq. zro)then
r2 = x1b(i)**2
else
prpq  = x1b(i) * sin(x2b(j))
pzpq  = x1b(i) * cos(x2b(j))
prpt  = X1C  * sin(X2C)
pzpt  = X1C  * cos(X2C)

r2   = ( pzpq - pzpt )**2
.           + ( prpq - prpt )**2
.           +   2.0 * prpq * prpt
.           * ( 1.0 - cos(x3b(k) - X3C) )
endif
endif

else
write(*,*) "wrong FLAG !"
stop

endif ! flag