! _______________________ START 3-D SIMULATIONS ! ! Klypin, August 1993 ! !------------------------------------------------------------- MODULE Randoma Contains !-------------------------------------- ! normal random numbers FUNCTION GAUSS (M) !-------------------------------------- X = 0. DO I = 1, 5 X = X + RANDd (M) EndDo X2 = 1.5491933 * (X - 2.5) GAUSS = X2 * (1. - 0.01 * (3. - X2 * X2) ) RETURN END FUNCTION GAUSS !------------------------------------------------ ! random number generator FUNCTION RANDd (M) !------------------------------------------------ DATA LC, AM, KI, K1, K2, K3, K4, L1, L2, L3, L4 / 453815927, & 2147483648., 2147483647, 536870912, 131072, 256, 16777216, 4, & 16384, 8388608, 128 / ML = M / K1 * K1 M1 = (M - ML) * L1 ML = M / K2 * K2 M2 = (M - ML) * L2 ML = M / K3 * K3 M3 = (M - ML) * L3 ML = M / K4 * K4 M4 = (M - ML) * L4 M5 = KI - M IF (M1.GE.M5) M1 = M1 - KI - 1 ML = M + M1 M5 = KI - ML IF (M2.GE.M5) M2 = M2 - KI - 1 ML = ML + M2 M5 = KI - ML IF (M3.GE.M5) M3 = M3 - KI - 1 ML = ML + M3 M5 = KI - ML IF (M4.GE.M5) M4 = M4 - KI - 1 ML = ML + M4 M5 = KI - ML IF (LC.GE.M5) ML = ML - KI - 1 M = ML + LC RANDd = M / AM RETURN END FUNCTION RANDd end module Randoma ! ! !------------------------------------------------------------------ MODULE setInitialConditions PUBLIC Character (LEN=12), DIMENSION(2) :: NamesDir=(/ & 'PMstartX.DAT',& 'PMstartY.DAT'/) INCLUDE 'PMinitialp.h' INTEGER, PARAMETER :: NPAGE = NROW**2, & ! # particles in a record NMAX = NGRID/2, & NRECL = NPAGE*6, & ! # particles in a record NARR = MAX(2*NROW+1,NGRID+1), & ! FFT NF67 = MAX(NROW,NGRID/2), & Nblocks=NROW/Lblock, & NSPEC = NROW/2, & ! No waves shorter than Ny marr = NROW+1, & ! Arrays for FFT mf67 = NROW/2 REAL, DIMENSION(10) :: timing =0. REAL :: alpha =0. COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras(100) COMMON / HEADDR/ HEADER CHARACTER*45 HEADER COMMON /FOURAR/Zf(NARR),Yf(NARR) COMMON /F67COM/ & IBC, IP, ISL, L1, N2, & N3, N4, N7, & SI(NF67), INDEX(NF67) COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE),& VX(NPAGE),VY(NPAGE),VZ(NPAGE) DIMENSION RECDAT(NRECL),wspecies(10),lspecies(10) EQUIVALENCE (RECDAT(1),XPAR(1)) & ,(wspecies(1),extras(1)), & (lspecies(1),extras(11)) COMMON / MASK/ iMask(Nblocks,Nblocks,Nblocks) ! Refinement Mask REAL, ALLOCATABLE, DIMENSION (:,:,:) :: GR COMMON / BUFFx / XPt(Nmaxpart) COMMON / BUFFv/ VXt(Nmaxpart) COMMON /SeedsPage/ NseedP(NROW) COMMON / LevelsM/ Max_lev_abs,Nmin_lev,Nmax_lev Contains !-------------------------------------------------- ! : Store and retrieve seeds for parallelization ! of random number generator !-------------------------------------------------- SUBROUTINE SetRandom(iCont) use Randoma INCLUDE 'luxuryp.h' Character (LEN=70) :: FileName Logical :: FileExists =.false., RNGluxury =.false. Logical :: iCont REAL :: t0=0,t1=0. ! counters for timing CALL CPU_TIME(t0) ! initialize time counter write(FileName,'(a,i1.1,".",i4.4,".",i8.8,".dat")') & 'Seeds.',INT(extras(80)),NROW,Nseed INQUIRE(file=FileName,EXIST=FileExists) If(abs(extras(80)-1.) < 1.e-5)RNGluxury =.true. If(FileExists)Then ! read the data write(*,'(10x,2a)')'Open Existing FileName=',FileName Open(62,file=FileName,form='unformatted') read(62)NROWdummy,Nseeddummy if(NROWdummy /= NROW .or. Nseeddummy /= Nseed)Then write(*,*) ' File with seeds does not match parameters' write(*,*) ' Read Nrow =',NROWdummy,' Internal=',NROW write(*,*) ' Read Nseed=',Nseeddummy,' Internal=',Nseed STOP endIf If(RNGluxury)Then ! read(62)i24A,j24A,in24A,kountA,carryA,gsetA,iFlagA read(62)SeedsPage Close(62) Else read(62)i24A EndIf Close(62) Else ! open file with seeds write(*,'(10x,2a)')'Create FileName=',FileName Open(62,file=FileName,form='unformatted') ! generate seeds If(RNGluxury)Then ! Detect random generator: luxury or nrandd iFlag = 0 gSet = 0. Do k=1,NROW ! luxury if(k/10*10==k)write(*,*) ' page number=',k i24A(k) = i24 j24A(k) = j24 in24A(k) = in24 kountA(k) = kount carryA(k) = carry gsetA(k) = gSet iFlagA(k) = iFlag Do m=1,24 SeedsPage(m,k) = seeds(m) EndDo Do j=1,NROW Do i=1,NROW ! If(i+j+k.ne.3)Then x = GAUSS3(gSet,iFlag) ! EndIf EndDo EndDo EndDo write(62)NROW,Nseed write(62)i24A,j24A,in24A,kountA,carryA,gsetA,iFlagA write(62)SeedsPage Close(62) Else write(*,*) ' Nseed=',Nseed NRAND = Nseed Do k=1,NROW ! gauss and randd if(k/10*10==k)write(*,*) ' page number=',k i24A(k) = NRAND Do j=1,NROW Do i=1,NROW If(i+j+k.ne.3)Then x = GAUSS(NRAND) EndIf EndDo EndDo EndDo write(62)NROW,Nseed write(62)i24A Close(62) EndIf ! end luxury/randd switch write(*,'(/2(20("-"),a/))')'File with seeds was created. STOP', & 'Do not change anything. Just restart the code to produce Initial Conditions' iCont = .false. EndIf ! file exists switch CALL CPU_TIME(t1) timing(2) = t1-t0 ! time for setting random numbers Return End SUBROUTINE SetRandom !-------------------------------------------------- ! sqrt(Power spectrum) ! k = (2pi/L) !-------------------------------------------------- FUNCTION TRUNF (WK) !------------------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE Real ns, k IF (WK.GE.FLOAT (NSPEC) ) THEN TRUNF = 0. RETURN ENDIF k = QSCALE * wk ! use approximations If (iFlag.gt.0) Then ! Holtzman approx If (Par (6) .ne.0.) Then sk = sqrt (k) TRUNF = k** (ns / 2.) / (1. + sk * (Par (2) + sk * (Par (3) & + sk * (Par (4) + sk * (Par (5) ) ) ) ) ) ** (Par (6) ) ! BBKS + Sugiyama approx Else ! Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om)) ! Q = wk/hsmall/Gammaeff Q = k * qqscaleb TRUNF = k** (ns / 2.) * LOG (1. + Par (1) * Q) / (Par (1) * Q) & / sqrt (sqrt (1. + Q * (Par (2) + Q * (Par (3) **2 + Q * (Par (4) & **3 + Q * (Par (5) **4) ) ) ) ) ) EndIf ! use table for P(k) Else TRUNF = sqrt (Ppk (k) ) EndIf RETURN END FUNCTION TRUNF !--------------------------------------- FUNCTION Ppk (x) ! interpolate table with p(k) ! x is in real 1/Mpc !--------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag ! slope is ns =-3 If (x.ge.xkt (Ntab) ) Then Ppk = Pkt (Ntab) / (x / xkt (Ntab) ) **3 Return EndIf If (x.lt.xkt (1) ) Then ! slope is ns=1 Ppk = Pkt (1) * (x / xkt (1) ) Return EndIf ind = INT ( (log10 (x) - alog0) / StepK) + 1 dk = xkt (ind+1) - xkt (ind) Ppk = (Pkt (ind) * (xkt (ind+1) - x) + Pkt (ind+1) * (x - xkt ( & ind) ) ) / dk ! write (*,'(5x," k=",g12.4," table=",2g12.4," P=",3g12.4)') ! & x,xkt(ind),xkt(ind+1), ! & Ppk,Pkt(ind),Pkt(ind+1) Return END FUNCTION Ppk !--------------------------------------- SUBROUTINE ReadTable (hsmall) ! interpolate table with p(k) !--------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag CHARACTER(80) Line ! W.Hu table If (iFlag.eq. - 1) Then Open (50, file = 'WHUpk.dat') twop = 2. * 3.14145926**2 Ntab = 0 10 READ (50, *, end = 30, err = 30) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx * hsmall ! scale to real Mpc Pkt (Ntab) = pp / xx**3 * twop ! k^3Pk -> Pk GOTO 10 30 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from WHUpk.dat table' Write (*,*) ' hsmall =', hsmall Goto 50 EndIf ! Eis.Hu table If (iFlag.eq. - 2) Then ! Open (50, file = 'pk_EHu_Om0.3_Omb0.043_s8=0.9.dat') ! write(*,*) ' Use pk_EHu_Om0.3_Omb0.043_s8=0.9.dat ' Open (50, file = 'pk_EHu_WMAP06.dat') write(*,*) ' Use pk_EHu_WMAP06.dat ' Else !Open (50, file = 'pk_EHu_WMAP06.dat') Open (50, file = 'pk_cambWMAP.dat') write(*,*) ' Use pk_cambWMAP.dat' EndIf ! read header Do i = 1, 9 READ (50, '(a)', end = 100, err = 100) Line EndDo Ntab = 0 12 READ (50, *, end = 32, err = 32) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx * hsmall ! scale to real Mpc Pkt (Ntab) = pp ! Pk GOTO 12 32 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from P(k) table ' Write (*,*) ' hsmall =', hsmall 50 If (Ntab.le.1) stop 'wrong table for p(k)' StepK = log10 (xkt (2) / xkt (1) ) alog0 = log10 (xkt (1) ) ! test that spacing is the same Do k = 2, Ntab - 1 ss = log10 (xkt (k + 1) / xkt (k) ) If (abs (ss / StepK - 1.) .gt.2.e-2) Then Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k) STOP EndIf EndDo CLOSE (50) Return 100 Write (*,*) ' Error reading table of P(k)' STOP END SUBROUTINE ReadTable !********************************************************************* ! Power spectrum for CDM. Wk- wavenumber(Mpc^-1, h=0.5) ! P = wk *T^2(k) ! QSCALE = 2pi/Box_size, QS =hubble**(-2) ! FUNCTION TRUNF(WK) ! COMMON / TRUNCOM/ QSCALE,g,gy ! COMMON / TRUNCDM/ QS,A1,A2,A3,A4,A5,A32,A43,A54 ! Q = QSCALE*QS*WK ! TRUNF = sqrt( WK* ! + (LOG(1.+A1*Q)/(A1*Q))**2/ ! + SQRT(1.+Q*(A2 +Q*(A32 +Q*(A43+Q*A54)))) ) ! ! RETURN ! END !********************************************************************* ! INITIALIZE CONSTANTS: ! Scalel =length of square in MPC ! AEXPN =expansion factor (current) ! ASTEP =increment of expansion factor: AEXPN =AEXP_0+N*ASTEP ! PARTW =weight of a particle: total 'mass' per cell ! even for LCDM or OCDM ! AMPLT =amplitude of perturbations inside the simulation box ! NBYTE =record length in bytes ! SUBROUTINE InitValues (NBYTE, SCALEL) !-------------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE COMMON / FERMI / Vscale Real ns Character Answer*1 Real INPUT DATA PI / 3.1415926535 / Do i = 1, 100 extras (i) = 0. enddo Write (*,*) 'Would you like to use a model provided in cdm.fit ' Write (*,*) 'or your own model?' Write (*,'(A,$)') ' Enter Yes for cdm.fit; No for your model =' Read (*,'(A)') Answer HEADER = 'N=128x256 L=20h-1CDMsig=0.239 z=30-----------' Write (*,*) '------ Enter Header for the run up to 45 characters' Write (*,*) ' Example:' Write (*,'(A)') HEADER READ (*,'(A)') HEADER Write (*,'(A)') HEADER AEXPN = INPUT (' Initial expansion parameter (0-1)=') ASTEP = INPUT (' Step of the expansion parameter =') AMPLT = INPUT (' Amplitude of density fluctuations=') ns = INPUT (' Slope of the Power spectrum n=') SCALEL= INPUT (' Box size (Mpc/h) =') Ocurv = INPUT (' Omega_curvature at present =') Nseed = INPUT (' Random seed (Integer 1-2^31) =') If (Nseed.le.0) Then lux = 2 ! level for luxury extras (80) = 1. extras (81) = lux If (Nseed.eq.0) Nseed = 121071 Nseed = ABS (Nseed) CALL rluxgo (lux, Nseed, 0, 0) ! initialize luxury Else ! use randd and gauss extras (80) = 0. EndIf ! extras(80) is used as switch for ! random numbers generators ! = 0 - old Randd+Gauss ! = 1 - luxury+Gauss3 ! extras(81) = lux - level of luxury ! store the box size extras (100) = SCALEL If (Answer.eq.'Y'.or.Answer.eq.'y') Then CALL MODEL (hsmall) hubble = hsmall Om0 = Om Oml0 = 1. - Om0 - Ocurv Else Write (*,*) ' You use your cosmological model.' Write (*,*) ' Be sure you provide routine TRUNF' hubble = INPUT (' The Hubble constant (h=H/100) =') Om0 = INPUT (' Omega_0 matter at present =') Oml0 = INPUT (' Omega_lambda at present =') EndIf ! scale it to real megaparsecs SCALEL = SCALEL / hubble NBYTE = NPAGE * 6 * 4 Write (*, '(A,$)') ' Enter Min and Max number of mass species =' READ ( *, * ) Nmin_lev, Nmax_lev W = (FLOAT (NGRID) / FLOAT (NROW) ) **3 PARTW = W ! multiple mass resolution If (Nmin_lev.lt.Nmax_lev) Then Write ( *, * ) Write (*,*) ' You use multiple mass resolution' ENDIF ! change the order of levels If(Nmin_lev.gt.Nmax_lev)Then i = Nmax_lev Nmax_lev = Nmin_lev Nmin_lev = i endif ISTEP = 0 TINTG = 0. AU0 = 0. AEU0 = 0. EKIN = 0. EKIN1 = 0. EKIN2 = 0. NROWC = NROW NGRIDC= NGRID QSCALE = 2.*PI/SCALEL QS = hubble**(-2) write (*,'(/3x,3(a,g12.3))') ' Qscale=',QSCALE, & ' Box=',SCALEL,' slope=',ns RETURN END SUBROUTINE InitValues !________________________________________Read parameters of the model ! Om = Omega_0 ! Omb = Omega_baryon ! Omnu = Omega_neutrino ! hsmall = hubble constant ! Par = fitting parameters ! Par(6) = 0 --- bbks-style+Hu&Sugiyama ! Par(6) ne 0 --- holtzman-style SUBROUTINE MODEL (hsmall) !--------------------------------------- COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag Real ns, INPUT Character Header1*79 Line1 = INPUT (' Enter Line Number in cdm.fit =') If (Line1.gt.0) Then ! use analytical approximation iFlag = Line1 OPEN (2, file = 'cdm.fit') Read (2, '(A)') Header1 If (Line1.gt.2) Then Do i = 1, Line1 - 2 Read (2, * ) a EndDo EndIf Read (2, * ) Om, Omb, Omc, Omnu, hsmall, Par CLOSE (2) ! = T_cmb/2.7 theta = 2.726 / 2.7 ! Hu&Sugiyama fitting Ob0 = Omb / Om Omh2 = Om * hsmall**2 a1 = (46.9 * Omh2) **0.670 * (1. + (32.1 * Omh2) ** ( - 0.532) ) a2 = (12.0 * Omh2) **0.424 * (1. + (45. * Omh2) ** ( - 0.582) ) alph = 1. / a1**Ob0 / a2** (Ob0**3) qqscaleb = theta**2 / (1. - Ob0) **0.60 / sqrt (alph) / Omh2 Write (*,20) Om, Omb, Omc, Omnu, hsmall, Par ! Write (16,20)Om,Omb,Omc,Omnu,hsmall,Par 20 Format (' Model: Om_0=', F5.3, ' Om_baryon=', F5.3, ' Om_cold=', & F5.3, ' O_nu=', F5.3, ' hsmall=', F4.2, / 8x, 'Parameters=', & (6G9.3) ) Else ! read spectrum from table iFlag = Line1 If (iFlag.lt. - 3) Stop ' Error: model Line # less than -3' ! W.Hu table If (iFlag.eq. - 1) Then hsmall = 0.7 Omb = 0.024 / hsmall**2 Om = 0.3 Ocurv = 0. Omc = Om EndIf ! Eis.Hu approximation If (iFlag.eq. - 2) Then hsmall = 0.73 Omb = 0.04 Om = 0.24 Ocurv = 0. Omc = Om EndIf ! Eis.Hu approximation: WMAP If (iFlag.eq. - 3) Then hsmall = 0.70 Omb = 0.046 Om = 0.27 Ocurv = 0. Omc = Om EndIf CALL ReadTable (hsmall) Write ( *, 22) Om, Omb, hsmall Write (1, 22) Om, Omb, hsmall 22 Format (' Model: Om_0=', F5.3, ' Om_baryon=', F5.3, ' hsmall=', & F4.2, ' <== Using table for P(k)') EndIf Return END SUBROUTINE MODEL !------------------------------------------------ ! Clear MASK for multiple mass resolution ! Initiat to level Nmin_lev SUBROUTINE ZeroMask !------------------------------------------------ DO MK3 = 1, Nblocks DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iMask (MK1, MK2, MK3) = Nmin_lev ENDDO ENDDO ENDDO RETURN END SUBROUTINE ZeroMask !------------------------------------------------ ! Set MASK for multiple mass resolution ! Read particles with lagrangian coordinat ! to set the mask ! Block = 1 = low resolution ! = 2 = higher resolution ! = 3 = even higher resolution ! After the mask is set, make maximum resolution Nmax_lev ! SUBROUTINE SetMask !------------------------------------------------ Dimension nLevel (10) Character *70 Line Do i = 1, 5 Read (30, '(A)') Line Write (*, '(A)') Line EnDDo Nlines = 0 10 Read (30, *, end = 50, err = 50) i, xw, yw, zw, vwx, vwy, vwz, im,& jm, km Nlines = Nlines + 1 ! write (*,*) ' Block=',im,jm,km If (im.lt.0.or.im.gt.Nblocks) Then Write (*,*) ' wrong x_block: ', im, jm, km, ' Line=', Nlines Stop EndIf If (jm.lt.0.or.jm.gt.Nblocks) Then Write (*,*) ' wrong y_block: ', im, jm, km, ' Line=', Nlines Stop EndIf If (km.lt.0.or.km.gt.Nblocks) Then Write (*,*) ' wrong z_block: ', im, jm, km, ' Line=', Nlines Stop EndIf ! If(iMask(im,jm,km).ne.0)Then ! write (*,*) ' Occupied block: ',im,jm,km,' Line=',Nlines ! write (*,*) ' Block=',iMask(im,jm,km) ! Stop ! EndIf ! mark high resolution block iMask (im, jm, km) = Nmax_lev Do ik = - 1, 1 k = km + ik If (k.le.0) k = k + Nblocks If (k.gt.Nblocks) k = k - Nblocks Do ij = - 1, 1 j = jm + ij If (j.le.0) j = j + Nblocks If (j.gt.Nblocks) j = j - Nblocks Do ii = - 1, 1 i = im + ii If (i.le.0) i = i + Nblocks If (i.gt.Nblocks) i = i - Nblocks ind = abs (ik) + abs (ij) + abs (ii) ! only seven blocks ! are marked out ! If(iMask(i,j,k).le.0.or.iMask(i,j,k).gt.Nspecies) ! + Then ! Stop ! Endif If (ind.le.1) Then iMask (i, j, k) = Nmax_lev ! change mask Endif ENDDO ENDDO ENDDO GoTo 10 50 Write (*,*) ' Total Lines read=', Nlines ! Change the mask: iteratevely add layers of lower resolution. ! Find blocks adjasent to filled blocks and ! mark them with kSpec =Nspecies-1, -2, ! ! DO kSpec = Nmax_lev - 1, 2, - 1 DO MK3 = 1, Nblocks DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iM = iMask (MK1, MK2, MK3) If (iM.eq.kSpec + 1) Then ! change to kSpec if not high res Do ik = - 1, 1 k = MK3 + ik If (k.le.0) k = k + Nblocks If (k.gt.Nblocks) k = k - Nblocks Do ij = - 1, 1 j = MK2 + ij If (j.le.0) j = j + Nblocks If (j.gt.Nblocks) j = j - Nblocks Do ii = - 1, 1 i = MK1 + ii If (i.le.0) i = i + Nblocks If (i.gt.Nblocks) i = i - Nblocks ! chang If (iMask (i, j, k) .lt.kSpec) iMask (i, j, k) = kSpec ENDDO ENDDO ENDDO EndIf ! end iM=kSpec+1 test ENDDO ! end MK1 ENDDO ! end MK2 ENDDO ! end MK3 Enddo ! end kSpec Do i = 1, Max_lev_abs ! count marked blocks nLevel (i) = 0 Enddo DO MK3 = 1, Nblocks DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iM = iMask (MK1, MK2, MK3) If (iM.le.Nmax_lev.and.iM.ge.Nmin_lev) Then nLevel (iM) = nLevel (iM) + 1 Else Write (*,*) 'wrong mask: ', iM, ' block=', MK1, MK2, MK3 Write (*, * ) ' it should be within limits:', Nmin_lev, Nmax_lev Stop EndIf ENDDO ENDDO ENDDO Ntot = 0 Do i = Max_lev_abs, 1, - 1 icurrent = nLevel (i) * 8** (i - 1) Write (*,*) ' Level=', i, ' Number of Blocks=', nLevel (i) , & ' N_Particles=', icurrent Ntot = Ntot + icurrent Enddo Write (*,*) ' Total expected particles =', Ntot RETURN END SUBROUTINE SetMask !------------------------------------------------ ! Make a realization of spectrum of perturbations ! ALPHA = normalization factor for displacement ! NRAND = seed for random numbers SUBROUTINE SPECTR (NRAND, iDirection) !------------------------------------------------ use Randoma INCLUDE 'luxuryp.h' REAL(8) :: SUMM=0., Wi3, Wj3, Wk3, WD, TS, TRX REAL :: t0=0,t1=0. ! counters for timing INTEGER, DIMENSION(NROW) :: mapz, map3 CALL CPU_TIME(t0) ! initialize time counter iRandom = INT (extras (80) + 1.e-5) ! set type of Random Numbers ! 0 - old GAUSS 1-GAUSS3(ranlux) SUMM = 0. Do j=1,NROW jsign = - 1 jz = j + NSPEC j3 = j - 1 If (j3.GT.NSPEC) j3 = j3 - NSPEC IF (jz.gt.NROW) THEN jz = jz - NROW jsign = 1 ENDIF mapz(j) = jz map3(j) = j3*jsign EndDo !.... Open_MP !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE(Mk3,Mk2,Mk1) DO Mk3 = 1, NROW DO Mk2 = 1, NROW DO Mk1 = 1, NROW GR (Mk1, Mk2, Mk3) = 0. ENDDO ENDDO ENDDO ! Set Spectrum !.... Open_MP !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE(k,ksign,kz,k3,Wk3,j,jsign,jz,j3,Wj3,i,isign,iz,i3,Wi3) & !$OMP PRIVATE(WD,Wk,TS,TRX,m,NRAND,gSet,iFlag) & !$OMP REDUCTION(+:SUMM) DO k = 1, NROW If (iRandom.eq.0) Then NRAND = i24A(k) Else i24 = i24A(k) j24 = j24A(k) in24 = in24A(k) kount = kountA(k) mkount = 0 carry = carryA(k) iflag = iFlagA(k) gset = gSetA(k) Do m=1,24 seeds(m) =SeedsPage(m,k) EndDo EndIf ! end iRandom Wk3 = map3(k)**2 DO j = 1, NROW Wj3 = map3(j)**2 DO i = 1, NROW Wi3 = map3(i)**2 If (iRandom.eq.0) Then TS = GAUSS (NRAND) Else TS = GAUSS3 (gSet,iFlag) EndIf IF (i + j + k.EQ.3) THEN GR (1, 1, 1) = 0. ELSE ! if(k.le.1)write (61,'(4g12.5,3i5)')TS,Wi3,Wj3,Wk3,i,j,k WD = Wi3 + Wj3 + Wk3 Wk = SQRT (WD) TS = TRUNF (Wk) * TS TRX = TS / WD ! start switch If (iDirection.eq.1) Then GR (mapz(i), j, k) = TRX * map3(i) ElseIF (iDirection.eq.2) Then GR (i, mapz(j), k) = TRX * map3(j) Else GR (i, j, mapz(k)) = TRX * map3(k) EndIf ! end iDirection SUMM = SUMM + TS**2 ENDIF ! end kx=ky=kz=0 switch ENDDO ! i ENDDO ! j ENDDO ! k IF (SUMM.LE.0.) Write (*, * ) ' Error!!! Summ over spectrum = 0' ALPHA = AMPLT / SQRT (SUMM) * sqrt (8.) Write (*,'(10x,a20,3g13.6)') 'SPECTR ==> SUMM=', SUMM, ALPHA CALL CPU_TIME(t1) timing(3) = timing(3) +t1-t0 ! time for setting SPECTR ! Do k=1,NROW,8 ! write (50,*) ' -- GR k=',k ! write (50,'(10g12.5)') ((GR(i,j,k),i=1,10),j=1,NROW,8) ! EndDo RETURN END SUBROUTINE SPECTR !------------------------------------------------ ! Define coordinates and velosities for ! particles with resolution = Indx ! Indx =1 = high resolution ! Indx =2 = medium resolution ! Icurrent = number of particles SUBROUTINE BLOCKS (XCONS, VCONS, Indx, Icurrent, iDirection) !------------------------------------------------ Dimension nLevel (10) Real(8) :: sDispl =0. REAL :: t0=0,t1=0. ! counters for timing CALL CPU_TIME(t0) ! initialize time counter Do i = 1, Max_lev_abs ! number of particles at nLevel (i) = 0 ! different levels Enddo Boff = Lblock / 2 + 0.5 sDispl = 0. shift = 0.125 QFACT = FLOAT (NGRID) / FLOAT (NROW) DO Mk3 = 1, Nblocks M3 = (Mk3 - 1) * Lblock ! part of offset DO Mk2 = 1, Nblocks M2 = (Mk2 - 1) * Lblock ! part of offset DO Mk1 = 1, Nblocks M1 = (Mk1 - 1) * Lblock ! part of offset iM = iMask (Mk1, Mk2, Mk3) ! -------------------------- ! M1 + Boff = lagrangian coordinate of 'mean' particle ! Lblock**3 = total number of averaged cells If (iM.eq.Indx) Then If (iM.eq.1) Then !-------------- Low resolution: one particle kp = M3 + (Lblock + 1) / 2 jp = M2 + (Lblock + 1) / 2 ip = M1 + (Lblock + 1) / 2 DX = GR (ip, jp, kp) Icurrent = Icurrent + 1 If (Icurrent.gt.Nmaxpart) Then Write (*,*) 'Too many particles:', Icurrent, ' block=', & Mk1,Mk2, Mk3 STOP EndIf nLevel (Indx) = nLevel (Indx) + 1 If (iDirection.eq.1) Then Q = QFACT * (M1 + Boff - 1.) + 1. ElseIF (iDirection.eq.2) Then Q = QFACT * (M2 + Boff - 1.) + 1. Else Q = QFACT * (M3 + Boff - 1.) + 1. EndIf XPt (Icurrent) = Q - XCONS * DX + shift VXt (Icurrent) = VCONS * DX sDispl = sDispl + DX**2 Else ! --------------- Highest resolution: If (iM.eq.Max_lev_abs) Then Do kp = M3 + 1, M3 + Lblock Do jp = M2 + 1, M2 + Lblock Do ip = M1 + 1, M1 + Lblock DX = GR (ip, jp, kp) ! find displacements Icurrent = Icurrent + 1 If (Icurrent.gt.Nmaxpart) Then Write (*,*) 'Too many particles:', Icurrent, ' block=', & Mk1,Mk2, Mk3 STOP EndIf nLevel (Indx) = nLevel (Indx) + 1 If (iDirection.eq.1) Then Q = QFACT * (ip - 1.) + 1. ElseIF (iDirection.eq.2) Then Q = QFACT * (jp - 1.) + 1. Else Q = QFACT * (kp - 1.) + 1. EndIf XPt (Icurrent) = Q - XCONS * DX + shift VXt (Icurrent) = VCONS * DX sDispl = sDispl + DX**2 EndDo ! end ip EndDo ! end jp EndDo ! end kp Else !------------------- Medium resolution: ! KBlock = 2**(iM-1) = number of particles in 1 ! Nbcells= Lblock/KBlock = number of cells per particl ! (MK-1)*Lblock +(i-1)*Nbcells = offset of cells in 1D ! (MK-1)*Lblock +(i-1)*Nbcells +Nbcells/2 +1/2 = lagra ! coordi ! M1 = (Mk1-1)*Lbloc KBlock = 2** (iM - 1) ! number of particles in 1D per block Nbcells = Lblock / KBlock ! number of cells per particle in 1D Xlag = Nbcells / 2 + 0.5 ! offset for lagrangian coordinates Nbcells3 = Nbcells**3 ! total number of cells per particle Do i = 1, KBlock ! loop over particles i1 = M1 + (i - 1) * Nbcells ! offset for cells in X Q1 = QFACT * (i1 + Xlag - 1.) + 1. ! scale it to box size Do j = 1, KBlock j1 = M2 + (j - 1) * Nbcells ! offset for cells in Y Q2 = QFACT * (j1 + Xlag - 1.) + 1. Do k = 1, KBlock k1 = M3 + (k - 1) * Nbcells ! offset for cells in Z Q3 = QFACT * (k1 + Xlag - 1.) + 1. kp = k1 + (Nbcells + 1) / 2 jp = j1 + (Nbcells + 1) / 2 ip = i1 + (Nbcells + 1) / 2 DX = GR (ip, jp, kp) Icurrent = Icurrent + 1 If (Icurrent.gt.Nmaxpart) Then Write (*,*) 'Too many particles:', Icurrent, ' block=', Mk1, & Mk2, Mk3 STOP EndIf nLevel (Indx) = nLevel (Indx) + 1 If (iDirection.eq.1) Then Q = Q1 ElseIF (iDirection.eq.2) Then Q = Q2 Else Q = Q3 EndIf XPt (Icurrent) = Q - XCONS * DX + shift VXt (Icurrent) = VCONS * DX sDispl = sDispl + DX**2 EndDo ! end k EndDo ! end j EndDo ! end i EndIf ! end iM=Max_lev_abs EndIf ! end iM = 1 EndIf ! end iM = Indx EndDo ! end Mk1 =1,Nblocks EndDo ! end Mk2 =1,Nblocks EndDo ! end Mk3 =1,Nblocks If (nLevel (indx) .gt.0) Then Nspecies = Nspecies + 1 lspecies (Nspecies) = Icurrent wspecies (Nspecies) = PARTW * Lblock**3 / 8** (Indx - 1) Write (*, '(10x,a20,2g13.5,i10)') ' RMS =',XCONS,sDispl,nLevel(indx) Write (*, '(10x,a20,g12.4)') ' RMS 3d diplacement=', & XCONS*sqrt (sDispl/nLevel (indx) ) write (*,'(10x,a20,i4,a,i10,a,6i10)') ' Species=',Nspecies, & ' Current Particles=',Icurrent,' on all levels:', & (nLevel(i),i=1,Max_lev_abs ) EndIf CALL CPU_TIME(t1) timing(4) = timing(4) +t1-t0 ! time for setting SPECTR RETURN END SUBROUTINE BLOCKS !------------------------------------------------ ! FFT of the spectru SUBROUTINE VECTOR (IFOUR) !------------------------------------------------ DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) REAL :: t0=0,t1=0. ! counters for timing CALL CPU_TIME(t0) jb1 = 4 jq = IFOUR CALL SETF67 (jb1, jq, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,& Uf, Vf) ! !.... Open_MP -1- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO K = 1, NROW DO J = 1, NROW DO I = 1, NROW Uf (I) = GR (I, J, K) ENDDO CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) DO I = 1, NROW GR (I, J, K) = Vf (I) ENDDO ENDDO ENDDO Write (*,'(10x,a)') 'Xfft is done' !.... Open_MP -2- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO K = 1, NROW DO I = 1, NROW DO J = 1, NROW Uf (J) = GR (I, J, K) ENDDO CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) DO J = 1, NROW GR (I, J, K) = Vf (J) ENDDO ENDDO ENDDO Write (*,'(10x,a)') 'Yfft is done' write (*,'(10x,a)') 'Swap k<->i' !.... Open_MP !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,aa) DO J=1,NROW DO K=1,NROW-1 DO I=K+1,NROW aa = GR(I,J,K) GR(I,J,K) =GR(K,J,I) GR(K,J,I) =aa ENDDO ENDDO ENDDO !.... Open_MP -3- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO K = 1, NROW DO J = 1, NROW DO I = 1, NROW Uf (I) = GR (I, J, K) ENDDO CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) DO I = 1, NROW GR (I, J, K) = Vf (I) / 8. ENDDO ENDDO ENDDO write (*,'(10x,a)') 'Swap i<->k' !.... Open_MP !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,aa) DO J=1,NROW DO K=1,NROW-1 DO I=K+1,NROW aa = GR(I,J,K) GR(I,J,K) =GR(K,J,I) GR(K,J,I) =aa ENDDO ENDDO ENDDO !!.... Open_MP -3- !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) ! DO J = 1, NROW ! DO I = 1, NROW ! DO K = 1, NROW ! Uf (K) = GR (I, J, K) ! ENDDO ! CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) ! DO K = 1, NROW ! GR (I, J, K) = Vf (K) / 8. ! ENDDO ! ENDDO ! ENDDO Write (*,'(10x,a)') ' FFT is done' CALL CPU_TIME(t1) timing(5) = timing(5) +t1-t0 ! time for VECTOR RETURN END SUBROUTINE VECTOR !--------------------------------------------- ! 1) Update SKINE and Wtotal ! 2) Write current data to disk/tape ! do not write files for 3rd direction SUBROUTINE WriteData (Vscale, AexpV, Wtotal,SKINE,iDirection) !---------------------------------------------- DOUBLE PRECISION Wtotal,SKINE,vvx,vx2 REAL :: t0=0,t1=0. ! counters for timing ! write(30,*) ' XXX' ! test data ! write(30,'(10g12.5)') (XPt(i),i=1,lspecies(Nspecies),5000) ! write(30,*) ' VVV' ! write(30,'(10g12.5)') (VXt(i),i=1,lspecies(Nspecies),5000) CALL CPU_TIME(t0) ! initialize time counter jstart = 1 Do j = 1, Nspecies vvx = 0. vx2 = 0. jend = lspecies (j) W = wspecies (j) If (jend.gt.jstart) Then Do i = jstart, jend vvx = vvx + VXt (i) vx2 = vx2 + VXt (i) **2 SKINE = SKINE+W * VXt (i) **2 Wtotal = Wtotal + w EndDo ! end current species npp = max (jend-jstart + 1, 1) v2 = sqrt ( vx2 / npp) / AexpV * Vscale vvx = vvx / npp / AexpV * Vscale Write ( *,'(22x,"V(km/s)=",2g11.3," Npart=",i10," Species=",i3)') & vvx, v2, jend-jstart + 1, j jstart = jend+1 Endif EndDo ! end all species If(iDirection == 3)Then CALL CPU_TIME(t1) timing(6) = timing(6) +t1-t0 ! time for VECTOR Return ! do not write the last page EndIf Ibuff = 0 ! write data page by page KROW = 0 Do i = 1, lspecies (Nspecies) Ibuff = Ibuff + 1 ! buffer counter XPAR (Ibuff) = XPt (i) VX (Ibuff) = VXt (i) If(Ibuff.ge.NPAGE)Then KROW = KROW +1 ! write (*,*) ' Write page=',KROW,' i=',i,Ibuff CALL WRIROW Ibuff =0 ! initialize buffer counter EndIf EndDo ! end lspecies loop If(Ibuff.ne.0)Then ! write leftover KROW = KROW +1 ! write (*,*) ' Write page=',KROW,' i_buff=',Ibuff CALL WRIROW Ibuff =0 EndIf RETURN CALL CPU_TIME(t1) timing(6) = timing(6) +t1-t0 ! time for VECTOR END SUBROUTINE WriteData !-------------------------------------------------- ! Write current page of particles (x,v) to disk ! SUBROUTINE WRIROW Write (1) XPAR Write (1) VX RETURN END SUBROUTINE WRIROW !-------------------------------------------------- ! Write current PAGE of particles (x,v) to disk ! NRECL - length of ROW block in words SUBROUTINE WRIROWd(IROW) WRITE (21,REC=IROW) RECDAT RETURN END SUBROUTINE WRIROWd !------------------------------------------------------- ! ! Print statistics of spicies Subroutine Species_Change (Wtotal) DOUBLE PRECISION Wtotal Write ( *, 300) Nspecies Write (16, 300) Nspecies 300 FORMAT ('--- New number of species=',i3,/ & 3x,'Specie',2x, & 'Nparticles Ntotal Weight/part Weight_Total ') 310 FORMAT (T5,i3,T13,i8,T22,i8,T30,2g10.3) Wtotal = 0. Do i = 1, Nspecies If (i.eq.1) Then Wtotal = wspecies (i) * lspecies (i) Write ( *, 310) i, lspecies (i), lspecies (i), wspecies (i), & Wtotal Write (16, 310) i, lspecies (i), lspecies (i), wspecies (i), & Wtotal Else ww = wspecies (i) * (lspecies (i) - lspecies (i - 1) ) Write ( *, 310) i, lspecies (i) - lspecies (i - 1), lspecies (i), & wspecies (i), ww Write (16, 310) i, lspecies (i) - lspecies (i - 1), lspecies (i), & wspecies (i), ww Wtotal = Wtotal + ww Endif Enddo Write (*,*) ' Total weight =', Wtotal Return END Subroutine Species_Change !--------------------------------------------------- ! Read current data from disk/tape, ! Open files ! Nrecl is the number of values in a re ! Npage is the number of particles in a SUBROUTINE FilesOpen Character (LEN=50) :: FileName Logical :: FileExists write(*,'(10x,a/,2(30x,a/))')'Create temporary Files====>',NamesDir ! Open file on a tape OPEN (UNIT = 9, FILE = 'PMcrd.DAT', form = 'unformatted') ! this clears old header file Write (9) CLOSE (9) OPEN (9, FILE = 'PMcrd.DAT', FORM = 'UNFORMATTED', STATUS = & 'UNKNOWN') OPEN (UNIT = 16, FILE = 'RESNEW.DAT') OPEN (60, file = 'pt.dat', form = 'unformatted') RETURN END SUBROUTINE FilesOpen !--------------------------------------------- ! Write current data to disk/tape ! SUBROUTINE FilesWrite !---------------------------------------------- ! write header and control data REAL :: t0=0,t1=0. ! counters for timing CALL CPU_TIME(t0) ! initialize time counter 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 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) OPEN(31,FILE=NamesDir(1),form='unformatted') ! x coord and veloc OPEN(32,FILE=NamesDir(2),form='unformatted') ! y XMAX = FLOAT (NGRID) + 1. XSHF = FLOAT (NGRID) Ibuff =0 KROW =0 Npages = (lspecies (Nspecies)-1)/NPAGE+1 If(Npages.eq.0)stop ' Zero number of pages requested. Stop' Do i=1,Npages If(i==Npages)Then NinPage = lspecies (Nspecies) -(i-1)*NPAGE ! # particles in the current page Else NinPage = NPAGE EndIf jfirst = (i-1)*NPAGE +1 jlast = jfirst + NinPage-1 ! write(*,'(10x,a,i5,a,4i9)')'Write page=',i,' Particles=',NinPage,jfirst,jlast read(31)XPAR read(31)VX read(32)YPAR read(32)VY Do j = jfirst,jlast ZPAR(j-jfirst+1) =XPt(j) VZ (j-jfirst+1) =Vxt(j) EndDo Do j=1,NinPage ! Impose Periodical boundary conditions IF (XPAR(j) .GT.XMAX) XPAR(j) = XPAR(j) - XSHF IF (XPAR(j) .LE.1.) XPAR(j) = XPAR(j) + XSHF If (XPAR (j) .lt.1.)write(*,*) ' Boundary Error X!!', i,j,XPAR(j) If (XPAR (j) .ge.XMAX)Then write (*,*)' Fixing Boundary:', j, XPAR(j) XPAR (j) = XPAR (j) - 1.e-4 ENDIF IF (YPAR(j) .GT.XMAX) YPAR(j) = YPAR(j) - XSHF IF (YPAR(j) .LE.1.) YPAR(j) = YPAR(j) + XSHF If (YPAR (j) .lt.1.)write(*,*) ' Boundary Error Y!!', i,j,YPAR(j) If (YPAR (j) .ge.XMAX)Then write (*,*)' Fixing Boundary:', j, YPAR(j) YPAR (j) = YPAR (j) - 1.e-4 ENDIF IF (ZPAR(j) .GT.XMAX) ZPAR(j) = ZPAR(j) - XSHF IF (ZPAR(j) .LE.1.) ZPAR(j) = ZPAR(j) + XSHF If (ZPAR (j) .lt.1.)write(*,*) ' Boundary Error Z!!', i,j,ZPAR(j) If (ZPAR (j) .ge.XMAX)Then write (*,*)' Fixing Boundary:', j, ZPAR(j) ZPAR (j) = ZPAR (j) - 1.e-4 ENDIF EndDo WRITE (21,REC=i) RECDAT EndDo ! end lspecies loop rewind(31) ! erase temporary files write(31)NROW rewind(32) write(32)NROW CALL CPU_TIME(t1) timing(7) = timing(7) +t1-t0 ! time for VECTOR RETURN END SUBROUTINE FilesWrite !-------------------------------------- ! normal random numbers FUNCTION GAUSS3 (gSet,iFlag) ! Uses ranlux for homogenous rand numbers ! Uses Box-Muller inversion for gaussian ! Excellent quality and speed ! N= 100000000 GAUSS3+ranlux ------ ! sigma= 0.99997 mean= 0.18692E-03 ! sigma Frac(>sig) Frac(<-sig) True n(>) n(<) ! 1.00 0.1587 0.1586 0.1587 15869095 15857724 ! 1.50 0.6682E-01 0.6680E-01 0.6681E-01 6681652 6679768 ! 2.00 0.2277E-01 0.2273E-01 0.2275E-01 2276651 2273197 ! 2.50 0.6222E-02 0.6206E-02 0.6210E-02 622190 620610 ! 3.00 0.1356E-02 0.1347E-02 0.1350E-02 135628 134674 ! 4.00 0.3242E-04 0.3145E-04 0.3170E-04 3242 3145 ! !-------------------------------------- DIMENSION RanNum (2) If (iFlag.eq.0) Then 1 CALL ranlux (RanNum, 2 ) x1 = 2. * RanNum (1) - 1. x2 = 2. * RanNum (2) - 1. R = x1**2 + x2**2 If (R.ge.1.) GoTo 1 F = sqrt ( - 2. * LOG (R) / R) gSet = x1 * F GAUSS3 = x2 * F iFlag = 1 Else GAUSS3 = gSet iFlag = 0 EndIf RETURN END FUNCTION GAUSS3 !-------------------------------------------------------------------- ! LUXURY LEVELS. ! level 0 (p=24): equivalent to the original RCARRY of Marsaglia ! and Zaman, very long period, but fails many tests. ! level 1 (p=48): considerable improvement in quality over level 0, ! now passes the gap test, but still fails spectral test. ! level 2 (p=97): passes all known tests, but theoretically still ! defective. ! level 3 (p=223): DEFAULT VALUE. Any theoretically possible ! correlations have very small chance of being observed. ! level 4 (p=389): highest possible luxury, all 24 bits chaotic. !--------------------------------------------- SUBROUTINE ranlux (rvec, lenv) ! returns a vector RVEC of LEN ! 32-bit random floating point numbers between ! zero (not included) and one (also not incl.). ! Next call to ranlux gives next LEN of random numbers !--------------------------------------------- INTEGER lenv REAL rvec (lenv) INTEGER iseeds (24) INCLUDE 'luxuryp.h' ! DATA ndskip / 0, 24, 73, 199, 365 / ! DATA i24, j24, luxlev / 24, 10, lxdflt / ! DATA notyet / .true. / ! DATA in24, kount, mkount, carry / 0, 0, 0, 0. / ! NOTYET is .TRUE. if no initialization has been performed yet. ! Default Initialization by Multiplicative Congruential IF (notyet) THEN stop ' Access to runlux before initialization STOP' ENDIF ! The Generator proper: "Subtract-with-borrow", ! as proposed by Marsaglia and Zaman, ! Florida State University, March, 1989 DO ivec = 1, lenv uni = seeds (j24) - seeds (i24) - carry IF (uni.LT.0.) THEN uni = uni + 1.0 carry = twom24 ELSE carry = 0. ENDIF seeds (i24) = uni i24 = next (i24) j24 = next (j24) rvec (ivec) = uni ! small numbers (with less than 12 "significant" bits) are "padde IF (uni.LT.twom12) THEN rvec (ivec) = rvec (ivec) + twom24 * seeds (j24) ! and zero is forbidden in case someone takes a logarithm IF (rvec (ivec) .EQ.0.) rvec (ivec) = twom24 * twom24 ENDIF ! Skipping to luxury. As proposed by Martin Luscher. in24 = in24 + 1 IF (in24.EQ.24) THEN in24 = 0 kount = kount + nskip DO isk = 1, nskip uni = seeds (j24) - seeds (i24) - carry IF (uni.LT.0.) THEN uni = uni + 1.0 carry = twom24 ELSE carry = 0. ENDIF seeds (i24) = uni i24 = next (i24) j24 = next (j24) ENDDO ENDIF ENDDO kount = kount + lenv IF (kount.GE.igiga) THEN mkount = mkount + 1 kount = kount - igiga ENDIF RETURN ! SUBROUTINE ranlux END SUBROUTINE ranlux !--------------------------------------------- ! Subroutine to initialize from one or three integers SUBROUTINE rluxgo (lux, ins, k1, k2) ! initializes the generator from ! one 32-bit integer INT and sets Luxury Level LUX ! which is integer between zero and MAXLEV, or if ! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ! should be set to zero unless restarting at a break ! point given by output of RLUXAT !--------------------------------------------- INCLUDE 'luxuryp.h' INTEGER iseeds(24) IF (lux.LT.0) THEN luxlev = lxdflt ELSEIF (lux.LE.maxlev) THEN luxlev = lux ELSEIF (lux.LT.24.OR.lux.GT.2000) THEN luxlev = maxlev Write (6, '(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ', lux ELSE luxlev = lux DO ilx = 0, maxlev IF (lux.EQ.ndskip (ilx) + 24) luxlev = ilx ENDDO ENDIF IF (luxlev.LE.maxlev) THEN nskip = ndskip (luxlev) Write (6, '(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', & luxlev, ' P=', nskip + 24 ELSE nskip = luxlev - 24 Write (6, '(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:', & luxlev ENDIF in24 = 0 IF (ins.LT.0) Write (6,'(A)')' Illegal initialization by RLUXGO ',& 'negative input seed' IF (ins.GT.0) THEN jseed = ins Write (6, '(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', & jseed, k1, k2 ELSE jseed = jsdflt Write (6, '(A)') ' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' ENDIF inseed = jseed notyet = .false. twom24 = 1. DO i = 1, 24 twom24 = twom24 * 0.5 k = jseed / 53668 jseed = 40014 * (jseed-k * 53668) - k * 12211 IF (jseed.LT.0) jseed = jseed+icons iseeds (i) = MOD (jseed, itwo24) ENDDO twom12 = twom24 * 4096. DO i = 1, 24 seeds (i) = REAL (iseeds (i) ) * twom24 next (i) = i - 1 ENDDO next (1) = 24 i24 = 24 j24 = 10 carry = 0. IF (seeds (24) .EQ.0.) carry = twom24 ! If restarting at a break point, skip K1 + IGIGA*K2 ! Note that this is the number of numbers delivered to ! the user PLUS the number skipped (if luxury .GT. 0). kount = k1 mkount = k2 IF (k1 + k2.NE.0) THEN DO iouter = 1, k2 + 1 inner = igiga IF (iouter.EQ.k2 + 1) inner = k1 DO isk = 1, inner uni = seeds (j24) - seeds (i24) - carry IF (uni.LT.0.) THEN uni = uni + 1.0 carry = twom24 ELSE carry = 0. ENDIF seeds (i24) = uni i24 = next (i24) j24 = next (j24) ENDDO ENDDO ! Get the right value of IN24 by direct calculation in24 = MOD (kount, nskip + 24) IF (mkount.GT.0) THEN izip = MOD (igiga, nskip + 24) izip2 = mkount * izip + in24 in24 = MOD (izip2, nskip + 24) ENDIF ! Now IN24 had better be between zero and 23 inclusive IF (in24.GT.23) THEN Write (6, '(A/A,3I11,A,I5)') ' Error in RESTARTING with RLUXGO:',& ' The values', ins, k1, k2, ' cannot occur at luxury level', luxlev in24 = 0 ENDIF ENDIF RETURN ! SUBROUTINE rluxgo END SUBROUTINE rluxgo !---------------------------------- Read in variables REAL FUNCTION INPUT (text) !------------------------------------------------ Character text * (*) Write (*,'(A,$)') text READ (*,*) x INPUT = x Return END FUNCTION INPUT !-------------------------------------- Fourier Transform SUBROUTINE SETF67 (JBC1, JQ1, jbc, jp, jsl, ll1, k2, k3, k4, k7, & Qi, jndx, Uf, Vf) ! ----------------------------------------------------- DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) PI = DATAN (1.D0) * 4. JBC = JBC1 JQ = JQ1 IF (JBC.LT.3) GOTO 101 JQ = JQ - 1 101 K3 = 2**JQ K7 = K3 / 2 N5 = K3 / 4 I = 1 JNDX (I) = N5 QI (I) = 0.5 * SQRT (2.) K = I I = I + 1 102 IL = I IF (I.EQ.K7) GOTO 104 103 K1 = JNDX (K) / 2 JNDX (I) = K1 QI (I) = SIN (PI * FLOAT (K1) / FLOAT (K3) ) K1 = K7 - K1 I = I + 1 JNDX (I) = K1 QI (I) = SIN (PI * FLOAT (K1) / FLOAT (K3) ) K = K + 1 I = I + 1 IF (K.EQ.IL) GOTO 102 GOTO 103 104 RETURN END SUBROUTINE SETF67 !----------------------------------------- SUBROUTINE TFOLD (IS, L, ZZZ, k2) DIMENSION ZZZ (MARR) IH2 = K2 / 2 - 1 DO 100 I = IS, IH2 I1 = I + L I2 = K2 - I + L A = ZZZ (I1) ZZZ (I1) = A - ZZZ (I2) ZZZ (I2) = A + ZZZ (I2) 100 END DO RETURN END SUBROUTINE TFOLD !------------------------------------ SUBROUTINE NEG (I1, I3, I2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi,& jndx, Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) DO K = I1, I3, I2 Vf (K) = - Vf (K) END DO RETURN END SUBROUTINE NEG !------------------------------------- SUBROUTINE REVNEG (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 DO I = 1, K7 J = K3 + 1 + I K = K4 + 1 - I A = Vf (J) Vf (J) = - Vf (K) Vf (K) = - A END DO RETURN END SUBROUTINE REVNEG !------------------------------------- SUBROUTINE ZEERO (L, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 DO I = 1, K2 LI = L + I Uf (LI - 1) = 0.0 END DO RETURN END SUBROUTINE ZEERO !------------------------------------------ SUBROUTINE TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 II = K2 - 1 DO I = 1, II I1 = JSL + I I2 = LL1 - I A = Uf (I1) Uf (I1) = A + Uf (I2) Uf (I2) = A - Uf (I2) END DO RETURN END SUBROUTINE TFOLD1 !------------------------------------- SUBROUTINE KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,& Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 JS1 = K2 I = 1 J5 = JSL + K2 IS1 = JSL IC1 = LL1 JS1 = JS1 / 2 IF (JS1.NE.1) GOTO 200 K1 = JNDX (I) SN = QI (I) IS0 = IS1 IS1 = IS1 + JS1 IC0 = IC1 IC1 = IC1 + JS1 ODD1 = SN * (Uf (IC1) - Uf (IS1) ) Vf (K1 + 1) = Uf (IC0) + ODD1 K3MK1 = K3 - K1 Vf (K3MK1 + 1) = Uf (IC0) - ODD1 IF (JBC.LT.3) GOTO 110 ODD2 = SN * (Uf (IC1) + Uf (IS1) ) K3PK1 = K3 + K1 Vf (K3PK1 + 1) = Uf (IS0) + ODD2 K4MK1 = K4 - K1 Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 110 RETURN 200 SN = QI (I) IS1 = IS1 + JS1 IC1 = IC1 + JS1 J3 = IS1 + JS1 210 IS0 = IS1 - JS1 IC0 = IC1 - JS1 ODD1 = SN * (Uf (IC1) - Uf (IS1) ) ODD2 = SN * (Uf (IC1) + Uf (IS1) ) Uf (IC1) = Uf (IC0) - ODD1 Uf (IS1) = - Uf (IS0) + ODD2 Uf (IC0) = Uf (IC0) + ODD1 Uf (IS0) = Uf (IS0) + ODD2 IS1 = IS1 + 1 IC1 = IC1 + 1 IF (IS1.NE.J3) GOTO 210 I = I + 1 300 IS1 = JSL IC1 = LL1 JS1 = JS1 / 2 IF (JS1.EQ.1) GOTO 400 310 SN = QI (I) I = I + 1 CS = QI (I) IS1 = IS1 + JS1 IC1 = IC1 + JS1 J3 = IS1 + JS1 320 IS0 = IS1 - JS1 IC0 = IC1 - JS1 ODD1 = CS * Uf (IC1) - SN * Uf (IS1) ODD2 = SN * Uf (IC1) + CS * Uf (IS1) Uf (IC1) = Uf (IC0) - ODD1 Uf (IC0) = Uf (IC0) + ODD1 Uf (IS1) = - Uf (IS0) + ODD2 Uf (IS0) = Uf (IS0) + ODD2 IS1 = IS1 + 1 IC1 = IC1 + 1 IF (IS1.NE.J3) GOTO 320 IS1 = IS1 + JS1 IC1 = IC1 + JS1 J3 = IS1 + JS1 330 IS0 = IS1 - JS1 IC0 = IC1 - JS1 ODD1 = SN * Uf (IC1) - CS * Uf (IS1) ODD2 = CS * Uf (IC1) + SN * Uf (IS1) Uf (IC1) = Uf (IC0) - ODD1 Uf (IC0) = Uf (IC0) + ODD1 Uf (IS1) = - Uf (IS0) + ODD2 Uf (IS0) = Uf (IS0) + ODD2 IS1 = IS1 + 1 IC1 = IC1 + 1 IF (IS1.NE.J3) GOTO 330 I = I + 1 IF (IS1.EQ.J5) GOTO 300 GOTO 310 400 K1 = JNDX (I) SN = QI (I) I = I + 1 CS = QI (I) IS0 = IS1 IS1 = IS1 + JS1 IC0 = IC1 IC1 = IC1 + JS1 ODD1 = CS * Uf (IC1) - SN * Uf (IS1) Vf (K1 + 1) = Uf (IC0) + ODD1 K3MK1 = K3 - K1 Vf (K3MK1 + 1) = Uf (IC0) - ODD1 IF (JBC.LT.3) GOTO 410 ODD2 = SN * Uf (IC1) + CS * Uf (IS1) K3PK1 = K3 + K1 Vf (K3PK1 + 1) = Uf (IS0) + ODD2 K4MK1 = K4 - K1 Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 410 IS1 = IS1 + 1 IC1 = IC1 + 1 K1 = JNDX (I) IS0 = IS1 IS1 = IS1 + JS1 IC0 = IC1 IC1 = IC1 + JS1 ODD1 = SN * Uf (IC1) - CS * Uf (IS1) Vf (K1 + 1) = Uf (IC0) + ODD1 K3MK1 = K3 - K1 Vf (K3MK1 + 1) = Uf (IC0) - ODD1 IF (JBC.LT.3) GOTO 420 ODD2 = CS * Uf (IC1) + SN * Uf (IS1) K3PK1 = K3 + K1 Vf (K3PK1 + 1) = Uf (IS0) + ODD2 K4MK1 = K4 - K1 Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 420 IS1 = IS1 + 1 IC1 = IC1 + 1 I = I + 1 IF (IS1.NE.J5) GOTO 400 RETURN END SUBROUTINE KFOLD ! ------------------------------------------------------ SUBROUTINE FOUR67 (JBC1, JQ1, jp1, jsl1, ll11, k21, k71, Qi, jndx,& Uf, Vf) ! ------------------------------------------------------ DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 JBC = JBC1 JQ = JQ1 jp = jp1 jsl = jsl1 ll1 = ll11 k2 = k21 k7 = k71 A5 = 0.5 * SQRT (2.0) K4 = 2**JQ K3 = K4 GOTO (103, 103, 101, 102), JBC 101 Uf (1) = Uf (1) / 2.0 Uf (K3 + 1) = Uf (1) K2 = K3 CALL TFOLD (0, 1, Uf, k2) K3 = K3 / 2 JQ = JQ - 1 GOTO 103 102 K3 = K3 / 2 JQ = JQ - 1 103 N5 = K3 / 4 K7 = K3 / 2 N11 = 3 * K7 K31 = K3 + 1 GOTO (300, 400, 500, 600), JBC 300 Uf (K31) = 0.0 Uf (1) = 0.0 K2 = K3 DO 301 I = 2, JQ CALL TFOLD (1, 1, Uf, k2) 301 K2 = K2 / 2 Vf (K7 + 1) = Uf (2) JF = N5 JSL = 1 DO 302 JP = 2, JQ LL1 = K2 + 1 CALL ZEERO (1, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,& Vf) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) I1 = 3 * JF + 1 I2 = 4 * JF I3 = I1 + (K2 / 2 - 1) * I2 CALL NEG (I1, I3, I2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, & jndx, Uf, Vf) K2 = K2 + K2 302 JF = JF / 2 RETURN 400 Uf (1) = Uf (1) / 2.0 Uf (K31) = Uf (K31) / 2.0 K2 = K3 DO 401 I = 2, JQ CALL TFOLD (0, K31 - K2, Uf, k2) 401 K2 = K2 / 2 LL1 = K31 - K2 A = Uf (LL1) + Uf (LL1 + 2) Vf (1) = - A - Uf (LL1 + 1) Vf (K31) = - A + Uf (LL1 + 1) Vf (K7 + 1) = Uf (LL1) - Uf (LL1 + 2) DO 402 JP = 2, JQ JSL = K31 - K2 LL1 = JSL - K2 idum = jsl CALL ZEERO (idum, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) 402 K2 = K2 + K2 CALL NEG (1, K31, 2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) RETURN 500 K2 = K3 L2 = K4 DO 501 JP = 2, JQ CALL TFOLD (1, 1, Uf, k2) CALL TFOLD (0, L2 - K2 + 1, Uf, k2) 501 K2 = K2 / 2 LL1 = L2 - K2 + 1 A = Uf (LL1) + Uf (LL1 + 2) Vf (K7 + 1) = 2.0 * ( - Uf (LL1) + Uf (LL1 + 2) ) Vf (1) = 2.0 * (A + Uf (LL1 + 1) ) Vf (K31) = 2.0 * (A - Uf (LL1 + 1) ) Vf (N11 + 1) = 2.0 * Uf (2) DO 502 JP = 2, JQ Uf (K2 + 1) = 2.0 * Uf (K2 + 1) JSL = K2 + 1 CALL TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) LL1 = LL1 - K2 Uf (LL1) = - 2.0 * Uf (LL1) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) 502 K2 = K2 + K2 Vf (1) = Vf (1) * A5 Vf (K31) = Vf (K31) * A5 RETURN 600 Uf (1) = Uf (1) * A5 Uf (K31) = Uf (K31) * A5 K2 = K3 DO 601 JP = 2, JQ CALL TFOLD (0, K31 - K2, Uf, k2) CALL TFOLD (1, K31, Uf, k2) 601 K2 = K2 / 2 LL1 = K31 - K2 A = Uf (LL1) + Uf (LL1 + 2) Vf (1) = 2.0 * (A + Uf (LL1 + 1) ) Vf (K31) = 2.0 * (A - Uf (LL1 + 1) ) Vf (K7 + 1) = 2.0 * ( - Uf (LL1) + Uf (LL1 + 2) ) Vf (N11 + 1) = 2.0 * Uf (K3 + 2) DO 602 JP = 2, JQ JSL = K31 + K2 Uf (JSL) = 2.0 * Uf (JSL) CALL TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) LL1 = LL1 - K2 Uf (LL1) = - 2.0 * Uf (LL1) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) 602 K2 = K2 + K2 CALL REVNEG (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, Vf) idum = k3 CALL NEG (2, idum, 2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,& Uf, Vf) K2 = K4 CALL TFOLD (1, 1, Vf, k2) Vf (K4 + 1) = 0.0 RETURN END SUBROUTINE FOUR67 end module setInitialConditions ! ! ! BLOCKDATA PARAMETER ( maxlev = 4, lxdflt = 3) LOGICAL notyet COMMON /LUX1/ & i24, j24, in24, kount, mkount, carry, seeds(24) COMMON /LUX2/ notyet, twom24, twom12, & nskip, ndskip(0:maxlev), next(24), inseed, luxlev DATA ndskip / 0, 24, 73, 199, 365 / DATA i24, j24, luxlev / 24, 10, lxdflt / DATA notyet / .true. / DATA in24, kount, mkount, carry / 0, 0, 0, 0. / End BLOCKDATA !------------------------------------------------------------------------- PROGRAM PMstartMp use setInitialConditions ! NROW = number of particles in 1 dimension ! NGRID= number of cells in 1 dimension ! Npage = number of particles per page = NROW**2 ! COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb,QSCALE DOUBLE PRECISION :: Wtotal=0., SKINE=0. REAL, PARAMETER :: PI=3.1415926535 REAL :: t0=0,t1=0. ! counters for timing Logical :: iCont =.true. ! continue (true) after setting seeds CALL CPU_TIME(t0) Open(30,file='TestDataParallel.dat') Wtotal = (float(NROW)**3*4+2.*Nmaxpart*4. & + 31.*NROW*4 +float(Nblocks)**3*4 & + 6. *NPAGE*4)/1024.**3 write (*,'(//10x,a,g13.4,a/)') ' Memory required=',Wtotal,'Gb' CALL InitValues (NBYTE, SCALEL) Write (*,*) ' InitValues is done. N bytes=', NBYTE Max_lev_abs = INT (ALOG (REAL (Lblock) ) / ALOG (2.) + 1.1) Write (*,'(10x,a,2i4)') 'Absolute Limits on mass resolution levels: ', 1,& Max_lev_abs Write (*,'(10x,2(a,i4))') 'Actual Max resolution Level=', Nmax_lev, & ' Min Level=', Nmin_lev Write (16,'(10x,2(a,i4))') 'Actual Max resolution Level=', Nmax_lev, & ' Min Level=', Nmin_lev If (Nmax_lev.gt.Max_lev_abs) Then Write (*,*) ' Wrong number of Levels for mass refinement (', & Nmax_lev, ') Block size =', Lblock, ' Use Absolute value=', & Max_lev_abs Nmax_lev = Max_lev_abs Endif CALL FilesOpen ! this opens files on disk Write (16, * ) ' Seed=', Nseed AEXP0 = AEXPN AEXPV = AEXPN - ASTEP / 2. CALL ZeroMASK ! set mask for multiple mass resolution If (Nmax_lev.ne.Nmin_lev) Then Open (30, file = 'SelCoords.DAT', status = 'unknown') CALL SetMASK EndIf ! set scaling constants for coordinates and velocities IFOUR = INT (ALOG (FLOAT (NROW) ) / ALOG (2.) + 0.5) Write (*,'(3x,a12,i4,a,i4)') ' Ifour =', IFOUR, ' Nrow=', NROW Fact = sqrt (Om0 + Oml0 * AEXPV**3 + Ocurv * AEXPV) QFACT = FLOAT (NGRID) / FLOAT (NROW) Vscale = ScaleL * 100. * hubble / NGRID ! get a realization of spectrum in /GRID/ NRAND = Nseed CALL CPU_TIME(t1) timing(1) = t1-t0 ! time for intialization CALL SetRandom(iCont) ! set seeds for ramdom numbers If(.NOT.iCont)Then write(*,'(2(10x,a,f10.2/))') & 'cpu_time for init =',timing(1), & 'cpu_time for seeds=',timing(2) Stop EndIf SKINE = 0. ! counter for kinetic energy Wtotal = 0. ! counter for total weight ALLOCATE(GR(NROW,NROW,NROW)) ! allocate memory for spectrum Do iDirection = 1, 3 ! ----------------- loop over x,y,z directions Write(*,'(//3x,a12,i3,15("-"))' )'- Direction=',iDirection If(iDirection < 3) & OPEN(1,FILE=NamesDir(iDirection), FORM='UNFORMATTED', & STATUS ='UNKNOWN') CALL SPECTR (NRAND, iDirection) ! get the displacement vector by FFT VCONS = - ALPHA/(2.*PI/NGRID)*(AEXPV/AEXP0)*SQRT(AEXPV)*Fact XCONS = ALPHA / (2. * PI / NGRID) * (AEXPN / AEXP0) Write(*,'(3x,a12,g12.4,a,g12.4)' ) 'Scaling:(x)=', XCONS, ' (v)=', VCONS CALL VECTOR (IFOUR) Write (*,'(3x,a12,i4,a,g12.4)' ) 'SPECTR done:', IFOUR, ' Alpha=', ALPHA write(*,'(2(a,f8.3))')' Spectr t=',timing(3),' FFT=',timing(5) ! find x,v and write to file, page by page Do i = 1, Max_lev_abs lspecies (i) = 0 Enddo Nspecies = 0 Icurrent = 0 ! current particle Do kSpec = Max_lev_abs, 1, - 1 CALL BLOCKS (XCONS, VCONS, kSpec, Icurrent, iDirection) Enddo CALL WriteData (Vscale, AexpV, Wtotal,SKINE,iDirection) CLOSE(1) EndDo DEALLOCATE (GR) CALL Species_Change (Wtotal) ! print info on species EKIN = 0.5 * SKINE / AEXPV**2 Write (*, '('' Ekin='',E12.4,'' Weght per cell='',g12.5,'' (must be 1)'')') EKIN,Wtotal/NGRID**3 Write (16, '('' Ekin='',E12.4,'' Weght per cell='',g12.5,'' (must be 1)'')') EKIN,Wtotal/NGRID**3 CALL FilesWrite ! write header and control data Write (60) astep ! write pt.dat file: time-step for particles CALL CPU_TIME(t1) write(*,'(10x,a/,(10x,a15,f10.2))') & 'Wall cpu_time (seconds) ', & 'Init =',timing(1), & 'Seeds =',timing(2), & 'Spectr =',timing(3), & 'Coords =',timing(4), & 'FFT =',timing(5), & 'Write temp =',timing(6), & 'Write final =',timing(7), & 'Total =',t1-t0 STOP END PROGRAM PMstartMp