C _______________________ START 3-D SIMULATIONS C C Klypin, August 1993 C INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' C NROW = number of particles in 1 dimension C NGRID= number of cells in 1 dimension C Npage = number of particles per page = NROW**2 C COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Par(6),ns,qqscaleb,QSCALE COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2 Character Hd*5,Tail*4 DATA PI /3.1415926535/ CALL InitValues(NBYTE,SCALEL) write (*,*) ' InitValues is done. N bytes=',NBYTE Max_lev_abs = INT(ALOG(REAL(Lblock))/ALOG(2.)+1.1) write (*,*) ' Absolute Limits on mass resolution levels: ', & 1,Max_lev_abs write (*,*) ' Actual Max resolution Level=',Nmax_lev, & ' Min Level=',Nmin_lev write (16,*) ' 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 OPEN(UNIT=9,FILE='PMcrd.DAT',form='unformatted') write (9) ! this clears old header file CLOSE (9) CALL RDTAPE ! this opens files on disk OPEN(UNIT=16,FILE='RESNEW.DAT') AEXP0 = AEXPN AEXPV = AEXPN - ASTEP/2. C set mask for multiple mass resolution CALL ZeroMASK If(Nmax_lev.ne.Nmin_lev)Then Open(30,file='SelCoords.DAT') CALL SetMASK EndIf C get a realization of spectrum in /GRID/ NRAND =Nseed CALL SPectr_CR(NRAND) CALL SPECTR(NRAND,Levels) ! set spectrum for higher levels CALL RELOCATE_SP ! remove harmonics included at lower levels C get the displacement vector by FFT CALL VECTOR write (*,*) ' SPECTR is done. Alpha=', & (ALPHA(i),i=1,Levels) Fact = sqrt(Om0+Oml0*AEXPV**3+Ocurv*AEXPV) VCONS = -1./(2.*PI/NGRID)*(AEXPV/AEXP0)* & SQRT(AEXPV)*Fact XCONS = 1./(2.*PI/NGRID)*(AEXPN/AEXP0) write (*,*) ' Cons=', XCONS ,VCONS QFACT = FLOAT(NGRID)/FLOAT(NROW) SKINE = 0. SX = 0. SY = 0. SZ = 0. SX2 = 0. SY2 = 0. SZ2 = 0. C find x,v and write to file, page by page Do i=1,Max_lev_abs lspecies(i) =0 Enddo Nspecies =0 C find x,v for multiple mass resolution code Icurrent =0 ! current particle Do kSpec =Max_lev_abs,1,-1 CALL BLOCKS(XCONS,VCONS,kSpec,Icurrent) write (*,*) ' Level=',Max_lev_abs-kSpec+1, & ' Nparticles=',lspecies(Max_lev_abs-kSpec+1) Enddo Vscale =ScaleL*100.*hubble/NGRID CALL WriteData(Vscale,AexpV,Wtotal) CALL Species_Change(Wtotal_o) ! rearrange 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 C CALL WRTAPE ! write header and control data C write pt.dat file: time-step for particles Nparticles =lspecies(Nspecies) c Do ic1 =1,Nparticles c XPt(ic1) =astep c Enddo open ( 60 , file = 'pt.dat' , form = 'unformatted' ) c write(60) (XPt(ic1),ic1=1,Nparticles) write(60) astep write (*,*) ' pt=', astep !(XPt(ic1),ic1=1,10) STOP END C-------------------------------------------------- C sqrt(Power spectrum) C k = (2pi/L) wk REAL FUNCTION TRUNF(wk) C-------------------------------------------------- INCLUDE 'PMparameters.h' PARAMETER (NSPEC =NROW/2) 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 If(iFlag.gt.0)Then ! use approximations If(Par(6).ne.0.)Then ! Holtzman approx sk= sqrt(k) TRUNF = k**(ns/2.) / . (1.+sk*(Par(2) . +sk*(Par(3) . +sk*(Par(4) . +sk*(Par(5) )))) )**(Par(6)) Else ! BBKS + Sugiyama approx c Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om)) c 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 Else ! use table for P(k) TRUNF = sqrt(Ppk(k)) EndIf RETURN END C--------------------------------------- FUNCTION Ppk(x) c interpolate table with p(k) c x is in real 1/Mpc C--------------------------------------- PARAMETER (NtabM = 100000) COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),Ntab,StepK,alog0,iFlag If(x.ge.xkt(Ntab))Then ! slope is ns =-3 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 c write (*,'(5x," k=",g12.4," table=",2g12.4," P=",3g12.4)') c & x,xkt(ind),xkt(ind+1), c & Ppk,Pkt(ind),Pkt(ind+1) Return End C--------------------------------------- SUBROUTINE ReadTable(hsmall) c interpolate table with p(k) C--------------------------------------- PARAMETER (NtabM = 100000) COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),Ntab,StepK,alog0,iFlag character*80 Line If(iFlag.eq.-1)Then ! W.Hu table 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 file WHUpk.dat' write (*,*) ' hsmall =',hsmall Goto 50 EndIf If(iFlag.eq.-2)Then ! Eis.Hu table Open(50,file='pk_EHu_Om0.3_Omb0.043_s8=0.9.dat') write (*,*)' Read data from file', & ' pk_EHu_Om0.3_Omb0.043_s8=0.9.dat' Else Open(50,file='pk_EHu_WMAP06.dat') write (*,*)' Read data from file', & ' pk_EHu_WMAP06.dat' EndIf Do i=1,9 ! read header 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 file ...' write (*,*) ' hsmall =',hsmall 50 If(Ntab.le.1)stop 'wrong table for p(k)' StepK =log10(xkt(2)/xkt(1)) alog0 = log10(xkt(1)) Do k =2,Ntab-1 ! test that spacing is the same 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 C********************************************************************* C Power spectrum for CDM. Wk- wavenumber(Mpc^-1, h=0.5) C P = wk *T^2(k) C QSCALE = 2pi/Box_size, QS =hubble**(-2) c FUNCTION TRUNF(WK) c COMMON / TRUNCOM/ QSCALE,g,gy c COMMON / TRUNCDM/ QS,A1,A2,A3,A4,A5,A32,A43,A54 c Q = QSCALE*QS*WK c TRUNF = sqrt( WK* c + (LOG(1.+A1*Q)/(A1*Q))**2/ c + SQRT(1.+Q*(A2 +Q*(A32 +Q*(A43+Q*A54)))) ) c c RETURN c END C********************************************************************* C INITIALIZE CONSTANTS: C Scalel =length of square in MPC C AEXPN =expansion factor (current) C ASTEP =increment of expansion factor: AEXPN =AEXP_0+N*ASTEP C PARTW =weight of a particle: total 'mass' per cell must be unity c even for LCDM or OCDM C AMPLT =amplitude of perturbations inside the simulation box C NBYTE =record length in bytes C SUBROUTINE InitValues(NBYTE,SCALEL) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' 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_CR(1)=INPUT(' Amplitude of density fluctuations CR =') AMPLT_CR(2)=INPUT(' Amplitude of density fluctuations noCR=') 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 extras(80) = 0. EndIf c extras(80) is used as switch for c random numbers generators c = 0 - old Randd+Gauss c = 1 - luxury+Gauss3 c extras(81) = lux - level of luxury extras(100) = SCALEL ! store the box size 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 SCALEL=SCALEL/hubble ! scale it to real megaparsecs 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 If(Nmin_lev.lt.Nmax_lev)Then ! multiple mass resolution write(*,*) write(*,*) ' You use multiple mass resolution' endif If(Nmin_lev.gt.Nmax_lev)Then ! change the order of levels 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 (*,*) ' Qscale=',QSCALE,SCALEL,ns RETURN END C________________________________________Read parameters of the model C Om = Omega_0 C Omb = Omega_baryon C Omnu = Omega_neutrino C hsmall = hubble constant C Par = fitting parameters C Par(6) = 0 --- bbks-style+Hu&Sugiyama C Par(6) ne 0 --- holtzman-style SUBROUTINE MODEL(hsmall) C--------------------------------------- INCLUDE 'PMparameters.h' 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 =') write (*,*) ' Line number=',Line1 If(Line1.gt.0)Then 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) theta = 2.726/2.7 ! = T_cmb/2.7 Ob0 =Omb/Om ! Hu&Sugiyama fitting 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)) alpha =1./a1**Ob0/a2**(Ob0**3) qqscaleb = theta**2/(1.-Ob0)**0.60/sqrt(alpha)/Omh2 Write (*,20)Om,Omb,Omc,Omnu,hsmall,Par c 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' If(iFlag.eq.-1)Then ! W.Hu table hsmall = 0.7 Omb = 0.024/hsmall**2 Om = 0.3 Ocurv = 0. Omc = Om EndIf If(iFlag.eq.-2)Then ! Eis.Hu approximation hsmall = 0.7 Omb = 0.04 Om = 0.3 Ocurv = 0. Omc = Om EndIf If(iFlag.eq.-3)Then ! Eis.Hu approximation: WMAP hsmall = 0.73 Omb = 0.04 Om = 0.24 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 C------------------------------------------------ C Clear MASK for multiple mass resolution C Initiat to level Nmin_lev SUBROUTINE ZeroMask C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' DO MK3 = 1,Nblocks DO MK2 = 1,Nblocks DO MK1 = 1,Nblocks iMask(MK1,MK2,MK3) =Nmin_lev ENDDO ENDDO ENDDO RETURN END C------------------------------------------------ C Set MASK for multiple mass resolution C Read particles with lagrangian coordinates C to set the mask C Block = 1 = low resolution C = 2 = higher resolution C = 3 = even higher resolution C After the mask is set, make maximum resolution Nmax_lev C SUBROUTINE SetMask C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' Dimension nLevel(10) Character*70 Line Do i=1,5 Read(30,'(A)') Line Write(*,'(A)') Line EnDDo Nlines =0 write (*,*) ' Nblocks=',Nblocks 10 Read(30,*,end=50,err=50) i,xw,yw,zw, + vwx,vwy,vwz,im,jm,km Nlines =Nlines +1 c 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 c If(iMask(im,jm,km).ne.0)Then c write (*,*) ' Occupied block: ',im,jm,km,' Line=',Nlines c write (*,*) ' Block=',iMask(im,jm,km) c Stop c EndIf iMask(im,jm,km) =Nmax_lev ! mark high resolution block 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 in crest ! are marked out c If(iMask(i,j,k).le.0.or.iMask(i,j,k).gt.Nspecies) c + Then c Stop c Endif If(ind.le.1)Then iMask(i,j,k)=Nmax_lev ! change mask Endif ENDDO ENDDO ENDDO GoTo 10 C Change the mask: iteratevely add layers of lower resolution 50 write(*,*) ' Total Lines read=',Nlines DO kSpec =Nmax_lev-1,2,-1 ! DO MK3 = 1,Nblocks ! find blocks adjasent to filled blocks DO MK2 = 1,Nblocks ! mark them with kSpec =Nspecies-1, -2, ..1 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 If(iMask(i,j,k).lt.kSpec)iMask(i,j,k)=kSpec ! change mask ENDDO ENDDO ENDDO EndIf ENDDO ENDDO ENDDO Enddo 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 C------------------------------------------------ C Make a realization of spectrum of perturbations: C ALPHA = normalization factor for displacements C NRAND = seed for random numbers C Lstart = start on level Lstart (=1 or 2) SUBROUTINE SPECTR(NRAND,Lstart) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' INCLUDE 'luxury.h' REAL*8 SUMM, SIJ,WI3,WJ3,WK3,WD,TS,TRX Dimension GG(NROW,NROW) ! temp array for random seeds C Set spectrum If(Lstart.gt.Levels)Then write (*,*) ' Wrong range of levels in SPECTR: start/total=', & Lstart,Levels stop Endif write (*,*) ' range of levels in SPECTR: start/total=', & Lstart,Levels iRandom =INT(extras(80)+1.e-5) ! set type of Random Numbers write (*,*)' Random numbers=',iRandom,' (0=Randd/1=lux)' C$omp parallel do default(shared) C$omp& private(MK3,MK2,MK1) DO MK3 = 1,NROW DO MK2 = 1,NROW DO MK1 = 1,NROW GRX(MK1,MK2,MK3) =0. GRY(MK1,MK2,MK3) =0. GRZ(MK1,MK2,MK3) =0. ENDDO ENDDO ENDDO C--------------------------- first level If(Lstart.eq.1)Then C$omp parallel do default(shared) C$omp& private(MK3,MK2,MK1) DO MK3 = 1,NROW2 DO MK2 = 1,NROW2 DO MK1 = 1,NROW2 GRX2(MK1,MK2,MK3) =0. GRY2(MK1,MK2,MK3) =0. GRZ2(MK1,MK2,MK3) =0. ENDDO ENDDO ENDDO SUMM = 0. ! set spectrum for first level DO k = 1,NROW2 ksign = -1 kz = k +NSPEC2 k3 = k - 1 If(k3.GT.NSPEC2)k3=k3-NSPEC2 IF(kz.gt.NROW2)THEN kz =kz -NROW2 ksign =1 ENDIF WK3= K3**2 c write (*,*) ' K=',k,' K3=',K3,' Nspec=',NSPEC2 DO j = 1,NROW2 jsign = -1 jz = j +NSPEC2 j3 = j - 1 If(j3.GT.NSPEC2)j3=j3-NSPEC2 IF(jz.gt.NROW2)THEN jz =jz -NROW2 jsign =1 ENDIF WJ3= j3**2 DO i = 1,NROW2 isign = -1 iz = i +NSPEC2 i3 = i - 1 If(i3.GT.NSPEC2)i3=i3-NSPEC2 IF(iz.gt.NROW2)THEN iz =iz -NROW2 isign =1 ENDIF WI3= i3**2 IF(i3+j3+k3.EQ.0)THEN GRX2(1,1,1) = 0. GRY2(1,1,1) = 0. GRZ2(1,1,1) = 0. ELSE WD = (WI3 + WJ3 + WK3) WK = SQRT(WD) If(iRandom.eq.0)Then TS = TRUNF(WK) * GAUSS(NRAND) Else TS = TRUNF(WK) * GAUSS3() EndIf TRX = TS / WD IF(i.EQ.NSPEC2+1.or. & j.EQ.NSPEC2+1.or. & k.EQ.NSPEC2+1)TRX =0 GRX2(iz, j, k) = i3 * TRX *isign GRY2(i,jz, k) = j3 * TRX *jsign GRZ2(i, j,kz) = k3 * TRX *ksign SUMM = SUMM + TS**2 ENDIF ENDDO ENDDO ENDDO IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0' write (*,*) ' ==== SPECTR: SUMM=',SUMM,' Level=',Lstart ALPHA(Lstart) = AMPLT_CR(Lstart) / SQRT(SUMM) *sqrt(8.) write (*,*) ' ALPHA(Lstart)=',ALPHA(Lstart),AMPLT_CR(Lstart) DO MK3 = 1,NROW2 DO MK2 = 1,NROW2 DO MK1 = 1,NROW2 GRX2(MK1,MK2,MK3) =GRX2(MK1,MK2,MK3)*ALPHA(Lstart) GRY2(MK1,MK2,MK3) =GRY2(MK1,MK2,MK3)*ALPHA(Lstart) GRZ2(MK1,MK2,MK3) =GRZ2(MK1,MK2,MK3)*ALPHA(Lstart) ENDDO ENDDO ENDDO EndIf ! end first level C------------------------- second level If(Levels.eq.1)Return SUMM = 0. DO k = 1,NROW ksign = -1 kz = k +NSPEC k3 = k - 1 If(k3.GT.NSPEC)k3=k3-NSPEC IF(kz.gt.NROW)THEN kz =kz -NROW ksign =1 ENDIF WK3= K3**2 DO j = 1,NROW DO i = 1,NROW If(iRandom.eq.0)Then GG(i,j) = GAUSS(NRAND) Else GG(i,j) = GAUSS3() c If(j.eq.1.and. i.le.20)write(30,'(2i4,2i10,3g12.4)') c & k,i,i24,j24,GG(i,j), c & seeds(i24),seeds(j24) EndIf ENDDO ENDDO SIJ =0. C$omp parallel do default(shared) C$omp& private(j,jsign,jz,j3,WJ3,i,isign,iz,i3,WI3,WD,WK,TS,TRX) C$omp& REDUCTION(+:SIJ) 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 WJ3= j3**2 DO i = 1,NROW isign = -1 iz = i +NSPEC i3 = i - 1 If(i3.GT.NSPEC)i3=i3-NSPEC IF(iz.gt.NROW)THEN iz =iz -NROW isign =1 ENDIF WI3= i3**2 IF(i3+j3+k3.EQ.0)THEN GRX(1,1,1) = 0. GRY(1,1,1) = 0. GRZ(1,1,1) = 0. ELSE WD = (WI3 + WJ3 + WK3) WK = SQRT(WD) TS = TRUNF(WK) * GG(i,j) c If(j.eq.1.and. i.le.2)write(30,'(i4,2i10,3g12.4)') c & k,i24,j24,GG(i,j), c & seeds(i24),seeds(j24) TRX = TS / WD IF(i.EQ.NSPEC+1.or. & j.EQ.NSPEC+1.or. & k.EQ.NSPEC+1)TRX =0 GRX(iz, j, k) = i3 * TRX *isign GRY( i,jz, k) = j3 * TRX *jsign GRZ( i, j,kz) = k3 * TRX *ksign SIJ = SIJ + TS**2 c write(*,'(4i5,3g12.4)')i,j,iz,jz,WD,TRX,TS/TRUNF(WK) c If(k.eq.2)STOP ENDIF ENDDO ENDDO SUMM = SUMM +SIJ ENDDO IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0' write (*,*) ' ==== SPECTR: SUMM=',SUMM,' Level=',Levels ALPHA(Levels) = AMPLT_CR(Levels)/ SQRT(SUMM) *sqrt(8.) DO MK3 = 1,NROW DO MK2 = 1,NROW DO MK1 = 1,NROW GRX(MK1,MK2,MK3) =GRX(MK1,MK2,MK3)*ALPHA(Levels) GRY(MK1,MK2,MK3) =GRY(MK1,MK2,MK3)*ALPHA(Levels) GRZ(MK1,MK2,MK3) =GRZ(MK1,MK2,MK3)*ALPHA(Levels) ENDDO ENDDO ENDDO c write(50,*) 'XXXX' c write(50,'(12g12.4)') c & ((GRx(i,1,k),i=1,24),k=1,2) RETURN END C------------------------------------------------ C Make spectrum for Constrained Realizations C C SUBROUTINE SPectr_CR(NRAND) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Par(6),ns,qqscaleb,QSCALE REAL*8 SUMM, WI3,WJ3,WK3,WD,TS,TRX Dimension DENS_CR(NROW2,NROW2,NROW2) C Set spectrum c R_f = 0.25/hubble ! correction for gaussian filter; R_f c RR_f = QSCALE*R_f ! dimensionless gaussian radius c R_cor = RR_f**2/2. c write (*,*) ' Apply correction for filter radius=',R_f,RR_f SUMM = 0. WD =0. open(40,file='dens_cr.dat',status='unknown') read(40,*) DENS_CR Do k=1,NROW2 Do j=1,NROW2 Do i=1,NROW2 SUMM = SUMM +dens_cr(i,j,k) WD = WD +dens_cr(i,j,k)**2 EndDo EndDo EndDo sigma = sqrt(MAX(1.d-10,WD/NROW2**3)) write (*,*) ' Average density=',SUMM/NROW2**3 write (*,*) ' RMS density=',sigma IFOUR = INT(ALOG(FLOAT(NROW2))/ALOG(2.)+0.5) Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW2 CALL VECTOR_dens(IFOUR,dens_cr) C$omp parallel do default(shared) C$omp& private(MK3,MK2,MK1) DO MK3 = 1,NROW2 DO MK2 = 1,NROW2 DO MK1 = 1,NROW2 GRX2(MK1,MK2,MK3) =0. GRY2(MK1,MK2,MK3) =0. GRZ2(MK1,MK2,MK3) =0. ENDDO ENDDO ENDDO DO k = 1,NROW2 ksign = -1 kz = k +NSPEC2 k3 = k - 1 If(k3.GT.NSPEC2)k3=k3-NSPEC2 IF(kz.gt.NROW2)THEN kz =kz -NROW2 ksign =1 ENDIF WK3= K3**2 c write (*,*) ' K=',k,' K3=',K3,' Nspec=',NSPEC2 DO j = 1,NROW2 jsign = -1 jz = j +NSPEC2 j3 = j - 1 If(j3.GT.NSPEC2)j3=j3-NSPEC2 IF(jz.gt.NROW2)THEN jz =jz -NROW2 jsign =1 ENDIF WJ3= j3**2 DO i = 1,NROW2 isign = -1 iz = i +NSPEC2 i3 = i - 1 If(i3.GT.NSPEC2)i3=i3-NSPEC2 IF(iz.gt.NROW2)THEN iz =iz -NROW2 isign =1 ENDIF WI3= i3**2 IF(i3+j3+k3.EQ.0)THEN GRX2(1,1,1) = 0. GRY2(1,1,1) = 0. GRZ2(1,1,1) = 0. ELSE WD = WI3 + WJ3 + WK3 c WK = SQRT(WD) c TS = -dens_cr(i,j,k) * exp(WD*R_cor) TS = -dens_cr(i,j,k) TRX = TS / WD IF(i.EQ.NSPEC2+1.or. & j.EQ.NSPEC2+1.or. & k.EQ.NSPEC2+1)TRX =0 GRX2(iz, j, k) = i3 * TRX *isign GRY2( i,jz, k) = j3 * TRX *jsign GRZ2( i, j,kz) = k3 *TRX *ksign SUMM = SUMM + TS**2 ENDIF ENDDO ENDDO ENDDO IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0' write (*,*) ' ==== SPECTR === SUMM=',SUMM ALPHA(1) = AMPLT_CR(1) / SQRT(SUMM) *sqrt(8.) *1.08 write (*,*) ' ALPHA(1) =',ALPHA(1) ,AMPLT_CR(1) DO MK3 = 1,NROW2 DO MK2 = 1,NROW2 DO MK1 = 1,NROW2 GRX2(MK1,MK2,MK3) =GRX2(MK1,MK2,MK3)*ALPHA(1) GRY2(MK1,MK2,MK3) =GRY2(MK1,MK2,MK3)*ALPHA(1) GRZ2(MK1,MK2,MK3) =GRZ2(MK1,MK2,MK3)*ALPHA(1) ENDDO ENDDO ENDDO RETURN END C------------------------------------------------ C Copy harmonics, which C were included at previous level C of resolution C Lstart = start on level Lstart (=1 or 2) SUBROUTINE RELOCATE_SP C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' c If(Levels.lt.2)Then c write (*,*) ' Wrong range of levels in RELOCATE_SP: =', c & Levels c Return c Endif J_cos = NROW/2 +1 ! limit for cosine harmonics J_sin_1 = NROW/2+2 ! limits for sin harmonics J2_cos = NROW2/2 +1 ! limit for cosine harmonics J2_sin_1 = NROW2/2+2 ! limits for sin harmonics C$omp parallel do default(shared) C$omp& private(K,J,I,Ishift,Jshift,Kshift) DO K = 1,NROW2 Kshift = 0 If(K.ge.J2_sin_1)Kshift =J_cos -J2_cos DO J = 1,NROW2 Jshift = 0 If(J.ge.J2_sin_1)Jshift =J_cos -J2_cos DO I = 1,NROW2 Ishift = 0 If(I.ge.J2_sin_1)Ishift =J_cos -J2_cos GRX(I+Ishift,J+Jshift,K+Kshift) =GRX2(I,J,K) GRY(I+Ishift,J+Jshift,K+Kshift) =GRY2(I,J,K) GRZ(I+Ishift,J+Jshift,K+Kshift) =GRZ2(I,J,K) ENDDO ENDDO ENDDO c do k=1,4 c write(30,*) ' k=',k c do j=1,10 c write(30,*)' j=',j c write(30,'((10g12.4))')(GRX(i,j,k),i=1,NROW) c write(30,*) '---- y' c write(30,'((10g12.4))')(GRY(i,j,k),i=1,NROW) c write(30,*) '---- z' c write(30,'((10g12.4))')(GRZ(i,j,k),i=1,NROW) c EndDo c EndDo RETURN END C------------------------------------------------ C Define coordinates and velosities for C particles with resolution = Indx C Indx =1 = high resolution C Indx =2 = medium resolution C Icurrent = number of particles SUBROUTINE BLOCKS(XCONS,VCONS,Indx,Icurrent) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' Dimension nLevel(10) Real*8 xDispl,yDispl,zDispl Do i=1, Max_lev_abs ! number of particles at nLevel(i)=0 ! different levels Enddo Boff = Lblock/2 +0.5 xDispl =0. yDispl =0. zDispl =0. c write (*,*) ' Blocks:', Nblocks,Lblock 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) C---------------------------------------- c M1 + Boff = lagrangian coordinate of 'mean' particle c Lblock**3 = total number of averaged cells If(iM.eq.Indx)Then If(iM.eq.1)Then ! Low resolution: one particle DX =0. DY =0. DZ =0. kp=(M3+1 +M3+Lblock)/2 jp=(M2+1 +M2+Lblock)/2 ip=(M1+1 +M1+Lblock)/2 DX = GRX(ip,jp,kp) DY = GRY(ip,jp,kp) DZ = GRZ(ip,jp,kp) Q3 = QFACT*(kp-1.) +1. ! scale them to box size Q2 = QFACT*(jp-1.) +1. Q1 = QFACT*(ip-1.) +1. Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & MK1,MK2,MK3 STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - XCONS*DX + 1.E-4 +0.5 YPt(Icurrent) = Q2 - XCONS*DY + 1.E-4 +0.5 ZPt(Icurrent) = Q3 - XCONS*DZ + 1.E-4 +0.5 VXt(Icurrent) = VCONS*DX VYt(Icurrent) = VCONS*DY VZt(Icurrent) = VCONS*DZ XDispl = xDispl +DX**2 yDispl = yDispl +DY**2 zDispl = zDispl +DZ**2 Else C------------------------------------------ If(iM.eq.Max_lev_abs)Then ! Highest resolution: Do kp=M3+1,M3+Lblock Do jp=M2+1,M2+Lblock Do ip=M1+1,M1+Lblock DX = GRX(ip,jp,kp) ! find displacements DY = GRY(ip,jp,kp) DZ = GRZ(ip,jp,kp) Q3 = QFACT*(kp-1.) +1. ! scale them to box size Q2 = QFACT*(jp-1.) +1. Q1 = QFACT*(ip-1.) +1. Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & MK1,MK2,MK3 STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - XCONS*DX + 1.E-4 +0.5 YPt(Icurrent) = Q2 - XCONS*DY + 1.E-4 +0.5 ZPt(Icurrent) = Q3 - XCONS*DZ + 1.E-4 +0.5 VXt(Icurrent) = VCONS*DX VYt(Icurrent) = VCONS*DY VZt(Icurrent) = VCONS*DZ XDispl = xDispl +DX**2 yDispl = yDispl +DY**2 zDispl = zDispl +DZ**2 EndDo EndDo EndDo Else C------------------------------------------- c Medium resolution: c KBlock = 2**(iM-1) = number of particles in 1D per block c Nbcells= Lblock/KBlock = number of cells per particle c (MK-1)*Lblock +(i-1)*Nbcells = offset of cells in 1D c (MK-1)*Lblock +(i-1)*Nbcells +Nbcells/2 +1/2 = lagrangian c coordinate of particle c 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 c write (*,*)' Lblock=',Lblock,' Nblocks=',Nblocks, c + ' Npart/block/1D=',KBlock 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 X Q2 = QFACT*(j1+Xlag-1.) +1. Do k=1,KBlock k1 =M3+(k-1)*Nbcells ! offset for cells in X Q3 = QFACT*(k1+Xlag-1.) +1. Nncells =0 kp=(k1+1 +k1+Nbcells)/2 ! jp=(j1+1 +j1+Nbcells)/2 ip=(i1+1 +i1+Nbcells)/2 Nncells =Nncells +1 If(ip.gt.NROW.or.ip.lt.1)Then write (*,*) ' error' endif If(jp.gt.NROW.or.jp.lt.1)Then write (*,*) ' error' endif If(kp.gt.NROW.or.kp.lt.1)Then write (*,*) ' error' endif DX = GRX(ip,jp,kp) DY = GRY(ip,jp,kp) DZ = GRZ(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 XPt(Icurrent) = Q1 - XCONS*DX + 1.E-4 +0.5 YPt(Icurrent) = Q2 - XCONS*DY + 1.E-4 +0.5 ZPt(Icurrent) = Q3 - XCONS*DZ + 1.E-4 +0.5 VXt(Icurrent) = VCONS*DX VYt(Icurrent) = VCONS*DY VZt(Icurrent) = VCONS*DZ XDispl = xDispl +DX**2 yDispl = yDispl +DY**2 zDispl = zDispl +DZ**2 EndDo EndDo EndDo EndIf EndIf EndIf ! end iM = Indx EndDo EndDo EndDo If(nLevel(indx).gt.0)Then Nspecies = Nspecies +1 lspecies(Nspecies) =Icurrent wspecies(Nspecies) =PARTW* & Lblock**3/8**(Indx-1) Write (*, '(10x,a20,3g13.5)') ' RMS 3d diplacement=', & XCONS*sqrt (xDispl/nLevel (indx) ), & XCONS*sqrt (yDispl/nLevel (indx) ), & XCONS*sqrt (zDispl/nLevel (indx) ) write(*,'(" Species=",i3," Particles=",i11," on levels=",10i11)') & Nspecies,Icurrent,(nLevel(i),i=1,Max_lev_abs) endif RETURN END C------------------------------------------------ C Define coordinates and velosities for C particles with resolution = Indx C Indx =1 = high resolution C Indx =2 = medium resolution C Icurrent = number of particles SUBROUTINE BLOCKSaver(XCONS,VCONS,Indx,Icurrent) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' Dimension nLevel(10) Do i=1, Max_lev_abs ! number of particles at nLevel(i)=0 ! different levels Enddo Boff = Lblock/2 +0.5 c write (*,*) ' Blocks:', Nblocks,Lblock 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) C---------------------------------------- c M1 + Boff = lagrangian coordinate of 'mean' particle c Lblock**3 = total number of averaged cells If(iM.eq.Indx)Then If(iM.eq.1)Then ! Low resolution: one particle DX =0. DY =0. DZ =0. Do kp=M3+1,M3+Lblock ! sum of displacements Do jp=M2+1,M2+Lblock Do ip=M1+1,M1+Lblock DX = DX +GRX(ip,jp,kp) DY = DY +GRY(ip,jp,kp) DZ = DZ +GRZ(ip,jp,kp) EndDo EndDo EndDo DX = DX/Lblock**3 ! mean displacement DY = DY/Lblock**3 DZ = DZ/Lblock**3 Q3 = QFACT*(M3+Boff-1.) +1. ! scale them to box size Q2 = QFACT*(M2+Boff-1.) +1. Q1 = QFACT*(M1+Boff-1.) +1. Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & MK1,MK2,MK3 STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - XCONS*DX + 1.E-4 +0.5 YPt(Icurrent) = Q2 - XCONS*DY + 1.E-4 +0.5 ZPt(Icurrent) = Q3 - XCONS*DZ + 1.E-4 +0.5 VXt(Icurrent) = VCONS*DX VYt(Icurrent) = VCONS*DY VZt(Icurrent) = VCONS*DZ Else C------------------------------------------ If(iM.eq.Max_lev_abs)Then ! Highest resolution: Do kp=M3+1,M3+Lblock Do jp=M2+1,M2+Lblock Do ip=M1+1,M1+Lblock DX = GRX(ip,jp,kp) ! find displacements DY = GRY(ip,jp,kp) DZ = GRZ(ip,jp,kp) Q3 = QFACT*(kp-1.) +1. ! scale them to box size Q2 = QFACT*(jp-1.) +1. Q1 = QFACT*(ip-1.) +1. Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & MK1,MK2,MK3 STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - XCONS*DX + 1.E-4 +0.5 YPt(Icurrent) = Q2 - XCONS*DY + 1.E-4 +0.5 ZPt(Icurrent) = Q3 - XCONS*DZ + 1.E-4 +0.5 VXt(Icurrent) = VCONS*DX VYt(Icurrent) = VCONS*DY VZt(Icurrent) = VCONS*DZ EndDo EndDo EndDo Else C------------------------------------------- c Medium resolution: c KBlock = 2**(iM-1) = number of particles in 1D per block c Nbcells= Lblock/KBlock = number of cells per particle c (MK-1)*Lblock +(i-1)*Nbcells = offset of cells in 1D c (MK-1)*Lblock +(i-1)*Nbcells +Nbcells/2 +1/2 = lagrangian c coordinate of particle c 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 c write (*,*)' Lblock=',Lblock,' Nblocks=',Nblocks, c + ' Npart/block/1D=',KBlock 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 X Q2 = QFACT*(j1+Xlag-1.) +1. Do k=1,KBlock k1 =M3+(k-1)*Nbcells ! offset for cells in X Q3 = QFACT*(k1+Xlag-1.) +1. DX =0. DY =0. DZ =0. Nncells =0 Do kp=k1+1,k1+Nbcells ! sum of displacements Do jp=j1+1,j1+Nbcells Do ip=i1+1,i1+Nbcells Nncells =Nncells +1 If(ip.gt.NROW.or.ip.lt.1)Then write (*,*) ' error' endif If(jp.gt.NROW.or.jp.lt.1)Then write (*,*) ' error' endif If(kp.gt.NROW.or.kp.lt.1)Then write (*,*) ' error' endif DX = DX +GRX(ip,jp,kp) DY = DY +GRY(ip,jp,kp) DZ = DZ +GRZ(ip,jp,kp) EndDo EndDo EndDo DX = DX/ Nbcells3 ! mean displacement DY = DY/ Nbcells3 DZ = DZ/ Nbcells3 If(Nncells.ne.Nbcells3)Then write (*,*) ' error: wrong number of cells=', + Nncells,Nbcells3 stop endif Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & MK1,MK2,MK3 STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - XCONS*DX + 1.E-4 +0.5 YPt(Icurrent) = Q2 - XCONS*DY + 1.E-4 +0.5 ZPt(Icurrent) = Q3 - XCONS*DZ + 1.E-4 +0.5 VXt(Icurrent) = VCONS*DX VYt(Icurrent) = VCONS*DY VZt(Icurrent) = VCONS*DZ EndDo EndDo EndDo EndIf EndIf EndIf ! end iM = Indx EndDo EndDo EndDo If(nLevel(indx).gt.0)Then Nspecies = Nspecies +1 lspecies(Nspecies) =Icurrent wspecies(Nspecies) =PARTW* & Lblock**3/8**(Indx-1) endif write(*,'(" Species=",i3," Particles=",i11," on levels=",10i11)') & Nspecies,Icurrent,(nLevel(i),i=1,Max_lev_abs) RETURN END C------------------------------------------------ SUBROUTINE SETtest(K,NRAND) C place particles with small velocities in a sphere C Radius = radius of the sphere in zero-level cell units C Xcentr = position of the sphere C C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2 QFACT = FLOAT(NGRID)/FLOAT(NROW) XMAX = FLOAT(NGRID) + 1. XSHF = FLOAT(NGRID) Xcentr = NGRID/2. Vss =1. Radius = 5. write (*,*) ' Xcentr=',Xcentr DO J = 1,NROW DO I = 1,NROW IJ = I +(J-1)*NROW DX = Radius*(2.*RANDd(NRAND)-1.) DY = Radius*(2.*RANDd(NRAND)-1.) DZ = Radius*(2.*RANDd(NRAND)-1.) c DX = 0.05*I c DY = 0.05*J c DZ = 0.05*K C the factor 0.5 in the following 3 lines places unperturbed c particle in the middle of a cell. This reduces the shot noise C Small number 1.E-4 is to avoid placing some particles C exactly in cell nodes. Does not make much differentce. XPAR(IJ) = Xcentr+ DX + 1.E-4 +0.5 YPAR(IJ) = Xcentr+ DY + 1.E-4 +0.5 ZPAR(IJ) = Xcentr+ DZ + 1.E-4 +0.5 VX(IJ) = Vss*GAUSS(NRAND) VY(IJ) = Vss*GAUSS(NRAND) VZ(IJ) = Vss*GAUSS(NRAND) C Periodical boundary conditions IF(XPAR(IJ).GT.XMAX) XPAR(IJ)=XPAR(IJ)-XSHF IF(XPAR(IJ).LE.1.) XPAR(IJ)=XPAR(IJ)+XSHF IF(YPAR(IJ).GT.XMAX) YPAR(IJ)=YPAR(IJ)-XSHF IF(YPAR(IJ).LE.1.) YPAR(IJ)=YPAR(IJ)+XSHF IF(ZPAR(IJ).GT.XMAX) ZPAR(IJ)=ZPAR(IJ)-XSHF IF(ZPAR(IJ).LE.1.) ZPAR(IJ)=ZPAR(IJ)+XSHF SX = SX + DX SY = SY + DY SZ = SZ + DZ SX2 = SX2 + DX**2 SY2 = SY2 + DY**2 SZ2 = SZ2 + DZ**2 SKINE = SKINE + VX(IJ)**2 +VY(IJ)**2 +VZ(IJ)**2 ENDDO ENDDO RETURN END C------------------------------------------------ C Define coordinates and velosities for C all particles of current row (given q2) C x = q - amplt*b(t)/b(t0)*S(q1,q2) C momentum = a**2*deriv x C C Twidled variables: C x= q - {1/(2 pi Ngrid)}{Alpha*a(t)/a(0)}*S C P= - {1/(2 pi Ngrid)}{Alpha*a(t)**1.5/a(0)}*S C C q(i) = Ngrid/Nrow*(i-1) + 1 C x = q - xcons*S C P = vcons*S SUBROUTINE SETXV(XCONS1,VCONS1,K) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2 QFACT = FLOAT(NGRID)/FLOAT(NROW) XMAX = FLOAT(NGRID) + 1. XSHF = FLOAT(NGRID) Q3 = QFACT*(K-1.) +1. DO J = 1,NROW Q2 = QFACT*(J-1.) +1. DO I = 1,NROW Q1 = QFACT*(I-1.) +1. IJ = I +(J-1)*NROW DX = XCONS1*GRX(I,J,K) DY = XCONS1*GRY(I,J,K) DZ = XCONS1*GRZ(I,J,K) C the factor 0.5 in the following 3 lines places unperturbed c particle in the middle of a cell. This reduces the shot noise C Small number 1.E-4 is to avoid placing some particles C exactly in cell nodes. Does not make much differentce. XPAR(IJ) = Q1 - DX + 1.E-4 +0.5 YPAR(IJ) = Q2 - DY + 1.E-4 +0.5 ZPAR(IJ) = Q3 - DZ + 1.E-4 +0.5 VX(IJ) = VCONS1*GRX(I,J,K) VY(IJ) = VCONS1*GRY(I,J,K) VZ(IJ) = VCONS1*GRZ(I,J,K) C Periodical boundary conditions IF(XPAR(IJ).GT.XMAX) XPAR(IJ)=XPAR(IJ)-XSHF IF(XPAR(IJ).LE.1.) XPAR(IJ)=XPAR(IJ)+XSHF IF(YPAR(IJ).GT.XMAX) YPAR(IJ)=YPAR(IJ)-XSHF IF(YPAR(IJ).LE.1.) YPAR(IJ)=YPAR(IJ)+XSHF IF(ZPAR(IJ).GT.XMAX) ZPAR(IJ)=ZPAR(IJ)-XSHF IF(ZPAR(IJ).LE.1.) ZPAR(IJ)=ZPAR(IJ)+XSHF SX = SX + DX SY = SY + DY SZ = SZ + DZ SX2 = SX2 + DX**2 SY2 = SY2 + DY**2 SZ2 = SZ2 + DZ**2 SKINE = SKINE + VX(IJ)**2 +VY(IJ)**2 +VZ(IJ)**2 ENDDO ENDDO RETURN END C------------------------------------------------ C FFT of the spectrum SUBROUTINE VECTOR_dens(IFOUR,dens_cr) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67),dens_cr(NROW2,NROW2,NROW2) jb1 =3 jq =IFOUR CALL SETF67(jb1,jq,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) write (*,*) ' SETF67 is done' C x-component C$omp parallel do default(shared) C$omp& private(k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO K=1,NROW2 DO J=1,NROW2 DO I=1,NROW2 Uf(I) =dens_cr(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW2 dens_cr(I,J,K) = Vf(I) ENDDO ENDDO ENDDO write (*,*) ' 1st is done' C$omp parallel do default(shared) C$omp& private(k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO K=1,NROW2 DO I=1,NROW2 DO J=1,NROW2 Uf(J) =dens_cr(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW2 dens_cr(I,J,K) = Vf(J) ENDDO ENDDO ENDDO C$omp parallel do default(shared) C$omp& private(k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO J=1,NROW2 DO I=1,NROW2 DO K=1,NROW2 Uf(K) =dens_cr(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW2 dens_cr(I,J,K) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' dens_cr is done' Return END C------------------------------------------------ C FFT of the spectrum SUBROUTINE VECTOR C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) c write(50,*) '-------' c write(50,'(12g12.4)') c & ((GRx(i,1,k),i=1,24),k=1,2) c write(50,'(12g12.4)') c & ((GRx(i,2,k),i=1,24),k=1,2) IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5) jb1 =4 jq =IFOUR CALL SETF67(jb1,jq,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) C x-component C$omp parallel do default(shared) C$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) =GRX(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW GRX(I,J,K) = Vf(I) ENDDO ENDDO ENDDO write (*,*) ' 1st is done Level=' C$omp parallel do default(shared) C$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) =GRX(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW GRX(I,J,K) = Vf(J) ENDDO ENDDO ENDDO c write(50,*) '======Ydone' c write(50,'(12g12.4)') c & ((GRX(i,j,1),i=1,24),j=1,4) C$omp parallel do default(shared) C$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) =GRX(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW GRX(I,J,K) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' x is done ' C Y-component C$omp parallel do default(shared) C$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) =GRY(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW GRY(I,J,K) = Vf(I) ENDDO ENDDO ENDDO C$omp parallel do default(shared) C$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) =GRY(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW GRY(I,J,K) = Vf(J) ENDDO ENDDO ENDDO C$omp parallel do default(shared) C$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) =GRY(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW GRY(I,J,K) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' y is done' C Z-component C$omp parallel do default(shared) C$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) =GRZ(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW GRZ(I,J,K) = Vf(I) ENDDO ENDDO ENDDO C$omp parallel do default(shared) C$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) =GRZ(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW GRZ(I,J,K) = Vf(J) ENDDO ENDDO ENDDO C$omp parallel do default(shared) C$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) =GRZ(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW GRZ(I,J,K) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' z is done' c write(50,*) '-------' c write(50,'(12g12.4)') c & ((GRx(i,1,k),i=1,24),k=1,2) c write(50,'(12g12.4)') c & ((GRx(i,2,k),i=1,24),k=1,2) RETURN END C------------------------------------------ c Rearange species so that the highest resolution c is fisrt, lowest resolution - last. c Print statistics of spicies Subroutine Species_Change(Wtotal) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' 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,i11,T26,i11,T40,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 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' Character Hd*5,Tail*4,Nm*40 Hd ='PMcrs' Tail='.DAT' 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=10,end=10) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, + NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 + ,Ocurv,extras Ocurv =0. 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=',I6,/ + 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 10 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT', + STATUS='UNKNOWN',RECL=NACCES) REWIND 9 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 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 cC-------------------------------------------------- cC Set Weights of particles for cC fast access c SUBROUTINE SetWeight cC-------------------------------------------------- c INCLUDE 'PMparameters.h' c If(Nspecies.eq.0)Then ! old constant weights c N_particles =lspecies(1) c Do i=1,N_particles c iWeight(i) =PARTW c EndDo c Else c N_particles =lspecies(Nspecies) c jstart =1 c Do j=1,Nspecies c jend =lspecies(j) c Do k=jstart ,jend c iWeight(k) =wspecies(j) c EndDo c jstart =jend c EndDo c EndIf c write (*,*) ' Set Weights for ',N_particles,' particles' c RETURN c 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------------------------------------------------ C random number generator FUNCTION RANDd(M) C------------------------------------------------ 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 C-------------------------------------- C normal random numbers FUNCTION GAUSS(M) C-------------------------------------- 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 C-------------------------------------- C normal random numbers FUNCTION GAUSS3() C Uses ranlux for homogenous rand numbers C Uses Box-Muller inversion for gaussian C Excellent quality and speed c N= 100000000 GAUSS3+ranlux ------ c sigma= 0.99997 mean= 0.18692E-03 c sigma Frac(>sig) Frac(<-sig) True n(>) n(<) c 1.00 0.1587 0.1586 0.1587 15869095 15857724 c 1.50 0.6682E-01 0.6680E-01 0.6681E-01 6681652 6679768 c 2.00 0.2277E-01 0.2273E-01 0.2275E-01 2276651 2273197 c 2.50 0.6222E-02 0.6206E-02 0.6210E-02 622190 620610 c 3.00 0.1356E-02 0.1347E-02 0.1350E-02 135628 134674 c 4.00 0.3242E-04 0.3145E-04 0.3170E-04 3242 3145 c C-------------------------------------- DIMENSION RanNum(2) DATA iFlag/0/ SAVE gSet 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 c-------------------------------------------------------------------- ! 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. c--------------------------------------------- SUBROUTINE ranlux(rvec, lenv) c returns a vector RVEC of LEN c 32-bit random floating point numbers between c zero (not included) and one (also not incl.). c Next call to ranlux gives next LEN of random numbers c--------------------------------------------- INTEGER lenv REAL rvec(lenv) INTEGER iseeds(24) INCLUDE 'luxury.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 notyet = .false. jseed = jsdflt inseed = jseed WRITE (6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ', jseed luxlev = lxdflt nskip = ndskip(luxlev) lp = nskip + 24 in24 = 0 kount = 0 mkount = 0 WRITE (6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL= ', luxlev, & ' p =', lp 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) END DO twom12 = twom24 * 4096. DO i = 1, 24 seeds(i) = REAL(iseeds(i)) * twom24 next(i) = i - 1 END DO next(1) = 24 i24 = 24 j24 = 10 carry = 0. IF (seeds(24).EQ.0.) carry = twom24 END IF ! 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. END IF seeds(i24) = uni i24 = next(i24) j24 = next(j24) rvec(ivec) = uni ! small numbers (with less than 12 "significant" bits) are "padded". 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 END IF ! 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. END IF seeds(i24) = uni i24 = next(i24) j24 = next(j24) END DO END IF END DO kount = kount + lenv IF (kount.GE.igiga) THEN mkount = mkount + 1 kount = kount - igiga END IF RETURN END ! SUBROUTINE ranlux c--------------------------------------------- c Subroutine to initialize from one or three integers SUBROUTINE rluxgo(lux, ins, k1, k2) c initializes the generator from c one 32-bit integer INT and sets Luxury Level LUX c which is integer between zero and MAXLEV, or if c LUX .GT. 24, it sets p=LUX directly. K1 and K2 c should be set to zero unless restarting at a break c point given by output of RLUXAT c--------------------------------------------- INCLUDE 'luxury.h' INTEGER iseeds(24) IF (lux.LT.0) THEN luxlev = lxdflt ELSE IF (lux.LE.maxlev) THEN luxlev = lux ELSE IF (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 END DO END IF 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 END IF 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' END IF 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) END DO twom12 = twom24 * 4096. DO i = 1, 24 seeds(i) = REAL(iseeds(i)) * twom24 next(i) = i - 1 END DO 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. END IF seeds(i24) = uni i24 = next(i24) j24 = next(j24) END DO END DO ! 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) END IF ! 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 END IF END IF RETURN END ! SUBROUTINE rluxgo C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End C-------------------------------------- Fourier Transform SUBROUTINE SETF67(JBC1,JQ1,jbc,jp,jsl,ll1,k2,k3,k4,k7, & Qi,jndx,Uf,Vf) c ----------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) PI=DATAN(1.D0)*4. JBC=JBC1 JQ=JQ1 IF(JBC.LT.3) GO TO 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) GO TO 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) GO TO 102 GO TO 103 104 RETURN END C----------------------------------------- SUBROUTINE TFOLD(IS,L,ZZZ,k2) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' 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 CONTINUE RETURN END C------------------------------------ SUBROUTINE NEG(I1,I3,I2,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) DO 100 K=I1,I3,I2 Vf(K)=-Vf(K) 100 CONTINUE RETURN END C------------------------------------- SUBROUTINE REVNEG(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) integer jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7 DO 100 I=1,K7 J=K3+1+I K=K4+1-I A=Vf(J) Vf(J)=-Vf(K) Vf(K)=-A 100 CONTINUE RETURN END C------------------------------------- SUBROUTINE ZEERO(L,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) integer jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7 DO 100 I=1,K2 LI=L+I Uf(LI-1)=0.0 100 CONTINUE RETURN END C------------------------------------------ SUBROUTINE TFOLD1(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) integer jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7 II=K2-1 DO 100 I=1,II I1=JSL+I I2=LL1-I A=Uf(I1) Uf(I1)=A+Uf(I2) Uf(I2)=A-Uf(I2) 100 CONTINUE RETURN END C------------------------------------- SUBROUTINE KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' 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) GO TO 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) GO TO 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) GO TO 210 I=I+1 300 IS1=JSL IC1=LL1 JS1=JS1/2 IF(JS1.EQ.1) GO TO 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) GO TO 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) GO TO 330 I=I+1 IF(IS1.EQ.J5) GO TO 300 GO TO 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) GO TO 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) GO TO 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) GO TO 400 RETURN END c ------------------------------------------------------ SUBROUTINE FOUR67(JBC1,JQ1,jp1,jsl1,ll11,k21,k71, & Qi,jndx,Uf,Vf) c ------------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' 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 GO TO (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 GO TO 103 102 K3=K3/2 JQ=JQ-1 103 N5=K3/4 K7=K3/2 N11=3*K7 K31=K3+1 GO TO(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 C--------------------------------------------- C Write current data to disk/tape C for PMstartM - multiple masses SUBROUTINE WriteData(Vscale,AexpV,Wtotal) C---------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR2.h' COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2 C 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) write(30,*) ' YYY' ! test data write(30,'(10g12.5)') (YPt(i),i=1,lspecies(Nspecies),5000) write(30,*) ' VVV' write(30,'(10g12.5)') (VYt(i),i=1,lspecies(Nspecies),5000) write(30,*) ' ZZZ' ! test data write(30,'(10g12.5)') (ZPt(i),i=1,lspecies(Nspecies),5000) write(30,*) ' VVV' write(30,'(10g12.5)') (VZt(i),i=1,lspecies(Nspecies),5000) XMAX = FLOAT(NGRID) + 1. XSHF = FLOAT(NGRID) SKINE =0. Wtotal =0. jstart =1 write (*,*) ' Nspecies=',Nspecies Do j =1,Nspecies vvx =0. vvy =0. vvz =0. vx2 =0. vy2 =0. vz2 =0. jend =lspecies(j) W =wspecies(j) If(jend.gt.jstart)Then write (*,*)' wrdata: jstart,jend=',jstart,jend Do i=jstart,jend vvx = vvx + VXt(i) vvy = vvy + VYt(i) vvz = vvz + VZt(i) vx2 = vx2 + VXt(i)**2 vy2 = vy2 + VYt(i)**2 vz2 = vz2 + VZt(i)**2 SKINE = SKINE + W*(VXt(i)**2 +VYt(i)**2 +VZt(i)**2) Wtotal =Wtotal +w EndDo npp =max(jend-jstart+1,1) v2 = sqrt((vx2 +vy2+vz2)/npp)/AexpV*Vscale vvx =vvx/npp/AexpV*Vscale vvy =vvy/npp/AexpV*Vscale vvz =vvz/npp/AexpV*Vscale write (*,350) vvx,vvy,vvz,v2,jend-jstart+1,j 350 format(' V(km/s)=',4g11.3, ' Npart=',i8,' Species=',i3) jstart = jend+1 Endif EndDo Ibuff =0 KROW =0 write (*,*) ' lspecies(Nspecies)=',lspecies(Nspecies) Do i=1,lspecies(Nspecies) C Periodical boundary conditions IF(XPt(i).GT.XMAX) XPt(i)=XPt(i)-XSHF IF(XPt(i).LE.1.) XPt(i)=XPt(i)+XSHF IF(YPt(i).GT.XMAX) YPt(i)=YPt(i)-XSHF IF(YPt(i).LE.1.) YPt(i)=YPt(i)+XSHF IF(ZPt(i).GT.XMAX) ZPt(i)=ZPt(i)-XSHF IF(ZPt(i).LE.1.) ZPt(i)=ZPt(i)+XSHF Ibuff = Ibuff +1 XPAR(Ibuff) = XPt(i) YPAR(Ibuff) = YPt(i) ZPAR(Ibuff) = ZPt(i) VX(Ibuff) = VXt(i) VY(Ibuff) = VYt(i) VZ(Ibuff) = VZt(i) If( XPAR(Ibuff).lt.1..or.YPAR(Ibuff).lt.1. & .or.ZPAR(Ibuff).lt.1.)write (*,*) ' Error!!!',i, & XPAR(Ibuff),YPAR(Ibuff),ZPAR(Ibuff) If( XPAR(Ibuff).ge.XMAX.or.YPAR(Ibuff).ge.XMAX & .or.ZPAR(Ibuff).ge.XMAX)write (*,*) ' Error!!!',i, & XPAR(Ibuff),YPAR(Ibuff),ZPAR(Ibuff) If( XPAR(Ibuff).ge.XMAX)Then XPAR(Ibuff) = XPAR(Ibuff) -1.e-4 write (*,500) XPAR(Ibuff) endif If( YPAR(Ibuff).ge.XMAX)Then YPAR(Ibuff) = YPAR(Ibuff) -1.e-4 write (*,500) YPAR(Ibuff) endif If( ZPAR(Ibuff).ge.XMAX)Then ZPAR(Ibuff) = ZPAR(Ibuff) -1.e-4 write (*,500) ZPAR(Ibuff) 500 format(' fixed boundary. It is now=',g14.7) endif If(Ibuff.ge.NPAGE)Then KROW = KROW +1 if(KROW/10*10.eq.KROW) & write (*,*) ' Write page=',KROW,' i=',i,Ibuff CALL WRIROW(KROW,1) Ibuff =0 EndIf EndDo If(Ibuff.ne.0)Then ! write last incomplete page KROW = KROW +1 write (*,*) ' Write page=',KROW,' i_buff=',Ibuff CALL WRIROW(KROW,1) Ibuff =0 EndIf RETURN END