C ______________________________ Profiles for Galaxy Simulations
C		             
C		                     
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMgalaxy.h'
      REAL         INPUT,moment
      Character    FileASCII*100,Analysis*3
      Character  File1*100,File2*100
      EQUIVALENCE (Ndisk,extras(90)),
     &               (Nbulge,extras(97)),(aMassDisk,extras(93)),
     &               (Box,extras(100))


      
c      aMassDisk =1.e10    ! "Bulge"
c      aMassDisk =5.71e8      ! "Dwarf"
c      aMassDisk =3.36e10/h      ! "MW_LARGE"
      Analysis ='AAB'

C...................................................................
C			Read data and open files
       write(*,*)' enter 6 digits of the snapshot Example: 0.6235'
       read(*,*) moment
       If(moment.lt.1.5)Then
         write(File1,'("FILES/PMcrda",f6.4,".DAT")')moment
         write(File2,'("FILES/PMcrs0a",f6.4,".DAT")')moment
       Else
          imoment =INT(moment+0.1)
c       read(*,*) imoment
         write(File1,'("FILES/PMcrda.",i4.4,".DAT")')imoment
         write(File2,'("FILES/PMcrs0a.",i4.4,".DAT")')imoment
       endIf
         write(*,*) File1
         write(*,*) File2

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

      CALL BLOCK_DATA
      CALL RDTAPE(File1,File2)

c      Nbulge = 33409
c      aMassDisk =4.28e10*hubble      ! MOVIE
c      Ndisk     =199996
      time = Age(AEXPN)-Age(0.60)
         write(FileASCII,'("RESULTS/Result",a,".DM.",f6.4,".dat")')
     &    Analysis,AEXPN
      OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') 
        write(FileASCII,'("RESULTS/Result",a,".SPH.",f6.4,".dat")')
     &    Analysis,AEXPN
      OPEN(39,FILE=FileASCII,STATUS='UNKNOWN') 

         write(FileASCII,'("RESULTS/Result",a,".CYL.",f6.4,".dat")')
     &    Analysis,AEXPN
      OPEN(41,FILE=FileASCII,STATUS='UNKNOWN') 

         write(FileASCII,'("RESULTS/Result",a,".BAR.",f6.4,".dat")')
     &    Analysis,AEXPN
      OPEN(43,FILE=FileASCII,STATUS='UNKNOWN') 
      OPEN(18,FILE='Stat.dat',STATUS='UNKNOWN',access='append')
         write(FileASCII,'("RESULTS/Result",a,".PRT.",f6.4,".dat")')
     &    Analysis,AEXPN
      OPEN(20,FILE=FileASCII,STATUS='UNKNOWN')
c      OPEN(30,FILE='particlesu.dat',STATUS='UNKNOWN',form='UNFORMATTED')

 
      time      = (age(AEXPN)-age(0.6))
      aMassDisk = aMassDisk/hubble
c      aMassDisk = 3.5e10/hubble

      write (*,*)  ' hubble =',hubble,' mass of disk=',aMassDisk
      write (*,*) ' RDTAPE is done Time=',time,' a=',AEXPN
          WRITE (*,100) HEADER,AEXPN,ASTEP,time,
     +                  ISTEP,Ndisk,Nbulge,aMassDisk,
     +                  NROWC,NGRID,bar_max*1.e3,dZ_slice*1.e3,
     +                  Om0,Oml0,hubble,Ocurv
      Do i=37,43,2
          WRITE (i,100) HEADER,AEXPN,ASTEP,time,
     +                  ISTEP,Ndisk,Nbulge,aMassDisk,
     +                  NROWC,NGRID,bar_max*1.e3,dZ_slice*1.e3,
     +                  Om0,Oml0,hubble,Ocurv
      EndDo 
100   FORMAT(1X,'Header=>',A45,/
     +    T2,'Aexp',T12,'dA',T22,'Time(Gyrs)',T35,'Step',
     &    T42,'Ndisk',T50,'Nbulge',T58'Mdisk(Msun)',/
     &    T1,f8.5,T10,f8.5,T22,g11.4,T35,i4,T40,i8,T48,i7,T58,g11.4,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Max Bar radius(kpc)=',
     &             f7.3,' Slice in Z for cylindr shells=',f8.3/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3)
c      Box  =INPUT(' Enter box size in comoving Mpc/h =')
c      Box =0.5 ! NG3109
c      Box = 2   ! 'bulge'

      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
      if(Ndisk.le.0)write (*,*) ' wrong number of disk particles:',Ndisk
      if(Ndisk.le.0)Stop
      ScaleM = aMassDisk/Ndisk
      Ndisk = Ndisk +Nbulge
      Nbulge = 0

      CALL ReadPnts(N,ScaleC,ScaleV,ScaleM)
      Call ChangePnts(N,Box,ScaleH)
c      Call Acceleration(N,Box)
c      Call DumpPnts(N,Box,ScaleM)

c      Call AngMom(N,30000,3.,Box,ScaleM)
c      CALL PV(N,Box,ScaleM)
c      CALL RandV(N,Box,ScaleM)
      Call Profile(N,Box,ScaleM)
      Call Statist(N,Box,ScaleM)
c      Call Projection(N,Box,0.)
c      Call Projection(N,Box,5.)

      END