C ______________________________ Profiles for 'Cooled' Simulations
C		             
C		                     
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMcool.h'
      REAL         AMprevious(Np),AMlast(Np),Daverg(0:500)
      Real*8       Xc,Yc,Zc
      DIMENSION    Iaverg(0:500)
      REAL         INPUT
      Character    FileASCII*100
      Character  File1*100,File2*100
      EQUIVALENCE (Ndisk,extras(90)),
     &               (Nbulge,extras(97)),(aMassDisk,extras(93)),
     &               (Box,extras(100))


c      aMassDisk =3.36e10/0.7     ! MOVIE mass of disk in real Msun
c      aMassDisk =1.e10    ! "Bulge"
c      aMassDisk =5.71e8      ! "Dwarf"
c      aMassDisk =4.8e10      ! "MW_LARGE"

C...................................................................
C			Read data and open files
       write(*,*)' enter 4 digits of the snapshot'
       read(*,*) moment
         write(File1,'("FILES/PMcrda.",i4.4,".DAT")')moment
         write(File2,'("FILES/PMcrs0a.",i4.4,".DAT")')moment
         write(*,*) File1
         write(*,*) File2

c      File1 ='PMcrd.DAT'
c      File2 ='PMcrs0.DAT'


      CALL RDTAPE(File1,File2)
         write(FileASCII,'("RESULTS/Result.DM.",i4.4,".dat")')
     &    INT(AEXPN*10000.)
      OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') 
 
        write(FileASCII,'("RESULTS/Result.SPH.",i4.4,".dat")')
     &    INT(AEXPN*10000.)
      OPEN(39,FILE=FileASCII,STATUS='UNKNOWN') 

         write(FileASCII,'("RESULTS/Result.CYL.",i4.4,".dat")')
     &    INT(AEXPN*10000.)
      OPEN(41,FILE=FileASCII,STATUS='UNKNOWN') 

      OPEN(18,FILE='Stat.dat',STATUS='UNKNOWN',access='append')
      OPEN(20,FILE='particles.dat',STATUS='UNKNOWN')
      OPEN(30,FILE='particlesu.dat',STATUS='UNKNOWN',form='UNFORMATTED')
      OPEN(40,FILE='FILES/cooling_part.dat',STATUS='UNKNOWN',
     &              form='UNFORMATTED')

      CALL ReadCool
      time      = (age(AEXPN)-age(0.6))/1.e9
      
      aMassOne= Om0*(Box/NGRID)**3 *wspecies(1)*2.746e+11/hubble ! mass scale
      write (*,'(6x,"Mass of smallest particle",/9x,
     +    "in units M_sun/h is   =",3x,g10.3)') 
     +      aMassOne
      ScaleM = aMassOne

      xfrac        = extras(51) ! fraction of particles to be cooled
      cool_frac    = extras(52) ! total cooling
      cfrac        = extras(53) ! cooling per step
      aexpn_s      = extras(54) ! start cooling
      aexpn_e      = extras(55) ! and cooling
      xm_part      = extras(59) ! particle weight (h^-1 M_s)
      i_cool_all    = INT(extras(58) + 0.1)  !total # of cooling steps
      

      write (*,*) ' RDTAPE is done Time=',time,' a=',AEXPN
      write (*,*)  ' mass particle(extr(59))=',xm_part
      Do i=37,41,2
          WRITE (i,100) HEADER,AEXPN,ASTEP,time,
     +                  ISTEP,
     +                  NROWC,NGRID,Om0,Oml0,hubble,
     +                  Ocurv
      write (i,104) xfrac
      write (i,105) cool_frac 
      write (i,106) cfrac 
      write (i,107) i_cool_all
      write (i,108) aexpn_s
      write (i,109) aexpn_e
      a_cool= cfrac
      b_cool= cfrac
      c_cool= 0.1*cfrac
      write (i,110) a_cool
      write (i,111) b_cool
      write (i,112) c_cool

      EndDo 
100   FORMAT(1X,'Header=>',A45,/
     +    T2,'Aexp',T12,'dA',T22,'Time(Gyrs)',T35,'Step',/
     &    T2,f8.4,T12,f8.5,T22,g11.3,T35,i4,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3)
 104  FORMAT ( ' frac. cooled particles  = ',F6.3 )
 105  FORMAT ( ' total amount of cooling = ',F6.3 )
 106  FORMAT ( ' cooling per step        = ',F6.3 )
 107  FORMAT ( ' total # of cooling steps= ',I6)
 108  FORMAT ( ' start cooling aexpn     = ',F6.3 )
 109  FORMAT ( ' end   cooling aexpn     = ',F6.3 )
 110  FORMAT ( ' a_cool                  = ',F6.3 )
 111  FORMAT ( ' b_cool                  = ',F6.3 )
 112  FORMAT ( ' c_cool                  = ',F6.3 )
      BoxV   = Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID      ! scale factor for Velocities (km/s)
      ScaleC = Box/NGRID*AEXPN/hubble  ! scale factor for Coordinates (real Mpc)
      ScaleH  = 100.*hubble*sqrt(Om0/AEXPN**3+Oml0) ! H(a)=km/s/Mpc

      i_halo = 1   ! read and analyze this halo
      kkF    = iph_cool(i_halo,2)  ! first particle
      kkL    = iph_cool(i_halo,3)  ! last particle
      CALL ReadPnts(N,ScaleC,ScaleV,ScaleM)
      Call ChangePnts(N,Box,ScaleH,Ndisk)
      Call DumpPnts(N,Box,ScaleM,Ndisk)
c      Call Statist(N,Box,ScaleM,Ndisk)
c      Call AngMom(N,30000,3.,Box,ScaleM)
      Call Profile(N,Box,ScaleM,Ndisk)
c      Call Projection(N,Box,0.)
c      Call Projection(N,Box,5.)

      END