 	implicit real*8 (a-h,o-z)
	parameter(nbodsmax=400000,nsphmax=200000,ndim=3)
	parameter(nghost=6)
	parameter(ngrid_without_ghostzones=256) ! ***** Should be divisible by 4 *******
	parameter(ngrid=ngrid_without_ghostzones+nghost)
	parameter(ngrid2=ngrid/2)
	parameter(nfra=1)
	parameter(gam=5.0/3.0)
	parameter(gam1=gam-1.0)
	parameter(con=3.0e0/3.141592e0/32.0e0)
        INTEGER, PARAMETER :: nlist_kind = SELECTED_INT_KIND(6)
	INTEGER, PARAMETER :: nnabomax = 800
        INTEGER(nlist_kind) n_list(nnabomax, nsphmax)
        INTEGER(nlist_kind) nnabo(nsphmax)
	character*80 ibodfile
	character*3 cfi,pjd
	logical projdir
        real*8 poscmold(ndim)
        real*8 jugas(3),justar(3),ju(3)
	real*8 ss(100),ss2(100),ppp(100),pp2(100),l,b,sumrr(100)
        real*8 x(nbodsmax),y(nbodsmax),z(nbodsmax),vx(nbodsmax)
        real*8 vy(nbodsmax),vz(nbodsmax)
        real*8 x2(nbodsmax),y2(nbodsmax),z2(nbodsmax),vx2(nbodsmax)
        real*8 vy2(nbodsmax),vz2(nbodsmax)
 	real*8 mass,tn(1000),jt(1000,4,10),posc(3),velc(3),jx,jy,jz,
     &    jtt(1000,4,10)
	real*8 po(nbodsmax,ndim),ve(nbodsmax,ndim),jtave(1000,4,10)
	real*8 pocm(nbodsmax,ndim),vecm(nbodsmax,ndim),jxc,jyc,jzc
	real*8 rm(nbodsmax),rr(nbodsmax)
	real*8 pot200,distkeep(nbodsmax),velkeep(nbodsmax),ang(nbodsmax) 
	real*8 lnh(nbodsmax),tl(nbodsmax),lnhmax,poscm(3),poscm1(3)
	real*8 velcm(3),jtot,j4(nfra,nfra),rh(1000,2),rrr(1000),m
        real*8 vx1(nbodsmax),vy1(nbodsmax),vz1(nbodsmax),rtot(nbodsmax)
	real*8 mdisk,mbulge,lsigma,lsigma0,mbd,mbulge0,mdisk0
        INTEGER p,pp,ubodsin,ubodsout,num(0:1000),pmax,jump(nbodsmax)
	dimension v(nbodsmax,ndim),rs(nbodsmax)
	dimension mass(nbodsmax),pos(nbodsmax,ndim),
     &	vel(nbodsmax,ndim),
     &	  ener(nbodsmax),phi(nbodsmax),rho(nbodsmax),
     &	smooth(nbodsmax)
	real*8 abundance(nbodsmax,9),tborn(nbodsmax)
        real*8 abundancesm(nbodsmax,9)
        real*8 abundancetot(nbodsmax),abundancetottot(nbodsmax)
	integer indx(nbodsmax),nbo(nfra) !,indx1(nbodsmax)
	dimension count(0:100),dens(0:100)
	real*8 lT(11),l_23(11),lTmin,lambdaT,lXbol,nH(nbodsmax),lXbolsm
	real*8 lT_tab(89),cf_metal_tab(89,7),oh_tab(7),lambdaM,lambdaMsm
	real*8 lXbo,lXbocum
	INTEGER nrtot(50),indxtot(50,nbodsmax),p1,ninow(1000)
	real*8 tnowm(1000),starmass(1000),tbornm(1000)
	real*8 om(0:100),mgm(0:100),fem(0:100),mghm,mgmcum,mgmcum1,mghmcum

	real readin(nbodsmax),tnowr

	real*8 logtlim,lognHlim
	integer indxhot(nsphmax),gridx,gridy,gridz,gridxx,gridyy,gridzz
	real*8 posg(3),wsumv(3),gr(0:10000),massir(10000)
	integer i_brute_force(2)

 	real*4 
     1   posgrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,3)
	real*4
     1   velgrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,3)
	real*4
     1   grgrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,3)
	real*4 
     1   energrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1)
	real*4 
     1   rhogrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1)
	real*4 
     1   smogrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1)
	real*4 
     1   sgrid(-ngrid2:ngrid2-1,-ngrid2:ngrid2-1,-ngrid2:ngrid2-1)

	real*4 readout(ngrid**3),dlr4,larr(-ngrid2:ngrid2-1)

	COMMON/tabs/lT_tab,cf_metal_tab,oh_tab

	common/temp/temp(0:12,-8:10,-2:10)

c       NCWH, 24/3/2004:
c       Put the big variables in a common block to force them on heap
c       otherwise they would go on the stack.  

        common/ncwh/posgrid,velgrid,grgrid, energrid, rhogrid,smogrid

	CALL window_setup
c
c	open ttable
c
	open(file='ttable.dat',unit=14,status='old')
c
	do nt=0,12
	do nr=-8,10
	read(14,*)(temp(nt,nr,nu),nu=-2,10)
	end do
	end do
c
	close(14)


c read in data from metal_cool_diff_table.dat
c----------------------------------------------
       OPEN(110,FILE='metal_cool_diff_table.dat',STATUS='OLD')
      
      do k=1,89
         read(110,*)lT_tab(k),cf_metal_tab(k,:)
      enddo

      CLOSE(110)

      oh_tab(1)=-2.5d0
      oh_tab(2)=-1.5d0
      oh_tab(3)=-1.0d0
      oh_tab(4)=-0.5d0
      oh_tab(5)=-0.25d0
      oh_tab(6)=0.0d0
      oh_tab(7)=0.5d0

c       Cooling function for primordial gas

        lT(1)=4.3
        lT(2)=5
        lT(3)=5.3
        lT(4)=5.6
        lT(5)=6
        lT(6)=6.3
        lT(7)=6.6
        lT(8)=7
        lT(9)=7.3
        lT(10)=7.6
	lT(11)=10.0

        l_23(1)=17.0
        l_23(2)=2.8
        l_23(3)=0.52
        l_23(4)=0.48
        l_23(5)=0.52 
        l_23(6)=0.55
        l_23(7)=0.62
        l_23(8)=0.83
        l_23(9)=1.1
        l_23(10)=1.45
	l_23(11)=23.0

	nTmax=11

	do i=1,nTmax
	   l_23(i)=log10(l_23(i))
	enddo


	grav=4.30d-6
	pi=acos(-1.0)
	tiny=1.0d-16

        rcrit=400000000 !340 !348 !1680 !4000
	rcritlow=0 !1714 !1870

	ncoldgasold=0
	nstarold=0
	sfrm=0
	tnowm=0
	starmassm=0
	ninow=0
c
c
C   Open input body data file.
C   --------------------------
	ubodsin=10
	ubodsout=16
c	print*,'infile?'
c	read(*,10)ibodfile
c	print*,' '
10	format(a)
12	format(f6.2)
c
c        OPEN(UNIT=ubodsin,FILE=ibodfile,STATUS='old',form='unformatted')     
        OPEN(UNIT=ubodsin,FILE='coma.12.3.stack',STATUS='old',
     1    form='unformatted')     
c	print*,'outfile?' 
c	read(*,10)ibodfile
c	print*,' '
c        OPEN(UNIT=ubodsout,FILE=ibodfile,STATUS='new') 
c
	read(ubodsin)nf
c	print*,nf
c
c	print*,'Frame#, rvir [Mpc] - please!'
c	read*,nff,rvir
c	rvir=rvir*1000.0
	nff=1
c	rvir=2800.0
c	projdir=.FALSE.
c	print*,'Projection direction (0 for x, 1 for y) - please'
c	read*,iprd
c	if(iprd==1) projdir=.TRUE.
c	if(projdir.EQV..FALSE.) print*,'x direction chosen'
c	if(projdir.EQV..TRUE.) print*,'y direction chosen'

	DO iff=1,nff
	print*,iff
	if=1
c 
C   Read in body data.
C   ------------------
        read(ubodsin)nbodies,nsph,ndm,nstars,ndimo,tnowr
	tnow=dble(tnowr)
	nbo(if)=nbodies
c	if(if.eq.1) then
c	print*,'nf,nbodies,ndimo: ',nf,nbodies,ndimo
c	print*,' '
c	endif
c	print*,'time'
c	print*,' ' 

	read(ubodsin)(readin(p),p=1,nbodies)
	mass=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	ener=dble(readin)

	read(ubodsin)(readin(p),p=1,nbodies)
	pos(:,1)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	pos(:,2)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	pos(:,3)=dble(readin)

	read(ubodsin)(readin(p),p=1,nbodies)
	vel(:,1)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	vel(:,2)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	vel(:,3)=dble(readin)

	read(ubodsin)(readin(p),p=1,nbodies)
	abundance(:,4)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	abundancesm(:,4)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	abundance(:,5)=dble(readin)
        read(ubodsin)(readin(p),p=1,nbodies)
        abundancesm(:,5)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	abundance(:,6)=dble(readin)
        read(ubodsin)(readin(p),p=1,nbodies)
        abundancesm(:,6)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	abundance(:,9)=dble(readin)
        read(ubodsin)(readin(p),p=1,nbodies)
        abundancesm(:,9)=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	abundancetot=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	abundancetottot=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	rho=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	smooth=dble(readin)
	read(ubodsin)(readin(p),p=1,nbodies)
	tborn=dble(readin)
c       
	if(iff<nff) CYCLE
        icrit=0

        do j=1,3
        poscmold(j)=0.0
        enddo
c

	m=mass(1)/grav

	nbodies=nbo(if)
c
14	n1=1
	n2=nbodies
c	
	do p=1,nbodies
	if(smooth(p).eq.0.)then
	n3=p-1
	n2=n3
	go to 300
	endif
	enddo
c	endif
300	continue
c
        nsph0=nsph+nstars
        istar0=nbodies-nsph0+n2
        nsphdm=istar0 
c
	nsph=n2-n1+1
	ndm=nsphdm-nsph
	nstar=nbodies-nsphdm
c
	print*,nsph,ndm,nstar	 

c
	time=tnow
	if(time.gt.6.) time=6.
c
	fac=.00722958
	deltae=0.01
	deltarho=1.0d-9
	lnhmax=-1.0d9
c
	pmax=0
c	
	sumf=0.0
	summ=0.0
c
	do j=1,3
	poscm1(j)=0.0
	poscm(j)=0.0
	posc(j)=0.0
	velcm(j)=0.0
	velc(j)=0.0
	end do
c
	do p=1,nbodies
        sum=0.0
        do j=1,3
        sum=sum+(pos(p,j)-poscmold(j))**2
        enddo
        r=sqrt(sum)
	if(p.gt.n2) goto 456
	do j=1,3
	poscm(j)=poscm(j)+mass(p)*pos(p,j)
	velcm(j)=velcm(j)+mass(p)*vel(p,j)
	end do
	summ=summ+mass(p)
c	if(if.eq.1) goto 567
	 call conuffe(time,log10(fac*rho(p)+deltarho)
     &   ,log10(ener(p)+deltae+0.1),templ)
        lnh(p)=log10(fac*rho(p)+deltarho)
	nH(p)=10.0**lnh(p)
	tl(p)=templ
c	write(17,*)p,ener(p),tl(p),lnh(p),iff
        if(icrit.eq.1.and.r>rcrit) goto 567 
        if(icrit.eq.1.and.r<rcritlow) goto 567 
	if(lnh(p).ge.lnhmax) then
	pmax=p
	lnhmax=lnh(p)
	end if
	goto 567
456	do j=1,3
	poscm(j)=poscm(j)+mass(p)*pos(p,j)
	velcm(j)=velcm(j)+mass(p)*vel(p,j)
	end do
	summ=summ+mass(p)
	sumf=sumf+fdm
567	end do
c
c	if(if.ne.1) then
         do j=1,3
         poscm1(j)=pos(pmax,j)
         end do
c	endif
c
	do j=1,3
	poscm(j)=poscm(j)/summ
	velcm(j)=velcm(j)/summ
	end do
c
	nit=0
c
234	ncm=0
c
	nit=nit+1
c
	do j=1,3
	posc(j)=0.0
	velc(j)=0.0
	end do
c
	summ=0.0
c
c	if(if.eq.1) goto 666
c
	rmax=500.0		!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
	do p1=1,nsph+nstar
	   if(p1<=n2) p=p1
	   if(p1>n2) p=p1+ndm
	sum=0.0
	do j=1,3
	a=pos(p,j)-poscm1(j)
	sum=sum+a**2
	end do
	rp=sqrt(sum)
c
	if(rp.ge.rmax) goto 123
c
	do j=1,3
	posc(j)=posc(j)+mass(p)*pos(p,j)
	velc(j)=velc(j)+mass(p)*vel(p,j)
	end do
	summ=summ+mass(p)
	ncm=ncm+1
123	end do
c
	do j=1,3
	posc(j)=posc(j)/summ
	velc(j)=velc(j)/summ
	end do
c
	sum=0.0
	do j=1,3
	sum=sum+(poscm1(j)-posc(j))**2
	poscm1(j)=posc(j)
	end do
c
	if(sqrt(sum).ge.0.0001) goto 234
c
	print*,'tnow,log(nH_max),ncm ',tnowr,lnhmax,ncm
c
c	print*,'number of iterations: ',nit
c
c666	continue
c **************************************
	rmax=100.0		!!!!!!!!!!!!!!!!!!!!!!!!
444	ncm=0
	summ=0.
	do p=nsphdm+1,nbodies
	  sum=0.
	  do j=1,3
	  a=pos(p,j)-poscm1(j)
	  sum=sum+a**2.
	  enddo
	 rp=sqrt(sum)
	 if(rp.ge.rmax) goto 333
	  do j=1,3
	  posc(j)=posc(j)+mass(p)*pos(p,j)
	  velc(j)=velc(j)+mass(p)*vel(p,j)
	  enddo
	summ=summ+mass(p)
	ncm=ncm+1
333	enddo
	
	do j=1,3
	posc(j)=posc(j)/summ
	velc(j)=velc(j)/summ
	enddo

	sum=0.
	do j=1,3
	sum=sum+(poscm1(j)-posc(j))**2
	poscm1(j)=posc(j)
	enddo
	if(sqrt(sum).ge.0.0001) goto 444
	print*,'tnow,log(nH_max),ncm ',tnowr,lnhmax,ncm
c****************************************
	rmax=5.0		!!!!!!!!!!!!!!!!!!!!!!!!
544	ncm=0
	summ=0.
	do p=nsphdm+1,nbodies
	  sum=0.
	  do j=1,3
	  a=pos(p,j)-poscm1(j)
	  sum=sum+a**2.
	  enddo
	 rp=sqrt(sum)
	 if(rp.ge.rmax) goto 533
	  do j=1,3
	  posc(j)=posc(j)+mass(p)*pos(p,j)
	  velc(j)=velc(j)+mass(p)*vel(p,j)
	  enddo
	summ=summ+mass(p)
	ncm=ncm+1
533	enddo
	
	do j=1,3
	posc(j)=posc(j)/summ
	velc(j)=velc(j)/summ
	enddo

	sum=0.
	do j=1,3
	sum=sum+(poscm1(j)-posc(j))**2
	poscm1(j)=posc(j)
	enddo
	if(sqrt(sum).ge.0.0001) goto 544
	print*,'tnow,log(nH_max),ncm ',tnowr,lnhmax,ncm
c****************************************
c        posc(1)=2040.586182 
c	 posc(2)=-5056.095703 
c        posc(3)=7777.654785
	do p=n1,nbodies
	do j=1,3
	pocm(p,j)=pos(p,j)-poscm(j)
	vecm(p,j)=vel(p,j)-velcm(j)	
	po(p,j)=pos(p,j)-posc(j)
	ve(p,j)=vel(p,j)-velc(j)
	end do
	end do

	sum=0.0
	do j=1,3
	sum=sum+(poscm(j)-posc(j))**2
        poscmold(j)=poscm(j)
	end do
	dcm=sqrt(sum)
	print*,'dcm= ',dcm
        icrit=1
        if(dcm.gt.rcrit.or.dcm<rcritlow) goto 14  
c
	if(n1.ne.1) then
	print*,'n1 problems ',1,n1
	stop
	endif
	print*,'posc',real(posc)
	print*,'velc',real(velc)
c
c*******************************************************************
c       Calculate gravity accelerations for the spherically averaged
c       mass distribution - zero gravity softenings assumed 
c*******************************************************************

	drg=1.0       ![kpc]
	rmaxg=3000.0  ![kpc]
	
	massir=0
	do i=1,nbodies
	   dist=sqrt(po(i,1)**2+po(i,2)**2+po(i,3)**2)
	   if(dist>=rmaxg) CYCLE
	   ir=int(dist/drg)+1
	   massir(ir)=massir(ir)+mass(i)
	enddo

	gr=0
	cummass=0
	do ir=1,nint(rmaxg/drg)
	   cummass=cummass+massir(ir)
	   gr(ir)=cummass/(dble(ir)*drg)**2   ! GM(<r)/r**2
	enddo


c*********************************************************
c       IDENTIFY HOT GAS PARTICLES
c*********************************************************
	logtlim=0 !6
	lognHlim=10 !-1.0
	print*,'rmax'
	read*,rmax
	nhot=0
	smoothmin=999999
	do i=1,nsph
	   sum=0
	   do j=1,3
	      sum=sum+po(i,j)**2
	   enddo
	   cdist=sqrt(sum)
	   if(cdist>rmax) CYCLE
	   if(tl(i)>logtlim.and.lnH(i)<lognHlim) then
	      nhot=nhot+1
	      indxhot(nhot)=i
	      if(smooth(i)<smoothmin) then
		 smoothmin=smooth(i)
		 lnhmax=lnH(i)
              endif
	   endif
	enddo

	print*,nhot,2*smoothmin,lnhmax

	nnabosum=0
	nnabomaxx=0
	nnabomin=9999999
	do ihot=1,nhot
	   i=indxhot(ihot)
	   w2=(4.0*smooth(i))**2
	   nn=0
	   do jhot=1,nhot
	      j=indxhot(jhot)
	      rij2=0.0
	      do jj=1,3
		 a=po(i,jj)-po(j,jj)
		 rij2=rij2+a**2
	      enddo
	      if(rij2<w2) then
		 nn=nn+1
		 n_list(nn,ihot)=jhot
	      endif
	   enddo
	   nnabo(ihot)=nn
	   if(nn>nnabomaxx) nnabomaxx=nn
	   if(nn<nnabomin) nnabomin=nn
	   nnabosum=nnabosum+nn
	enddo

	do ihot=nhot-1,nhot
	   print*,nnabo(ihot)
	   do jhot=1,nnabo(ihot)
	      write(26,*)jhot,n_list(jhot,ihot),indxhot(n_list(jhot,ihot))
	   enddo
	enddo

	print*,nnabosum/nhot,nnabomin,nnabomaxx

c
c       Start interpolating onto a grid
c       
c	dl=2.5            ! Cell-length in kpc (2.5 for testing at 32^3)
	dl=25./2.0            ! Cell-length in kpc
	N=ngrid-nghost    ! Grid-size (Number divisible by 4)
	N2=N/2
	N4=N/4
	rhomax=138.9/10.0 ! max allowed density log(n_H)=-1 
                          ! Temp is changed to preserve pressure

c       Find the nearest (hot) particle to gridpoint (-N/2,-N/2,-N/2)
c
	posg(1)=-dble(N2)*dl+0.5*dl
	posg(2)=-dble(N2)*dl+0.5*dl
	posg(3)=-dble(N2)*dl+0.5*dl

	distmin=99999999
	do ihot=1,nhot
	   i=indxhot(ihot)
	   rig2=0.0
	   do jj=1,3
	      a=po(i,jj)-posg(jj)
	      rig2=rig2+a**2
	   enddo
	   dist=sqrt(rig2)
	   if(dist<distmin) then
	      distmin=dist
	      ihotmin=ihot
	   endif
	enddo

c       Check whether 2*smooth>dl
c       Check whether 2*smooth>distmin
c
	i=indxhot(ihotmin)
	if(dl>2.0*smooth(i)) then
	   print*,'dl > 2.0*smooth'
	   stop
	endif
	if(distmin>2.0*smooth(i)) then
	   print*,'distmin > 2.0*smooth'
	   stop
	endif
	print*,posg
	print*,ihotmin,i,distmin,2.0*smooth(i)

c       Calculate spline-kernal smoothed gas properties at (-N/2,-N/2,-N/2)
c



	print*,'start main grid loop'
c
c       Now loop over entire grid by "drawing" a space filling curve
c
	do gridz=-N4,N4-1
	   print*,gridz
	   do gridy=-N4,N4-1
	      do gridx=-N2,N2-1
c	if(gridz==-N4.and.gridy==-N4.and.gridx==-N2) CYCLE !Done already
	
		 posg(1)=dble(gridx)*dl+0.5*dl
		 posg(2)=dble(2*gridy)*dl+0.5*dl
		 posg(3)=dble(2*gridz)*dl+0.5*dl
		 
		 posgrid(gridx,2*gridy,2*gridz,:)=posg(:)

c       Search for the nearest particle in n_list of ihotmin_i
		 distmin=99999999
		 do j=1,nnabo(ihotmin)
		    i=indxhot(n_list(j,ihotmin))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin1=n_list(j,ihotmin)
		     endif
		  enddo

c       Search for the nearest particle in n_list of ihotmin_i+1
		 distmin=99999999
		 do j=1,nnabo(ihotmin1)
		    i=indxhot(n_list(j,ihotmin1))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin2=n_list(j,ihotmin1)
		     endif
		  enddo

		  i=indxhot(ihotmin2)

c       Check whether 2*smooth_i+1>dl
c       Check whether 2*smooth_i+1>distmin
c
	i_brute_force=0
	if(dl>2.0*smooth(i)) then
	   i_brute_force(1)=1
	endif
	if(distmin>2.0*smooth(i)) then
	   i_brute_force(2)=1
	endif

	if(i_brute_force(1)/=0.or.i_brute_force(2)/=0) then 

c       Find the nearest (hot) particle to the gridpoint (Brute force) 
c
	distmin=99999999
	do ihot=1,nhot
	   i=indxhot(ihot)
	   rig2=0.0
	   do jj=1,3
	      a=po(i,jj)-posg(jj)
	      rig2=rig2+a**2
	   enddo
	   dist=sqrt(rig2)
	   if(dist<distmin) then
	      distmin=dist
	      ihotmin2=ihot 
	   endif
	enddo
	i=indxhot(ihotmin2)
		  write(27,*)real(posg),i_brute_force
		  write(27,*)ihotmin2,indxhot(ihotmin2),real(distmin),
     1     real(2.0*smooth(i)),real(rho(i)),real(ener(i))


	endif

	ihotmin=ihotmin2

	smoothnew=smooth(indxhot(ihotmin))

         wsumv(:)=0.0d0
         wsumgasmass=0.0d0
         wsumgasener=0.0d0
         wsumgasrho=0.0d0
         wsumgassmooth=0.0d0

	do j=1,nnabo(ihotmin)
	   pp=indxhot(n_list(j,ihotmin))
	   sum=0
	   do jj=1,3
	      a=po(pp,jj)-posg(jj)
	      sum=sum+a**2
	   enddo
	   dist=sqrt(sum)
          wsumgasmass=wsumgasmass+mass(pp)*windowheat(dist,
     1              smoothnew)
          wsumv(:)=wsumv(:)+mass(pp)*windowheat(dist,
     1              smoothnew)*vel(pp,:)
          wsumgasener=wsumgasener+mass(pp)*windowheat(dist,
     1              smoothnew)*ener(pp)	   
          wsumgasrho=wsumgasrho+mass(pp)*windowheat(dist,
     1              smoothnew)*rho(pp)
          wsumgassmooth=wsumgassmooth+mass(pp)*windowheat(dist,
     1              smoothnew)*smooth(pp)
	enddo

	   
	energrid(gridx,2*gridy,2*gridz)=wsumgasener/wsumgasmass
	rhogrid(gridx,2*gridy,2*gridz)=wsumgasrho/wsumgasmass
	rhofac=wsumgasrho/wsumgasmass/rhomax
	if(rhofac>1.0) then
	   rhogrid(gridx,2*gridy,2*gridz)=rhomax
	   energrid(gridx,2*gridy,2*gridz)=
     1       energrid(gridx,2*gridy,2*gridz)*rhofac
	endif
	smogrid(gridx,2*gridy,2*gridz)=wsumgassmooth/wsumgasmass
	velgrid(gridx,2*gridy,2*gridz,:)=wsumv(:)/wsumgasmass


	dist=sqrt(posg(1)**2+posg(2)**2+posg(3)**2)
	if(dist>=rmaxg) stop
	ir=int(dist/drg)
	distir=ir*drg
	del=(dist-distir)/drg

	gravacc=gr(ir)+(gr(ir+1)-gr(ir))*del

	grgrid(gridx,2*gridy,2*gridz,:)=-gravacc*posg(:)/dist



	      enddo

c	      
c       ----------------------------------------------------------------
c

	      do gridxx=-(N2-1),N2
		 gridx=-gridxx
		 
		 posg(1)=dble(gridx)*dl+0.5*dl
		 posg(2)=dble(2*gridy+1)*dl+0.5*dl
		 posg(3)=dble(2*gridz)*dl+0.5*dl

		 posgrid(gridx,2*gridy+1,2*gridz,:)=posg(:)

c       Search for the nearest particle in n_list of ihotmin_i
		 distmin=99999999
		 do j=1,nnabo(ihotmin)
		    i=indxhot(n_list(j,ihotmin))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin1=n_list(j,ihotmin)
		     endif
		  enddo

c       Search for the nearest particle in n_list of ihotmin_i+1
		 distmin=99999999
		 do j=1,nnabo(ihotmin1)
		    i=indxhot(n_list(j,ihotmin1))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin2=n_list(j,ihotmin1)
		     endif
		  enddo

		  i=indxhot(ihotmin2)

c       Check whether 2*smooth_i+1>dl
c       Check whether 2*smooth_i+1>distmin
c
	i_brute_force=0
	if(dl>2.0*smooth(i)) then
	   i_brute_force(1)=1
	endif
	if(distmin>2.0*smooth(i)) then
	   i_brute_force(2)=1
	endif

	if(i_brute_force(1)/=0.or.i_brute_force(2)/=0) then 

c       Find the nearest (hot) particle to the gridpoint (Brute force) 
c
	distmin=99999999
	do ihot=1,nhot
	   i=indxhot(ihot)
	   rig2=0.0
	   do jj=1,3
	      a=po(i,jj)-posg(jj)
	      rig2=rig2+a**2
	   enddo
	   dist=sqrt(rig2)
	   if(dist<distmin) then
	      distmin=dist
	      ihotmin2=ihot 
	   endif
	enddo
	i=indxhot(ihotmin2)
		  write(27,*)real(posg),i_brute_force
		  write(27,*)ihotmin2,indxhot(ihotmin2),real(distmin),
     1     real(2.0*smooth(i)),real(rho(i)),real(ener(i))


	endif

	ihotmin=ihotmin2

	smoothnew=smooth(indxhot(ihotmin))

         wsumv(:)=0.0d0
         wsumgasmass=0.0d0
         wsumgasener=0.0d0
         wsumgasrho=0.0d0
         wsumgassmooth=0.0d0

	do j=1,nnabo(ihotmin)
	   pp=indxhot(n_list(j,ihotmin))
	   sum=0
	   do jj=1,3
	      a=po(pp,jj)-posg(jj)
	      sum=sum+a**2
	   enddo
	   dist=sqrt(sum)
          wsumgasmass=wsumgasmass+mass(pp)*windowheat(dist,
     1              smoothnew)
          wsumv(:)=wsumv(:)+mass(pp)*windowheat(dist,
     1              smoothnew)*vel(pp,:)
          wsumgasener=wsumgasener+mass(pp)*windowheat(dist,
     1              smoothnew)*ener(pp)	   
          wsumgasrho=wsumgasrho+mass(pp)*windowheat(dist,
     1              smoothnew)*rho(pp)
          wsumgassmooth=wsumgassmooth+mass(pp)*windowheat(dist,
     1              smoothnew)*smooth(pp)
	enddo

	   
	energrid(gridx,2*gridy+1,2*gridz)=wsumgasener/wsumgasmass
	rhogrid(gridx,2*gridy+1,2*gridz)=wsumgasrho/wsumgasmass
	rhofac=wsumgasrho/wsumgasmass/rhomax
	if(rhofac>1.0) then
	   rhogrid(gridx,2*gridy+1,2*gridz)=rhomax
	   energrid(gridx,2*gridy+1,2*gridz)=
     1       energrid(gridx,2*gridy+1,2*gridz)*rhofac
	endif
	smogrid(gridx,2*gridy+1,2*gridz)=wsumgassmooth/wsumgasmass
	velgrid(gridx,2*gridy+1,2*gridz,:)=wsumv(:)/wsumgasmass

	dist=sqrt(posg(1)**2+posg(2)**2+posg(3)**2)
	if(dist>=rmaxg) stop
	ir=int(dist/drg)
	distir=ir*drg
	del=(dist-distir)/drg

	gravacc=gr(ir)+(gr(ir+1)-gr(ir))*del

	grgrid(gridx,2*gridy+1,2*gridz,:)=-gravacc*posg(:)/dist

	      enddo
	      
	   enddo

c
c       --------------------------------------------------------
c

	   do gridyy=-(N4-1),N4
	      gridy=-gridyy
	      do gridx=-N2,N2-1
		
		 posg(1)=dble(gridx)*dl+0.5*dl
		 posg(2)=dble(2*gridy+1)*dl+0.5*dl
		 posg(3)=dble(2*gridz+1)*dl+0.5*dl

		 posgrid(gridx,2*gridy+1,2*gridz+1,:)=posg(:)

c       Search for the nearest particle in n_list of ihotmin_i
		 distmin=99999999
		 do j=1,nnabo(ihotmin)
		    i=indxhot(n_list(j,ihotmin))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin1=n_list(j,ihotmin)
		     endif
		  enddo

c       Search for the nearest particle in n_list of ihotmin_i+1
		 distmin=99999999
		 do j=1,nnabo(ihotmin1)
		    i=indxhot(n_list(j,ihotmin1))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin2=n_list(j,ihotmin1)
		     endif
		  enddo

		  i=indxhot(ihotmin2)

c       Check whether 2*smooth_i+1>dl
c       Check whether 2*smooth_i+1>distmin
c
	i_brute_force=0
	if(dl>2.0*smooth(i)) then
	   i_brute_force(1)=1
	endif
	if(distmin>2.0*smooth(i)) then
	   i_brute_force(2)=1
	endif

	if(i_brute_force(1)/=0.or.i_brute_force(2)/=0) then 

c       Find the nearest (hot) particle to the gridpoint (Brute force) 
c
	distmin=99999999
	do ihot=1,nhot
	   i=indxhot(ihot)
	   rig2=0.0
	   do jj=1,3
	      a=po(i,jj)-posg(jj)
	      rig2=rig2+a**2
	   enddo
	   dist=sqrt(rig2)
	   if(dist<distmin) then
	      distmin=dist
	      ihotmin2=ihot 
	   endif
	enddo
	i=indxhot(ihotmin2)
		  write(27,*)real(posg),i_brute_force
		  write(27,*)ihotmin2,indxhot(ihotmin2),real(distmin),
     1     real(2.0*smooth(i)),real(rho(i)),real(ener(i))


	endif

	ihotmin=ihotmin2

	smoothnew=smooth(indxhot(ihotmin))

         wsumv(:)=0.0d0
         wsumgasmass=0.0d0
         wsumgasener=0.0d0
         wsumgasrho=0.0d0
         wsumgassmooth=0.0d0

	do j=1,nnabo(ihotmin)
	   pp=indxhot(n_list(j,ihotmin))
	   sum=0
	   do jj=1,3
	      a=po(pp,jj)-posg(jj)
	      sum=sum+a**2
	   enddo
	   dist=sqrt(sum)
          wsumgasmass=wsumgasmass+mass(pp)*windowheat(dist,
     1              smoothnew)
          wsumv(:)=wsumv(:)+mass(pp)*windowheat(dist,
     1              smoothnew)*vel(pp,:)
          wsumgasener=wsumgasener+mass(pp)*windowheat(dist,
     1              smoothnew)*ener(pp)	   
          wsumgasrho=wsumgasrho+mass(pp)*windowheat(dist,
     1              smoothnew)*rho(pp)
          wsumgassmooth=wsumgassmooth+mass(pp)*windowheat(dist,
     1              smoothnew)*smooth(pp)
	enddo

	   
	energrid(gridx,2*gridy+1,2*gridz+1)=wsumgasener/wsumgasmass
	rhogrid(gridx,2*gridy+1,2*gridz+1)=wsumgasrho/wsumgasmass
	rhofac=wsumgasrho/wsumgasmass/rhomax
	if(rhofac>1.0) then
	   rhogrid(gridx,2*gridy+1,2*gridz+1)=rhomax
	   energrid(gridx,2*gridy+1,2*gridz+1)=
     1       energrid(gridx,2*gridy+1,2*gridz+1)*rhofac
	endif
	smogrid(gridx,2*gridy+1,2*gridz+1)=wsumgassmooth/wsumgasmass
	velgrid(gridx,2*gridy+1,2*gridz+1,:)=wsumv(:)/wsumgasmass

	dist=sqrt(posg(1)**2+posg(2)**2+posg(3)**2)
	if(dist>=rmaxg) stop
	ir=int(dist/drg)
	distir=ir*drg
	del=(dist-distir)/drg

	gravacc=gr(ir)+(gr(ir+1)-gr(ir))*del

	grgrid(gridx,2*gridy+1,2*gridz+1,:)=-gravacc*posg(:)/dist


	      enddo
c 
c       -------------------------------------------------------------
c      
	      do gridxx=-(N2-1),N2
		 gridx=-gridxx
		
		 posg(1)=dble(gridx)*dl+0.5*dl
		 posg(2)=dble(2*gridy)*dl+0.5*dl
		 posg(3)=dble(2*gridz+1)*dl+0.5*dl

		 posgrid(gridx,2*gridy,2*gridz+1,:)=posg(:)

c       Search for the nearest particle in n_list of ihotmin_i
		 distmin=99999999
		 do j=1,nnabo(ihotmin)
		    i=indxhot(n_list(j,ihotmin))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin1=n_list(j,ihotmin)
		     endif
		  enddo

c       Search for the nearest particle in n_list of ihotmin_i+1
		 distmin=99999999
		 do j=1,nnabo(ihotmin1)
		    i=indxhot(n_list(j,ihotmin1))
		    sum=0
		    do jj=1,3
		       a=po(i,jj)-posg(jj)
		       sum=sum+a**2
		     enddo
		     dist=sqrt(sum)
		     if(dist<distmin) then
			distmin=dist
			ihotmin2=n_list(j,ihotmin1)
		     endif
		  enddo

		  i=indxhot(ihotmin2)

c       Check whether 2*smooth_i+1>dl
c       Check whether 2*smooth_i+1>distmin
c
	i_brute_force=0
	if(dl>2.0*smooth(i)) then
	   i_brute_force(1)=1
	endif
	if(distmin>2.0*smooth(i)) then
	   i_brute_force(2)=1
	endif

	if(i_brute_force(1)/=0.or.i_brute_force(2)/=0) then 

c       Find the nearest (hot) particle to the gridpoint (Brute force) 
c
	distmin=99999999
	do ihot=1,nhot
	   i=indxhot(ihot)
	   rig2=0.0
	   do jj=1,3
	      a=po(i,jj)-posg(jj)
	      rig2=rig2+a**2
	   enddo
	   dist=sqrt(rig2)
	   if(dist<distmin) then
	      distmin=dist
	      ihotmin2=ihot 
	   endif
	enddo
	i=indxhot(ihotmin2)
		  write(27,*)real(posg),i_brute_force
		  write(27,*)ihotmin2,indxhot(ihotmin2),real(distmin),
     1     real(2.0*smooth(i)),real(rho(i)),real(ener(i))


	endif

	ihotmin=ihotmin2

	smoothnew=smooth(indxhot(ihotmin))

         wsumv(:)=0.0d0
         wsumgasmass=0.0d0
         wsumgasener=0.0d0
         wsumgasrho=0.0d0
         wsumgassmooth=0.0d0

	do j=1,nnabo(ihotmin)
	   pp=indxhot(n_list(j,ihotmin))
	   sum=0
	   do jj=1,3
	      a=po(pp,jj)-posg(jj)
	      sum=sum+a**2
	   enddo
	   dist=sqrt(sum)
          wsumgasmass=wsumgasmass+mass(pp)*windowheat(dist,
     1              smoothnew)
          wsumv(:)=wsumv(:)+mass(pp)*windowheat(dist,
     1              smoothnew)*vel(pp,:)
          wsumgasener=wsumgasener+mass(pp)*windowheat(dist,
     1              smoothnew)*ener(pp)	   
          wsumgasrho=wsumgasrho+mass(pp)*windowheat(dist,
     1              smoothnew)*rho(pp)
          wsumgassmooth=wsumgassmooth+mass(pp)*windowheat(dist,
     1              smoothnew)*smooth(pp)
	enddo

	   
	energrid(gridx,2*gridy,2*gridz+1)=wsumgasener/wsumgasmass
	rhogrid(gridx,2*gridy,2*gridz+1)=wsumgasrho/wsumgasmass
	rhofac=wsumgasrho/wsumgasmass/rhomax
	if(rhofac>1.0) then
	   rhogrid(gridx,2*gridy,2*gridz+1)=rhomax
	   energrid(gridx,2*gridy,2*gridz+1)=
     1       energrid(gridx,2*gridy,2*gridz+1)*rhofac
	endif
	smogrid(gridx,2*gridy,2*gridz+1)=wsumgassmooth/wsumgasmass
	velgrid(gridx,2*gridy,2*gridz+1,:)=wsumv(:)/wsumgasmass

	dist=sqrt(posg(1)**2+posg(2)**2+posg(3)**2)
	if(dist>=rmaxg) stop
	ir=int(dist/drg)
	distir=ir*drg
	del=(dist-distir)/drg

	gravacc=gr(ir)+(gr(ir+1)-gr(ir))*del

	grgrid(gridx,2*gridy,2*gridz+1,:)=-gravacc*posg(:)/dist


	      enddo
	   enddo
	enddo
c
c  needs to come before explorer output (for pencil code)
c
	dlr4=real(dl)

	goto 747 !Skip Explorer style writeout

c----------------------------------------------------------------
c	Test writeout - Explorer style
c-----------------------------------------------------------------

	nbodies=N**3
	write(ubodsout)1
	write(ubodsout)nbodies,ndim,real(0)

	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=1000.0
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)

	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=energrid(gridx,gridy,gridz)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=posgrid(gridx,gridy,gridz,1)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=posgrid(gridx,gridy,gridz,2)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=posgrid(gridx,gridy,gridz,3)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=velgrid(gridx,gridy,gridz,1)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=velgrid(gridx,gridy,gridz,2)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=velgrid(gridx,gridy,gridz,3)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)

	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=1000.0
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=rhogrid(gridx,gridy,gridz)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)
	
	n=0
	do gridx=-N2,N2-1
	   do gridy=-N2,N2-1
	      do gridz=-N2,N2-1
		 n=n+1
		 readout(n)=smogrid(gridx,gridy,gridz)
	      enddo
	   enddo
	enddo

	write(ubodsout)(readout(p),p=1,nbodies)

	do p=1,nbodies
	   write(28,*)p,readout(p)
	enddo


 747	Continue

c----------------------------------------------------------------
c	Output for the Pencil Code
c-----------------------------------------------------------------

	rhogrid=rhogrid/grav  ! [Msun/kpc**3]

	sgrid=log(gam*gam1*energrid/rhogrid**gam1)/gam

	do gridx=-N2,N2-1
	   larr(gridx)=(real(gridx)+real(0.5))*dlr4
	enddo

        write(99) velgrid,alog(rhogrid),sgrid,grgrid
        write(99) tnowr,larr,larr,larr,dlr4,dlr4,dlr4
      print*,'t=',tnowr
      print*,'dlr4=',dlr4
      print*,'larr=',larr


 991	format(i4,8f9.2,i4)
 992	format(i4,4f9.2)
 993	format(i4,4f9.2,2f9.3,i6)
 994	format(i4,2f9.4,i6)	
 995	format(2f7.2,8f7.3,i6)
 996	format(f8.3,3f9.2,5f9.6)
 997	format(i4,f8.2,i5)
 998	format(i4,4f9.2,i6)
 999	format(3f9.2,i6)

	enddo

	close(ubodsin)
	close(ubodsout)
	end        
c
c
c
	subroutine conuffe(time,rhol,ul,T)
c
c     Conversion from specific energy to temperature, when the
c     time (i.e. UV-flux) and density are known.
c
	implicit real*8 (a-h,o-z)
	real*8 mymin,mymax
c
	common/temp/temp(0:12,-8:10,-2:10)
c
        mymin=.588
	mymax=1.22
c
	if(ul.ge.5.0) then
	T=log10(80.8*mymin*10**ul)
	return
	end if
c
	if(rhol.le.-4.0) then
	T=log10(80.8*mymin*10**ul)
	return
	end if
c
	minul=(2.*ul)
	maxul=minul+1
c
	minrhol=(2.*rhol)
	maxrhol=minrhol+1
c
	mintime=(2*time)
	maxtime=mintime+1
c
	a11=(2.*ul-real(minul))/(real(maxul)-real(minul))*
     1  (temp(mintime,minrhol,maxul)-temp(mintime,minrhol,minul))+
     2  temp(mintime,minrhol,minul)
	a12=(2.*ul-real(minul))/(real(maxul)-real(minul))*
     1  (temp(mintime,maxrhol,maxul)-temp(mintime,maxrhol,minul))+
     2  temp(mintime,maxrhol,minul)
	a21=(2.*ul-real(minul))/(real(maxul)-real(minul))*
     1  (temp(maxtime,minrhol,maxul)-temp(maxtime,minrhol,minul))+
     2  temp(maxtime,minrhol,minul)
	a22=(2.*ul-real(minul))/(real(maxul)-real(minul))*
     1  (temp(maxtime,maxrhol,maxul)-temp(maxtime,maxrhol,minul))+
     2  temp(maxtime,maxrhol,minul)
c
	a1=(2.*rhol-real(minrhol))/(real(maxrhol)-real(minrhol))*
     1  (a12-a11)+a11
	a2=(2.*rhol-real(minrhol))/(real(maxrhol)-real(minrhol))*
     1  (a22-a21)+a21
c
	T=(2.*time-real(mintime))/(real(maxtime)-real(mintime))*
     1  (a2-a1)+a1
c
	return
	end





      FUNCTION cf_metal(lT, Z_O)
      IMPLICIT NONE
c********************************************
c 
c The "differential metal cooling function"
c===================================================

      REAL*8 :: cf_metal       ! Result
      REAL*8 :: lT
      REAL*8 :: Z_O

	REAL*8 za1,za2,za,frac,frac1,X,oh,Z,ten
	REAL*8 lT_tab(89),cf_metal_tab(89,7),oh_tab(7)

      INTEGER i,j

      COMMON/tabs/lT_tab,cf_metal_tab,oh_tab

	ten=10.0d0

      if(lT<=lT_tab(1).or.lT>=lT_tab(89)) then
         cf_metal = -99.0
         RETURN
      endif

      if(Z_O<=0.0d0) then

         cf_metal = -99.0
         RETURN
      endif

        Z=2.0*Z_O  ! approximation

	X=0.76-(Z/0.02)*0.06 ! approximation
	if(X<0.1) X=0.1

      oh=12.0d0+log10(Z_O/X/16.0d0)-8.87d0

      i=INT(lT/0.05d0)-81

      frac=(lT-lT_tab(i))/(lT_tab(i+1)-lT_tab(i))

      if(oh<oh_tab(1)) then

         frac1=ten**oh/ten**oh_tab(1)

         za1=ten**cf_metal_tab(i,1)
         za2=ten**cf_metal_tab(i+1,1)

         za=za1+frac*(za2-za1)

         za=frac1*za

         cf_metal = log10(za+1.0d-99)
         RETURN
         
      endif

      do j=2,7
         if(oh<oh_tab(j)) goto 1
      enddo
      
      j=7
         
      goto 2

 1    j=j-1

 2    za1=cf_metal_tab(i,j)+frac*(cf_metal_tab(i+1,j)-cf_metal_tab(i,j))

      if(j<7) then

         frac1=(oh-oh_tab(j))/(oh_tab(j+1)-oh_tab(j))

         za2=cf_metal_tab(i,j+1)+frac*(cf_metal_tab(i+1,j+1)-
     1     cf_metal_tab(i,j+1))

         za=za1+frac1*(za2-za1)

      else

         za=za1

      endif

      cf_metal = za

      RETURN
      END



      FUNCTION windowheat(dist,smoothl)
      IMPLICIT NONE
c***********************************************************************
c     calculates the value of the smoothing kernel at distance dist
c     from the center for smoothinglength smoothl
c     
c=======================================================================

      INCLUDE 'windows.inc'

      REAL*8 :: windowheat 
      
      INTEGER ix
      REAL*8 :: dist, smoothl, x, y

      x = MIN(dist/smoothl,xm)
      y = x*rdwx
      ix = INT(y)
      windowheat = 2.0*(win(ix)+wind(ix)*(y-ix))/smoothl**3

      END


      SUBROUTINE window_setup
      IMPLICIT NONE
c***********************************************************888
c     setup the four arrays used for calculating the smoothing kernel
c     and its derivative by interpolation rather than directly.
c     

c     NOTE NOTE: CURRENT VERSION USES THE WINDOWS GIVEN BELOW C
C     WINDOW: SPLINE FUNCTION APPROXIMATING GAUSSIAN (HK STYLE)

C     GRAD_WINDOW: THOMAS AND COUCHMAN IMPROVEMENT OF GRAD OF ABOVE
*     SPLINE

c===============================================================

      INCLUDE 'windows.inc'

      REAL*8 x,win_gaus,grad_win_gaus
      INTEGER j

      DO j=0,nwin
         x=dwx*j
         win(j) = win_gaus(x)
      END DO
      win(nwin) = 0

c     multiply by 0.5 in preparation for HK symmetrization

      DO j=0,nwin
         win(j) = win(j)/2
      END DO

      DO j=0,nwin-1
         wind(j) = win(j+1)-win(j)
      END DO
      wind(nwin) = 0

c     now setup grad_window:

      DO j=0,nwin
         x=dwx*j
c        gwin(j) = grad_win_thomas(x)
         gwin(j) = grad_win_gaus(x)
      END DO
      gwin(nwin) = 0

c     Multiply by 0.5 for HK symmetrization

      DO j=0,nwin
         gwin(j)=gwin(j)/2
      END DO

      DO j=0,nwin-1
         gwind(j)=gwin(j+1)-gwin(j)
      END DO
      gwind(nwin) = 0

      END


      FUNCTION win_gaus(x)
      IMPLICIT NONE
c****************************************
c     this gives the HK (and other) spline function approximation
c     to a gaussian. By having compact support it is superior to
c     using a real gaussian.

      REAL*8, PARAMETER :: pi = 3.14159265358979323846
      
      REAL*8 :: win_gaus
      REAL*8 :: x
      
      IF (x.gt.2) THEN
         win_gaus = 0
         RETURN
      END IF

      IF (x.lt.1) THEN
         win_gaus=4-6*x**2+3*x**3
      ELSE
         win_gaus=(2-x)**3
      END IF

      win_gaus=win_gaus/(4*pi)

      END

      FUNCTION grad_win_gaus(x)
      IMPLICIT NONE
c****************************************
c     this gives the gradient of the usual gaussian approximating
c     spline function used by e.g. HK

      REAL*8, PARAMETER :: pi = 3.14159265358979323846
      
      REAL*8 :: grad_win_gaus
      REAL*8 :: x

      IF (x.gt.2) THEN
         grad_win_gaus = 0
         RETURN
      END IF

      IF (x.lt.1) THEN
         grad_win_gaus=3*x*(4-3*x)
      ELSE
         grad_win_gaus=3*(2-x)**2
      END IF

      grad_win_gaus = -grad_win_gaus/(4*pi)

      END




c	if(dl>2.0*smooth(i)) then
c	   print*,'dl > 2.0*smooth'
c	   i_brute_force(3)=1
c	endif
c	if(distmin>2.0*smooth(i)) then
c	   print*,'distmin > 2.0*smooth'
c	   i_brute_force(4)=1
c	endif

c	if(i_brute_force(3)/=0.or.i_brute_force(4)/=0) then 
c	   print*,i_brute_force
c	endif



