diff --git a/source/merger_trees.construct.build.masses.read.F90 b/source/merger_trees.construct.build.masses.read.F90 index 9a6202c86..ec1b3f4c3 100644 --- a/source/merger_trees.construct.build.masses.read.F90 +++ b/source/merger_trees.construct.build.masses.read.F90 @@ -37,7 +37,7 @@ contains !![ - + !!] procedure :: construct => readConstruct @@ -58,13 +58,16 @@ subroutine readConstruct(self,time,mass,massMinimum,massMaximum,weight) !!{ Construct a set of merger tree masses by reading from a file. !!} - use :: Sorting , only : sort + use :: Sorting, only : sort implicit none class (mergerTreeBuildMassesRead), intent(inout) :: self double precision , intent(in ) :: time - double precision , intent( out), allocatable, dimension(:) :: mass , weight , & - & massMinimum , massMaximum - integer :: i + double precision , intent( out), allocatable, dimension(:) :: mass , weight , & + & massMinimum , massMaximum + double precision :: massPrevious + integer :: iStart , iEnd , & + & iStartPrevious, iEndPrevious, & + & i !$GLC attributes unused :: time call self%read(mass,weight) @@ -74,19 +77,34 @@ subroutine readConstruct(self,time,mass,massMinimum,massMaximum,weight) call sort(mass ) allocate(massMinimum,mold=mass) allocate(massMaximum,mold=mass) - do i=1,size(mass) - massMinimum(i)=+mass(i)/sqrt(1.0d0+self%massIntervalFractional) - massMaximum(i)=+mass(i)*sqrt(1.0d0+self%massIntervalFractional) - if (i > 1 .and. massMinimum(i) < massMaximum(i-1)) then - massMinimum(i )=sqrt( & - & +mass(i-1) & - & *mass(i ) & - & ) - massMaximum(i-1)=sqrt( & - & +mass(i-1) & - & *mass(i ) & - & ) + massPrevious =-huge(0.0d0) + iStart =-huge(0 ) + iEnd =-huge(0 ) + iStartPrevious=-huge(0 ) + iEndPrevious =-huge(0 ) + i =+ 0 + do while (i <= size(mass)) + i=i+1 + if (i == size(mass)+1 .or. mass(i) /= massPrevious) then + ! Process a block of trees of the same mass. + if (massPrevious > 0.0d0) then + massMinimum(iStart:iEnd)=+mass(iStart)/sqrt(1.0d0+self%massIntervalFractional) + massMaximum(iStart:iEnd)=+mass(iStart)*sqrt(1.0d0+self%massIntervalFractional) + if (iStart > 1 .and. massMinimum(iStart) < massMaximum(iStartPrevious)) then + massMinimum(iStart :iEnd )=+sqrt( & + & +mass (iStartPrevious) & + & *mass (iStart ) & + & ) + massMaximum(iStartPrevious:iEndPrevious)=+ massMinimum(iStart ) + end if + end if + ! Update to the start of the next block. + iStartPrevious= iStart + iEndPrevious = iEnd + iStart = i + massPrevious =mass(i ) end if + iEnd=i end do end if return diff --git a/source/merger_trees.construct.build.masses.read.XML.F90 b/source/merger_trees.construct.build.masses.read.XML.F90 index 74cbbd3eb..dd1abe443 100644 --- a/source/merger_trees.construct.build.masses.read.XML.F90 +++ b/source/merger_trees.construct.build.masses.read.XML.F90 @@ -120,8 +120,8 @@ subroutine readXMLRead(self,mass,weight) use :: File_Utilities, only : File_Name_Expand implicit none class (mergerTreeBuildMassesReadXML), intent(inout) :: self - double precision , intent( out), allocatable, dimension(:) :: mass, weight - type (node ), pointer :: doc , rootNode + double precision , intent( out), allocatable, dimension(:) :: mass , weight + type (node ), pointer :: doc , rootNode integer :: ioErr !$omp critical (FoX_DOM_Access) @@ -131,7 +131,10 @@ subroutine readXMLRead(self,mass,weight) ! Read all tree masses. call XML_Array_Read(doc,"treeRootMass",mass) ! Extract tree weights if available. - if (XML_Path_Exists(rootNode,"treeWeight")) call XML_Array_Read_Static(doc,"treeWeight",weight) + if (XML_Path_Exists(rootNode,"treeWeight")) then + call XML_Array_Read(doc,"treeWeight",weight) + if (size(weight) /= size(mass)) call Error_Report('`mass` and `weight` arrays must have the same size'//{introspection:location}) + end if ! Finished - destroy the XML document. call destroy(doc) !$omp end critical (FoX_DOM_Access) diff --git a/source/merger_trees.construct.build.masses.replicate.F90 b/source/merger_trees.construct.build.masses.replicate.F90 index b38280768..dad7dd43f 100644 --- a/source/merger_trees.construct.build.masses.replicate.F90 +++ b/source/merger_trees.construct.build.masses.replicate.F90 @@ -105,8 +105,8 @@ subroutine replicateConstruct(self,time,mass,massMinimum,massMaximum,weight) !!{ Construct a set of merger tree masses by sampling from a distribution. !!} - use :: Error , only : Error_Report - use, intrinsic :: ISO_C_Binding , only : c_size_t + use :: Error , only : Error_Report + use, intrinsic :: ISO_C_Binding, only : c_size_t implicit none class (mergerTreeBuildMassesReplicate), intent(inout) :: self double precision , intent(in ) :: time @@ -117,17 +117,32 @@ subroutine replicateConstruct(self,time,mass,massMinimum,massMaximum,weight) integer (c_size_t ) :: i call self%mergerTreeBuildMasses_%construct(time,massTmp,massMinimumTmp,massMaximumTmp,weightTmp) - if (allocated(massTmp).and.allocated(weightTmp)) then - allocate(mass (size(massTmp )*self%replicationCount)) - allocate(weight(size(weightTmp)*self%replicationCount)) + if (allocated(massTmp)) then + allocate (mass (size(massTmp )*self%replicationCount)) + if (allocated(massMinimumTmp)) & + & allocate(massMinimum(size(massMinimumTmp)*self%replicationCount)) + if (allocated(massMaximumTmp)) & + & allocate(massMaximum(size(massMaximumTmp)*self%replicationCount)) + if (allocated(weightTmp)) & + & allocate(weight (size(weightTmp )*self%replicationCount)) do i=1,size(massTmp) - mass ((i-1)*self%replicationCount+1:i*self%replicationCount)=massTmp (i) - weight((i-1)*self%replicationCount+1:i*self%replicationCount)=weightTmp(i)/dble(self%replicationCount) + mass ((i-1)*self%replicationCount+1:i*self%replicationCount)=massTmp (i) + if (allocated(massMinimumTmp)) & + & massMinimum((i-1)*self%replicationCount+1:i*self%replicationCount)=massMinimumTmp(i) + if (allocated(massMaximumTmp)) & + & massMaximum((i-1)*self%replicationCount+1:i*self%replicationCount)=massMaximumTmp(i) + if (allocated(weightTmp)) & + & weight ((i-1)*self%replicationCount+1:i*self%replicationCount)=weightTmp (i)/dble(self%replicationCount) end do - deallocate(massTmp ) - deallocate(weightTmp) + deallocate (massTmp ) + if (allocated(massMinimumTmp)) & + & deallocate(massMinimumTmp) + if (allocated(massMaximumTmp)) & + & deallocate(massMaximumTmp) + if (allocated(weightTmp )) & + & deallocate(weightTmp ) else - call Error_Report('masses and weights are required'//{introspection:location}) + call Error_Report('masses are required'//{introspection:location}) end if return end subroutine replicateConstruct