!!!----------------------------------------------------------
!!!  sumCompProp
!!!----------------------------------------------------------
 !> This routine may be called to determine the distribution of bulk composition
 !> on the nodes from the bulk c of particles close to each node. Each particle
 !> is assumed to be made up of a mixture of either hazburgite and lherzolite (C<0.20)
 !> or lherzolite and basalt (0.20<C<1.0) or basalt and primordial (1.0<C<2.0.
 !> The proportions for each particle are interpolated to the nodes, weighted by particle mass.
 !> The array 'comp_prop' is filled with the composition proportions and may be
 !> output as desired. This subroutine was based off sumxfield. As in sumxfield, we need to use an
 !> extendded array (c_prop1, cprop1p) to allow for comminucation of particle
 !> information at the edge of local subdomains.
 !
 !> \param comp_prop    Array storing compositional proportion of hzb,lhz,bas,prim 

 subroutine sumCompProp(comp_prop)

   use, non_intrinsic :: meltingParam, only : initBulkComp

   real(dp), dimension(0:nt,nt+1,nd,nr+1,4),   intent(inout) :: comp_prop


   real(dp) :: attr0, attrp0
   common /ptc0/ attr0(nt,nt,nd,na,np,nr), attrp0(1,1,nd/5,na,np,nr)
   integer(i32) :: nptl, nptlp
   common /ptc2/ nptl(nt,nt,nd,nr), nptlp(1,1,nd/5,nr)
   real(dp) :: rr, dr
   common /rtrc/ rr(nr+1), dr(nr)
   real(dp) :: xr, xrp
   common /xrtn/ xr(0:nt,nt+1,nd,3,3), xrp(0:1,2,nd,3,3)
   real(dp) :: xnbr, xnbrp, qnbr, qnbrp, rdet, rdetp
   common /ngbr/ xnbr(nt,nt,6,2), xnbrp(1,1,6,2),&
        qnbr(nt,nt,6,2), qnbrp(1,1,6,2),&
        rdet(nt,nt,6),   rdetp(1,1,6)
   real(dp) :: coor(nt,nt,nd,np,2), coorp(1,1,nd/5,np,2)
   integer(i32) :: itri(nt,nt,nd,np), itrip(1,1,nd/5,np)
   real(dp) :: tin, tout, rad, rup, rdn, alph,beta,gamm, pmas, amas
   integer(i32) :: ir, id, i2,i1, ip, m1,m2, j1,j2, k1,k2, ia, ipole
   integer(i32) :: idur,idul,idlr,idll, i, isize, jd
   integer(i32) :: ii, idx, count_max, idx_max, freq_plate(0:5)

   real(dp), dimension(0:nt+1,0:nt+1,nd,5,nr+1)  :: c_prop1
   real(dp), dimension(0:1,1,nd/5,5,nr+1)       :: c_prop1p
   !real     :: c_harz=0.0, c_lhz=0.25, c_bas=1.0
   real(dp) :: c_harz=initBulkComp(1), c_lhz=initBulkComp(2)
   real(dp) :: c_bas=initBulkComp(3), c_prim=2.0
   real(dp) :: prop_harz, prop_lhz, prop_bas, prop_prim, frac

   ! Check size of flds array:
   isize = ((nt+2)**2*nd*(na-3)*(nr+1)) + ((nd/5)*(na-3)*(nr+1))

   !re-set arrays
   c_prop1  = 0D0
   c_prop1p = 0D0
   comp_prop = 0D0

 ! First sum all the active particle fields onto the nodes.
   do ir=1,nr
 
      ! Compute the barycentric coordinates of the particles.
      call ptclbary(attr0(1,1,1,1,1,ir),coor,itri,nptl(1,1,1,ir),&
           xr,xnbr,qnbr,rdet,nd,nt)
 
      if( time_terra ) call mytime(tin)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  TO DO: Add steps for is iptcl==3?!
!!!!       Add steps for nproc == 1  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         
      ! Take the weighted average / sum:
      do id=1,nd
        do i2=1,nt
          do i1=1,nt
            do ip=1,nptl(i1,i2,id,ir)
              if(attr0(i1,i2,id,4,ip,ir) > 0.) then

                m1 = itri(i1,i2,id,ip)
                m2 = mod(m1, 6) + 1
                
                rad = sqrt( attr0(i1,i2,id,1,ip,ir)**2 + &
                     attr0(i1,i2,id,2,ip,ir)**2 + &
                     attr0(i1,i2,id,3,ip,ir)**2 )
                rup = (rad - rr(ir))/dr(ir)
                rdn = 1. - rup

                ! Account for the fact that the m=6 neighbor is non-existent
                ! for the left corner node of each icosahedral diamond.

                if(nt*mod(mynum, mproc)==mt-nt .and. &
                     i1==nt .and. i2==1 .and. m1==5) m2 = 1

                j1 = i1 + ispoke(m1)
                j2 = i2 + ispoke(m1+4)
                k1 = i1 + ispoke(m2)
                k2 = i2 + ispoke(m2+4)

                alph = coor(i1,i2,id,ip,1)
                beta = coor(i1,i2,id,ip,2)
                gamm = 1. - alph - beta
                
                ! Note: pmas is the particle mass.
                pmas = attr0(i1,i2,id,4,ip,ir)
                
                
                !Splitting bulk C into 4 components (harz=0.0, lhz=0.20, bas=1.0,)
                !and so representing C as a mechanical mixture of these endmembers.
                !Store prop_lhz and prop_bas (prop_harz = 1 - prop_lhz - prop_bas)
                !We assume that all particles are composed of a mixture of two of
                !these three end members, that is particles with c_hzb<C<c_lhz are a mixture
                !of harz and lhz, and particles with c_lhz<C<c_bas are a mxture of lhz and
                !basalt.
                if(attr0(i1,i2,id,attr_index%Comp,ip,ir)<c_lhz) then
                  !Calculate fraction harz
                  frac=(c_lhz-attr0(i1,i2,id,attr_index%Comp,ip,ir))/(c_lhz-c_harz)
                  !weighted frac harz
                  prop_harz=frac*pmas
                  !weighted frac lhz (dont need harz)
                  prop_lhz=(1.-frac)*pmas

                  !Distribute harz
                  c_prop1(i1,i2,id,2,nr+2-ir) = &
                    c_prop1(i1,i2,id,2,nr+2-ir) + prop_harz*rdn*gamm
                  c_prop1(j1,j2,id,2,nr+2-ir) = &
                    c_prop1(j1,j2,id,2,nr+2-ir) + prop_harz*rdn*alph
                  c_prop1(k1,k2,id,2,nr+2-ir) = &
                    c_prop1(k1,k2,id,2,nr+2-ir) + prop_harz*rdn*beta
                  c_prop1(i1,i2,id,2,nr+1-ir) = &
                    c_prop1(i1,i2,id,2,nr+1-ir) + prop_harz*rup*gamm
                  c_prop1(j1,j2,id,2,nr+1-ir) = &
                    c_prop1(j1,j2,id,2,nr+1-ir) + prop_harz*rup*alph
                  c_prop1(k1,k2,id,2,nr+1-ir) = &
                    c_prop1(k1,k2,id,2,nr+1-ir) + prop_harz*rup*beta

                  !Distribute lhz
                  c_prop1(i1,i2,id,3,nr+2-ir) = &
                    c_prop1(i1,i2,id,3,nr+2-ir) + prop_lhz*rdn*gamm
                  c_prop1(j1,j2,id,3,nr+2-ir) = &
                    c_prop1(j1,j2,id,3,nr+2-ir) + prop_lhz*rdn*alph
                  c_prop1(k1,k2,id,3,nr+2-ir) = &
                    c_prop1(k1,k2,id,3,nr+2-ir) + prop_lhz*rdn*beta
                  c_prop1(i1,i2,id,3,nr+1-ir) = &
                    c_prop1(i1,i2,id,3,nr+1-ir) + prop_lhz*rup*gamm
                  c_prop1(j1,j2,id,3,nr+1-ir) = &
                    c_prop1(j1,j2,id,3,nr+1-ir) + prop_lhz*rup*alph
                  c_prop1(k1,k2,id,3,nr+1-ir) = &
                    c_prop1(k1,k2,id,3,nr+1-ir) + prop_lhz*rup*beta


                else if(attr0(i1,i2,id,attr_index%Comp,ip,ir)<c_bas)then

                  !Calculate fraction lhz
                  frac=(c_bas-attr0(i1,i2,id,attr_index%Comp,ip,ir))/(c_bas-c_lhz)
                  !weighted frac lhz
                  prop_lhz=frac*pmas
                  !weighted frac bas
                  prop_bas=(1.-frac)*pmas

                  !Distribute lhz
                  c_prop1(i1,i2,id,3,nr+2-ir) = &
                    c_prop1(i1,i2,id,3,nr+2-ir) + prop_lhz*rdn*gamm
                  c_prop1(j1,j2,id,3,nr+2-ir) = &
                    c_prop1(j1,j2,id,3,nr+2-ir) + prop_lhz*rdn*alph
                  c_prop1(k1,k2,id,3,nr+2-ir) = &
                    c_prop1(k1,k2,id,3,nr+2-ir) + prop_lhz*rdn*beta
                  c_prop1(i1,i2,id,3,nr+1-ir) = &
                    c_prop1(i1,i2,id,3,nr+1-ir) + prop_lhz*rup*gamm
                  c_prop1(j1,j2,id,3,nr+1-ir) = &
                    c_prop1(j1,j2,id,3,nr+1-ir) + prop_lhz*rup*alph
                  c_prop1(k1,k2,id,3,nr+1-ir) = &
                    c_prop1(k1,k2,id,3,nr+1-ir) + prop_lhz*rup*beta

                  !Distribute bas
                  c_prop1(i1,i2,id,4,nr+2-ir) = &
                    c_prop1(i1,i2,id,4,nr+2-ir) + prop_bas*rdn*gamm
                  c_prop1(j1,j2,id,4,nr+2-ir) = &
                    c_prop1(j1,j2,id,4,nr+2-ir) + prop_bas*rdn*alph
                  c_prop1(k1,k2,id,4,nr+2-ir) = &
                    c_prop1(k1,k2,id,4,nr+2-ir) + prop_bas*rdn*beta
                  c_prop1(i1,i2,id,4,nr+1-ir) = &
                    c_prop1(i1,i2,id,4,nr+1-ir) + prop_bas*rup*gamm
                  c_prop1(j1,j2,id,4,nr+1-ir) = &
                    c_prop1(j1,j2,id,4,nr+1-ir) + prop_bas*rup*alph
                  c_prop1(k1,k2,id,4,nr+1-ir) = &
                    c_prop1(k1,k2,id,4,nr+1-ir) + prop_bas*rup*beta

                else if(attr0(i1,i2,id,attr_index%Comp,ip,ir).le.c_prim)then

                  !Calculate fraction bas
                  frac=(c_prim-attr0(i1,i2,id,attr_index%Comp,ip,ir))/(c_prim-c_bas)
                  !weighted frac bas
                  prop_bas=frac*pmas
                  !weighted frac prim
                  prop_prim=(1.-frac)*pmas

                  !Distribute bas
                  c_prop1(i1,i2,id,4,nr+2-ir) = &
                    c_prop1(i1,i2,id,4,nr+2-ir) + prop_bas*rdn*gamm
                  c_prop1(j1,j2,id,4,nr+2-ir) = &
                    c_prop1(j1,j2,id,4,nr+2-ir) + prop_bas*rdn*alph
                  c_prop1(k1,k2,id,4,nr+2-ir) = &
                    c_prop1(k1,k2,id,4,nr+2-ir) + prop_bas*rdn*beta
                  c_prop1(i1,i2,id,4,nr+1-ir) = &
                    c_prop1(i1,i2,id,4,nr+1-ir) + prop_bas*rup*gamm
                  c_prop1(j1,j2,id,4,nr+1-ir) = &
                    c_prop1(j1,j2,id,4,nr+1-ir) + prop_bas*rup*alph
                  c_prop1(k1,k2,id,4,nr+1-ir) = &
                    c_prop1(k1,k2,id,4,nr+1-ir) + prop_bas*rup*beta

                  !Distribute prim
                  c_prop1(i1,i2,id,5,nr+2-ir) = &
                    c_prop1(i1,i2,id,5,nr+2-ir) + prop_prim*rdn*gamm
                  c_prop1(j1,j2,id,5,nr+2-ir) = &
                    c_prop1(j1,j2,id,5,nr+2-ir) + prop_prim*rdn*alph
                  c_prop1(k1,k2,id,5,nr+2-ir) = &
                    c_prop1(k1,k2,id,5,nr+2-ir) + prop_prim*rdn*beta
                  c_prop1(i1,i2,id,5,nr+1-ir) = &
                    c_prop1(i1,i2,id,5,nr+1-ir) + prop_prim*rup*gamm
                  c_prop1(j1,j2,id,5,nr+1-ir) = &
                    c_prop1(j1,j2,id,5,nr+1-ir) + prop_prim*rup*alph
                  c_prop1(k1,k2,id,5,nr+1-ir) = &
                    c_prop1(k1,k2,id,5,nr+1-ir) + prop_prim*rup*beta
                end if !c<lhz
              end if !mass>0.
            end do
          end do
        end do
      end do
         

      if(  time_terra ) then
         call mytime(tout)
         sec(28) = sec(28) + tout - tin
      endif
 
   end do
 
   if( time_terra ) call mytime(tin)
 
   ! Repeat for the particles at the pole.
   
   if(mod(mynum, mproc) == 0) then
     do ir=1,nr
        call ptclbary(attrp0(1,1,1,1,1,ir),coorp,itrip,nptlp(1,1,1,ir), &
              xrp,xnbrp,qnbrp,rdetp,nd/5,1)

 
        do id=1,nd/5
           do ip=1,nptlp(1,1,id,ir)
              if(attrp0(1,1,id,4,ip,ir) > 0.) then
                 
                 m1 = itrip(1,1,id,ip)
                 m2 = mod(m1, 6) + 1

                 ! Account for fact that m=6 neighbor does not exist.
                 if(m1 == 5) m2 = 1
                 
                 m1 = m1 + (id-1)*(nd-5)
                 m2 = m2 + (id-1)*(nd-5)
              
                 rad = sqrt( attrp0(1,1,id,1,ip,ir)**2 + &
                             attrp0(1,1,id,2,ip,ir)**2 + &
                             attrp0(1,1,id,3,ip,ir)**2 )
                 rup = (rad - rr(ir))/dr(ir)
                 rdn = 1. - rup

                 alph = coorp(1,1,id,ip,1)
                 beta = coorp(1,1,id,ip,2)
                 gamm = 1. - alph - beta
                 
                 pmas = attrp0(1,1,id,4,ip,ir)
                 
                 !Splitting bulk C into 3 components (harz=0.0, lhz=0.25, bas=1.0)
                 !and so representing C as a mechanical mixture of these endmembers.
                 !Store prop_lhz and prop_bas (prop_harz = 1 - prop_lhz - prop_bas)
                 !We assume that all particles are composed of a mixture of two of
                 !these three end members, that is particles with 0<c<0.25 are a mixture
                 !of harz and lhz, and particles with 0.25<c<1.0 are a mxture of lhz and
                 !basalt.
                 if(attr0(1,1,id,attr_index%Comp,ip,ir)<c_lhz) then
                   
                   !Calculate fraction harz
                   frac=(c_lhz-attr0(1,1,id,attr_index%Comp,ip,ir))/(c_lhz-c_harz)
                   !weighted frac harz
                   prop_harz=frac*pmas
                   !weighted frac lhz
                   prop_lhz=(1.-frac)*pmas
                   !Distribute harz
                   c_prop1p(1,1,id,2,nr+2-ir) = &
                     c_prop1p(1,1,id,2,nr+2-ir) + prop_harz*rdn*gamm
                   c_prop1(1,1,m1,2,nr+2-ir) = &
                     c_prop1(1,1,m1,2,nr+2-ir) + prop_harz*rdn*alph
                   c_prop1(1,1,m2,2,nr+2-ir) = &
                     c_prop1(1,1,m2,2,nr+2-ir) + prop_harz*rdn*beta
                   c_prop1p(1,1,id,2,nr+1-ir) = &
                     c_prop1p(1,1,id,2,nr+1-ir) + prop_harz*rup*gamm
                   c_prop1(1,1,m1,2,nr+1-ir) = &
                     c_prop1(1,1,m1,2,nr+1-ir) + prop_harz*rup*alph
                   c_prop1(1,1,m2,2,nr+1-ir) = &
                     c_prop1(1,1,m2,2,nr+1-ir) + prop_harz*rup*beta

                   !Distribute lhz
                   c_prop1p(1,1,id,3,nr+2-ir) = &
                     c_prop1p(1,1,id,3,nr+2-ir) + prop_lhz*rdn*gamm
                   c_prop1(1,1,m1,3,nr+2-ir) = &
                     c_prop1(1,1,m1,3,nr+2-ir) + prop_lhz*rdn*alph
                   c_prop1(1,1,m2,3,nr+2-ir) = &
                     c_prop1(1,1,m2,3,nr+2-ir) + prop_lhz*rdn*beta
                   c_prop1p(1,1,id,3,nr+1-ir) = &
                     c_prop1p(1,1,id,3,nr+1-ir) + prop_lhz*rup*gamm
                   c_prop1(1,1,m1,3,nr+1-ir) = &
                     c_prop1(1,1,m1,3,nr+1-ir) + prop_lhz*rup*alph
                   c_prop1(1,1,m2,3,nr+1-ir) = &
                     c_prop1(1,1,m2,3,nr+1-ir) + prop_lhz*rup*beta

                 else if(attr0(1,1,id,attr_index%Comp,ip,ir)<c_bas)then
                   !Calculate fraction lhz
                   frac=(c_bas-attr0(1,1,id,attr_index%Comp,ip,ir))/(c_bas-c_lhz)
                   !weighted frac lhz
                   prop_lhz=frac*pmas
                   !weighted frac bas
                   prop_bas=(1.-frac)*pmas

                   !Distribute lhz
                   c_prop1p(1,1,id,3,nr+2-ir) = &
                     c_prop1p(1,1,id,3,nr+2-ir) + prop_lhz*rdn*gamm
                   c_prop1(1,1,m1,3,nr+2-ir) = &
                     c_prop1(1,1,m1,3,nr+2-ir) + prop_lhz*rdn*alph
                   c_prop1(1,1,m2,3,nr+2-ir) = &
                     c_prop1(1,1,m2,3,nr+2-ir) + prop_lhz*rdn*beta
                   c_prop1p(1,1,id,3,nr+1-ir) = &
                     c_prop1p(1,1,id,3,nr+1-ir) + prop_lhz*rup*gamm
                   c_prop1(1,1,m1,3,nr+1-ir) = &
                     c_prop1(1,1,m1,3,nr+1-ir) + prop_lhz*rup*alph
                   c_prop1(1,1,m2,3,nr+1-ir) = &
                     c_prop1(1,1,m2,3,nr+1-ir) + prop_lhz*rup*beta

                   !Distribute bas
                   c_prop1p(1,1,id,4,nr+2-ir) = &
                     c_prop1p(1,1,id,4,nr+2-ir) + prop_bas*rdn*gamm
                   c_prop1(1,1,m1,4,nr+2-ir) = &
                     c_prop1(1,1,m1,4,nr+2-ir) + prop_bas*rdn*alph
                   c_prop1(1,1,m2,4,nr+2-ir) = &
                     c_prop1(1,1,m2,4,nr+2-ir) + prop_bas*rdn*beta
                   c_prop1p(1,1,id,4,nr+1-ir) = &
                     c_prop1p(1,1,id,4,nr+1-ir) + prop_bas*rup*gamm
                   c_prop1(1,1,m1,4,nr+1-ir) = &
                     c_prop1(1,1,m1,4,nr+1-ir) + prop_bas*rup*alph
                   c_prop1(1,1,m2,4,nr+1-ir) = &
                     c_prop1(1,1,m2,4,nr+1-ir) + prop_bas*rup*beta

                 else if(attr0(1,1,id,attr_index%Comp,ip,ir).le.c_prim)then
                   !Calculate fraction bas
                   frac=(c_prim-attr0(1,1,id,attr_index%Comp,ip,ir))/(c_prim-c_bas)
                   !weighted frac bas
                   prop_bas=frac*pmas
                   !weighted frac prim
                   prop_prim=(1.-frac)*pmas

                   !Distribute bas
                   c_prop1p(1,1,id,4,nr+2-ir) = &
                     c_prop1p(1,1,id,4,nr+2-ir) + prop_bas*rdn*gamm
                   c_prop1(1,1,m1,4,nr+2-ir) = &
                     c_prop1(1,1,m1,4,nr+2-ir) + prop_bas*rdn*alph
                   c_prop1(1,1,m2,4,nr+2-ir) = &
                     c_prop1(1,1,m2,4,nr+2-ir) + prop_bas*rdn*beta
                   c_prop1p(1,1,id,4,nr+1-ir) = &
                     c_prop1p(1,1,id,4,nr+1-ir) + prop_bas*rup*gamm
                   c_prop1(1,1,m1,4,nr+1-ir) = &
                     c_prop1(1,1,m1,4,nr+1-ir) + prop_bas*rup*alph
                   c_prop1(1,1,m2,4,nr+1-ir) = &
                     c_prop1(1,1,m2,4,nr+1-ir) + prop_bas*rup*beta

                   !Distribute prim
                   c_prop1p(1,1,id,5,nr+2-ir) = &
                     c_prop1p(1,1,id,5,nr+2-ir) + prop_prim*rdn*gamm
                   c_prop1(1,1,m1,5,nr+2-ir) = &
                     c_prop1(1,1,m1,5,nr+2-ir) + prop_prim*rdn*alph
                   c_prop1(1,1,m2,5,nr+2-ir) = &
                     c_prop1(1,1,m2,5,nr+2-ir) + prop_prim*rdn*beta
                   c_prop1p(1,1,id,5,nr+1-ir) = &
                     c_prop1p(1,1,id,5,nr+1-ir) + prop_prim*rup*gamm
                   c_prop1(1,1,m1,5,nr+1-ir) = &
                     c_prop1(1,1,m1,5,nr+1-ir) + prop_prim*rup*alph
                   c_prop1(1,1,m2,5,nr+1-ir) = &
                     c_prop1(1,1,m2,5,nr+1-ir) + prop_prim*rup*beta

                 end if
              end if ! attrp0 > 0
           end do ! ip
        end do ! id
      end do ! ir
      
   end if ! poles
 
   ! Communicate and sum data at points outside subdomain:
 
   if(nt*mod(mynum, mproc) == mt-nt) then
      do ir=1,nr+1
         do id=1,nd
            !same for proportioned composition
            do ia=2,5
              c_prop1(nt+1,1,id,ia,ir) = c_prop1(nt+1,1,id,ia,ir)  &
                   + c_prop1(nt+1,0,id,ia,ir)
              c_prop1(nt+1,0,id,ia,ir) = 0.
            end do
         end do
      end do
   end if
 
   if(mod(mynum, mproc) == 0) then
      do ir=1,nr+1
         do id=1,nd
            ipole = 1 + id/6
            !same for proportioned composition
            do ia=2,5
              c_prop1p(1,1,ipole,ia,ir) = c_prop1p(1,1,ipole,ia,ir) &
                    + c_prop1p(0,1,ipole,ia,ir)
              c_prop1p(0,1,ipole,ia,ir) = 0.
            end do
         end do
      end do
   end if

 
   ! The following is for calculation on a single processor.
   if(nproc == 1) then
     call pstop(mynum,"Cannot use comp_prop on single process")
   else
      ! The following is for calculation on multiple processors.
      call pghostadd(c_prop1,5*(nr+1),1)
   end if


   !Convert data structure for proportion of compositons
   do ir=1,nr+1
     do id=1,nd
       do i2=1,nt
         do i1=1,nt
           c_prop1(i1,i2,id,1,ir)=sum(c_prop1(i1,i2,id,2:5,ir))
           do ia=1,4
           comp_prop(i1,i2,id,ir,ia) = c_prop1(i1,i2,id,ia+1,ir) &
                   / (c_prop1(i1,i2,id,1,ir) + 1.e-10  )
           end do
         end do
       end do
     end do
   end do

   if(mod(mynum, mproc) == 0) then
     do ir=1,nr+1
       do id=1,nd/5
         c_prop1p(1,1,id,1,ir)=sum(c_prop1p(1,1,id,2:5,ir))
         do ia=1,4
           jd = id + (id-1)*(nd-5)

           comp_prop(0,1,jd,ir,ia) = c_prop1p(1,1,id,ia+1,ir) &
                  / (c_prop1p(1,1,id,1,ir) + 1.e-10 )
         end do
       end do
     end do
   end if
 

   call comm3s(comp_prop,nr,nt,4)
   
 
   if( time_terra ) then
      call mytime(tout)
      sec(28) = sec(28) + tout - tin
   endif
   

 end subroutine sumCompProp



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Below are code fragments implemented for converting from simulation !!
!! properties to seismic properties assuming a mechanical mixture.     !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


...
 else if( params%isctabs .ge. 2) then !new seismic conversion files - from Bob Myhill

      funit1 = io_get_free_unit()
      ! Open first file for Harzburgite material
      open(funit1,file='hzb_anelastic_properties.dat',status='old')
      call log_write('Loading harzburgite data file for conversion')
      do i = 1, isize
         read(funit1,*) (wrk(i,1,j), j=1,8) ! Harz stored in wrk(x,1,x)
      end do
      close(funit1)

      funit1 = io_get_free_unit()
      ! Open first file for Lherzolite material
      open(funit1,file='lhz_anelastic_properties.dat',status='old')
      call log_write('Loading lherzolite data file for conversion')
      do i = 1, isize
         read(funit1,*) (wrk(i,2,j), j=1,8) ! lherz stored in wrk(x,2,x)
      end do
      close(funit1)


      funit1 = io_get_free_unit()
      ! Open first file for basaltic material
      open(funit1,file='crust_anelastic_properties.dat',status='old')
      call log_write('Loading basalt data file for conversion')
      do i = 1, isize
         read(funit1,*) (wrk(i,3,j), j=1,8) ! Crust stored in wrk(x,3,x)
      end do
      close(funit1)

      if(params%isctabs.ge.3) then
        funit1 = io_get_free_unit()
        ! Open first file for primordial material
        open(funit1,file='ceb_anelastic_properties.dat',status='old')
        call log_write('Loading chondrite-enriched basalt (CEB) data file for conversion')
        do i = 1, isize
           read(funit1,*) (wrk(i,4,j), j=1,8) ! Prim stored in wrk(x,4,x)
        end do
        close(funit1)
      end if
...

...
else if(params%isctabs .ge. 2) then
               !Using new mineral physics tables and new method:
               !We now store the fraction of Harz/Lherz/Bas at each node,
               !so we take a mean in order to determine properties of a mechanical mixture
               
               do j = 3, 8

                  !first perform a linear interpolation
                  call int_linear(wrk(ip1t1,1,j),wrk(ip1t2,1,j),&
                       wrk(ip2t1,1,j),wrk(ip2t2,1,j),tx1,py1,vj)

                  call int_linear(wrk(ip1t1c,2,j),wrk(ip1t2c,2,j),&
                          wrk(ip2t1c,2,j),wrk(ip2t2c,2,j),tx2,py2,vj2)

                  call int_linear(wrk(ip1t1s,3,j),wrk(ip1t2s,3,j),&
                            wrk(ip2t1s,3,j),wrk(ip2t2s,3,j),tx3,py3,vj3)

                  if(params%isctabs.ge.3)then
                    call int_linear(wrk(ip1t1f,4,j),wrk(ip1t2f,4,j),&
                            wrk(ip2t1f,4,j),wrk(ip2t2f,4,j),tx4,py4,vj4)
                  end if

                  !arithmetic mean
!                  seis(ii,id,ir,j-2) = (vj*comp_prop(ii,id,ir,1)+ &
!                  vj2*comp_prop(ii,id,ir,2) + vj3*comp_prop(ii,id,ir,3))

                  !harmonic mean
                  if(params%isctabs==2)then
                    !Use only the crust table.
                    !If C>1.0 material exists then it is treated as oceanic crust (basalt)
                    seis(ii,id,ir,j-2) = sum(comp_prop(ii,id,ir,1:4)) / &
                                        ((1/vj)*comp_prop(ii,id,ir,1)+  &
                                        (1/vj2)*comp_prop(ii,id,ir,2)+  &
                                        (1/vj3)*sum(comp_prop(ii,id,ir,3:4)))
                  else if(params%isctabs.ge.3)then
                    seis(ii,id,ir,j-2) = sum(comp_prop(ii,id,ir,1:4)) / &
                                         ((1/vj)*comp_prop(ii,id,ir,1)+ &
                                        (1/vj2)*comp_prop(ii,id,ir,2) + &
                                        (1/vj3)*comp_prop(ii,id,ir,3) + &
                                        (1/vj4)*comp_prop(ii,id,ir,4))
                  end if
  
               end do
...




