C ______________________________ Convert PM from one ng to another ng C A. Klypin (aklypin@nmsu.edu) C C Scales internal PM coordinates and velocities C Produces new set of basic files C NROW = number of particles in 1D C NGRID= number of cells in 1D INCLUDE 'PMparameters.h' REAL INPUT logical inside Character FileASCII*50 C................................................................... C Read data and open files c WRITE (*,'(A,$)') 'Enter Name for output file = ' c READ (*,'(A)') FileASCII c OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') CALL RDTAPE Close (9) OPEN(UNIT=9,FILE='PMcrd.NEW.DAT', + FORM='UNFORMATTED', STATUS='UNKNOWN') OPEN(UNIT=22,FILE='PMcrs0.NEW.DAT',ACCESS='DIRECT', + STATUS='UNKNOWN',RECL=NACCES) write (*,*) ' RDTAPE is done' Scale =INPUT(' Enter scaling factor for new box/old box=') ! rescale NGRID; make scaling factors If(Scale.gt.1.)Then Scale = INT(Scale+0.1) iScale = INT(Scale+0.1) if((iScale.ne.2).and.(iScale.ne.4).and.(iScale.ne.8)) & STOP ' wrong scaling factor' NGRIDC = NGRID*iScale Else Scale = 1./Scale Scale = INT(Scale+0.1) iScale = INT(Scale+0.1) if((iScale.ne.2).and.(iScale.ne.4).and.(iScale.ne.8)) & STOP ' wrong scaling factor' Scale =1./Scale NGRIDC = NGRID/iScale EndIf NROWC = NROW ! do not change NROW ! rescale weights PARTW = PARTW *Scale**3 Do j=1,Nspecies wspecies(j) =wspecies(j)*Scale**3 EndDo CALL WRTAPE ! save control files write (*,*) ' Scaling =',Scale,iscale WRITE (*,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRIDC,NRECL,Om0,Oml0,hubble, + Ocurv 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3) C Loop over particles Icount = 0 DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last write (*,*)' Read page=',IROW,' N=',In_page CALL GETROW(IROW,1) DO IN=1,In_page XPAR(IN) =Scale* (XPAR(IN)-1.)+1. YPAR(IN) =Scale* (YPAR(IN)-1.)+1. ZPAR(IN) =Scale* (ZPAR(IN)-1.)+1. If(XPAR(IN).ge.NGRIDC+1)write(*,*)' ErrorX',XPAR(IN),IN If(YPAR(IN).ge.NGRIDC+1)write(*,*)' ErrorY',YPAR(IN),IN If(ZPAR(IN).ge.NGRIDC+1)write(*,*)' ErrorZ',ZPAR(IN),IN If(XPAR(IN).lt.1.)write(*,*)' ErrorX',XPAR(IN),IN If(YPAR(IN).lt.1.)write(*,*)' ErrorY',YPAR(IN),IN If(ZPAR(IN).lt.1.)write(*,*)' ErrorZ',ZPAR(IN),IN Vxs=Scale* VX(IN) Vys=Scale* VY(IN) Vzs=Scale* VZ(IN) Icount =Icount +1 xmax =MAX(xmax,XPAR(IN),YPAR(IN),ZPAR(IN)) xmin =MIN(xmin,XPAR(IN),YPAR(IN),ZPAR(IN)) ENDDO CALL WRIROW(IROW,2) ENDDO write (*,*) ' Scaled Coordinates were written' write (*,*) ' Number of particles =',Icount write (*,*) ' Min/Max of coordinates(Mpc/h)=',xmin,xmax END C--------------------------------------------------- C Read current data from disk/tape, C Open files C Nrecl is the number of values in a record C Npage is the number of particles in a page SUBROUTINE RDTAPE C--------------------------------------------------- INCLUDE 'PMparameters.h' C Open file on a tape OPEN(UNIT=9,FILE='PMcrd.DAT', + FORM='UNFORMATTED', STATUS='UNKNOWN') C Read control information C and see whether it has proper C structure READ (9,err=20,end=20) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, + NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 + ,Ocurv,extras WRITE (*,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN,EKIN1,EKIN2, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv WRITE (16,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN,EKIN1,EKIN2, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3) IF(NROW.NE.NROWC) THEN WRITE (*,*) + ' NROW in PARAMETER and in TAPE-FILE are different' ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) + ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF C Open work file on disk write (*,*) ' start openning PMcrs0' 10 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT', + STATUS='UNKNOWN',RECL=NACCES) REWIND 9 write (*,*) ' done openning PMcrs0' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END C-------------------------------------------------- C Write current PAGE of particles (x,v) to disk C NRECL - length of ROW block in words SUBROUTINE WRIROW(IROW,Ifile) C-------------------------------------------------- INCLUDE 'PMparameters.h' WRITE (20+Ifile,REC=IROW) RECDAT RETURN END C--------------------------------------------- C Write current data to disk/tape C SUBROUTINE WRTAPE C---------------------------------------------- INCLUDE 'PMparameters.h' C write header and control data WRITE (9) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, + NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5, + Ocurv,extras REWIND 9 RETURN END C-------------------------------------------- C Read current PAGE of particles from disk C NRECL - length of ROW block in words SUBROUTINE GETROW(IROW,Ifile) C-------------------------------------------- INCLUDE 'PMparameters.h' READ (20+Ifile,REC=IROW) RECDAT RETURN END C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End