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
      REAL(rl) ::  radii2 
      INTEGER  ::  i, j, k

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

      call cluster_read  !galactic "profile" parameters
     .          (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)
      call star_read
#endif /* STAREV || STARFORM */

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

#ifdef    RADCOOL
      call ism_rad_read
#endif /* RADCOOL */

CC read the artificial limits
#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  
               r2a(i,j,k) = radii2("a",i,j,k)
               r2b(i,j,k) = radii2("b",i,j,k)
            enddo
         enddo
      enddo

      return 
      end

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

      subroutine cluster_read(sig,rc,d0,beta,rmin,ff,mstar,reff)
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)      
         read  (1,cluster)
         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
      end ! of subroutine cluster_read

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 

C     Assuming relevant r is scale radius radius
      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

      function radii2(flag,i,j,k)

      use real_prec
      use config
      use param
      use grid

      implicit none 
      REAL(rl) ::  pzpt, pzpq, prpt, prpq
      REAL(rl) ::  r2, radii2
      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

      radii2 = r2

      return
      end ! of function radii2