../patch/hydro/merger/condinit.f90 !================================================================================== !=== Hydro ic for Toomre/exponential radial density profile galaxies. === !=== Sets up a galactic binary merger. === !=== Author : D.Chapon === !=== Date : 2010/09/01 === !================================================================================== module merger_parameters!{{{ use amr_commons !--------------------! ! Galactic merger IC ! !--------------------! ! Gas disk masses, given in GMsun in namelist, ! then converted in user unit. ! !!!!!! The galaxy #1 must be the heaviest !!!!! real(dp)::Mgas_disk1 = 1.0D2 real(dp)::Mgas_disk2 = 1.0D2 ! Set to 0.0 for isolated galaxy runs ! Galactic centers, 0-centered, given in user unit real(dp), dimension(3)::gal_center1 = 0.0D0 ! Set gal_center2 to values larger than 2*Lbox for isolated galaxy runs real(dp), dimension(3)::gal_center2 = 0.0D0 ! Galactic disks rotation axis real(dp), dimension(3)::gal_axis1= (/ 0.0D0, 0.0D0, 1.0D0 /) real(dp), dimension(3)::gal_axis2= (/ 0.0D0, 0.0D0, 1.0D0 /) ! Particle ic ascii file names for the galaxies. !~~~~~~~~~~~~~~~~~~~~~~~WARNING~~~~~~~~~~~~~~~~~~~~~~~! ! Assumed to be J=z-axis, 0-centered galaxy ic files. ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! character(len=256)::ic_part_file_gal1='ic_part1' character(len=256)::ic_part_file_gal2='ic_part2' ! Rotation curve files for the galaxies !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ WARNING ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! Those circular velocity files must contain : ! ! - Column #1 : radius (in pc) ! ! - Column #2 : circular velocity (in km/s) ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! character(len=256)::Vcirc_dat_file1='Vcirc1.dat' character(len=256)::Vcirc_dat_file2='Vcirc2.dat' ! Galactic global velocities, given in km/s in namelist, ! then converted in user unit. real(dp), dimension(3)::Vgal1 = 0.0D0 real(dp), dimension(3)::Vgal2 = 0.0D0 ! Gas disk typical/truncation radius/height, given in kpc in namelist, ! then converted in user unit. real(dp)::typ_radius1 = 3.0D0 real(dp)::typ_radius2 = 3.0D0 real(dp)::cut_radius1 = 10.0D0 real(dp)::cut_radius2 = 10.0D0 real(dp)::typ_height1 = 1.5D-1 real(dp)::typ_height2 = 1.5D-1 real(dp)::cut_height1 = 4.0D-1 real(dp)::cut_height2 = 4.0D-1 ! Gas density profile : radial =>'exponential' (default) or 'Toomre' ! vertical => 'exponential' or 'gaussian' character(len=16)::rad_profile='exponential' character(len=16)::z_profile='gaussian' ! Inter galactic gas density contrast factor real(dp)::IG_density_factor = 1.0D-5 !~~~~~~~~~~~~~~~~~~~~~~~~ NOTE ~~~~~~~~~~~~~~~~~~~~~~~~! ! For isolated galaxy runs : ! ! -------------------------- ! ! - set 'Mgas_disk2' to 0.0 ! ! - set 'gal_center2' to values larger than Lbox*2 ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! end module merger_parameters!}}} module merger_commons use merger_parameters real(dp), dimension(:,:), allocatable::Vcirc_dat1, Vcirc_dat2 end module merger_commons subroutine read_merger_params! {{{ use merger_parameters implicit none #ifndef WITHOUTMPI include 'mpif.h' #endif logical::nml_ok=.true. character(LEN=80)::infile !-------------------------------------------------- ! Local variables !-------------------------------------------------- real(dp)::norm_u logical::vcirc_file1_exists, vcirc_file2_exists logical::ic_part_file1_exists, ic_part_file2_exists real(dp)::scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2 !-------------------------------------------------- ! Namelist definitions !-------------------------------------------------- namelist/merger_params/ IG_density_factor & & ,gal_center1, gal_center2, Mgas_disk1, Mgas_disk2 & & ,typ_radius1, typ_radius2, cut_radius1, cut_radius2 & & ,typ_height1, typ_height2, cut_height1, cut_height2 & & ,rad_profile, z_profile, Vcirc_dat_file1, Vcirc_dat_file2 & & , ic_part_file_gal1, ic_part_file_gal2 & & ,gal_axis1, gal_axis2, Vgal1, Vgal2 CALL getarg(1,infile) open(1,file=infile) rewind(1) read(1,NML=merger_params,END=106) goto 107 106 write(*,*)' You need to set up namelist &MERGER_PARAMS in parameter file' call clean_stop 107 continue close(1) call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) !------------------------------------------------- ! This section deals with the galactic merger initial conditions !------------------------------------------------- ! Gaseous disks masses : GMsun => user unit Mgas_disk1 = Mgas_disk1 * 1.0D9 * 1.9891D33 / (scale_d * scale_l**3) Mgas_disk2 = Mgas_disk2 * 1.0D9 * 1.9891D33 / (scale_d * scale_l**3) ! Galaxy global velocities : km/s => user unit Vgal1 = Vgal1 * 1.0D5 / scale_v Vgal2 = Vgal2 * 1.0D5 / scale_v ! Gas disk typical radii (a) : kpc => user unit typ_radius1 = typ_radius1 * 3.085677581282D21 / scale_l typ_radius2 = typ_radius2 * 3.085677581282D21 / scale_l ! Gas disk max radii : kpc => user unit cut_radius1 = cut_radius1 * 3.085677581282D21 / scale_l cut_radius2 = cut_radius2 * 3.085677581282D21 / scale_l ! Gas disk typical thicknesses (h) : kpc => user unit typ_height1 = typ_height1 * 3.085677581282D21 / scale_l typ_height2 = typ_height2 * 3.085677581282D21 / scale_l ! Gas disk max thicknesses(zmax) : kpc => user unit cut_height1 = cut_height1 * 3.085677581282D21 / scale_l cut_height2 = cut_height2 * 3.085677581282D21 / scale_l select case (rad_profile) case ('Toomre') if(myid==1) write(*,*) 'Chosen hydro radial density profile :'Toomre'' case ('exponential') if(myid==1) write(*,*) 'Chosen hydro radial density profile :'exponential'' case default if(myid==1) write(*,*) 'Chosen hydro radial density profile :'exponential'' rad_profile='exponential' end select select case (z_profile) case ('gaussian') if(myid==1) write(*,*) 'Chosen hydro vertical density profile :'gaussian'' case ('exponential') if(myid==1) write(*,*) 'Chosen hydro vertical density profile :'exponential'' case default if(myid==1) write(*,*) 'Chosen hydro vertical density profile :'exponential'' z_profile='exponential' end select if(Mgas_disk2 .GT. Mgas_disk1)then if(myid==1)write(*,*)'Error: The galaxy #1 must be bigger than #2' nml_ok=.false. endif ! Check whether velocity radial profile files exist or not. inquire(file=trim(Vcirc_dat_file1), exist=vcirc_file1_exists) if(.NOT. vcirc_file1_exists) then if(myid==1)write(*,*)'Error: Vcirc_dat_file1 ''', trim(Vcirc_dat_file1), ''' doesn''t exist ' nml_ok=.false. end if inquire(file=trim(Vcirc_dat_file2), exist=vcirc_file2_exists) if(.NOT. vcirc_file2_exists) then if(myid==1)write(*,*)'Error: Vcirc_dat_file2 ''', trim(Vcirc_dat_file2), ''' doesn''t exist ' nml_ok=.false. end if if ((vcirc_file1_exists) .AND. (vcirc_file2_exists)) then call read_vcirc_files() end if ! Check whether ic_part files exist or not. inquire(file=trim(initfile(levelmin))//'/'//trim(ic_part_file_gal1), exist=ic_part_file1_exists) if(.NOT. ic_part_file1_exists) then if(myid==1)write(*,*)'Error: ic_part_file1 ''', trim(ic_part_file_gal1), ''' doesn''t exist in '''& & , trim(initfile(levelmin)) nml_ok=.false. end if inquire(file=trim(initfile(levelmin))//'/'//trim(ic_part_file_gal2), exist=ic_part_file2_exists) if(.NOT. ic_part_file2_exists) then if(myid==1)write(*,*)'Error: ic_part_file2 ''', trim(ic_part_file_gal2), ''' doesn''t exist in '''& & , trim(initfile(levelmin)) nml_ok=.false. end if ! galactic rotation axis norm_u = sqrt(gal_axis1(1)**2 + gal_axis1(2)**2 + gal_axis1(3)**2) if(norm_u .EQ. 0.0D0) then if(myid==1)write(*,*)'Error: Galactic axis(1) is zero ' nml_ok=.false. else if(norm_u .NE. 1.0D0) then gal_axis1 = gal_axis1 / norm_u end if end if norm_u = sqrt(gal_axis2(1)**2 + gal_axis2(2)**2 + gal_axis2(3)**2) if(norm_u .EQ. 0.0D0) then if(myid==1)write(*,*)'Error: Galactic axis(2) is zero ' nml_ok=.false. else if(norm_u .NE. 1.0D0) then gal_axis2 = gal_axis2 / norm_u end if end if if(.not. nml_ok)then if(myid==1)write(*,*)'Too many errors in the namelist' if(myid==1)write(*,*)'Aborting...' call clean_stop end if end subroutine read_merger_params ! }}} subroutine condinit(x,u,dx,nn) use amr_parameters use hydro_parameters use merger_commons implicit none integer ::nn ! Number of cells real(dp)::dx ! Cell size real(dp),dimension(1:nvector,1:nvar)::u ! Conservative variables real(dp),dimension(1:nvector,1:ndim)::x ! Cell center position. !================================================================ ! This routine generates initial conditions for RAMSES. ! Positions are in user units: ! x(i,1:3) are in [0,boxlen]**ndim. ! U is the conservative variable vector. Conventions are here: ! U(i,1): d, U(i,2:ndim+1): d.u,d.v,d.w and U(i,ndim+2): E. ! Q is the primitive variable vector. Conventions are here: ! Q(i,1): d, Q(i,2:ndim+1):u,v,w and Q(i,ndim+2): P. ! If nvar >= ndim+3, remaining variables are treated as passive ! scalars in the hydro solver. ! U(:,:) and Q(:,:) are in user units. !================================================================ integer::ivar,i, ind_gal real(dp),dimension(1:nvector,1:nvar),save::q ! Primitive variables real(dp)::v,M,rho,dzz,zint, HH, rdisk, dpdr,dmax real(dp)::r, rr, rr1, rr2, abs_z real(dp), dimension(3)::vgal, axe_rot, xx1, xx2, xx, xx_rad real(dp)::rgal, sum,sum2,dmin,dmin1,dmin2,zmin,zmax,pi,tol real(dp)::scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2 real(dp)::M_b,az,eps real(dp)::rmin,rmax,a2,aa,Vcirc, HH_max real(dp)::rho_0_1, rho_0_2, rho_0, weight, da1, Vrot logical, save:: init_nml=.false. call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Read user-defined merger parameters from the namelist if (.not. init_nml) then call read_merger_params init_nml = .true. end if if(maxval(abs(gal_center1)) .GT. (boxlen/2.0D0))then write(*,*)'Error: galactic center (1) coordinates must be in the box [', & & (-boxlen/2.0D0), ';', (boxlen/2.0D0), ']^3' call clean_stop endif if((maxval(abs(gal_center2)) .GT. (boxlen/2.0D0)) .AND. (Mgas_disk2 .NE. 0.0D0))then write(*,*)'Error: galactic center (2) coordinates must be in the box [', & & (-boxlen/2.0D0), ';', (boxlen/2.0D0), ']^3' call clean_stop endif a2=1d4 / 1.2d0 / scale_T2 ! sound speed squared aa=sqrt(a2) ! Galactic central gas densities pi=dacos(-1.0D0) select case (rad_profile) case ('exponential') rho_0_1 = 1.0D0 - exp(-cut_radius1 / typ_radius1) * (1.0D0 + cut_radius1 / typ_radius1) rho_0_2 = 1.0D0 - exp(-cut_radius2 / typ_radius2) * (1.0D0 + cut_radius2 / typ_radius2) case ('Toomre') rho_0_1 = sqrt(1.0D0 + cut_radius1**2/typ_radius1**2) - 1.0D0 rho_0_2 = sqrt(1.0D0 + cut_radius2**2/typ_radius2**2) - 1.0D0 end select select case (z_profile) case ('exponential') rho_0_1 = rho_0_1 * (1.0D0 - exp(-cut_height1 / typ_height1)) rho_0_2 = rho_0_2 * (1.0D0 - exp(-cut_height2 / typ_height2)) case ('gaussian') rho_0_1 = rho_0_1 * (dsqrt(pi/2.0D0) * erf(cut_height1 / (dsqrt(2.0D0)*typ_height1))) rho_0_2 = rho_0_2 * (dsqrt(pi/2.0D0) * erf(cut_height2 / (dsqrt(2.0D0)*typ_height2))) end select rho_0_1 = Mgas_disk1 / (4.0D0 * pi * typ_radius1**2 * typ_height1 * rho_0_1) rho_0_2 = Mgas_disk2 / (4.0D0 * pi * typ_radius2**2 * typ_height2 * rho_0_2) ! Intergalactic gas density select case (rad_profile) case ('exponential') dmin1 = rho_0_1 * exp(-cut_radius1 / typ_radius1) dmin2 = rho_0_2 * exp(-cut_radius2 / typ_radius2) case ('Toomre') dmin1 = rho_0_1 / sqrt(1.0D0 + cut_radius1**2/typ_radius1**2) dmin2 = rho_0_2 / sqrt(1.0D0 + cut_radius2**2/typ_radius2**2) end select select case (z_profile) case ('exponential') dmin1 = dmin1 * exp(-cut_height1 / typ_height1) dmin2 = dmin2 * exp(-cut_height2 / typ_height2) case ('gaussian') dmin1 = dmin1 * exp(-0.5D0 * (cut_height1 / typ_height1)**2) dmin2 = dmin2 * exp(-0.5D0 * (cut_height2 / typ_height2)**2) end select if(Mgas_disk2 .NE. 0.0D0) then dmin = min(dmin1, dmin2) else dmin = dmin1 end if dmin = IG_density_factor * dmin ! Loop over cells do i=1,nn xx1=x(i,:)-(gal_center1(:)+boxlen/2.0D0) xx2=x(i,:)-(gal_center2(:)+boxlen/2.0D0) ! Compute angular velocity ! Distance between cell and both galactic centers rr1 = norm2(xx1) rr2 = norm2(xx2) ! Projected cell position over galactic centers axis da1 = dot_product(gal_center2 - gal_center1, xx1) / norm2(gal_center2 - gal_center1)**2 if(da1 .LT. (Mgas_disk1 / (Mgas_disk1 + Mgas_disk2))) then ! cell belongs to galaxy #1 ind_gal = 1 rr = rr1 xx = xx1 axe_rot = gal_axis1 vgal = Vgal1 rgal = typ_radius1 rdisk = cut_radius1 HH = typ_height1 HH_max = cut_height1 rho_0 = rho_0_1 else ! cell belongs to galaxy #2 ind_gal = 2 rr = rr2 xx = xx2 axe_rot = gal_axis2 vgal = Vgal2 rgal = typ_radius2 rdisk = cut_radius2 HH = typ_height2 HH_max = cut_height2 rho_0 = rho_0_2 end if ! Cylindric radius : distance between the cell and the galactic rotation axis xx_rad = xx - dot_product(xx,axe_rot) * axe_rot r = norm2(xx_rad) ! vertical position absolute value abs_z = sqrt(rr**2 - r**2) if(((r-dx/2.0D0).lt.rdisk) .and. ((abs_z-dx/2.0D0) .lt. HH_max)) then ! Cell in the disk : analytical density profile + rotation velocity weight = (min(r+dx/2.0D0,rdisk)-(r-dx/2.0D0))/dx if (weight .NE. 1.0D0) then r = r + (weight-1.0D0)*dx/2.0D0 end if ! Circular velocity Vcirc= find_Vcirc(r, ind_gal) weight = weight*(min(abs_z+dx/2.0D0,HH_max)-(abs_z-dx/2.0D0))/dx ! Density select case (rad_profile) case ('exponential') q(i,1)= rho_0 * exp(-r / rgal) case ('Toomre') q(i,1)= rho_0 / sqrt(1.0D0 + r**2/rgal**2) end select select case (z_profile) case ('exponential') q(i,1)= q(i,1) * exp(-abs_z / HH) case ('gaussian') q(i,1)= q(i,1) * exp(-0.5D0 * (abs_z / HH)**2) end select q(i,1) = max(weight * q(i,1), dmin) ! P=rho*cs^2 q(i,ndim+2)=a2*q(i,1) ! V = Vrot * (u_rot^xx_rad)/r + Vx_gal ! -> Vrot = sqrt(Vcirc**2 - 3*Cs^2 + r/rho * grad(rho) * Cs) select case (rad_profile) case ('exponential') Vrot = sqrt(max(Vcirc**2 - 3.0D0*a2 - r/rgal * a2,0.0D0)) case ('Toomre') Vrot = sqrt(max(Vcirc**2 - 3.0D0*a2 - r**2/(r**2+rgal**2)*a2,0.0D0)) end select Vrot = weight * Vrot q(i,ndim-1:ndim+1) = Vrot * vect_prod(axe_rot,xx_rad)/r + vgal if(metal)q(i,6)=z_ave*0.02 ! z_ave is in solar units ! if(metal)q(i,6)=z_ave*0.02*10.**(0.5-r/cut_radius1) ! z_ave is in solar units else ! Cell out of the gaseous disk : density = peanut and velocity = V_gal q(i,1)=1d-6/scale_nH q(i,ndim+2)=1d6/scale_T2*q(i,1) ! V = Vgal q(i,ndim-1:ndim+1)= vgal if(metal)q(i,6)=0.0 endif enddo ! Convert primitive to conservative variables ! density -> density u(1:nn,1)=q(1:nn,1) ! velocity -> momentum : Omega = rho * V u(1:nn,2)=q(1:nn,1)*q(1:nn,2) #if NDIM>1 u(1:nn,3)=q(1:nn,1)*q(1:nn,3) #endif #if NDIM>2 u(1:nn,4)=q(1:nn,1)*q(1:nn,4) #endif ! kinetic energy ! Total system global velocity : 0 u(1:nn,ndim+2)=0.0D0 u(1:nn,ndim+2)=u(1:nn,ndim+2)+0.5D0*q(1:nn,1)*q(1:nn,2)**2 #if NDIM>1 u(1:nn,ndim+2)=u(1:nn,ndim+2)+0.5D0*q(1:nn,1)*q(1:nn,3)**2 #endif #if NDIM>2 u(1:nn,ndim+2)=u(1:nn,ndim+2)+0.5D0*q(1:nn,1)*q(1:nn,4)**2 #endif ! pressure -> total fluid energy ! E = Ec + P / (gamma - 1) u(1:nn,ndim+2)=u(1:nn,ndim+2)+q(1:nn,ndim+2)/(gamma-1.0d0) ! passive scalars do ivar=ndim+3,nvar u(1:nn,ivar)=q(1:nn,1)*q(1:nn,ivar) end do contains !{{{ function find_Vcirc(rayon, indice) implicit none real(dp), intent(in) :: rayon integer, intent(in) :: indice real(dp) :: find_Vcirc real(dp) :: vitesse, rayon_bin, vitesse_old, rayon_bin_old integer :: k, indmax k=2 if (indice .EQ. 1) then indmax = size(Vcirc_dat1,1) rayon_bin = Vcirc_dat1(k,1) rayon_bin_old = Vcirc_dat1(k-1,1) vitesse = Vcirc_dat1(k,2) vitesse_old = Vcirc_dat1(k-1,2) else indmax = size(Vcirc_dat2,1) rayon_bin = Vcirc_dat2(k,1) rayon_bin_old = Vcirc_dat2(k-1,1) vitesse = Vcirc_dat2(k,2) vitesse_old = Vcirc_dat2(k-1,2) end if do while (rayon .GT. rayon_bin) if(k .GE. indmax) then write(*,*) 'Hydro IC error : Radius out of rotation curve !!!' call clean_stop end if k = k + 1 if (indice .EQ. 1) then vitesse_old = vitesse vitesse = Vcirc_dat1(k,2) rayon_bin_old = rayon_bin rayon_bin = Vcirc_dat1(k,1) else vitesse_old = vitesse vitesse = Vcirc_dat2(k,2) rayon_bin_old = rayon_bin rayon_bin = Vcirc_dat2(k,1) end if end do find_Vcirc = vitesse_old + (rayon - rayon_bin_old) * (vitesse - vitesse_old) / (rayon_bin - rayon_bin_old) return end function find_Vcirc function vect_prod(a,b) implicit none real(dp), dimension(3), intent(in)::a,b real(dp), dimension(3)::vect_prod vect_prod(1) = a(2) * b(3) - a(3) * b(2) vect_prod(2) = a(3) * b(1) - a(1) * b(3) vect_prod(3) = a(1) * b(2) - a(2) * b(1) end function vect_prod function norm2(x) implicit none real(dp), dimension(3), intent(in)::x real(dp) :: norm2 norm2 = sqrt(dot_product(x,x)) end function norm2 !}}} end subroutine condinit !------------------------------------------------------------------------------------- ! Circular velocity files reading ! {{{ subroutine read_vcirc_files use merger_commons implicit none integer:: nvitesses, ierr, i real(dp)::scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2 call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Galaxy #1 nvitesses = 0 open(unit=123, file=trim(Vcirc_dat_file1), iostat=ierr) do while (ierr==0) read(123,*,iostat=ierr) if(ierr==0) then nvitesses = nvitesses + 1 ! Number of samples end if end do allocate(Vcirc_dat1(nvitesses,2)) Vcirc_dat1 = 0.0D0 rewind(123) do i=1,nvitesses read(123,*) Vcirc_dat1(i,:) end do close(123) ! Unit conversion kpc -> code unit and km/s -> code unit Vcirc_dat1(:,1) = Vcirc_dat1(:,1) * 3.085677581282D21 / scale_l Vcirc_dat1(:,2) = Vcirc_dat1(:,2) * 1.0D5 / scale_v ! Galaxy #2 nvitesses = 0 open(unit=123, file=trim(Vcirc_dat_file2), iostat=ierr) do while (ierr==0) read(123,*,iostat=ierr) if(ierr==0) then nvitesses = nvitesses + 1 ! Number of samples end if end do allocate(Vcirc_dat2(nvitesses,2)) Vcirc_dat2 = 0.0D0 rewind(123) do i=1,nvitesses read(123,*) Vcirc_dat2(i,:) end do close(123) ! Unit conversion kpc -> code unit and km/s -> code unit Vcirc_dat2(:,1) = Vcirc_dat2(:,1) * 3.085677581282D21 / scale_l Vcirc_dat2(:,2) = Vcirc_dat2(:,2) * 1.0D5 / scale_v end subroutine read_vcirc_files ! }}} !-------------------------------------------------------------------------------------- ../patch/hydro/merger/init_part.f90 subroutine init_part use amr_commons use merger_commons use pm_commons use clfind_commons #ifdef RT use rt_parameters,only: convert_birth_times #endif implicit none #ifndef WITHOUTMPI include 'mpif.h' #endif !------------------------------------------------------------ ! Allocate particle-based arrays. ! Read particles positions and velocities from grafic files !------------------------------------------------------------ real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v,scale_m integer::npart2,ndim2,ncpu2 integer::ipart,jpart,ipart_old,ilevel,idim integer::i,igrid,ncache,ngrid,iskip,isink integer::ind,ix,iy,iz,ilun,info,icpu,nx_loc integer::i1,i2,i3,i1_min,i1_max,i2_min,i2_max,i3_min,i3_max integer::buf_count,indglob,npart_new real(dp)::dx,xx1,xx2,xx3,vv1,vv2,vv3,mm1,ll1,ll2,ll3 real(dp)::scale,dx_loc,rr,rmax,dx_min integer::ncode,bit_length,temp real(kind=8)::bscale real(dp),dimension(1:twotondim,1:3)::xc integer ,dimension(1:nvector)::ind_grid,ind_cell,cc,ii integer ,dimension(1:ncpu)::npart_cpu,npart_all real(dp),allocatable,dimension(:)::xdp integer ,allocatable,dimension(:)::isp real(kind=4),allocatable,dimension(:,:)::init_plane,init_plane_x real(dp),allocatable,dimension(:,:,:)::init_array,init_array_x real(kind=8),dimension(1:nvector,1:3)::xx,vv,xs real(dp),dimension(1:nvector,1:3)::xx_dp integer,dimension(1:nvector)::ixx,iyy,izz real(qdp),dimension(1:nvector)::order real(kind=8),dimension(1:nvector)::mm real(kind=8)::dispmax=0.0,dispall real(dp),dimension(1:3)::skip_loc real(dp),dimension(1:3)::centerofmass real(dp),dimension(1:3, 1:3)::Rot_gal, Prot, Qrot, Irot real(dp),dimension(1:3)::rot_axis real(dp)::cangle, sangle integer::ibuf,tag=101,tagf=102,tagu=102 integer::countsend,countrecv #ifndef WITHOUTMPI integer,dimension(MPI_STATUS_SIZE,2*ncpu)::statuses integer,dimension(2*ncpu)::reqsend,reqrecv integer,dimension(ncpu)::sendbuf,recvbuf #endif logical::error,keep_part,eof,jumped,ic_sink=.false.,read_pos=.false.,ok character(LEN=80)::filename,filename_x character(LEN=80)::fileloc character(LEN=20)::filetype_loc character(LEN=5)::nchar if(verbose)write(*,*)'Entering init_part' if(allocated(xp))then if(verbose)write(*,*)'Initial conditions already set' return end if ! Allocate particle variables allocate(xp (npartmax,ndim)) allocate(vp (npartmax,ndim)) allocate(mp (npartmax)) allocate(nextp (npartmax)) allocate(prevp (npartmax)) allocate(levelp(npartmax)) allocate(idp (npartmax)) xp=0.0; vp=0.0; mp=0.0; levelp=0; idp=0 if(star.or.sink)then allocate(tp(npartmax)) tp=0.0 if(metal)then allocate(zp(npartmax)) zp=0.0 end if end if !-------------------- ! Read part.tmp file !-------------------- if(nrestart>0)then ilun=2*ncpu+myid+10 call title(nrestart,nchar) fileloc='output_'//TRIM(nchar)//'/part_'//TRIM(nchar)//'.out' call title(myid,nchar) fileloc=TRIM(fileloc)//TRIM(nchar) open(unit=ilun,file=fileloc,form='unformatted') rewind(ilun) read(ilun)ncpu2 read(ilun)ndim2 read(ilun)npart2 read(ilun)localseed read(ilun)nstar_tot read(ilun)mstar_tot read(ilun)mstar_lost read(ilun)nsink if(ncpu2.ne.ncpu.or.ndim2.ne.ndim.or.npart2.gt.npartmax)then write(*,*)'File part.tmp not compatible' write(*,*)'Found =',ncpu2,ndim2,npart2 write(*,*)'Expected=',ncpu,ndim,npartmax call clean_stop end if ! Read position allocate(xdp(1:npart2)) do idim=1,ndim read(ilun)xdp xp(1:npart2,idim)=xdp end do ! Read velocity do idim=1,ndim read(ilun)xdp vp(1:npart2,idim)=xdp end do ! Read mass read(ilun)xdp mp(1:npart2)=xdp deallocate(xdp) ! Read identity allocate(isp(1:npart2)) read(ilun)isp idp(1:npart2)=isp ! Read level read(ilun)isp levelp(1:npart2)=isp deallocate(isp) if(star.or.sink)then ! Read birth epoch allocate(xdp(1:npart2)) read(ilun)xdp tp(1:npart2)=xdp #ifdef RT if(convert_birth_times) then do i = 1, npart2 ! Convert birth time to proper for RT postpr. call getProperTime(tp(i),tp(i)) enddo endif #endif if(metal)then ! Read metallicity read(ilun)xdp zp(1:npart2)=xdp end if deallocate(xdp) end if close(ilun) if(debug)write(*,*)'part.tmp read for processor ',myid npart=npart2 else filetype_loc=filetype if(.not. cosmo)filetype_loc='ascii' select case (filetype_loc) case ('grafic') !---------------------------------------------------- ! Reading initial conditions GRAFIC2 multigrid arrays !---------------------------------------------------- ipart=0 ! Loop over initial condition levels do ilevel=levelmin,nlevelmax if(initfile(ilevel)==' ')cycle ! Mesh size at level ilevel in coarse cell units dx=0.5D0**ilevel ! Set position of cell centers relative to grid center do ind=1,twotondim iz=(ind-1)/4 iy=(ind-1-4*iz)/2 ix=(ind-1-2*iy-4*iz) if(ndim>0)xc(ind,1)=(dble(ix)-0.5D0)*dx if(ndim>1)xc(ind,2)=(dble(iy)-0.5D0)*dx if(ndim>2)xc(ind,3)=(dble(iz)-0.5D0)*dx end do !-------------------------------------------------------------- ! First step: compute level boundaries and particle positions !-------------------------------------------------------------- i1_min=n1(ilevel)+1; i1_max=0 i2_min=n2(ilevel)+1; i2_max=0 i3_min=n3(ilevel)+1; i3_max=0 ipart_old=ipart ! Loop over grids by vector sweeps ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do do i=1,ngrid xx1=xg(ind_grid(i),1)+xc(ind,1) xx1=(xx1*(dxini(ilevel)/dx)-xoff1(ilevel))/dxini(ilevel) xx2=xg(ind_grid(i),2)+xc(ind,2) xx2=(xx2*(dxini(ilevel)/dx)-xoff2(ilevel))/dxini(ilevel) xx3=xg(ind_grid(i),3)+xc(ind,3) xx3=(xx3*(dxini(ilevel)/dx)-xoff3(ilevel))/dxini(ilevel) i1_min=MIN(i1_min,int(xx1)+1) i1_max=MAX(i1_max,int(xx1)+1) i2_min=MIN(i2_min,int(xx2)+1) i2_max=MAX(i2_max,int(xx2)+1) i3_min=MIN(i3_min,int(xx3)+1) i3_max=MAX(i3_max,int(xx3)+1) keep_part=son(ind_cell(i))==0 if(keep_part)then ipart=ipart+1 if(ipart>npartmax)then write(*,*)'Maximum number of particles incorrect' write(*,*)'npartmax should be greater than',ipart call clean_stop endif if(ndim>0)xp(ipart,1)=xg(ind_grid(i),1)+xc(ind,1) if(ndim>1)xp(ipart,2)=xg(ind_grid(i),2)+xc(ind,2) if(ndim>2)xp(ipart,3)=xg(ind_grid(i),3)+xc(ind,3) mp(ipart)=0.5d0**(3*ilevel)*(1.0d0-omega_b/omega_m) end if end do end do ! End loop over cells end do ! End loop over grids ! Check that all grids are within initial condition region error=.false. if(active(ilevel)%ngrid>0)then if(i1_min<1.or.i1_max>n1(ilevel))error=.true. if(i2_min<1.or.i2_max>n2(ilevel))error=.true. if(i3_min<1.or.i3_max>n3(ilevel))error=.true. end if if(error) then write(*,*)'Some grid are outside initial conditions sub-volume' write(*,*)'for ilevel=',ilevel write(*,*)i1_min,i1_max write(*,*)i2_min,i2_max write(*,*)i3_min,i3_max write(*,*)n1(ilevel),n2(ilevel),n3(ilevel) call clean_stop end if if(debug)then write(*,*)myid,i1_min,i1_max,i2_min,i2_max,i3_min,i3_max endif !--------------------------------------------------------------------- ! Second step: read initial condition file and set particle velocities !--------------------------------------------------------------------- ! Allocate initial conditions array if(active(ilevel)%ngrid>0)then allocate(init_array(i1_min:i1_max,i2_min:i2_max,i3_min:i3_max)) allocate(init_array_x(i1_min:i1_max,i2_min:i2_max,i3_min:i3_max)) init_array=0d0 init_array_x=0d0 end if allocate(init_plane(1:n1(ilevel),1:n2(ilevel))) allocate(init_plane_x(1:n1(ilevel),1:n2(ilevel))) ! Loop over input variables do idim=1,ndim ! Read dark matter initial displacement field if(multiple)then call title(myid,nchar) if(idim==1)filename=TRIM(initfile(ilevel))//'/dir_velcx/ic_velcx.'//TRIM(nchar) if(idim==2)filename=TRIM(initfile(ilevel))//'/dir_velcy/ic_velcy.'//TRIM(nchar) if(idim==3)filename=TRIM(initfile(ilevel))//'/dir_velcz/ic_velcz.'//TRIM(nchar) else if(idim==1)filename=TRIM(initfile(ilevel))//'/ic_velcx' if(idim==2)filename=TRIM(initfile(ilevel))//'/ic_velcy' if(idim==3)filename=TRIM(initfile(ilevel))//'/ic_velcz' if(idim==1)filename_x=TRIM(initfile(ilevel))//'/ic_poscx' if(idim==2)filename_x=TRIM(initfile(ilevel))//'/ic_poscy' if(idim==3)filename_x=TRIM(initfile(ilevel))//'/ic_poscz' INQUIRE(file=filename_x,exist=ok) if(.not.ok)then read_pos = .false. else read_pos = .true. if(myid==1)write(*,*)'Reading file '//TRIM(filename_x) end if endif if(myid==1)write(*,*)'Reading file '//TRIM(filename) if(multiple)then ilun=myid+10 open(ilun,file=filename,form='unformatted') rewind ilun read(ilun) ! skip first line do i3=1,n3(ilevel) read(ilun)((init_plane(i1,i2),i1=1,n1(ilevel)),i2=1,n2(ilevel)) if(active(ilevel)%ngrid>0)then if(i3.ge.i3_min.and.i3.le.i3_max)then init_array(i1_min:i1_max,i2_min:i2_max,i3) = & & init_plane(i1_min:i1_max,i2_min:i2_max) end if endif end do close(ilun) else if(myid==1)then open(10,file=filename,form='unformatted') rewind 10 read(10) ! skip first line end if do i3=1,n3(ilevel) if(myid==1)then if(debug.and.mod(i3,10)==0)write(*,*)'Reading plane ',i3 read(10)((init_plane(i1,i2),i1=1,n1(ilevel)),i2=1,n2(ilevel)) else init_plane=0.0 endif buf_count=n1(ilevel)*n2(ilevel) #ifndef WITHOUTMPI call MPI_BCAST(init_plane,buf_count,MPI_REAL,0,MPI_COMM_WORLD,info) #endif if(active(ilevel)%ngrid>0)then if(i3.ge.i3_min.and.i3.le.i3_max)then init_array(i1_min:i1_max,i2_min:i2_max,i3) = & & init_plane(i1_min:i1_max,i2_min:i2_max) end if endif end do if(myid==1)close(10) if(read_pos) then if(myid==1)then open(10,file=filename_x,form='unformatted') rewind 10 read(10) ! skip first line end if do i3=1,n3(ilevel) if(myid==1)then if(debug.and.mod(i3,10)==0)write(*,*)'Reading plane ',i3 read(10)((init_plane_x(i1,i2),i1=1,n1(ilevel)),i2=1,n2(ilevel)) else init_plane_x=0.0 endif buf_count=n1(ilevel)*n2(ilevel) #ifndef WITHOUTMPI call MPI_BCAST(init_plane_x,buf_count,MPI_REAL,0,MPI_COMM_WORLD,info) #endif if(active(ilevel)%ngrid>0)then if(i3.ge.i3_min.and.i3.le.i3_max)then init_array_x(i1_min:i1_max,i2_min:i2_max,i3) = & & init_plane_x(i1_min:i1_max,i2_min:i2_max) end if endif end do if(myid==1)close(10) end if endif if(active(ilevel)%ngrid>0)then ! Rescale initial displacement field to code units init_array=dfact(ilevel)*dx/dxini(ilevel)*init_array/vfact(ilevel) if(read_pos)then init_array_x = init_array_x/boxlen_ini endif ! Loop over grids by vector sweeps ipart=ipart_old ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do do i=1,ngrid xx1=xg(ind_grid(i),1)+xc(ind,1) xx1=(xx1*(dxini(ilevel)/dx)-xoff1(ilevel))/dxini(ilevel) xx2=xg(ind_grid(i),2)+xc(ind,2) xx2=(xx2*(dxini(ilevel)/dx)-xoff2(ilevel))/dxini(ilevel) xx3=xg(ind_grid(i),3)+xc(ind,3) xx3=(xx3*(dxini(ilevel)/dx)-xoff3(ilevel))/dxini(ilevel) i1=int(xx1)+1 i1=int(xx1)+1 i2=int(xx2)+1 i2=int(xx2)+1 i3=int(xx3)+1 i3=int(xx3)+1 keep_part=son(ind_cell(i))==0 if(keep_part)then ipart=ipart+1 vp(ipart,idim)=init_array(i1,i2,i3) if(.not. read_pos)then dispmax=max(dispmax,abs(init_array(i1,i2,i3)/dx)) else xp(ipart,idim)=xg(ind_grid(i),idim)+xc(ind,idim)+init_array_x(i1,i2,i3) dispmax=max(dispmax,abs(init_array_x(i1,i2,i3)/dx)) endif end if end do end do ! End loop over cells end do ! End loop over grids endif end do ! End loop over input variables ! Deallocate initial conditions array if(active(ilevel)%ngrid>0)then deallocate(init_array,init_array_x) end if deallocate(init_plane,init_plane_x) if(debug)write(*,*)'npart=',ipart,'/',npartmax,' for PE=',myid end do ! End loop over levels ! Initial particle number npart=ipart ! Move particle according to Zeldovich approximation if(.not. read_pos)then xp(1:npart,1:ndim)=xp(1:npart,1:ndim)+vp(1:npart,1:ndim) endif ! Scale displacement to velocity vp(1:npart,1:ndim)=vfact(1)*vp(1:npart,1:ndim) ! Periodic box do ipart=1,npart #if NDIM>0 if(xp(ipart,1)< 0.0d0 )xp(ipart,1)=xp(ipart,1)+dble(nx) if(xp(ipart,1)>=dble(nx))xp(ipart,1)=xp(ipart,1)-dble(nx) #endif #if NDIM>1 if(xp(ipart,2)< 0.0d0 )xp(ipart,2)=xp(ipart,2)+dble(ny) if(xp(ipart,2)>=dble(ny))xp(ipart,2)=xp(ipart,2)-dble(ny) #endif #if NDIM>2 if(xp(ipart,3)< 0.0d0 )xp(ipart,3)=xp(ipart,3)+dble(nz) if(xp(ipart,3)>=dble(nz))xp(ipart,3)=xp(ipart,3)-dble(nz) #endif end do #ifndef WITHOUTMPI ! Compute particle Hilbert ordering sendbuf=0 do ipart=1,npart xx(1,1:3)=xp(ipart,1:3) xx_dp(1,1:3)=xx(1,1:3) call cmp_cpumap(xx_dp,cc,1) if(cc(1).ne.myid)sendbuf(cc(1))=sendbuf(cc(1))+1 end do ! Allocate communication buffer in emission do icpu=1,ncpu ncache=sendbuf(icpu) if(ncache>0)then allocate(emission(icpu,1)%up(1:ncache,1:twondim+1)) end if end do ! Fill communicators jpart=0 sendbuf=0 do ipart=1,npart xx(1,1:3)=xp(ipart,1:3) xx_dp(1,1:3)=xx(1,1:3) call cmp_cpumap(xx_dp,cc,1) if(cc(1).ne.myid)then icpu=cc(1) sendbuf(icpu)=sendbuf(icpu)+1 ibuf=sendbuf(icpu) emission(icpu,1)%up(ibuf,1)=xp(ipart,1) emission(icpu,1)%up(ibuf,2)=xp(ipart,2) emission(icpu,1)%up(ibuf,3)=xp(ipart,3) emission(icpu,1)%up(ibuf,4)=vp(ipart,1) emission(icpu,1)%up(ibuf,5)=vp(ipart,2) emission(icpu,1)%up(ibuf,6)=vp(ipart,3) emission(icpu,1)%up(ibuf,7)=mp(ipart) else jpart=jpart+1 xp(jpart,1:3)=xp(ipart,1:3) vp(jpart,1:3)=vp(ipart,1:3) mp(jpart) =mp(ipart) endif end do ! Communicate virtual particle number to parent cpu call MPI_ALLTOALL(sendbuf,1,MPI_INTEGER,recvbuf,1,MPI_INTEGER,MPI_COMM_WORLD,info) ! Compute total number of newly created particles npart_new=0 do icpu=1,ncpu npart_new=npart_new+recvbuf(icpu) end do if(jpart+npart_new.gt.npartmax)then write(*,*)'No more free memory for particles' write(*,*)'Increase npartmax' write(*,*)myid write(*,*)jpart,npart_new write(*,*)bound_key call MPI_ABORT(MPI_COMM_WORLD,1,info) end if ! Allocate communication buffer in reception do icpu=1,ncpu ncache=recvbuf(icpu) if(ncache>0)then allocate(reception(icpu,1)%up(1:ncache,1:twondim+1)) end if end do ! Receive particles countrecv=0 do icpu=1,ncpu ncache=recvbuf(icpu) if(ncache>0)then buf_count=ncache*(twondim+1) countrecv=countrecv+1 call MPI_IRECV(reception(icpu,1)%up,buf_count, & & MPI_DOUBLE_PRECISION,icpu-1,& & tagu,MPI_COMM_WORLD,reqrecv(countrecv),info) end if end do ! Send particles countsend=0 do icpu=1,ncpu ncache=sendbuf(icpu) if(ncache>0)then buf_count=ncache*(twondim+1) countsend=countsend+1 call MPI_ISEND(emission(icpu,1)%up,buf_count, & & MPI_DOUBLE_PRECISION,icpu-1,& & tagu,MPI_COMM_WORLD,reqsend(countsend),info) end if end do ! Wait for full completion of receives call MPI_WAITALL(countrecv,reqrecv,statuses,info) ! Wait for full completion of sends call MPI_WAITALL(countsend,reqsend,statuses,info) ! Create new particles do icpu=1,ncpu do ibuf=1,recvbuf(icpu) jpart=jpart+1 xp(jpart,1)=reception(icpu,1)%up(ibuf,1) xp(jpart,2)=reception(icpu,1)%up(ibuf,2) xp(jpart,3)=reception(icpu,1)%up(ibuf,3) vp(jpart,1)=reception(icpu,1)%up(ibuf,4) vp(jpart,2)=reception(icpu,1)%up(ibuf,5) vp(jpart,3)=reception(icpu,1)%up(ibuf,6) mp(jpart) =reception(icpu,1)%up(ibuf,7) end do end do ! Erase old particles do ipart=jpart+1,npart xp(ipart,1)=0d0 xp(ipart,2)=0d0 xp(ipart,3)=0d0 vp(ipart,1)=0d0 vp(ipart,2)=0d0 vp(ipart,3)=0d0 mp(ipart) =0d0 end do npart=jpart ! Deallocate communicators do icpu=1,ncpu if(sendbuf(icpu)>0)deallocate(emission(icpu,1)%up) if(recvbuf(icpu)>0)deallocate(reception(icpu,1)%up) end do write(*,*)'npart=',ipart,'/',npartmax,' for PE=',myid #endif ! Compute particle initial level do ipart=1,npart levelp(ipart)=levelmin end do ! Compute particle initial age and metallicity if(star.or.sink)then do ipart=1,npart tp(ipart)=0d0 if(metal)then zp(ipart)=0d0 end if end do end if ! Compute particle initial identity npart_cpu=0; npart_all=0 npart_cpu(myid)=npart #ifndef WITHOUTMPI call MPI_ALLREDUCE(npart_cpu,npart_all,ncpu,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,info) npart_cpu(1)=npart_all(1) #endif do icpu=2,ncpu npart_cpu(icpu)=npart_cpu(icpu-1)+npart_all(icpu) end do if(myid==1)then do ipart=1,npart idp(ipart)=ipart end do else do ipart=1,npart idp(ipart)=npart_cpu(myid-1)+ipart end do end if case ('ascii') ! Conversion factor from user units to cgs units call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) scale_m=scale_d*scale_l**3 ! Local particle count ipart=0 if(TRIM(initfile(levelmin)).NE.' ')then if(myid==1) write(*,*) 'Init particles : rotated/shifted/kicked galaxies' !################### Galaxy #1 particles ############################ filename=TRIM(initfile(levelmin))//'/'//trim(ic_part_file_gal1) if(myid==1)then open(10,file=filename,form='formatted') indglob=0 end if eof=.false. ! Identity matrix Irot = 0.0D0 Irot(1,1)=1.0D0 Irot(2,2)=1.0D0 Irot(3,3)=1.0D0 if ((gal_axis1(1)==0.0D0) .AND. (gal_axis1(2)==0.0D0) .AND. (gal_axis1(3)==1.0D0)) then Rot_gal = Irot else if ((gal_axis1(1)==0.0D0) .AND. (gal_axis1(2)==0.0D0) .AND. (gal_axis1(3)==-1.0D0)) then rot_axis = (/ 1.0D0, 0.0D0, 0.0D0 /) ! x-axis cangle = -1.0D0 sangle = 0.0D0 else ! rotation axis is vect(z-axis, gal_axis1) rot_axis(1)=-gal_axis1(2) rot_axis(2)=gal_axis1(1) rot_axis(3)=0.0D0 ! sinus of rotation angle is the norm of vect(z-axis, gal_axis1) sangle = dsqrt(gal_axis1(1)**2+gal_axis1(2)**2) ! Normalisation of rotation axis rot_axis = rot_axis / sangle ! Cosine of the rotation angle cangle = gal_axis1(3) end if ! theta = datan2(sangle, cangle) ! P = rot_axis * transpose(rot_axis) => Projection along the ! rotation axis Prot(:,1) = rot_axis * rot_axis(1) Prot(:,2) = rot_axis * rot_axis(2) Prot(:,3) = rot_axis * rot_axis(3) ! Q => antisymetrical representation of the rotation axis Qrot = 0.0D0 Qrot(1,2) = -rot_axis(3) Qrot(2,1) = rot_axis(3) Qrot(1,3) = rot_axis(2) Qrot(3,1) = -rot_axis(2) Qrot(2,3) = -rot_axis(1) Qrot(3,2) = rot_axis(1) Rot_gal = Prot + (Irot-Prot)*cangle + Qrot * sangle end if ! File units are supposed to be in kpc, km/s and 10^9 Msol do while (.not.eof) xx=0.0 if(myid==1)then jpart=0 do i=1,nvector read(10,*,end=100)xx1,xx2,xx3,vv1,vv2,vv3,mm1 xx1=xx1*3.085677581282D21/scale_l xx2=xx2*3.085677581282D21/scale_l xx3=xx3*3.085677581282D21/scale_l vv1=vv1*1D5/scale_v vv2=vv2*1D5/scale_v vv3=vv3*1D5/scale_v mm1=mm1*1D9*1.9891D33/scale_m jpart=jpart+1 indglob=indglob+1 xx(i,:) = matmul(Rot_gal, (/ xx1, xx2, xx3 /)) vv(i,:) = matmul(Rot_gal, (/ vv1, vv2, vv3 /)) + Vgal1 mm(i ) = mm1 ii(i ) = indglob end do 100 continue if(jpart Projection along the ! rotation axis Prot(:,1) = rot_axis * rot_axis(1) Prot(:,2) = rot_axis * rot_axis(2) Prot(:,3) = rot_axis * rot_axis(3) ! Q => antisymetrical representation of the rotation axis Qrot = 0.0D0 Qrot(1,2) = -rot_axis(3) Qrot(2,1) = rot_axis(3) Qrot(1,3) = rot_axis(2) Qrot(3,1) = -rot_axis(2) Qrot(2,3) = -rot_axis(1) Qrot(3,2) = rot_axis(1) Rot_gal = Prot + (Irot-Prot)*cangle + Qrot * sangle end if do while (.not.eof) xx=0.0 if(myid==1)then jpart=0 do i=1,nvector read(10,*,end=103)xx1,xx2,xx3,vv1,vv2,vv3,mm1 xx1=xx1*3.085677581282D21/scale_l xx2=xx2*3.085677581282D21/scale_l xx3=xx3*3.085677581282D21/scale_l vv1=vv1*1D5/scale_v vv2=vv2*1D5/scale_v vv3=vv3*1D5/scale_v mm1=mm1*1D9*1.9891D33/scale_m jpart=jpart+1 indglob=indglob+1 xx(i,:) = matmul(Rot_gal, (/ xx1, xx2, xx3 /)) vv(i,:) = matmul(Rot_gal, (/ vv1, vv2, vv3 /)) + Vgal2 mm(i ) = mm1 ii(i ) = indglob end do 103 continue if(jpart