!--------------------------------------------------------- ! Module Density ! !--------------------------------------------------------- ! use Tools Contains !------------------------------------------------------------- ! ! Make density contrast FI(NGRID,NGRID,NGRID) ! for Nparticles XPAR,YPAR,ZPAR(Nparticles) in 1-Ngrid ! SUBROUTINE DENSIT !------------------------------------------------------------- real*8 :: XN,YN, & X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W integer*8 :: IN integer*4 :: iSplit = 8 Call TimingMain(3,-1) XN =Real(NGRID)+1.-1.E-7 YN =Real(NGRID) Wpar = YN**3/REAL(Nparticles,8) ! Subtract mean density !FI(:,:,:) = -1. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID FI(M1,M2,M3) = -1. END DO END DO END DO write(*,*) ' init finished: ',NGRID,seconds() !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W) & !$OMP PRIVATE (I,J,K,I1,J1,K1) DO IN=1,Nparticles ! loop over particles X=XPAR(IN) Y=YPAR(IN) Z=ZPAR(IN) I=INT(X) J=INT(Y) K=INT(Z) !if(I<1.or.J<1.or.K<1)write(*,'(a,i12,3x,3es13.5,3x,3i8)') ' Error: out off bounds: ',IN,X,Y,Z,I,J,K !if(I>NGRID.or.J>NGRID.or.K>NGRID)write(*,'(a,i12,3x,3f12.4,3x,3i8)') ' Error: out off bounds: ',IN,X,Y,Z,I,J,K D1=X-REAL(I) D2=Y-REAL(J) D3=Z-REAL(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*WPAR D2W =D2*WPAR I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 !$OMP ATOMIC FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W !$OMP ATOMIC FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W !$OMP ATOMIC FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W !$OMP ATOMIC FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W !$OMP ATOMIC FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W !$OMP ATOMIC FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W !$OMP ATOMIC FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W !$OMP ATOMIC FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W ENDDO Call TimingMain(3,1) D1 = 0. ; D2 =0. write(*,*) ' finishing density ',seconds() !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (M1,M2,M3) REDUCTION(+:D1,D2) ! DO M3=1,NGRID ! DO M2=1,NGRID ! DO M1=1,NGRID ! D1 = D1 + FI(M1,M2,M3) ! D2 = D2 + FI(M1,M2,M3)**2 ! END DO ! END DO ! END DO ! D1 = D1/(Real(NGRID))**3 ! D2 = sqrt(D2/(Real(NGRID))**3) ! write(*,*) ' Finished density: aver/rms=',D1,D2 END SUBROUTINE DENSIT !------------------------------------------------------------- ! ! Make three projections of density field ! ! SUBROUTINE ProjectDENSIT ! !------------------------------------------------------------- Real*4, Allocatable, dimension(:,:) :: DenXY,DenXZ,DenYZ Real*8 :: D1,D2,D3 character*120 :: Name Call Timing(5,-1) ALLOCATE(DenXY(NGRID,NGRID),DenXZ(NGRID,NGRID),DenYZ(NGRID,NGRID)) moment = ISTEP t0 = seconds() !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2) DO M2=1,NGRID DO M1=1,NGRID DenXY(M1,M2) = 0. DenXZ(M1,M2) = 0. DenYZ(M1,M2) = 0. END DO END DO t1 = seconds() write(*,*) ' init finished time = ',t1-t0 write(*,*) ' Ngrid = ',NGRID !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3,D1,D2,D3) DO M3=1,NGRID DO M2=1,NGRID D1 = 0. D2 = 0. D3 = 0. DO M1=1,NGRID D1 = D1 + FI(M1,M2,M3)+1. D2 = D2 + FI(M2,M1,M3)+1. D3 = D3 + FI(M2,M3,M1)+1. end DO DenYZ(M2,M3) = DenYZ(M2,M3) +D1 DenXZ(M2,M3) = DenXZ(M2,M3) +D2 DenXY(M2,M3) = DenXY(M2,M3) +D3 END DO END DO t2 = seconds() !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) DO M3=1,NGRID DO M2=1,NGRID DenYZ(M2,M3) = DenYZ(M2,M3)/NGRID DenXZ(M2,M3) = DenXZ(M2,M3)/NGRID DenXY(M2,M3) = DenXY(M2,M3)/NGRID end DO end DO write(*,*) ' finished time = ',t2-t1 write(*,'(a,12es12.4)') ' DenXY(1..2,1..2) : ',DenXY(1,1),DenXY(2,1),DenXY(1,2),DenXY(2,2) write(Name,'(2(a,i4.4),3(a,i3.3))')'ProjectionXY.',moment,'.',Nrealization, & '.dat' OPEN(19,FILE=TRIM(Name),STATUS='UNKNOWN',form='unformatted') write(Name,'(2(a,i4.4),3(a,i3.3))')'ProjectionXZ.',moment,'.',Nrealization, & '.dat' OPEN(20,FILE=TRIM(Name),STATUS='UNKNOWN',form='unformatted') write(Name,'(2(a,i4.4),3(a,i3.3))')'ProjectionYZ.',moment,'.',Nrealization, & '.dat' OPEN(21,FILE=TRIM(Name),STATUS='UNKNOWN',form='unformatted') write(*,*) AEXPN,ISTEP,NGRID write(19) AEXPN,ISTEP,NGRID write(19) DenXY write(20) AEXPN,ISTEP,NGRID write(20) DenXZ write(21) AEXPN,ISTEP,NGRID write(21) DenYZ close(19) close(20) close(21) Call Timing(5,1) end SUBROUTINE ProjectDENSIT !--------------------------------------- ! generate a vector of gaussian numbers SUBROUTINE getGauss(Gg,jp,kp,N) ! !--------------------------------------- use Random use LUXURY Integer*8 :: kp,jp Real*4 :: Gg(N) Ns =SeedsPage(jp,kp) lux = 2 Call rluxgo(lux,Ns,0,0) Do i=1,N Gg(i) = GAUSS3(gSet,iFlag) end Do end SUBROUTINE getGauss ! !-------------------------------------------------- ! : Store and retrieve seeds for parallelization ! of random number generator !-------------------------------------------------- SUBROUTINE SetRandom use Tools use LUXURY use Random write(*,*) ' Inside SetRandom' write(*,*) NROW,NGRID ALLOCATE(SeedsPage(NROW,NROW)) Ns = Nseed Do j=1,NROW Do i=1,NROW SeedsPage(i,j) = Ns dummy = RANDd(Ns) End Do End Do Nseed = Ns ! new seed !write(*,'(a,6i12)') ' Initialized random seeds: ',(SeedsPage(i,1),i=1,6) End SUBROUTINE SetRandom ! !--------------------------------------- ! DM density at particle position ! SUBROUTINE DENSPARTall ! ! !--------------------------------------- use Tools Real*8 :: ss,X,Y,Z,D1,D2,D3,T1,T2,T3 Integer*8 :: IN Integer*4,parameter :: Npd = 100 Integer*8 :: iPDF(-30:Npd) Real*4 :: dLogR= 0.05 !Integer*4, parameter :: NpeakM = 300000 Integer*4 :: Npeak Real*4 :: MassOne,MassOver character*120 :: Name Call Timing(5,-1) ! start timing write(*,'(a,i5,3x,3es12.4)') ' Inside DensPart: Ngrid =',NGRID,Box,Om !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,I,J,K,I1,J1,K1,D1,D2,D3,T1,T2,T3) Do IN =1,Nparticles X = XPAR(IN) !--- find 8 nodes Y = YPAR(IN) Z = ZPAR(IN) I = INT(X) J = INT(Y) K = INT(Z) I1 = I + 1 J1 = J + 1 K1 = K +1 If(I1>NGRID)I1=1 If(J1>NGRID)J1=1 If(K1>NGRID)K1=1 D1=X-REAL(I) D2=Y-REAL(J) D3=Z-REAL(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 dens(IN)= FI(I ,J ,K )*T3*T1*T2 + & ! density at particle position FI(I1,J ,K )*T3*D1*T2 + & FI(I ,J1,K )*T3*T1*D2 + & FI(I1,J1,K )*T3*D1*D2 + & FI(I ,J ,K1)*D3*T1*T2 + & FI(I1,J ,K1)*D3*D1*T2 + & FI(I ,J1,K1)*D3*T1*D2 + & FI(I1,J1,K1)*D3*D1*D2 + 1. end Do Call Timing(5,1) END SUBROUTINE DENSPARTall !----------------------------------------------------- ! Compute mean density and rms SUBROUTINE DENTES(DELR) !----------------------------------------------------- real*8 :: SUM1,SUM2 real*8, parameter :: dlog =0.05 Call Timing(3,-1) ! start reading time SUM1 = 0. SUM2= 0. Nn = 0 Am = 0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i) & !$OMP REDUCTION(+:SUM1,SUM2) DO K=1,NGRID DO J=1,NGRID DO I=1,NGRID SUM1 = SUM1 + FI(I,J,K) SUM2= SUM2 + (FI(I,J,K))**2 ENDDO ENDDO ENDDO Total =(REAL(NGRID))**3 DENM = SUM1/Total DELR = DSQRT(SUM2/Total-DENM**2) WRITE (*,150) DELR,DENM WRITE (17,150) DELR,DENM 150 format(20x,'Density is in units of average density', & ' in the Box',/20x, & ' RMS Delta Rho/Rho =',G11.4,/20x, & ' Mean Delta Rho/Rho =',G11.4) Call Timing(3,1) END SUBROUTINE DENTES !----------------------------------------------------- ! Compute statistics of Density SUBROUTINE DensDistr !----------------------------------------------------- real*8, parameter :: dlog =0.05 real*8 :: SUM1,SUM2 Integer*8, Allocatable,DIMENSION(:) :: nCells Real*8, Allocatable,DIMENSION(:) :: den Integer*8, Allocatable,DIMENSION(:,:) :: nCellsT Real*8, Allocatable,DIMENSION(:,:) :: denT Integer*4 :: OMP_GET_MAX_THREADS,OMP_GET_THREAD_NUM DensMax = 1.e5 !-- assumed density maximum DensMin = 0.001 !-- assumed density minimum iMin = log10(DensMin)/dlog-1 iMax = log10(DensMax)/dlog+1 iThreads = OMP_GET_MAX_THREADS() allocate(nCells(iMin:iMax),den(iMin:iMax)) allocate(nCellsT(iMin:iMax,iThreads),denT(iMin:iMax,iThreads)) nCells(:) =0 den(:) =0. nCellsT(:,:) =0 denT(:,:) =0. SUM1 = 0. SUM2= 0. Nn = 0 Am = 0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i) & !$OMP REDUCTION(+:SUM1,SUM2) DO K=1,NGRID DO J=1,NGRID DO I=1,NGRID SUM1 = SUM1 + FI(I,J,K)+1. SUM2= SUM2 + (FI(I,J,K)+1.)**2 ENDDO ENDDO ENDDO Total =(REAL(NGRID))**3 DENM = SUM1/Total DELR = DSQRT(SUM2/Total-DENM**2) WRITE (*,150) DELR,DENM,NGRID WRITE (18,150) DELR,DENM,NGRID 150 format(20x,'Density is in units of average density', & ' in the Box',/20x, & ' RMS Delta Rho/Rho =',G11.4,/20x, & ' Mean Delta Rho =',G11.4,/20x, & ' Number grid points =',i5) !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( ind,k,j,i,iOMP) Do K=1,NGRID iOMP = OMP_GET_THREAD_NUM()+1 DO J=1,NGRID DO I=1,NGRID ind = Floor(log10(FI(I,J,K)+1.)/dlog) ind = MIN(MAX(ind,iMin),iMax) nCellsT(ind,iOMP) = nCellsT(ind,iOMP) +1 denT(ind,iOMP) = denT(ind,iOMP) +FI(I,J,K)+1. END DO END DO END Do DO i=iMin,iMax ! sum results of all threads DO iP=1,iThreads nCells(i) = nCells(i) + nCellsT(i,iP) den(i) = den(i) + denT(i,iP) end DO END DO dNorm = Real(Ngrid)**3 write(18,*)' dens_left dens_right density dens*dN/d(dens)/Ncells cells' DO i=iMin,iMax ! print results d1 = 10.**(i*dlog) if(i==iMin)d1 =0. d2 = 10.**((i+1)*dlog) den(i) = den(i)/max(nCells(i),1) if(nCells(i)>0)write(18,'(3es12.4,3x,es13.5,i12)') d1,d2,den(i), & nCells(i)/(d2-d1)/dnorm*den(i),nCells(i) END DO deallocate(nCells,den) deallocate(nCellsT,denT) END SUBROUTINE DensDistr !------------------------------------------------------------- ! ! ! SUBROUTINE SelectParticles !------------------------------------------------------------- use Tools use Random real*8 :: X,Y,Z,D1,D2,D3 Integer*8 :: IN, ip fraction = BiasPars(1) threshold= BiasPars(2) Ngalaxies = Nparticles*fraction*1.05 ALLOCATE(dens(Ngalaxies)) Allocate (Xb(Ngalaxies),Yb(Ngalaxies),Zb(Ngalaxies)) Allocate (VXb(Ngalaxies),VYb(Ngalaxies),VZb(Ngalaxies)) XscaleBox = Box/NGRID Vscale = 100.*XscaleBox/AEXPN ! Scale for velocities Ng = 0 Nseed = 9845619 t0 = seconds() print '(a,i11,a,2es12.4)',' Ngalaxies= ',Ngalaxies,' Scales = ',XscaleBox,Vscale DO IN=1,Nparticles ! loop over particles X= XPAR(IN) ! grid units Y= YPAR(IN) Z= ZPAR(IN) I=INT(X) J=INT(Y) K=INT(Z) D1=X-REAL(I) !--- find nearest cell center If(D1>0.5)I = I+1 D2=Y-REAL(J) If(D2>0.5)J = J+1 D3=Z-REAL(K) If(D3>0.5)K = K+1 IF(I.GT.NGRID)I=1 IF(J.GT.NGRID)J=1 IF(K.GT.NGRID)K=1 IF(FI(I,J,K) +1.>threshold)Then dd = Randd(Nseed) IF(ddNgalaxies)Stop ' Error in number of galaxies in SelectGalaxies' Xb(Ng) = (XPAR(IN)-1.)*XscaleBox Yb(Ng) = (YPAR(IN)-1.)*XscaleBox Zb(Ng) = (ZPAR(IN)-1.)*XscaleBox VXb(Ng) = VX(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K) VYb(Ng) = VY(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K) VZb(Ng) = VZ(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K) dens(Ng) = FI(I,J,K) +1. !print '(2i11,3x,3f8.2,f8.1)',IN,Ng,X,Y,Z,FI(I,J,K) end IF end IF end DO t1 = seconds()-t0 print '(a,i9,a,f8.4)',' Number of particles selected = ',Ng,' t= ',t1 Ngalaxies = Ng end SUBROUTINE SelectParticles !------------------------------------------------------------- ! ! ! SUBROUTINE WriteDataGalaxies(iSample) !------------------------------------------------------------- use Tools character :: Fname*120, txt1*5 write(Fname,'(3(a,i4.4),a)')'Galaxies.',ISTEP,'.',iSample,'.',Nrealization,'.dat' open(50,file=TRIM(Fname),form='unformatted') !write(50,'(3a)') ' X Y Z ', & ! ' Vx Vy Vz Density' !Do ip = 1, Ngalaxies ! write(50,'(3f10.3,4f10.2)') Xgal(ip), Ygal(ip), Zgal(ip), & ! VXgal(ip), VYgal(ip), VZgal(ip), Wgal(ip) !end Do fraction = BiasPars(1) threshold= BiasPars(2) write(50) AEXPN,ISTEP,Box,Ngalaxies,NGRID,fraction,threshold write(50) (Xb(i),i=1,Ngalaxies) write(50) (Yb(i),i=1,Ngalaxies) write(50) (Zb(i),i=1,Ngalaxies) write(50) (VXb(i),i=1,Ngalaxies) write(50) (VYb(i),i=1,Ngalaxies) write(50) (VZb(i),i=1,Ngalaxies) write(50) (dens(i),i=1,Ngalaxies) close(50) write(Fname,'(3(a,i4.4),a)')'GalaxiesShort.',ISTEP,'.',iSample,'.',Nrealization,'.dat' open(50,file=TRIM(Fname)) write(50,'(a,f9.4)') 'Aexpn = ',AEXPN write(50,'(a,i4)') 'ISTEP = ',ISTEP write(50,'(a,f8.2)') 'BOX = ',BOX write(50,'(a,i10)') 'Ngalx = ',Ngalaxies write(50,'(a,i4)') 'Mesh = ',NGRID write(50,'(a,f8.4)') 'fraction = ',fraction write(50,'(a,f8.4)') 'threshold = ',threshold write(50,'(3a)') ' X Y Z ', & ' Vx Vy Vz Density' Do ip = 1, min(100,Ngalaxies) write(50,'(3f10.3,4f10.2)') Xb(ip), Yb(ip), Zb(ip), & VXb(ip), VYb(ip), VZb(ip), dens(ip) end Do close(50) DEALLOCATE(Xb,Yb,Zb,VXb,VYb,VZb,dens) end SUBROUTINE WriteDataGalaxies !----------------------------------------------------- ! Density Field SUBROUTINE DENSITrsd(iSwitch) ! ! iSwitch = 0 - real space ! = 1,2,3 - redshift space !---------------------------------------------------- use Tools use Random Real*8 :: ss,X,Y,Z,D1,D2,D3,T1,T2,T3 !Integer*8, parameter :: Nslip= 160481183, Nf=4538127 !Integer*8 :: IN, ip,jp,kp,jpp,kpp !Integer*4, save :: Nseed =198239321 !Real*4, allocatable :: Gg(:) Call Timing(2,-1) ! start reading time XN =FLOAT(NGRID)+1.-1.E-7 YN =FLOAT(NGRID) Xscale = NGRID/Box Vscale = sqrt(AEXPN/(Om+AEXPN**3*OmL))/AEXPN ! Vscale = sqrt(AEXPN/(Om + Oml/a**(3*(w0+wa))/exp(3*wa*(1-a))))/100. write(*,*) ' Inside Densit: Ngrid =',NGRID,iSwitch print *,' Box = ',Box print *,' AEXPN= ',AEXPN print *,' w0a = ',w0,wa print *,' Om = ',Om,Oml Nthreads = OMP_GET_MAX_THREADS() ! Subtract mean density !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID FI(M1,M2,M3) = -1. END DO END DO END DO W = FLOAT(NGRID)**3/float(Nparticles) ! factor = sqrt(AEXPN/(Om0+AEXPN**3*Oml0))/100. ! V(km/s)->space ! Vscale = (Box/Ngrid)100/Aexpn ! km/s PARTW = W write(*,'(a,es12.4,a,i5,a,8es12.4)') ' DENSITrsd:Particle weight = ',PARTW, & ' Switch=',iSwitch, & ' Factors= ',Xscale,Vscale xmin =1.e+5 xmax =-1.e+5 ymin =1.e+5 ymax =-1.e+5 zmin =1.e+5 zmax =-1.e+5 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,I,J,K,D1,D2,D3,T1,T2,T3,T2W,D2W,I1,J1,K1, & !$OMP vv,ip,jpp,kpp,jp,kp) Do IN=1,Nparticles If(iSwitch == 0)Then X = (XPAR(IN)-1.)*Xscale+1. Y = (YPAR(IN)-1.)*Xscale+1. Z = (ZPAR(IN)-1.)*Xscale+1. Else SELECT CASE (iSwitch) CASE (1) Z = (XPAR(IN)-1.)*Xscale+1. + Vx(IN)*Vscale Y = (YPAR(IN)-1.)*Xscale+1. X = (ZPAR(IN)-1.)*Xscale+1. CASE (2) X = (XPAR(IN)-1.)*Xscale+1. Z = (YPAR(IN)-1.)*Xscale+1. + Vy(IN)*Vscale Y = (ZPAR(IN)-1.)*Xscale+1. CASE (3) X = (XPAR(IN)-1.)*Xscale+1. Y = (YPAR(IN)-1.)*Xscale+1. Z = (ZPAR(IN)-1.)*Xscale+1. + Vz(IN)*Vscale End SELECT end If IF(X.ge.NGRID+1.)X=X-NGRID IF(Y.ge.NGRID+1.)Y=Y-NGRID IF(Z.ge.NGRID+1.)Z=Z-NGRID IF(X.lt.1.)X=X+NGRID IF(Y.lt.1.)Y=Y+NGRID IF(Z.lt.1.)Z=Z+NGRID I=INT(X) J=INT(Y) K=INT(Z) If(I.le.0)write (*,*) ' X:',X,Y,Z,' Irow=',IROW,IN If(J.le.0)write (*,*) ' Y:',X,Y,Z,' Irow=',IROW,IN If(K.le.0)write (*,*) ' Z:',X,Y,Z,' Irow=',IROW,IN If(I.gt.NGRID+1)write (*,*) ' X:',X,Y,Z,' Irow=',IROW,IN If(J.gt.NGRID+1)write (*,*) ' Y:',X,Y,Z,' Irow=',IROW,IN If(K.gt.NGRID+1)write (*,*) ' Z:',X,Y,Z,' Irow=',IROW,IN !---------------------------------------- CIC D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*W D2W =D2*W I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 !$OMP Atomic FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W !$OMP Atomic FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W !$OMP Atomic FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W !$OMP Atomic FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W !$OMP Atomic FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W !$OMP Atomic FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W !$OMP Atomic FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W !$OMP Atomic FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W !------------------------------------------ NGP ! FI(I ,J ,K ) =FI(I ,J ,K ) + W ENDDO END SUBROUTINE DENSITrsd !------------------------------------------------------------ end Module Density