C _______________________ START 3-D SIMULATIONS C C Klypin, August 1993 C INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.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/ Hd ='PMcrs' Tail='.DAT' CALL InitValues(NBYTE,SCALEL) write (*,*) ' InitValues is done. N bytes=',NBYTE NspecM = Levels ! Nspecies -1 Nmin_lev = Levels ! Nspecies write (*,*) ' Max resolution Level=',NspecM, & ' Min Level=',Nmin_lev write (16,*) ' Max resolution Level=',NspecM, & ' Min Level=',Nmin_lev If(Nspecies.ne.0.and.2**(Nspecies-1).ne.Lblock)Then write (*,*) ' Wrong number of Levels for mass refinement (', & Nspecies,') Block size =',Lblock, & ' 2**(Levels-1) must = Block size' 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') write (16,*) ' Seed=',Nseed AEXP0 = AEXPN AEXPV = AEXPN - ASTEP/2. C set mask for multiple mass resolution If(Nspecies.ne.0)Then CALL ZeroMASK(Nmin_lev) c If(Nspecies.ge.2)Then c Open(30,file='SelCoords.DAT') c CALL SetMASK(NspecM) c EndIf EndIf C get a realization of spectrum in /GRID/ NRAND =Nseed c CALL SPectr_CR(NRAND) c If(Levels.gt.1)Then CALL SPECTR(NRAND,1) ! set spectrum for higher levels CALL ZERO_SP ! remove harmonics included at lower levels c Endif C get the displacement vector by FFT IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5) Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW CALL VECTOR(IFOUR) write (*,*) ' SPECTR is done. IFOUR = ',IFOUR,' Alpha=', & (ALPHA(i),i=1,Levels) Fact = sqrt(Om0+Oml0*AEXPV**3+Ocurv*AEXPV) Do Lv =1,Levels VCONS(Lv) = -ALPHA(Lv)/(2.*PI/NGRID)*(AEXPV/AEXP0)* & SQRT(AEXPV)*Fact*2**(Lv-1) XCONS(Lv) = ALPHA(Lv)/(2.*PI/NGRID)*(AEXPN/AEXP0) & *2**(Lv-1) Enddo 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,Nspecies lspecies(i) =0 Enddo If(Nspecies.eq.0)Then ! old code with one set DO KROW = 1, NROW CALL SETXV(XCONS(1),VCONS(1),KROW) c CALL SETtest(XCONS,VCONS,KROW,NRAND) CALL WRIROW(KROW,1) ENDDO lspecies(1) =NROW**3 wspecies(1) =PARTW Wtotal = wspecies(1)*lspecies(1) EKIN = 0.5*SKINE/AEXPV**2*PARTW DP = NROW**3 SX = SX/DP SY = SY/DP SZ = SZ/DP SX2 = SQRT(SX2/DP) SY2 = SQRT(SY2/DP) SZ2 = SQRT(SZ2/DP) WRITE (*,'('' Displacement of a particle: MEAN'',3F9.4)') . SX,SY,SZ WRITE (*,'('' in cell units RMS '',3F9.4)') . SX2,SY2,SZ2 WRITE (16,'('' Displacement of a particle: RMS'',3F9.4)') . SX2,SY2,SZ2 Veloc=sqrt( SKINE/AEXPV**2/NROW**3)*SCALEL/NGRID* . 100.*hubble WRITE (*,'('' Ekin='',E12.4,'' Veloc_cold(km/s)='',F8.2, . '' Weght per cell='',f8.2,'' (must be 1)'')') . EKIN,Veloc,Wtotal/NGRID**3 WRITE (16,'('' Ekin='',E12.4,'' Veloc_cold(km/s)='',F8.2, . '' Weght per cell='',f8.2,'' (must be 1)'')') . EKIN,Veloc,Wtotal/NGRID**3 Else C find x,v for multiple mass resolution code Icurrent =0 ! current particle Iprevious =0 Jspecies = 0 Do Lv =Levels,1,-1 write (*,*) ' Go to BLOCKS. Level =',Lv CALL BLOCKS(Lv,Icurrent) If(Icurrent.gt.Iprevious)Then Jspecies =Jspecies +1 lspecies(Jspecies) =Icurrent wspecies(Jspecies) =PARTW *8/8**Lv Iprevious =Icurrent Endif Enddo Nspecies = Jspecies Vscale =ScaleL*100.*hubble/NGRID CALL WriteData(Vscale,AexpV,Wtotal) CALL Species_Print(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 Endif C CALL WRTAPE ! write header and control data C write pt.dat file: time-step for particles If(Nspecies.eq.0)Then Nparticles =lspecies(1) Else Nparticles =lspecies(Nspecies) EndIf Do ic1 =1,Nparticles XPt(ic1) =astep Enddo open ( 60 , file = 'pt.dat' , form = 'unformatted' ) write(60) (XPt(ic1),ic1=1,Nparticles) write (*,*) ' pt=', (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) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Par(6),ns,qqscaleb,QSCALE Real ns,k c IF (WK.GE.FLOAT(NSPEC)) THEN c TRUNF =0. c RETURN c ENDIF If(Par(6).ne.0.)Then ! Holtzman approx k = QSCALE*wk 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 k = QSCALE*wk 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 RETURN 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 'PMinitialCR.h' 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/ 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 =') Do i=1,Levels AMPLT(i) =INPUT(' Amplitude of density fluctuations=') Enddo 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(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 Nspecies =INPUT(' Number of mass species = ') W = (FLOAT(NGRID)/FLOAT(NROW))**3 PARTW = W If(Nspecies.eq.0)Then ! old constant mass resolution write(*,*) ' You use PM code with one particle set ' Else ! multiple mass resolution write(*,*) write(*,*) ' You use multiple mass resolution' 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 Real ns,INPUT Character Header1*79 OPEN(2,file='cdm.fit') Line1 =INPUT(' Enter Line Number in 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)) Return End C------------------------------------------------ C Clear MASK for multiple mass resolution C Initiat to level Nmin_lev SUBROUTINE ZeroMask(Nmin_lev) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' C$omp parallel do default(shared) C$omp& private(MK3,MK2,MK1) 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 NspecM C SUBROUTINE SetMask(NspecM) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' 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,ww,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 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) =NspecM ! 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 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)=NspecM ! 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 =Nspecies-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 If(NspecM.lt.Nspecies)Then ! change maximum resolution to NspecM C$omp parallel do default(shared) C$omp& private(MK3,MK2,MK1) DO MK3 = 1,Nblocks DO MK2 = 1,Nblocks DO MK1 = 1,Nblocks iM = iMask(MK1,MK2,MK3) If(iM.gt.NspecM)iMask(MK1,MK2,MK3)=NspecM ENDDO ENDDO ENDDO Endif Do i=1,Nspecies ! 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.Nspecies.and.iM.ge.1)Then nLevel(iM) = nLevel(iM) +1 Else write (*,*) 'wrong mask: ',iM, ' block=',MK1,MK2,MK3 Stop EndIf ENDDO ENDDO ENDDO Ntot =0 Do i=Nspecies,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 'PMinitialCR.h' REAL*8 SUMM, WI3,WJ3,WK3,WD,TS,TRX c Dimension SummA(100) C Set spectrum If(Lstart.gt.Levels)Then write (*,*) ' Wrong range of levels in SPECTR: start/total=', & Lstart,Levels Return Endif DO ML =Lstart,Levels c.... Exemplar c$DIR LOOP_PARALLEL c$DIR LOOP_PRIVATE(MK3,MK2,MK1) c.... Power Challenge cc$DOACROSS LOCAL (MK3,MK2,MK1) 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,ML) =0. GRY(MK1,MK2,MK3,ML) =0. GRZ(MK1,MK2,MK3,ML) =0. ENDDO ENDDO ENDDO ENDDO DO Lv =Lstart,Levels ! set spectrum for each level c DO Lv =2,2 ! set spectrum for each level 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 c write (*,*) ' K=',k,' K3=',K3,' Nspec=',NSPEC 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,Lv) = 0. GRY(1,1,1,Lv) = 0. GRZ(1,1,1,Lv) = 0. ELSE WD = (WI3 + WJ3 + WK3)*2**(2*Lv-2) WK = SQRT(WD) TS = TRUNF(WK) * GAUSS(NRAND) TRX = TS / WD GRX(iz, j, k,Lv) = i3 * TRX *isign GRY( i,jz, k,Lv) = j3 * TRX *jsign GRZ( i, j,kz,Lv) = 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=',Lv ALPHA(Lv) = AMPLT(Lv) / SQRT(SUMM) *sqrt(8.)*1.087 ENDDO RETURN END C------------------------------------------------ C Make spectrum for Constrained Realizations C C SUBROUTINE SPectr_CR(NRAND) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' REAL*8 SUMM, WI3,WJ3,WK3,WD,TS,TRX Dimension DENS_CR(NROW,NROW,NROW) C Set spectrum SUMM = 0. WD =0. open(40,file='dens_cr.dat',status='unknown') read(40,*) DENS_CR Do k=1,NROW Do j=1,NROW Do i=1,NROW 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/NROW**3)) write (*,*) ' Average density=',SUMM/NROW**3 write (*,*) ' RMS density=',sigma IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5) Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW CALL VECTOR_dens(IFOUR,dens_cr) 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,1) =0. GRY(MK1,MK2,MK3,1) =0. GRZ(MK1,MK2,MK3,1) =0. ENDDO ENDDO ENDDO 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 write (*,*) ' K=',k,' K3=',K3,' Nspec=',NSPEC 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,1) = 0. GRY(1,1,1,1) = 0. GRZ(1,1,1,1) = 0. ELSE WD = WI3 + WJ3 + WK3 WK = SQRT(WD) TS = -dens_cr(i,j,k) TRX = TS / WD GRX(iz, j, k,1) = i3 * TRX *isign GRY( i,jz, k,1) = j3 * TRX *jsign GRZ( i, j,kz,1) = 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(1) / SQRT(SUMM) *sqrt(8.) RETURN END C------------------------------------------------ C Set to zero harmonics, which C were included at previous level C of resolution C Lstart = start on level Lstart (=1 or 2) SUBROUTINE ZERO_SP C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' Logical in_x,in_y,in_z,inside If(Levels.lt.2)Then write (*,*) ' Wrong range of levels in ZERO_SP: =', & Levels Return Endif J_cos = NROW/4 +1 ! limit for cosine harmonics J_sin_1 = NROW/2+2 ! limits for sin harmonics J_sin_2 = NROW/4*3 C$omp parallel do default(shared) C$omp& private(K,J,I,in_z,in_y,in_x,inside,Lv) DO K = 1,NROW If(K.le.J_cos.or.(K.ge.J_sin_1.and.K.le.J_sin_2))Then in_z =.true. else in_z =.false. Endif DO J = 1,NROW If(J.le.J_cos.or.(J.ge.J_sin_1.and.J.le.J_sin_2))Then in_y =.true. else in_y =.false. Endif DO I = 1,NROW If(I.le.J_cos.or.(I.ge.J_sin_1.and.I.le.J_sin_2))Then in_x =.true. else in_x =.false. Endif inside = (in_x.and.in_y.and.in_z) If(inside)Then c write (*,*) I,J,K, in_x,in_y,in_z DO Lv =2,Levels GRX(I,J,K,Lv) =0. GRY(I,J,K,Lv) =0. GRZ(I,J,K,Lv) =0. ENDDO Endif ENDDO ENDDO 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(Indx,Icurrent) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' Dimension nLevel(10) Do i=1,Nspecies ! number of particles at nLevel(i)=0 ! different levels Enddo If(Indx.eq.Levels)Then Do i =1,Levels Delta(i) = Float(NGRID)/NROW *2./2**i Delta2(i) = Delta(i)/2. iNROW(i) =NROW *2**i/2 Enddo Endif QFACT = FLOAT(NGRID)/FLOAT(NROW) DO K = 1,Nblocks qq3 = Delta(1)*(K-1.) +1. DO J = 1,Nblocks qq2 = Delta(1)*(J-1.) +1. DO I = 1,Nblocks qq1 = Delta(1)*(I-1.) +1. iM =iMask(I,J,K) C---------------------------------------- c If(iM.eq.Indx)Then If(iM.eq.1)Then ! Low resolution: one particle DX = GRX(I,J,K,1) DY = GRY(I,J,K,1) DZ = GRZ(I,J,K,1) Q3 = qq3 +Delta2(1) ! scale them to box size Q2 = qq2 +Delta2(1) Q1 = qq1 +Delta2(1) Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & I,J,K, ' Level=',iM STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - XCONS(1)*DX + 1.E-3 +0.5 YPt(Icurrent) = Q2 - XCONS(1)*DY + 1.E-3 +0.5 ZPt(Icurrent) = Q3 - XCONS(1)*DZ + 1.E-3 +0.5 VXt(Icurrent) = VCONS(1)*DX VYt(Icurrent) = VCONS(1)*DY VZt(Icurrent) = VCONS(1)*DZ Else C------------------------------- High resolution: C iM = current level ( > 1) Nlev =2**(iM-1) ! number of particles in one dimenstion per Block Do kp =1,Nlev Q3 = qq3+Delta(iM)*(kp-1.) +Delta2(iM) ! scale them to box size kl = INT((Q3-1.-Delta2(iM))/Delta(iM)) +1 ! inside large box if(kl.lt.1.or.kl.gt.NROW*2**(iM-1)) & write (*,*) ' error kl=',kl,' Lev=',iM,' k=',kp,K kll = kl -INT((kl-1)/NROW)*NROW ! inside small box of GRZ if(kll.lt.1.or.kll.gt.NROW) & write (*,*) ' error kll=',kll,' Lev=',iM,' kl=',kl Do jp =1,Nlev Q2 = qq2+Delta(iM)*(jp-1.) +Delta2(iM) jl = INT((Q2-1.-Delta2(iM))/Delta(iM)) +1 ! inside large box if(jl.lt.1.or.jl.gt.NROW*2**(iM-1)) & write (*,*) ' error jl=',jl,' Lev=',iM,' j=',jp,J jll = jl -INT((jl-1)/NROW)*NROW ! inside small box of GRY if(jll.lt.1.or.jll.gt.NROW) & write (*,*) ' error jll=',jll,' Lev=',iM,' jl=',jl Do ip =1,Nlev Q1 = qq1+Delta(iM)*(ip-1.) +Delta2(iM) il = INT((Q1-1.-Delta2(iM))/Delta(iM)) +1 ! inside large box if(il.lt.1.or.il.gt.NROW*2**(iM-1)) & write (*,*) ' error il=',il,' Lev=',iM,' j=',ip,I ill = il -INT((il-1)/NROW)*NROW ! inside small box of GRX if(ill.lt.1.or.ill.gt.NROW) & write (*,*) ' error ill=',ill,' Lev=',iM,' il=',il DX = XCONS(iM)*GRX(ill,jll,kll,iM) DY = XCONS(iM)*GRY(ill,jll,kll,iM) DZ = XCONS(iM)*GRZ(ill,jll,kll,iM) DvX = VCONS(iM)*GRX(ill,jll,kll,iM) DvY = VCONS(iM)*GRY(ill,jll,kll,iM) DvZ = VCONS(iM)*GRZ(ill,jll,kll,iM) Do Lv =iM-1,1,-1 ! take contributions from lower levels ii =INT((Q1-1.)/Delta(Lv)+1000.) -999 ! inside large box jm = INT((Q2-1.)/Delta(Lv)+1000.)-999 ! inside large box km = INT((Q3-1.)/Delta(Lv)+1000.)-999 ! inside large box c ii =INT((Q1-1.-Delta2(Lv))/Delta(Lv)+1000.) -999 ! inside large box c ii1 = ii +1 c dQ1 = Q1 - (Delta(Lv)*(ii-1)+Delta2(Lv)+1.) c if(ii.lt.1)ii =ii+iNROW(Lv) c if(ii.gt.iNROW(Lv))ii =ii - iNROW(Lv) c if(ii.lt.1.or.ii.gt.iNROW(Lv)) c & write (*,*) ' error ii=',ii,' Lev=',Lv,' j=',ip,I iim = ii -INT((ii -1)/NROW)*NROW ! inside small box of GRX c iim1 = ii1-INT((ii1-1)/NROW)*NROW if(iim.lt.1.or.iim.gt.NROW) & write (*,*) ' error iim=',iim,' Lev=',Lv,' ii=',ii c jm = INT((Q2-1.-Delta2(Lv))/Delta(Lv)+1000.)-999 ! inside large box c jm1 = jm +1 c dQ2= Q2 - (Delta(Lv)*(jm-1)+Delta2(Lv)+1.) c if(jm.lt.1)jm =jm+iNROW(Lv) c if(jm.gt.iNROW(Lv))jm =jm - iNROW(Lv) c if(jm.lt.1.or.jm.gt.iNROW(Lv)) c & write (*,*) ' error jm=',jm,' Lev=',Lv,' j=',jp,J jmm = jm -INT((jm -1)/NROW)*NROW ! inside small box of GRY c jmm1 = jm1-INT((jm1-1)/NROW)*NROW if(jmm.lt.1.or.jmm.gt.NROW) & write (*,*) ' error jmm=',jmm,' Lev=',Lv,' jm=',jm c km = INT((Q3-1.-Delta2(Lv))/Delta(Lv)+1000.)-999 ! inside large box c km1 = km +1 c dQ3 = Q3 - (Delta(Lv)*(km-1)+Delta2(Lv)+1.) c if(km.lt.1)km =km+iNROW(Lv) c if(km.gt.iNROW(Lv))km =km - iNROW(Lv) c if(km.lt.1.or.km.gt.iNROW(Lv)) c & write (*,*) ' error km=',km,' Lev=',Lv,' j=',jp,J kmm = km -INT((km -1)/NROW)*NROW ! inside small box of GRY c kmm1 = km1-INT((km1-1)/NROW)*NROW if(kmm.lt.1.or.kmm.gt.NROW) & write (*,*) ' error kmm=',kmm,' Lev=',Lv,' km=',km c write (*,300) I,J,K,ip,jp,kp,Q1,Q2,Q3, c & ill,jll,kll,ii,jm,km,iim,jmm,kmm c 300 format(3x,' IJK=',3i3,' p=',3i3,' Q=',3f6.1, c & ' ill=',3i4,' i=',3i4,3i4) uX = GRX(iim ,jmm ,kmm ,Lv) uY = GRY(iim ,jmm ,kmm ,Lv) uZ = GRZ(iim ,jmm ,kmm ,Lv) c dQ1 =dQ1/Delta(Lv) ! normalize weights c dQ2 =dQ2/Delta(Lv) c dQ3 =dQ3/Delta(Lv) c w111 = (1.-dQ1) *(1.-dQ2) *(1.-dQ3) c w211 = dQ1 *(1.-dQ2) *(1.-dQ3) c w121 = (1.-dQ1) * dQ2 *(1.-dQ3) c w221 = dQ1 * dQ2 *(1.-dQ3) c w112 = (1.-dQ1) *(1.-dQ2) * dQ3 c w212 = dQ1 *(1.-dQ2) * dQ3 c w122 = (1.-dQ1) * dQ2 * dQ3 c w222 = dQ1 * dQ2 * dQ3 c Interpolate displacements c uX = GRX(iim ,jmm ,kmm ,Lv)*w111 c & +GRX(iim1,jmm ,kmm ,Lv)*w211 c & +GRX(iim ,jmm1,kmm ,Lv)*w121 c & +GRX(iim1,jmm1,kmm ,Lv)*w221 c & +GRX(iim ,jmm ,kmm1,Lv)*w112 c & +GRX(iim1,jmm ,kmm1,Lv)*w212 c & +GRX(iim ,jmm1,kmm1,Lv)*w122 c & +GRX(iim1,jmm1,kmm1,Lv)*w222 c uY = GRY(iim ,jmm ,kmm ,Lv)*w111 c & +GRY(iim1,jmm ,kmm ,Lv)*w211 c & +GRY(iim ,jmm1,kmm ,Lv)*w121 c & +GRY(iim1,jmm1,kmm ,Lv)*w221 c & +GRY(iim ,jmm ,kmm1,Lv)*w112 c & +GRY(iim1,jmm ,kmm1,Lv)*w212 c & +GRY(iim ,jmm1,kmm1,Lv)*w122 c & +GRY(iim1,jmm1,kmm1,Lv)*w222 c uZ = GRZ(iim ,jmm ,kmm ,Lv)*w111 c & +GRZ(iim1,jmm ,kmm ,Lv)*w211 c & +GRZ(iim ,jmm1,kmm ,Lv)*w121 c & +GRZ(iim1,jmm1,kmm ,Lv)*w221 c & +GRZ(iim ,jmm ,kmm1,Lv)*w112 c & +GRZ(iim1,jmm ,kmm1,Lv)*w212 c & +GRZ(iim ,jmm1,kmm1,Lv)*w122 c & +GRZ(iim1,jmm1,kmm1,Lv)*w222 DX = DX +XCONS(Lv)*uX ! add contrubutions DY = DY +XCONS(Lv)*uY DZ = DZ +XCONS(Lv)*uZ DvX = DvX +VCONS(Lv)*uX DvY = DvY +VCONS(Lv)*uY DvZ = DvZ +VCONS(Lv)*uZ Enddo ! end Level Icurrent = Icurrent +1 If(Icurrent.gt.Nmaxpart)Then write (*,*)'Too many particles:',Icurrent,' block=', & I,J,K,' Levels=',iM STOP EndIf nLevel(Indx) =nLevel(Indx) +1 XPt(Icurrent) = Q1 - DX - 1.5999E-3 +0.5 YPt(Icurrent) = Q2 - DY - 1.5999E-3 +0.5 ZPt(Icurrent) = Q3 - DZ - 1.5999E-3 +0.5 VXt(Icurrent) = DvX VYt(Icurrent) = DvY VZt(Icurrent) = DvZ EndDo EndDo EndDo Endif ! end iM =1 EndIf ! end iM = Indx EndDo EndDo EndDo write (*,*) ' Particles=',Icurrent,' on levels:', & (nLevel(i),i=1,Nspecies) 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 'PMinitialCR.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 'PMinitialCR.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,1) DY = XCONS1*GRY(I,J,K,1) DZ = XCONS1*GRZ(I,J,K,1) 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,1) VY(IJ) = VCONS1*GRY(I,J,K,1) VZ(IJ) = VCONS1*GRZ(I,J,K,1) 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 'PMinitialCR.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67),dens_cr(NROW,NROW,NROW) jb1 =3 jq =IFOUR CALL SETF67(jb1,jq,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) C x-component C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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) =dens_cr(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW dens_cr(I,J,K) = Vf(I) ENDDO ENDDO ENDDO write (*,*) ' 1st is done' C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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) =dens_cr(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW dens_cr(I,J,K) = Vf(J) ENDDO ENDDO ENDDO C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(i,k,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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) =dens_cr(I,J,K) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW dens_cr(I,J,K) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' x is done' Return END C------------------------------------------------ C FFT of the spectrum SUBROUTINE VECTOR(IFOUR) C------------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' dimension Uf(marr) , Vf(marr) dimension Qi(mf67) , jndx(mf67) jb1 =4 jq =IFOUR CALL SETF67(jb1,jq,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf) DO Lv=1,Levels C x-component C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW GRX(I,J,K,Lv) = Vf(I) ENDDO ENDDO ENDDO write (*,*) ' 1st is done Level=',Lv C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW GRX(I,J,K,Lv) = Vf(J) ENDDO ENDDO ENDDO C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(i,k,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW GRX(I,J,K,Lv) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' x is done ',Lv C Y-component C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW GRY(I,J,K,Lv) = Vf(I) ENDDO ENDDO ENDDO C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW GRY(I,J,K,Lv) = Vf(J) ENDDO ENDDO ENDDO C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(k,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW GRY(I,J,K,Lv) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' y is done' C Z-component C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO I=1,NROW GRZ(I,J,K,Lv) = Vf(I) ENDDO ENDDO ENDDO C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO J=1,NROW GRZ(I,J,K,Lv) = Vf(J) ENDDO ENDDO ENDDO C.... Exemplar C$DIR LOOP_PARALLEL C$DIR LOOP_PRIVATE(k,i,Uf,Vf,jp,jsl,ll1,k2,k4) C.... Power Challenge Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) 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,Lv) ENDDO CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf) DO K=1,NROW GRZ(I,J,K,Lv) = Vf(K)/8. ENDDO ENDDO ENDDO write (*,*) ' z is done. Level=',Lv ENDDO RETURN END C------------------------------------------ c c Print statistics of spicies Subroutine Species_Print INCLUDE 'PMparameters.h' INCLUDE 'PMinitialCR.h' write (*,300)Nspecies write (16,300)Nspecies 300 format('--- 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 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 Nm =Hd//Char(48)//Tail iun=21 write (*,*) ' File unit=',iun,' Name=',Nm OPEN(UNIT=iun,FILE=Nm,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---------------------------------- 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 'PMinitialCR.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 'PMinitialCR.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 'PMinitialCR.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 'PMinitialCR.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 'PMinitialCR.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 'PMinitialCR.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 'PMinitialCR.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 'PMinitialCR.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' COMMON / BUFF / XPt(Nmaxpart),YPt(Nmaxpart),ZPt(Nmaxpart), + VXt(Nmaxpart),VYt(Nmaxpart),VZt(Nmaxpart) COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2 C 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