text:
!=================================================================================
!
!  bfragments.f90
!
!  routines to break up Lanczos vectors into "fragments"
!
!============================================================
!
!  The next module is used for
!  breaking up the Lanczos vectors
!
!=========================================================
!
!  FRAGMENTS
!
!   contiguous (sub) sectors of Lanczos vectors grouped together
!
! started 8/2011 by CWJ
!

module fragments
   use precisions
   use operation_stats
   implicit none

   logical :: fragment_vectors   ! flag to break up vectors into fragments
   integer :: nfragments         ! number of fragments
   logical :: useNewReorthog1 = .true.  ! Use new reorthog even for single fragment
   logical :: useNewReorthog  ! set to (nfragments > 1) .or. useNewReorthog1

   type frag
        integer :: nsectors  ! # of subsectors
        integer(8) :: ssectorstart,ssectorend  ! start, end of sectors; these are defined by "proton" sectors
		  integer(8) :: csectorstart,csectorend  ! start, end of conjugate sectors, defined by "neutron" sectors, must check contiguity
        integer(8):: basisstart,basisend  ! where the basis starts and stops for this sector
        integer(8):: localdim      ! local dimension = basisend-basisstart+1
        real(8)    :: totops    ! tot# of operations to/from this fragment

   end type frag

   type (frag), allocatable :: fragmentlist(:)

!   integer, allocatable :: sectormap2frag(:)  ! maps a sector to a fragment  SCHEDULED FOR OBSOLESCENCE

!--------- THIS ALLOWS US TO STORE JUST FRAGMENTS OF A VECTOR --

   integer (kind=basis_prec), allocatable :: basestart(:), basestop(:)

   logical :: test_fragments = .false.  ! added in 7.4.6, used to test new fragment scheme on a single processor

!------------- JUMP STORAGE ACROSS FRAGMENTS (f2f matvec blocks)-------------- added in 7.6.2 -----------------

   type f2fjump
	   integer(8) :: nX1bjump(2),nXX2bjump(2),nXXX3bjump(2)
	   real(8)    :: totjumpstorage
	   integer(4) :: minprocs     ! minimal assigned MPI procs/nodes
	   integer(4) :: nprocs       ! # of actual assigned procs in the end
   end type f2fjump

   type (f2fjump), allocatable,target :: f2fjumpstat(:,:)

!------- INFORMATION ON WHAT FRAGMENTS ON WHAT MPI PROCESSES
!        previously in module tribution, moved in 7.6.7

	type whatsonanode
	      integer :: ifragment,ffragment   ! which fragments are here
	      logical :: ifirst ! first node in group servicing ifragment

!	      integer :: sfragment  ! which fragment is stored here; later might want to store more than one
!	      real :: frac_ops   ! what fraction of possible ops are stored here
	end type whatsonanode

	type (whatsonanode), allocatable :: nodal(:)

contains
!======================================================================
!
!  SUBROUTINES CONTAINED IN THIS FILE
!    fragmenter:  master routine to break up Lanczos vectors into fragments
!        which calls:
!    vectorstorageceiling
!    defaultfragments:  routine to set fragment information when only 1 fragment
!    count_create_fragments_blocks: updated fragmentation routine which breaks fragments along haiku blocks
!    reconstruct_sector: when sectors are subdivided, reconstruct proton sectors
!    fragmentvectorsize:  sets base start/stop for a given fragment
!       also in this file (called by distro_opbundles_over_fragments )
!    divide_procs_over_fragments: routine to divide up available nodes over fragments
!
!=======================================================================
!
!  Notes 7.4.6 (March 2015) regarding "new" fragmentation based upon proton haiku blocks
!  To test/debug fragmentation, run on a single processor;
!  in order to do this, we must modify routines which divide up work over MPI processes
!
!  key routines include:
!          setnodaltribution  (in bparallel_lib1.f90)  assigns initial, final fragments to MPI procs/compute nodes
!          setup_localvectors (in blanczoslib1.f90)  allocates memory for (fragmented) lanczos vectors
!  If one simulates a fragmented case on a single processor, the above have to be modified
!

!======================================================================
!  subroutine set_nfragments
! Sets nfragments and side effects of assignment
subroutine set_nfragments(nf)
 !  use fragments
   use nodeinfo
   implicit none
   integer :: nf

   nfragments = nf
   useNewReorthog = useNewReorthog1 .or. nfragments > 1
   if(iproc == 0) print *, "useNewReorthog=", useNewReorthog
end subroutine set_nfragments

!======================================================================
!  subroutine fragmenter
!
!  master routine to break up Lanczos vectors into "fragments"
!  when necessary it breaks up sectors into "subsectors"
!  by breaking along the proton basis Slater Determinants
!
!  The "non-optimized" version simply breaks lanczos vectors by sectors,
!  because that is easiest to implement.
!
!  The "block" version breaks up sectors by proton haiku blocks, which is a natural division
!  and preserves most of the downstream routines intact; the divided sectors must be then redefined
!
!  first count up the number of break points
!
!  CALLED BY:
!    main routine
!
!  SUBROUTINES CALLED:
!    vectorstorageceiling   checks to see if vectors get broken up, and how
!    count_create_fragments_noopt
!        [this creates non-optimized fragments by simply bundling together sectors;
!         later will use optimized fragments with subroutine count_create_fragments ]
!       OR
!    count_create_fragments_blocks
!    reconstruct_sector
!    [map_sectors2noopt_fragments]
!    defaultfragments      only one fragment
!    fragmentvectorsize    sets start and stop for basis in fragments
!    fragmentconjugatesectors  checks which neutron sectors belong (added 7.6.2)
!
subroutine fragmenter
   use flagger
!   use fragments
   use nodeinfo
   use sectors
   use menu_choices

   implicit none

   call vectorstorageceiling  !  check to see if vectors get broken up into fragments
   if(break_vectors .or. test_fragments)then
!	   call test_sector_breaks(83)  ! only for testing
       call count_create_fragments_blocks(.false.)
       call count_create_fragments_blocks(.true.)
       call reconstruct_sectors
!   call test_sector_breaks(84)  ! only for testing
   else
      call defaultfragments
   endif
   call fragmentvectorsize
   call fragmentconjugatesectors
   return
end subroutine fragmenter

!===============================================================

!
!  VECTORSTORAGECEILING
!
!  computes a ceiling on storing pieces of vectors, if any
!
!  initiated 8/2011 by CWJ @ SDSU
!
!  CALLED BY:
!     fragmenter
!
!  CALLS:
!    fragmentsurveysays

subroutine vectorstorageceiling
   use nodeinfo
   use system_parameters
   use basis
   use menu_choices
   use flagger
   use nodeinfo
   use bmpi_mod
!   use fragments
   use io
   implicit none
   integer(8) :: minmaxfragsize  ! largest block found by splitting along haikus
   integer :: ierr
   character :: dummchar
   logical   :: interrflag  ! error in reading an integer

!...... parameters ask2break_vectors and maxfragmentsize_default found in module flagger..
!... MODIFIED IN 7.6.5 to ask for dummy fragmentsize even if not in MPI...

   if(ask2break_vectors)then  ! should always default to this...
	   if(nproc > 1 .or. menu_char=='m' .or. test_fragments)then
         if(iproc==0)then
		  call fragmentsurveysays(minmaxfragsize)
          print*,' '
          print*,' Enter desired limit on fragment size for breaking Lanczos vectors'
		  print*,' (Largest un-splittable block = ',minmaxfragsize,')'
		  if(minmaxfragsize < maxfragmentsize_default)then
             print*,'( Default = ',maxfragmentsize_default,')'
 		  else
			 print*,'( Default = ',minmaxfragsize,')'
!... PUT IN HINT IF Z /= N.....
             if( np(1)/= np(2))print*,' You might be able to reduce this by swapping (valence) protons and neutrons '
  		  end if
          print*,'( Enter 0 to use default )'
          read(5,*,err=1111)maxfragmentsize
          if(maxfragmentsize < 1 .and. minmaxfragsize <= maxfragmentsize_default)maxfragmentsize = maxfragmentsize_default
          if(maxfragmentsize < 1 .and. minmaxfragsize > maxfragmentsize_default)maxfragmentsize = minmaxfragsize

          print*,' Using target Lanczos fragment dimension of ',maxfragmentsize
          print*,' '
		  write(autoinputfile,'(i12,"    !  FRAGMENT SIZE (0 = use default)")')maxfragmentsize

         end if
         call BMPI_Bcast(maxfragmentsize,1,0,icomm,ierr)

       else
!..... ADDED 7.6.5......NEED TO LEAVE ROOM FOR FRAGMENTSIZE IN AUTOINPUT
         if(auto_input)then
	         read(autoinputfile,'(a)')dummchar
		 else
			 write(autoinputfile,'("0     !  LANCZOS FRAGMENT SIZE (0 = use default)")')
		 end if
         maxfragmentsize=maxfragmentsize_default

       end if

   else
      maxfragmentsize=maxfragmentsize_default

   end if

   if(break_vectors_enabled .and. dimbasis > maxfragmentsize)then
       if(menu_char/='m' .and. nproc ==1  .and. .not.test_fragments)then
          break_vectors=.false.
!		  if(.not.test_fragments)then
             print*,' Only one process; cannot break vector '
             print*,' You may run out of memory '
!		 end if
       else
          break_vectors = .true.
          if(iproc==0)then
             print*,' '
             print*,' Basis dimension greater than limit of ',maxfragmentsize
             print*,' Breaking up into fragments '
             print*,' (NOTE: You can set parameter maxfragmentsize in file bmodule_flags.f90 ) '
             print*,' '
          end if
        end if
   else
       break_vectors = .false.
   end if

   return
1111 continue
   if(iproc==0)then
	   print*,' ERROR you did not set the fragmentsize ERROR'
	   write(resultfile,*)' ERROR you did not set the fragmentsize ERROR'
	   write(logfile,*)' ERROR you did not set the fragmentsize ERROR'
	   close(resultfile)
	   close(logfile)
   endif
   call BMPI_ABORT(icomm,101,ierr)
   stop
end subroutine vectorstorageceiling
!===============================================================
!
! routine to survery fragment sizes
! returns the largest haiku-limited fragment
!  added 7.5.0
!
! called by vectorstorageceiling
!
subroutine fragmentsurveysays(maxfragsize)
    use sectors
!    use fragments
    use flagger
    use nodeinfo
    use basis
    use butil_mod

	implicit none
    integer(8) :: maxfragsize
    integer isector
    integer(8) :: npsd, nnsd
    integer(8) :: localpsd  ! how many proton SDs are described by a (right, left) pair of haiku blocks
    integer ic,ics
    integer i
    integer :: nphblocks  ! # of proton haiku blocks

	maxfragsize=0
    do isector = 1,nsectors(1)  ! loop over proton sectors
       nphblocks = xsd(1)%sector(isector)%nhblocks  ! # of proton haiku blocks in this sector
 !.............. COUNT UP HOW MANY CONJUGATE NEUTRON SDs THERE ARE
       nnsd = 0
       do ic = 1,xsd(1)%sector(isector)%ncsectors   ! LOOP OVER CONJUGATE NEUTRON SECTORS
          ics = xsd(1)%sector(isector)%csector(ic)
          nnsd = nnsd + xsd(2)%sector(ics)%nxsd
       end do   ! ic
 !.............. COUNT UP ALONG PROTON HAIKU BLOCKS..................

       do i = 1,nphblocks
          localpsd = xsd(1)%sector(isector)%blockend(i)-xsd(1)%sector(isector)%blockstart(i)+1
		  maxfragsize=bmax(maxfragsize,localpsd*nnsd)
	  end do ! i
  end do ! isector
  return

end subroutine fragmentsurveysays

!===============================================================
!
!    defaultfragments
! routine to set fragment information when only 1 fragment
!
! CALLED BY
!    fragmenter
!
subroutine defaultfragments

   use sectors
!   use fragments
   use basis
   use nodeinfo
   implicit none
   integer :: aerr

   if(iproc==0)print*,' defaulting to 1 Lanczos vector fragment (no break-up) '
   call set_nfragments(1)
   if(allocated(fragmentlist))deallocate(fragmentlist)
   allocate(fragmentlist(1), stat=aerr)
   if(aerr /= 0) call memerror("defaultfragments")
   fragmentlist(1)%nsectors = nsectors(1)
   fragmentlist(1)%ssectorstart = 1
   fragmentlist(1)%ssectorend  = nsectors(1)
   fragmentlist(1)%basisstart = 1
   fragmentlist(1)%basisend   = dimbasis
   fragmentlist(1)%localdim   = dimbasis

   return
end subroutine defaultfragments

!======================================================================

!
!  this version optimizes the distribution by counting up  proton haiku blocks
!  started 7.4.1
!
!  This approach will minimize the changes needed in other routines downstream
!  Note: any sector is constructed by pairs of "left" and "right" "blocks" of "haikus"
!    -- a left haiku is a slater determinant constructed only by s.p. states with jz < 0
!    -- a right haiku is a slater determinant constructed only by s.p. states with jz >= 0
!    -- a block is a set of haikus with the same abelian quantum numbers
!
!  Why this is useful, briefly (see detailed routines in bjumplib1,2,3....)
!  Jumps are created from geneologies, which are created from descendents, which descend from blocks
!  By splitting along blocks, all the jump-creation technology will continue to be used unchanged
!
!  INPUT:
!    create: logical flag to either count (create=F) or create (T)
!
!  CALLED BY:
!     fragmenter (above)
!
subroutine count_create_fragments_blocks(create)
   use sectors
!   use fragments
   use flagger
   use nodeinfo
   use basis
   use butil_mod
   use basis

   implicit none
   logical create
   integer isector
   integer(8) :: localbasiscount
   integer(8) :: nnsd
   integer(8) :: localpsd  ! how many proton SDs are described by a (right, left) pair of haiku blocks
   integer(8) :: ic,ics
   integer(8) :: psd_splitsize
   integer :: nfrag        ! # of fragments
   integer :: nss          ! # of subsectors
   integer(8) :: minfragsize,maxfragsize
   integer(8) :: sumbasis
   integer :: nsplit       ! temporary # of splits in a sector
   integer :: i, nphblocks  ! # of proton haiku blocks
   integer :: aerr
   integer :: ii
   integer :: udbg

   logical, parameter :: writefrag = .true.

   udbg = 0
   if(writefrag .and. iproc==0 .and. create) udbg = 124;
   if(udbg .ne. 0) then
      open(unit=udbg,file='fraginfo.bigstick',status='unknown')
      if(create)write(udbg,*)' # of fragments ',nfragments
      if(create)write(udbg,*)' fragment label     /       dim '
   end if

   nfrag = 1
   nss = 0
   localbasiscount = 0
   maxfragsize = 0
   sumbasis = 0
   minfragsize = maxfragmentsize

   if(iproc==0 .and. .not.create)print*,' Improved fragments (splitting via proton haiku blocks) '

   do isector = 1,nsectors(1)  ! loop over proton sectors
      nphblocks = xsd(1)%sector(isector)%nhblocks  ! # of proton haiku blocks in this sector
      if(create)then
         nsplit =  xsd(1)%sector(isector)%ncuts
         allocate(xsd(1)%sector(isector)%cut(nsplit), stat=aerr)
         if(aerr /= 0) call memerror("count_create_fragments_blocks 1")
      endif
      nsplit = 0
      nss = nss+1   ! start a new sector
!.............. COUNT UP HOW MANY CONJUGATE NEUTRON SDs THERE ARE
      nnsd = 0
      do ic = 1,xsd(1)%sector(isector)%ncsectors   ! LOOP OVER CONJUGATE NEUTRON SECTORS
         ics = xsd(1)%sector(isector)%csector(ic)
         nnsd = nnsd + xsd(2)%sector(ics)%nxsd
      end do   ! ic
!.............. COUNT UP ALONG PROTON HAIKU BLOCKS..................

      do i = 1,nphblocks
         localpsd = xsd(1)%sector(isector)%blockend(i)-xsd(1)%sector(isector)%blockstart(i)+1
         if(localbasiscount + nnsd*localpsd > maxfragmentsize .and. i < nphblocks)then  ! reset
                                                            ! create a new fragment, split off a new sector
            if(i==1) then   ! close out fragment at prior sector
               if(localbasiscount==0)then
                  print*,' SOME PROBLEM WITH FRAGMENTING '
                  print*,' FAILING IN FIRST SECTOR '
                  print*,localbasiscount,nnsd*localpsd,maxfragmentsize
                  print*,isector,nfrag
                  stop
  				   end if
               if(create)then  ! fill current fragment
					   if(nfrag==1)then
                     fragmentlist(nfrag)%basisstart = 1
                     fragmentlist(nfrag)%basisend = localbasiscount
                     fragmentlist(nfrag)%ssectorstart = 1
                  else
                     fragmentlist(nfrag)%basisstart = fragmentlist(nfrag-1)%basisend+1
                     fragmentlist(nfrag)%basisend = localbasiscount+fragmentlist(nfrag-1)%basisend
                     fragmentlist(nfrag)%ssectorstart = 1+fragmentlist(nfrag-1)%ssectorend
                  end if
                  fragmentlist(nfrag)%ssectorend   = nss - 1 ! prev sector in this case  bug  fix KSM 7.6.2
               end if
            else
               nsplit = nsplit+1
               if(create)then  ! fill current fragment
                  fragmentlist(nfrag)%localdim = localbasiscount
                  xsd(1)%sector(isector)%cut(nsplit)= i-1
                  if(nfrag==1)then
                      fragmentlist(nfrag)%basisstart = 1
                      fragmentlist(nfrag)%basisend = localbasiscount
                      fragmentlist(nfrag)%ssectorstart = 1
                  else
                      fragmentlist(nfrag)%basisstart = fragmentlist(nfrag-1)%basisend+1
                      fragmentlist(nfrag)%basisend = localbasiscount+fragmentlist(nfrag-1)%basisend
                      fragmentlist(nfrag)%ssectorstart = 1+fragmentlist(nfrag-1)%ssectorend
                  end if
                  fragmentlist(nfrag)%ssectorend   = nss
               end if
               nss = nss+1      ! start new subsector
            end if
            maxfragsize = bmax(maxfragsize, localbasiscount)
            minfragsize = bmin(minfragsize, localbasiscount)
            nfrag = nfrag +1  ! start new fragment
            sumbasis = sumbasis+localbasiscount    ! update basis so far
            if(udbg /= 0)then
               write(udbg,*)nfrag-1,localbasiscount
            end if
            localbasiscount = 0
         end if
         localbasiscount = localbasiscount+nnsd*localpsd
         xsd(1)%sector(isector)%ncuts = nsplit
      end do

      ! If we just finished the last sector
      if(isector==nsectors(1))then
         maxfragsize = bmax(maxfragsize, localbasiscount)
         minfragsize = bmin(minfragsize, localbasiscount)
         sumbasis = sumbasis+localbasiscount
         if(create)then
            if(nfrag==1)then
               fragmentlist(nfrag)%basisstart = 1
               fragmentlist(nfrag)%basisend = localbasiscount
               fragmentlist(nfrag)%ssectorstart = 1
            else
               fragmentlist(nfrag)%basisstart = fragmentlist(nfrag-1)%basisend+1
               fragmentlist(nfrag)%basisend = localbasiscount+fragmentlist(nfrag-1)%basisend
               fragmentlist(nfrag)%ssectorstart = 1+fragmentlist(nfrag-1)%ssectorend
            end if
            fragmentlist(nfrag)%ssectorend   = nss
         end if
      end if
   end do ! isector
   if(create)then
      do ii=1,nfrag
         fragmentlist(ii)%localdim = fragmentlist(ii)%basisend - fragmentlist(ii)%basisstart + 1
      end do
   end if
   if(iproc==0 .and. .not. create)then
      print*,' There are a total of ',nfrag,' Lanczos fragments'
      print*,' from a new total of ',nss,' sectors '
      print*,' min/max fragment dimension = ',minfragsize,maxfragsize
   end if
   if(udbg /= 0) then ! KSM
      write(udbg, *) "DEBUG - fragment dump1"
      do ii = 1, nfragments
         write(udbg, "(I3,A,I12,A,I12,A)") ii, ":(", fragmentlist(ii)%basisstart, ",", fragmentlist(ii)%basisend, ")"
         write(udbg, "(A,I0,A,I0)") "   sectorstart=", fragmentlist(ii)%ssectorstart, ", sectorend=", fragmentlist(ii)%ssectorend
         write(udbg, "(A,I0)")      "   nsectors=", fragmentlist(ii)%nsectors
      end do
   end if
   if(udbg /= 0) then
      close(udbg)
      udbg = 0
   end if
   if(sumbasis/=dimbasis)then
      print*,' Mismatch in basis dimension '
      stop
   end if

   if(.not.create)then
      call set_nfragments(nfrag)
      if(.not.allocated(fragmentlist)) then
         allocate(fragmentlist(nfragments) , stat=aerr)
         if(aerr /= 0) call memerror("count_create_fragments_blocks 10")
      end if
   end if
   return
end subroutine count_create_fragments_blocks

!======================================================================
!
! when sectors are subdivided, reconstruct proton sectors
! added 7.4.1; used after count_create_fragments_blocks
!
!  CALLED BY:
!    fragmenter

subroutine reconstruct_sectors
   use sectors
   use basis
   use nodeinfo
   implicit none
   type(mastersect) :: draftpsd
   integer nnewsectors
   integer :: isector,ics,jsector,jcs,ncs
   integer :: icut
   integer :: nblocks,iblock,jblock
   integer :: aerr
   integer(8) :: npsd

!......... COUNT UP NEW SECTORS
   nnewsectors = 0
   do isector = 1,nsectors(1)
      nnewsectors = nnewsectors + xsd(1)%sector(isector)%ncuts + 1
   end do
   if(nnewsectors == nsectors(1)) return   ! nothing to do

!........ CREATE TEMPORARY DRAFT OF SECTORS......
   allocate( draftpsd%sector(nnewsectors) , stat=aerr)
   if(aerr /= 0) call memerror("reconstruct_sectors 1")

!............ ISECTOR labels old sectors, while NNEWSECTORS labels new sectors......

   nnewsectors = 0
   do isector = 1,nsectors(1)
      nnewsectors = nnewsectors+1   ! INCREMENT # OF SECTORS

      draftpsd%sector(nnewsectors)%xsdstart = xsd(1)%sector(isector)%xsdstart   ! SET BASIC QUANTUM NUMBERS
      draftpsd%sector(nnewsectors)%jzX      = xsd(1)%sector(isector)%jzX
      draftpsd%sector(nnewsectors)%parX     = xsd(1)%sector(isector)%parX
      draftpsd%sector(nnewsectors)%WX       = xsd(1)%sector(isector)%wX
!      draftpsd%sector(nnewsectors)%parentsector = isector
!................ COPY OVER CONJUGATE SECTORS...................
      draftpsd%sector(nnewsectors)%ncsectors= xsd(1)%sector(isector)%ncsectors
      allocate ( draftpsd%sector(nnewsectors)%csector( draftpsd%sector(nnewsectors)%ncsectors ), stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 10")
      do ics = 1, xsd(1)%sector(isector)%ncsectors
         draftpsd%sector(nnewsectors)%csector(ics) = xsd(1)%sector(isector)%csector(ics)
      end do

      if( xsd(1)%sector(isector)%ncuts > 0)then
         do icut = 1, xsd(1)%sector(isector)%ncuts
!............. BEFORE STARTING NEW "SECTOR" FINISH DATA FOR CURRENT ONE........
            nblocks  = xsd(1)%sector(isector)%cut(icut)
            if(icut/=1)then
               nblocks = nblocks - xsd(1)%sector(isector)%cut(icut-1)
            end if
            draftpsd%sector(nnewsectors)%nhblocks  = nblocks
            draftpsd%sector(nnewsectors)%xsdend   = xsd(1)%sector(isector)%blockend( xsd(1)%sector(isector)%cut(icut)  )
            draftpsd%sector(nnewsectors)%nxsd = draftpsd%sector(nnewsectors)%xsdend - draftpsd%sector(nnewsectors)%xsdstart+1
            allocate (draftpsd%sector(nnewsectors)%rhblock(nblocks), draftpsd%sector(nnewsectors)%lhblock(nblocks), stat=aerr)
            if(aerr /= 0) call memerror("reconstruct_sectors 20")
            allocate (draftpsd%sector(nnewsectors)%blockstart(nblocks), draftpsd%sector(nnewsectors)%blockend(nblocks), stat=aerr)
            if(aerr /= 0) call memerror("reconstruct_sectors 21")

            do iblock = 1,nblocks
               jblock = iblock
               if(icut > 1)jblock = jblock + xsd(1)%sector(isector)%cut(icut-1)
               draftpsd%sector(nnewsectors)%rhblock(iblock) = xsd(1)%sector(isector)%rhblock(jblock)
               draftpsd%sector(nnewsectors)%lhblock(iblock) = xsd(1)%sector(isector)%lhblock(jblock)
               draftpsd%sector(nnewsectors)%blockstart(iblock) = xsd(1)%sector(isector)%blockstart(jblock)
               draftpsd%sector(nnewsectors)%blockend(iblock) = xsd(1)%sector(isector)%blockend(jblock)
            end do
            draftpsd%sector(nnewsectors)%nxsd = draftpsd%sector(nnewsectors)%blockend(nblocks)- &
                                                draftpsd%sector(nnewsectors)%blockstart(1)+1
            draftpsd%sector(nnewsectors)%xsdend = draftpsd%sector(nnewsectors)%xsdstart + draftpsd%sector(nnewsectors)%nxsd -1
!................... START NEW SECTOR..........................
            nnewsectors = nnewsectors+1
			draftpsd%sector(nnewsectors)%xsdstart= draftpsd%sector(nnewsectors-1)%xsdend +1
            draftpsd%sector(nnewsectors)%jzX      = xsd(1)%sector(isector)%jzX    ! SET BASIC QUANTUM NUMBERS
            draftpsd%sector(nnewsectors)%parX     = xsd(1)%sector(isector)%parX
            draftpsd%sector(nnewsectors)%WX       = xsd(1)%sector(isector)%wX
!           draftpsd%sector(nnewsectors)%parentsector = isector
!................ COPY OVER CONJUGATE SECTORS...................
            draftpsd%sector(nnewsectors)%ncsectors= xsd(1)%sector(isector)%ncsectors
            allocate ( draftpsd%sector(nnewsectors)%csector( draftpsd%sector(nnewsectors)%ncsectors ), stat=aerr)
            if(aerr /= 0) call memerror("reconstruct_sectors 30")
            do ics = 1, xsd(1)%sector(isector)%ncsectors
                draftpsd%sector(nnewsectors)%csector(ics) = xsd(1)%sector(isector)%csector(ics)
            end do
         end do
      end if
!............. FINISH DATA FOR CURRENT "SECTOR"...............

      draftpsd%sector(nnewsectors)%xsdend        = xsd(1)%sector(isector)%xsdend
      draftpsd%sector(nnewsectors)%nxsd = draftpsd%sector(nnewsectors)%xsdend - draftpsd%sector(nnewsectors)%xsdstart+1
      nblocks = xsd(1)%sector(isector)%nhblocks
      if(xsd(1)%sector(isector)%ncuts /=0)then
         nblocks = nblocks - xsd(1)%sector(isector)%cut(xsd(1)%sector(isector)%ncuts)
      end if
      draftpsd%sector(nnewsectors)%nhblocks  = nblocks
      allocate (draftpsd%sector(nnewsectors)%rhblock(nblocks), draftpsd%sector(nnewsectors)%lhblock(nblocks), stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 40")
      allocate (draftpsd%sector(nnewsectors)%blockstart(nblocks), draftpsd%sector(nnewsectors)%blockend(nblocks), stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 41")
      do iblock = 1,nblocks
               jblock = iblock
               if(xsd(1)%sector(isector)%ncuts > 0)jblock = jblock + xsd(1)%sector(isector)%cut(xsd(1)%sector(isector)%ncuts)
               draftpsd%sector(nnewsectors)%rhblock(iblock) = xsd(1)%sector(isector)%rhblock(jblock)
               draftpsd%sector(nnewsectors)%lhblock(iblock) = xsd(1)%sector(isector)%lhblock(jblock)
               draftpsd%sector(nnewsectors)%blockstart(iblock) = xsd(1)%sector(isector)%blockstart(jblock)
               draftpsd%sector(nnewsectors)%blockend(iblock) = xsd(1)%sector(isector)%blockend(jblock)
      end do

   end do  ! isector
!...... TEST THE SDs ARE ALIGNED............
   npsd = 0
   do isector = 1,nnewsectors
      npsd = npsd+ draftpsd%sector(isector)%xsdend - draftpsd%sector(isector)%xsdstart+1
	  if(draftpsd%sector(isector)%xsdend - draftpsd%sector(isector)%xsdstart+1 /= draftpsd%sector(isector)%nxsd )then
	     if(iproc==0)print*,' mismatch in new sectors ', isector
		 stop
	  end if
	  if(isector > 1)then
		  if( draftpsd%sector(isector)%xsdstart /= draftpsd%sector(isector-1)%xsdend+1)then
			  if(iproc==0)then
				  print*,' Mismatch in start/stop of proton sds ',isector
				  print*,draftpsd%sector(isector)%xsdstart, draftpsd%sector(isector-1)%xsdend+1
			  end if
			  stop
		  end if
	  end if

   end do
   if(npsd /= nxsd(1))then
	   if(iproc==0)print*,' Mismatch in new # of proton SDs ',npsd,nxsd(1)
	   stop
   end if
!........ CLEAR OUT OLD SECTORS
   do isector = 1,nsectors(1)
      nullify(xsd(1)%sector(isector)%csector)
      nullify(xsd(1)%sector(isector)%rhblock,xsd(1)%sector(isector)%lhblock )
      nullify(xsd(1)%sector(isector)%blockstart,xsd(1)%sector(isector)%blockend )
      nullify(xsd(1)%sector(isector)%cut)
   end do
   nullify(xsd(1)%sector)

!........ COPY DRAFT BACK INTO SECTORS
   nsectors(1) = nnewsectors
   allocate( xsd(1)%sector(nnewsectors) , stat=aerr)
   if(aerr /= 0) call memerror("reconstruct_sectors 50")
   do isector = 1,nsectors(1)
      xsd(1)%sector(isector)%xsdstart = draftpsd%sector(isector)%xsdstart
      xsd(1)%sector(isector)%xsdend   = draftpsd%sector(isector)%xsdend
      xsd(1)%sector(isector)%nxsd     = draftpsd%sector(isector)%nxsd
      xsd(1)%sector(isector)%jzX      = draftpsd%sector(isector)%jzX
      xsd(1)%sector(isector)%parX     = draftpsd%sector(isector)%parX
      xsd(1)%sector(isector)%WX       = draftpsd%sector(isector)%wX
      xsd(1)%sector(isector)%ncsectors= draftpsd%sector(isector)%ncsectors
      xsd(1)%sector(isector)%nhblocks = draftpsd%sector(isector)%nhblocks
!      xsd(1)%sector(isector)%parentsector = draftpsd%sector(isector)%parentsector
      allocate (xsd(1)%sector(isector)%csector(  draftpsd%sector(isector)%ncsectors ) , stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 60")
      allocate (xsd(1)%sector(isector)%lhblock(  draftpsd%sector(isector)%nhblocks ) , stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 61")
      allocate (xsd(1)%sector(isector)%rhblock(  draftpsd%sector(isector)%nhblocks ) , stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 62")
      allocate (xsd(1)%sector(isector)%blockstart(  draftpsd%sector(isector)%nhblocks ) , stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 63")
      allocate (xsd(1)%sector(isector)%blockend(  draftpsd%sector(isector)%nhblocks ) , stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 64")
      do ics = 1, draftpsd%sector(isector)%ncsectors
          xsd(1)%sector(isector)%csector(ics) = draftpsd%sector(isector)%csector(ics)
      end do
      do ics = 1, draftpsd%sector(isector)%nhblocks
          xsd(1)%sector(isector)%rhblock(ics) = draftpsd%sector(isector)%rhblock(ics)
          xsd(1)%sector(isector)%lhblock(ics) = draftpsd%sector(isector)%lhblock(ics)
          xsd(1)%sector(isector)%blockstart(ics) = draftpsd%sector(isector)%blockstart(ics)
          xsd(1)%sector(isector)%blockend(ics) = draftpsd%sector(isector)%blockend(ics)
      end do  ! ics
   end do  ! isector

!........ ASSIGN NEW CONJUGATE SECTORS TO NEUTRONS

   do jsector = 1,nsectors(2)
      ncs = 0
      do isector = 1,nsectors(1)
         do ics = 1,xsd(1)%sector(isector)%ncsectors
            if(xsd(1)%sector(isector)%csector(ics) == jsector)then
               ncs = ncs+1
               exit
            end if
         end do ! ics
      end do
      xsd(2)%sector(jsector)%ncsectors = ncs
      if(ncs ==0)cycle
      nullify( xsd(2)%sector(jsector)%csector )
      allocate ( xsd(2)%sector(jsector)%csector(ncs) , stat=aerr)
      if(aerr /= 0) call memerror("reconstruct_sectors 70")
      ncs = 0
      do isector = 1,nsectors(1)
         do ics = 1,xsd(1)%sector(isector)%ncsectors
            if(xsd(1)%sector(isector)%csector(ics) == jsector)then
                ncs = ncs+1
                xsd(2)%sector(jsector)%csector(ncs) = isector
                exit
            end if
         end do ! ics
      end do

   end do  ! jsector
   return
end subroutine reconstruct_sectors
!======================================================================
!
!  check new (7.4.1) break up of sectors for fragments is okay
!  writes to file fort.unitout; run two different ways to compare
!
subroutine test_sector_breaks(unitout)
   use sectors
   use basis

   implicit none
   integer unitout
   integer(8) :: basiscount
   integer    :: ips,ins
   integer    :: pblock,nblock
   integer    :: ncs   ! # of conjugate sectors
   integer    :: ics
   integer(8) :: pblocksize, nblocksize,nlocal
   basiscount = 0
   if(unitout ==83)then
       open(unit=83,file='sector.in',status='unknown')
   else
       open(unit=84,file='sector.out',status='unknown')

   endif
   do ips = 1,nsectors(1)  ! loop over proton sectors
      ncs = xsd(1)%sector(ips)%ncsectors   ! fetch the number of conjugate neutron sectors
	  nlocal = 0
      do pblock = 1,xsd(1)%sector(ips)%nhblocks  ! within this sector, loop over the pairs of proton haiku blocks
          pblocksize = int(xsd(1)%sector(ips)%blockend(pblock) - xsd(1)%sector(ips)%blockstart(pblock)+1,8)
                 ! fetch the proton contribution
		  nlocal = nlocal+pblocksize
		  write(unitout,*)pblock,pblocksize
          do ics = 1,ncs
              ins = xsd(1)%sector(ips)%csector(ics)
              do nblock = 1,xsd(2)%sector(ins)%nhblocks
                 nblocksize = int(xsd(2)%sector(ins)%blockend(nblock) - xsd(2)%sector(ins)%blockstart(nblock)+1,8)
                 basiscount = basiscount + pblocksize*nblocksize
!				 nlocal = nlocal + pblocksize*nblocksize
!                 write(unitout,*)basiscount

              end do ! nblock

          end do !ics

      end do ! pblock
	  write(unitout,*)ips,' proton sector SDs ',nlocal, xsd(1)%sector(ips)%nxsd

   end do
   write(unitout,*)basiscount,dimbasis
   if(basiscount /= dimbasis)then  ! error flag
      print*, '  v v v v v v v v v v v v v v v v v v v v v '
      print*,' Oh noes! Mismatch in basis ! '
      print*,dimbasis,basiscount
      print*, '  ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ '
      stop
   end if
   close(unitout)
   return

end subroutine test_sector_breaks
!======================================================================

!
!  fragmentvectorsize
!     sets base start/stop for a given fragment
!
! CALLED BY:
!  fragmenter
!
subroutine fragmentvectorsize
!   use fragments
   implicit none
   integer :: ifrag
   integer :: aerr

   allocate(basestart(nfragments), basestop(nfragments) , stat=aerr)
   if(aerr /= 0) call memerror("fragmentvectorsize")

   do ifrag = 1,nfragments
       basestart(ifrag) = fragmentlist(ifrag)%basisstart
       basestop(ifrag)  = fragmentlist(ifrag)%basisend

   end do
   return

end subroutine fragmentvectorsize

!======================================================================

!
!  printopfragstat
!     sets base start/stop for a given fragment
!
! CALLED BY:
!  fragmenter
!
subroutine printopfragstat
   use nodeinfo
   use operation_stats
!   use fragments
   implicit none
   integer :: ifrag, jfrag

   if(iproc==0)then
	   do ifrag = 1,nfragments
		   do jfrag=1,nfragments
			   print*,ifrag,jfrag,' nodes: ', opfragstat(ifrag,jfrag)%nnodes
		   end do
	   end do
   end if
end subroutine printopfragstat
!===============================================================
!
!  initiated July 2014 @ LLNL by CWJ
!
!  routines to guide distributing work across and between
!  Lanczos vectors "fragments."
!
!==========================================================================================
!   divide_procs_over_fragments
!
! routine to divide up available nodes over fragments
!
! main goal is to set up opfragstat(if,ff)%nnodes
! which assign nodes to initial, final fragments if,ff
!
! INPUT:
!    nnodesx  : sample # of MPI processes
!                This might differ the actual # of nodes available,
!                if one needs to reserve nodes
!    assignnodes: flag to actually assign nodes
!
! CALLED BY:
!   distro_opbundles_over_fragments
!
subroutine divide_procs_over_fragments(nnodesx,assignnodes)
!   use fragments
   use operation_stats
   use nodeinfo
!   use tribution
   use verbosity
   use butil_mod
   implicit none
!.......INPUT........
   integer :: nnodesx  !  # of nodes made available
   logical :: assignnodes  ! flag to assign nodes

!....... INTERNAL............
   real(8) :: opfrac      ! frac of operations
   integer :: ifrag, jfrag,i,j
   real(8) :: optotall,opavg,optotsingle
   integer :: nodesum, nodediff,nodesingle,nodeadd
   real(8) :: diffcf
   ! KSM  - just use ifrag, jfrag for position:   integer :: inodecf
   integer :: si, sj  ! used for search for biggest discrepency
!.... ESTIMATE EFFICIENCY.....
   real(8) :: maxnodeops, avgnodeops

!...... KEEP TRACK OF WHICH FRAGMENT BLOCKS HAVE BIGGEST DISCREPANCIES....
   ! KSM:  real(8) :: opdisc(nfragments*nfragments)
   real(8) :: opdisca(nfragments, nfragments)
   logical :: localverbose  = .false.   ! used for debugging

   real :: assignratio = 0.9

!........... DEFAULT IF NO FRAGMENTS.........

   opdisca = 0.0  ! initialize

!............. COMPUTE TOTAL # OF OPS...............
   optotall = 0.0d0

   do ifrag = 1,nfragments
      do jfrag = 1,nfragments
         ! shan: using new weight op
         if (iproc == 0 .and. assignnodes) write(*, 111) "FRAGOP ", ifrag, jfrag, " stateop ", opfragstat(ifrag,jfrag)%optot, &
              " fragop ", frag2fragops(ifrag, jfrag)

         opfragstat(ifrag,jfrag)%optot = frag2fragops(ifrag, jfrag)
          optotall = optotall + opfragstat(ifrag,jfrag)%optot
   end do
   end do

111 format(a7, i6, i6, a9, f16.2, a8, f16.2)

   opavg = optotall/real(nnodesx,8)
!!............ IF ALREADY ASSIGNED.................

   if(allocated(nodal))then
!	   if(iproc==0)print*,' nodes already assigned '
!	   if(iproc==0)print*,opfragstat(:,:)%nnodes
	   return

   end if
!   print*,' operation totals for TESTING ',optotall,opavg
!........... COUNT THE NUMBER OF SINGLE-NODE FRAGMENT BLOCKS.......
   nodesingle= 0
   optotsingle = 0.0d0
   do ifrag = 1,nfragments
      do jfrag = 1,nfragments
          if(opfragstat(ifrag,jfrag)%optot < opavg .and. opfragstat(ifrag,jfrag)%optot> 0.0d0)then
              opfragstat(ifrag,jfrag)%nnodes=1
              ! KSM - already set to 0:   opdisc(nfragments*(ifrag-1)+jfrag) = 0.0d0
              nodesingle = nodesingle+1
              optotsingle=optotsingle + opfragstat(ifrag,jfrag)%optot
          else
            opfragstat(ifrag,jfrag)%nnodes=0
          end if
      end do
   end do

   if (nnodesx * assignratio -nodesingle < 0) then
      assignratio = 1.0
   endif

   if(iproc==0 .and. verbose_fragments)then
        print*,nodesingle,' fragment blocks on single MPI nodes '
        print*,' fraction of operations on single MPI nodes = ',optotsingle, optotall, optotsingle/optotall, assignratio
   end if

!........... FIRST ATTEMPT AT A DISTRIBUTION........
   opavg = (optotall-optotsingle)/real(nnodesx-nodesingle,8)
   nodesum = 0
   do ifrag = 1,nfragments
      do jfrag = 1,nfragments
          if(opfragstat(ifrag,jfrag)%nnodes==1)then
              nodesum = nodesum+1
              opfragstat(ifrag,jfrag)%nnodes=1
              ! KSM:  already 0:  opdisc(nfragments*(ifrag-1)+jfrag) = 0.0d0
              cycle
          endif
          if(opfragstat(ifrag,jfrag)%optot== 0.0d0)then
              opfragstat(ifrag,jfrag)%nnodes=0
              ! KSM:  already 0:  opdisc(nfragments*(ifrag-1)+jfrag) = 0.0d0
              cycle
          end if
          opfrac = opfragstat(ifrag,jfrag)%optot/(optotall-optotsingle)
          opfragstat(ifrag,jfrag)%nnodes = nint(opfrac*real(nnodesx * assignratio -nodesingle,8))

          if(opfragstat(ifrag,jfrag)%nnodes==0 .and. opfragstat(ifrag,jfrag)%optot>0)then
              opfragstat(ifrag,jfrag)%nnodes=1
              ! KSM:  already 0:  opdisc(nfragments*(ifrag-1)+jfrag) = 0.0d0
          else
              ! opdisc(nfragments*(ifrag-1)+jfrag) = opfragstat(ifrag,jfrag)%optot/opavg-opfragstat(ifrag,jfrag)%nnodes
             !opdisca(ifrag,jfrag) = opfragstat(ifrag,jfrag)%optot/opavg-opfragstat(ifrag,jfrag)%nnodes
             opdisca(ifrag,jfrag) = (opfragstat(ifrag,jfrag)%optot/opavg- &
			                        opfragstat(ifrag,jfrag)%nnodes)/opfragstat(ifrag,jfrag)%nnodes
          end if
          nodesum = nodesum + opfragstat(ifrag,jfrag)%nnodes
      end do
   end do

!........ TOO FEW/TOO MANY NODES PASSED OUT........
   nodediff = nnodesx - nodesum
! -- CHECK DISTRIBUTION ---- ADDED in 7.4.9.....
   if(localverbose) then
      call printopfragstat
      if(iproc == 0) print*,' The difference is ',nodediff
   end if
!......... NOW ADD/SUBTRACT NODES FROM FRAGMENT BLOCKS.........................
!
!  because we won't have exactly used up all the nodes, find out how many extra
!  (or missing) nodes there are; and add or subtract from fragment blocks
!
! ADDED 7.4.9: cannot subtract nodes if only one node assigned

   if (iproc == 0) print *, "NODEDIFF ", nnodesx, nodesum, assignratio

   if(nodediff == 0)return
   if(nodediff < 0)then
		   nodeadd = -1
   else
		   nodeadd = +1
   end if

   do i= 1,abs(nodediff)
	!....... SEARCH FOR LARGEST DIFFERENCE........................
		  diffcf = 0.00
		  ifrag = 0
		  jfrag = 0
		  do si = 1, nfragments
			 do sj = 1, nfragments
				if( opfragstat(si,sj)%nnodes  < 2 .and. nodeadd < 0)cycle  !don't subtract if only 1 nodes assigned
				if(opdisca(si, sj)*nodeadd > diffcf) then
				   diffcf = opdisca(si, sj) * nodeadd
				   ! save position
				   ifrag = si
				   jfrag = sj
				end if
			 end do
		  end do
		  if(jfrag .eq. 0) then
			  print *, "iproc=", iproc, ", Can't find suitable fragment pair to  add/remove nodes"
			  print*,nfragments,' fragments '
			  print*,nodediff,' nodes to be reassigned out of ',nprocs,' with ',nprocs-nfragments**2,' available '
			  stop
		  end if
		  opfragstat(ifrag,jfrag)%nnodes= opfragstat(ifrag,jfrag)%nnodes+nodeadd
		  ! KSM:  opdisc(nfragments*(ifrag-1)+jfrag) = opfragstat(ifrag,jfrag)%optot/opavg-opfragstat(ifrag,jfrag)%nnodes
		  ! opdisca(ifrag, jfrag) = opfragstat(ifrag,jfrag)%optot/opavg-opfragstat(ifrag,jfrag)%nnodes
                  opdisca(ifrag, jfrag) = (opfragstat(ifrag,jfrag)%optot/opavg- &
				  opfragstat(ifrag,jfrag)%nnodes)/opfragstat(ifrag,jfrag)%nnodes
   end do
   if(iproc==0)write(24,*)'     init frag / final frag  / # nodes assigned'
	   do ifrag = 1,nfragments
              do jfrag = 1,nfragments
        if(iproc==0)write(24,*)ifrag,jfrag,opfragstat(ifrag,jfrag)%nnodes
              end do
	   end do
	   if(.not.assignnodes)return
	!....... NEXT STEP (AT LEAST WHEN MODELING) ASSIGN FRAGMENTS TO NODES.....

	   i = -1
   if(iproc==0)write(24,*)'       node /    init frag  / final frag  / operations load '
	   maxnodeops = 0.0
	   avgnodeops = 0.0
	   do ifrag = 1,nfragments
		  do jfrag = 1,nfragments
			  avgnodeops = avgnodeops + opfragstat(ifrag,jfrag)%optot
			  do j = 1, opfragstat(ifrag,jfrag)%nnodes
				   i = i+1
				   maxnodeops = bmax(maxnodeops, opfragstat(ifrag,jfrag)%optot/real(opfragstat(ifrag,jfrag)%nnodes,8))

               if(iproc==0)write(24,*)i,ifrag,jfrag,opfragstat(ifrag,jfrag)%optot/opfragstat(ifrag,jfrag)%nnodes
			  end do
		  end do ! jfrag
	   end do  ! ifrag
   avgnodeops = avgnodeops/float(nprocs)
!   if(iproc==0)print*,' estimated efficiency = ',avgnodeops/maxnodeops
   close(24)
   return
end subroutine divide_procs_over_fragments
!==========================================================================================
!
!  sets variable "ibelong" in xsd(it)%sector to denote if a particular sector
!  belongs to fragments residing on a MPI process
!
!  added in 7.5.9
!
subroutine setsectorpossession(iprocs)
	use nodeinfo
	use sectors
!	use fragments
	implicit none
	integer :: iprocs   ! dummy label for MPI processes
	integer :: it   ! label of species
	integer(8) :: is,cs,ics  ! label of sector
	integer :: ifrag,ffrag  ! fragment labels

	do it =1,2
		if(nfragments==1)then    ! ALL SECTORS BELONG
			do is = 1,nsectors(it)
				xsd(it)%sector(is)%ibelong=.true.
			end do ! is
			cycle
		end if

		do is = 1,nsectors(it)        ! DEFAULT
			xsd(it)%sector(is)%ibelong=.false.
        end do

	end do ! it
!-------------- "POSSESS" proton sectors

	ifrag = nodal(iprocs)%ifragment
	ffrag = nodal(iprocs)%ffragment
	do is = 1,fragmentlist(ifrag)%ssectorstart,fragmentlist(ifrag)%ssectorend
		xsd(1)%sector(is)%ibelong=.true.
	end do ! is
	do is = 1,fragmentlist(ffrag)%ssectorstart,fragmentlist(ffrag)%ssectorend
		xsd(1)%sector(is)%ibelong=.true.
	end do ! is

!......... NOW POSSESS NEUTRON FRAGMENTS-----------
    do is = 1,nsectors(1)
		if(xsd(1)%sector(is)%ibelong)then
			do ics = 1,xsd(1)%sector(is)%ncsectors
				cs = xsd(1)%sector(is)%csector(ics)
				xsd(2)%sector(cs)%ibelong=.true.
			end do ! ics
		end if
	end do
	return
end subroutine setsectorpossession

!======================================================================
!
! added to 7.6.2
!
! determines which conjugate ("neutron") sectors belongs to a particular fragment
! and confirms contiguity
!
!  CALLED BY fragmenter

subroutine fragmentconjugatesectors
	use sectors
!	use fragments
    implicit none

    integer ifrag
	logical, allocatable :: sectornocc(:)   ! logical flag to denote occupied neutron sectors
	integer(8) :: is,cs,ics,firstsector,lastsector

	allocate(sectornocc( nsectors(2)))

	do ifrag = 1,nfragments
		sectornocc = .false.
		do is = 1,fragmentlist(ifrag)%ssectorstart,fragmentlist(ifrag)%ssectorend
			do ics = 1,xsd(1)%sector(is)%ncsectors
				cs = xsd(1)%sector(is)%csector(ics)
                sectornocc(cs)=.true.
			end do ! ics
		end do ! is
!----------- FIND LIMITS
        firstsector = nsectors(2)
		lastsector  = 1
		do is = 1,nsectors(2)
			if(sectornocc(is))then
				firstsector = min(firstsector,is)
				lastsector  = max(lastsector, is)
			end if
		end do

!--------- CHECK FOR CONTIGUITY--------------------

        do is = firstsector,lastsector
	        if(.not.sectornocc(is))then
				print*,' OOPS! Non continuous neutron sectors in fragment ',ifrag
				print*,' first, last sectors = ',firstsector,lastsector
				print*,' but neutrons sector ',is,' not included '
				stop
			end if
        end do
!----------- SET LIMITS ---------------
        fragmentlist(ifrag)%csectorstart = firstsector
        fragmentlist(ifrag)%csectorend =   lastsector

	end do  ! ifrag

    deallocate(sectornocc)
    return
end	subroutine fragmentconjugatesectors

!===============================================
! debug print of tribution data structure.
! should print on all processes
subroutine printtribution
!   use fragments
   use nodeinfo
   use io
   integer :: i

   if(iproc /= 0) return
   write(logfile,*) "nprocs = ", nprocs
   write(logfile,*)"nodeal(proc) = {ifragment, ffragment }"
   do i = 0, nprocs-1
      write(logfile,*) "nodeal(", i, ") =", "{", nodal(i)%ifragment, ", ", nodal(i)%ffragment, "}"
   end do
   return
end subroutine printtribution

!===============================================================
!   setnodal!
!  assigns MPI processes / compute nodes to initial and final fragments
!  sets up "nodal" and fills data
!
!  IMPORTANT: opfragstat(ifrag,ffrag)%nnodes must be previous set
!            done in distro_opbundles_over_fragments and
!             and in divide_procs_over_fragments
!
!  CALLED BY:
!    master_para_distribute
!    density1b_from_oldwfn
!    overlap
!
subroutine setnodaltribution
!   use fragments
   use nodeinfo
   use io
   use reporter
   implicit none
   integer ifrag,jfrag
   integer i,j,ii
   integer :: aerr

   allocate(nodal(0:nprocs-1), stat=aerr )
   if(aerr /= 0) call memerror("setnodaltribution")
   if(iproc==0 .and. nfragments > 1) then
      open(unit=fragfile,file='fraginfo.bigstick',status='unknown')
      write(fragfile,*)' # of fragments ',nfragments
      write(fragfile,*)' fragment label     /       dim '
      write(logfile,*)' # of fragments ',nfragments
      write(logfile,*)' fragment label     /       dim '
	  do ii = 1,nfragments
		  write(fragfile,*)ii,fragmentlist(ii)%localdim
		  write(logfile,*)ii,fragmentlist(ii)%localdim

	  end do
      do ii = 1, nfragments
         write(fragfile, "(I3,A,I12,A,I12,A)") ii, ":(", fragmentlist(ii)%basisstart, ",", fragmentlist(ii)%basisend, ")"
         write(fragfile, "(A,I0,A,I0)") "   sectorstart=", fragmentlist(ii)%ssectorstart, ", sectorend=", &
		                                                          fragmentlist(ii)%ssectorend
         write(fragfile, "(A,I0)")      "   nsectors=", fragmentlist(ii)%nsectors
      end do
   end if

   i = -1
   if(nprocs==1)then  ! KSM - change from nproc.  Was clearing MPI data during modeling
      opfragstat(:,:)%nnodes = 1
   end if
   do ifrag = 1,nfragments
      do jfrag = 1,nfragments
          if (iproc == 0) then
            write(logfile,*) "opfragstat(", ifrag, ", ", jfrag, ")%nnodes = ", opfragstat(ifrag,jfrag)%nnodes
          end if
		  if(iproc==0 .and.opfragstat(ifrag,jfrag)%nnodes >0)then
			  write(fragfile,*)' frags ',ifrag,'->',jfrag,' on procs ',i+1,i+opfragstat(ifrag,jfrag)%nnodes
			  write(logfile,*)' frags ',ifrag,'->',jfrag,' on procs ',i+1,i+opfragstat(ifrag,jfrag)%nnodes

		  end if
          do j = 1, opfragstat(ifrag,jfrag)%nnodes
            i = i+1
            if(i >= nprocs) then    ! nprocs works even for modeling
               if(iproc == 0) then
                  print *, "setnodaltribution: nprocs=", nprocs, ", is too small"
                  write(logfile,*) "setnodaltribution: nprocs=", nprocs, ", is too small"
                  call printstacktrace
               end if
               stop 1
            end if
            nodal(i)%ifirst = j .eq. 1 .and. jfrag .eq. 1 ! marks one node which reads ifrag
            nodal(i)%ifragment=ifrag
            nodal(i)%ffragment=jfrag
!            nodal(i)%sfragment=0  ! KSM: init for debug checking
!            nodal(i)%frac_ops=0.0  ! KSM: init for debug checking
            ! don't overwrite ifirst here, look up 5 lines!
          end do
       end do
   end do

   call checknodaltribution
   if(nfragments > 1)close(fragfile)
   return
end subroutine setnodaltribution

!===============================================================
!
!  checks assignment of MPI processes/nodes to fragments against assignment of bundles
!  added 7.6.8
!
subroutine checknodaltribution
    use nodeinfo
    use io
	use opbundles
	use bmpi_mod
    implicit none
    integer ifrag,ffrag
	integer :: jprocs
    integer ibundle
	logical :: muststop

	integer :: ierr

    muststop = .false.
	do jprocs = 0,nprocs-1
		ifrag = nodal(jprocs)%ifragment
		ffrag = nodal(jprocs)%ffragment
        if(opbundlestart(jprocs) < 1 .or. opbundleend(jprocs) < 1)cycle
		do ibundle = opbundlestart(jprocs),opbundleend(jprocs)
			if(opbundle(ibundle)%ifragment/=ifrag .or. opbundle(ibundle)%ffragment/=ffrag)then
				muststop = .true.
				if(iproc==0)then
					print*,' Mismatch of fragments on MPI process/node ',jprocs
					print*,' Variable NODAL expects initial, final ',ifrag,ffrag
					print*,' but assigned bundle ',ibundle,' expects ', opbundle(ibundle)%ifragment,opbundle(ibundle)%ffragment

				end if
			end if

		end do

	end do
	if(muststop)then
		call BMPI_ABORT(icomm,101,ierr)
		stop
	end if
	return

end subroutine checknodaltribution
!===============================================================
!
! alternate routine for assigning MPI processes
! based upon OPBUNDLESTART/END
!
subroutine setnodaltribution_alt

!    use fragments
    use nodeinfo
    use io
	use opbundles
    implicit none
    integer ifrag,ffrag
	integer :: bstart,bend
	integer :: jprocs,kprocs
    integer ibundle
	integer :: aerr

    allocate(nodal(0:nprocs-1), stat=aerr )
    if(aerr /= 0) call memerror("setnodaltribution_alt")

    if(nprocs==1)then  ! KSM - change from nproc.  Was clearing MPI data during modeling
       opfragstat(:,:)%nnodes = 1
    end if

    if(iproc==0)print*,' Alternate nodal distribution routine '
	do jprocs = 0,nprocs-1
		bstart = opbundlestart(jprocs)
		bend   = opbundleend(jprocs)

		if(bstart==0 .or. bend ==0)then
			nodal(jprocs)%ifragment = nfragments  ! assign to last fragment
			nodal(jprocs)%ffragment = nfragments  ! assign to last fragment
			cycle
		end if
!............ CHECK ALL ASSIGNED OPBUNDLES HAVE SAME FRAGMENTS .............

		ifrag = opbundle(bstart)%ifragment
		ffrag = opbundle(bend)%ffragment
		do ibundle = bstart,bend
			if(opbundle(ibundle)%ifragment /= ifrag .or. opbundle(ibundle)%ffragment/=ffrag)then
				if(iproc==0)then
					print*,' Bad choice of fragments on opbundles for process ',jprocs
					print*,' opbundle start/end = ',bstart,bend
					print*,' Expect initial/final fragments = ',ifrag,ffrag
					print*,' but bundle ',ibundle,' has ',opbundle(ibundle)%ifragment,opbundle(ibundle)%ffragment
				end if
				stop
			end if

		end do
		nodal(jprocs)%ifragment =ifrag
		nodal(jprocs)%ffragment =ffrag

!......... FIND 'FIRST' NODE WITH THESE FRAGMENTS ........
        nodal(jprocs)%ifirst = .true.

		do 	kprocs = 0,jprocs-1
			if(nodal(kprocs)%ifragment ==ifrag .and. nodal(kprocs)%ffragment==ffrag)then
		        nodal(jprocs)%ifirst = .false.
				exit
			end if
		end do

	end do
	return

end subroutine setnodaltribution_alt
!======================================================================
!
! Debugging routine to dump the vec1/vec2 bounds for a given rank
!
!  CALLED BY:
!     check_opbundle_bounds
!
subroutine print_node_bounds(fn, jnode)
   use precisions
   use nodeinfo
!   use fragments
   use basis
   use butil_mod
   implicit none
   integer, intent(in) :: fn  ! file number
   integer, intent(in) :: jnode
   integer :: frag1, frag2 ! fragments for vec1 and vec2
   integer(kind=basis_prec) :: f1s, f1e, f2s, f2e  ! fragment bounds

   frag1 = nodal(jnode)%ifragment
   frag2 = nodal(jnode)%ffragment
   ! get bounds for vec1 and vec2
   f1s = basestart(frag1)
   f1e = basestop(frag1)
   f2s = basestart(frag2)
   f2e = basestop(frag2)
   write(fn, "(A,I0,A,I0,A,I0,A,I0,A,I0,A)") "node=", jnode, &
      ", vec1 bounds=(", f1s, ":", f1e, "), vec2 bounds=(", f2s, ":", f2e, ")"
   flush(fn)
end subroutine print_node_bounds
!===============================================================

end module fragments

max_stars_repo_name: cwjsdsu/BigstickPublick

max_issues_repo_name: cwjsdsu/BigstickPublick

max_forks_repo_name: cwjsdsu/BigstickPublick


subroutine addCountxy(Count,NGrid,thisx,thisy,XX,YY,r_cutoff)
	implicit none
	integer, intent(in) :: NGrid
!f2p intent(in) ::  NGrid
	real(8), intent(inout), dimension(0:NGrid-1, 0:NGrid-1) :: Count, XX, YY
!f2p intent(in,out) :: Count, XX, YY
	real(8), intent(in) :: thisx, thisy, r_cutoff
!f2p intent(in) :: thisx, thisy, r_cutoff
	integer :: i,j,xndx(1),yndx(1)
	integer :: xlowndx, ylowndx, xhighndx, yhighndx
	real(8) :: dist, dx, dy
	real(8), dimension(0:NGrid-1) :: xvec, yvec

	dx = XX(0,1) - XX(0,0)  !assumes linspace used to generate spacing
	dy = YY(1,0) - YY(0,0)
	xvec = XX(0,:)
	yvec = YY(:,0)
	xndx = minloc(abs(xvec-thisx))
	yndx = minloc(abs(yvec-thisy))
	!write(*,*) xndx,yndx,thisx,thisy
	xlowndx = nint(xndx(1) - (r_cutoff / dx * 2.5))
	xhighndx = nint(xndx(1) + (r_cutoff / dx * 2.5))
	ylowndx = nint(yndx(1) - (r_cutoff / dy * 2.5))
	yhighndx = nint(yndx(1) + (r_cutoff / dy * 2.5))
	if (xlowndx < 0) xlowndx = 0
	if (ylowndx < 0) ylowndx = 0
	if (xhighndx > NGrid-1) xhighndx = NGrid - 1
	if (yhighndx > NGrid-1) yhighndx = Ngrid - 1
	do j = xlowndx, xhighndx
		do i = ylowndx, yhighndx
			dist = (thisx-XX(i,j))**2 + (thisy-YY(i,j))**2
			dist = sqrt(dist)
			if (dist .LT. r_cutoff) then
				Count(i,j) = Count(i,j) + 1
			end if
			!write(*,*) thisx,thisy, XX(i,j), YY(i,j), dist
		end do
	end do

end subroutine

subroutine addCountxz(Count,NGrid,thisx,thisz,XX,ZZ,r_cutoff)
	implicit none
	integer, intent(in) :: NGrid
!f2p intent(in) ::  NGrid
	real(8), intent(inout), dimension(0:NGrid-1, 0:NGrid-1) :: Count, XX, ZZ
!f2p intent(in,out) :: Count, XX, ZZ
	real(8), intent(in) :: thisx, thisz, r_cutoff
!f2p intent(in) :: thisx, thisz, r_cutoff
	integer :: i,j,xndx(1),zndx(1)
	integer :: xlowndx, zlowndx, xhighndx, zhighndx
	real(8) :: dist, dx, dz
	real(8), dimension(0:NGrid-1) :: xvec, zvec

	dx = XX(0,1) - XX(0,0)  !assumes linspace used to generate spacing
	dz = ZZ(1,0) - ZZ(0,0)
	xvec = XX(0,:)
	zvec = ZZ(:,0)
	xndx = minloc(abs(xvec-thisx))
	zndx = minloc(abs(zvec-thisz))
	!write(*,*) xndx,yndx,thisx,thisy
	xlowndx = nint(xndx(1) - (r_cutoff / dx * 2.5))
	xhighndx = nint(xndx(1) + (r_cutoff / dx * 2.5))
	zlowndx = nint(zndx(1) - (r_cutoff / dz * 2.5))
	zhighndx = nint(zndx(1) + (r_cutoff / dz * 2.5))
	if (xlowndx < 0) xlowndx = 0
	if (zlowndx < 0) zlowndx = 0
	if (xhighndx > NGrid-1) xhighndx = NGrid - 1
	if (zhighndx > NGrid-1) zhighndx = Ngrid - 1
	do j = xlowndx, xhighndx
		do i = zlowndx, zhighndx
			dist = (thisx-XX(i,j))**2 + (thisz-ZZ(i,j))**2
			dist = sqrt(dist)
			if (dist .LT. r_cutoff) then
				Count(i,j) = Count(i,j) + 1
			end if
			!write(*,*) thisx,thisy, XX(i,j), YY(i,j), dist
		end do
	end do

end subroutine

subroutine addCountyz(Count,NGrid,thisy,thisz,YY,ZZ,r_cutoff)
	implicit none
	integer, intent(in) :: NGrid
!f2p intent(in) ::  NGrid
	real(8), intent(inout), dimension(0:NGrid-1, 0:NGrid-1) :: Count, YY, ZZ
!f2p intent(in,out) :: Count, YY, ZZ
	real(8), intent(in) :: thisy, thisz, r_cutoff
!f2p intent(in) :: thisy, thisz, r_cutoff
	integer :: i,j,yndx(1),zndx(1)
	integer :: ylowndx, zlowndx, yhighndx, zhighndx
	real(8) :: dist, dy, dz
	real(8), dimension(0:NGrid-1) :: yvec, zvec

	dy = YY(0,1) - YY(0,0)  !assumes linspace used to generate spacing
	dz = ZZ(1,0) - ZZ(0,0)
	yvec = YY(0,:)
	zvec = ZZ(:,0)
	yndx = minloc(abs(yvec-thisy))
	zndx = minloc(abs(zvec-thisz))
	!write(*,*) yndx,zndx,thisy,thisz
	ylowndx = nint(yndx(1) - (r_cutoff / dy * 2.5))
	yhighndx = nint(yndx(1) + (r_cutoff / dy * 2.5))
	zlowndx = nint(zndx(1) - (r_cutoff / dz * 2.5))
	zhighndx = nint(zndx(1) + (r_cutoff / dz * 2.5))
	if (ylowndx < 0) ylowndx = 0
	if (zlowndx < 0) zlowndx = 0
	if (yhighndx > NGrid-1) yhighndx = NGrid - 1
	if (zhighndx > NGrid-1) zhighndx = Ngrid - 1
	do j = ylowndx, yhighndx
		do i = zlowndx, zhighndx
			dist = (thisy-YY(i,j))**2 + (thisz-ZZ(i,j))**2
			dist = sqrt(dist)
			if (dist .LT. r_cutoff) then
				Count(i,j) = Count(i,j) + 1
			end if
			!write(*,*) thisy,thisz, YY(i,j), ZZ(i,j), dist
		end do
	end do

end subroutine


text:
      SUBROUTINE PMAT(V,W)
C
      PARAMETER(NFAC=49)
      COMMON/OACOF/AA(NFAC,NFAC)
      DIMENSION V(1),W(1),C(NFAC)
C
      DO 10,I=1,NFAC
      DO 20,J=1,NFAC
   20 C(J)=AA(I,J)
C
      W(I)=SDOT(NFAC,C,1,V,1)
   10 CONTINUE
C
C
      RETURN
      END
       FUNCTION SDOT(N,X,INCX,Y,INCY)
       REAL*4 X(1),Y(1),SDOT
       SDOT=0.0
       IX=1
       IY=1
       DO 10 I=1,N
          SDOT=SDOT+X(IX)*Y(IY)
          IX=IX+INCX
          IY=IY+INCY
  10   CONTINUE
       RETURN
       END

ext: f


text:
!
!  Written by Leandro Martínez, 2009-2011.
!  Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez,
!  Ernesto G. Birgin.
!

!
! Arrays required by the flashsort package. Used only in heuristics, but
! defined here to be allocated dynamically
!

module flashsort

  use sizes
  implicit none
  integer, allocatable :: indflash(:) ! (ntotat)
  integer, allocatable :: lflash(:) ! (ntotat)
  integer :: mflash

end module flashsort

max_stars_repo_name: donshen/packmol

max_issues_repo_name: donshen/packmol

max_forks_repo_name: donshen/packmol


text:
c     comment section
c
c     fm037
c
c         this routine tests arithmetic assignment statements of the
c     form
c              integer variable = arithmetic expression
c     where the arithmetic expression is formed with three integer
c     constants and the arithmetic operator /.  both positive and nega-
c     tive constants are used in the arithmetic expression.
c
c         there are tests which require no truncation of the result
c     and tests where the result must be truncated before being stored
c     in the resultant integer variable.  the standard states 'the value
c     of an integer factor or term is the nearest integer whose magni-
c     tude does not exceed the magnitude of the mathematical value
c     represented by that factor or term.  the associative and commuta-
c     tive laws do not apply in the evaluation of integer terms con-
c     taining division, hence the evaluation of such terms must effec-
c     tively proceed from left to right.'
c
c         there are tests where the arithmetic expression contains
c             (1)  integer constant/integer constant/integer constant
c                      no truncation required
c             (2)  integer constant/integer constant/integer constant
c                      truncation required
c
c      references
c        american national standard programming language fortran,
c              x3.9-1978
c
c        section 4.3, integer type
c        section 4.3.1, integer constant
c        section 6.1, arithmetic expressions
c        section 6.6, evaluation of expressions
c        section 10.1, arithmetic assignment statement
c
c      **********************************************************
c
c         a compiler validation system for the fortran language
c     based on specifications as defined in american national standard
c     programming language fortran x3.9-1978, has been developed by the
c     federal cobol compiler testing service.  the fortran compiler
c     validation system (fcvs) consists of audit routines, their related
c     data, and an executive system.  each audit routine is a fortran
c     program, subprogram or function which includes tests of specific
c     language elements and supporting procedures indicating the result
c     of executing these tests.
c
c         this particular program/subprogram/function contains features
c     found only in the subset as defined in x3.9-1978.
c
c         suggestions and comments should be forwarded to -
c
c                  department of the navy
c                  federal cobol compiler testing service
c                  washington, d.c.  20376
c
c      **********************************************************
c
c
c
c     initialization section
c
c     initialize constants
c      **************
c     i01 contains the logical unit number for the card reader.
      i01 = 5
c     i02 contains the logical unit number for the printer.
      i02 = 6
c     system environment section
c
cx010    this card is replaced by contents of fexec x-010 control card.
c     the cx010 card is for overriding the program default i01 = 5
c     (unit number for card reader).
cx011    this card is replaced by contents of fexec x-011 control card.
c     the cx011 card is for systems which require additional
c     fortran statements for files associated with cx010 above.
c
cx020    this card is replaced by contents of fexec x-020 control card.
c     the cx020 card is for overriding the program default i02 = 6
c     (unit number for printer).
cx021    this card is replaced by contents of fexec x-021 control card.
c     the cx021 card is for systems which require additional
c     fortran statements for files associated with cx020 above.
c
      ivpass=0
      ivfail=0
      ivdele=0
      iczero=0
c
c     write page headers
      write (i02,90000)
      write (i02,90001)
      write (i02,90002)
      write (i02, 90002)
      write (i02,90003)
      write (i02,90002)
      write (i02,90004)
      write (i02,90002)
      write (i02,90011)
      write (i02,90002)
      write (i02,90002)
      write (i02,90005)
      write (i02,90006)
      write (i02,90002)
c
c     test section
c
c         arithmetic assignment statement
c
c     test 491 through test 519 contain three integer constants and
c     operator / in an arithmetic expression.  the form tested is
c         integer variable = integer constant/integer constant/int.con.
c
c
c     test 491 through test 496 - positive integer constants
c                       no truncation required
c
 4911 continue
      ivtnum = 491
c
c      ****  test 491  ****
c
      if (iczero) 34910, 4910, 34910
 4910 continue
      ivcomp = 24/3/4
      go to 44910
34910 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44910, 4921, 44910
44910 if (ivcomp - 2) 24910,14910,24910
14910 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4921
24910 ivfail = ivfail + 1
      ivcorr = 2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4921 continue
      ivtnum = 492
c
c      ****  test 492  ****
c
      if (iczero) 34920, 4920, 34920
 4920 continue
      ivcomp = 330/3/2
      go to 44920
34920 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44920, 4931, 44920
44920 if (ivcomp - 55) 24920,14920,24920
14920 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4931
24920 ivfail = ivfail + 1
      ivcorr = 55
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4931 continue
      ivtnum = 493
c
c      ****  test 493  ****
c
      if (iczero) 34930, 4930, 34930
 4930 continue
      ivcomp = 15249/13/51
      go to 44930
34930 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44930, 4941, 44930
44930 if (ivcomp - 23) 24930,14930,24930
14930 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4941
24930 ivfail = ivfail + 1
      ivcorr = 23
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4941 continue
      ivtnum = 494
c
c      ****  test 494  ****
c
      if (iczero) 34940, 4940, 34940
 4940 continue
      ivcomp = 7150/2/25
      go to 44940
34940 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44940, 4951, 44940
44940 if (ivcomp - 143) 24940,14940,24940
14940 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4951
24940 ivfail = ivfail + 1
      ivcorr = 143
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4951 continue
      ivtnum = 495
c
c      ****  test 495  ****
c
      if (iczero) 34950, 4950, 34950
 4950 continue
      ivcomp = 32766/2/3
      go to 44950
34950 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44950, 4961, 44950
44950 if (ivcomp - 5461) 24950,14950,24950
14950 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4961
24950 ivfail = ivfail + 1
      ivcorr = 5461
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4961 continue
      ivtnum = 496
c
c      ****  test 496  ****
c
      if (iczero) 34960, 4960, 34960
 4960 continue
      ivcomp = 32766/1/1
      go to 44960
34960 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44960, 4971, 44960
44960 if (ivcomp - 32766) 24960,14960,24960
14960 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4971
24960 ivfail = ivfail + 1
      ivcorr = 32766
      write (i02,80004) ivtnum, ivcomp ,ivcorr
c
c     test 497 through test 502 - positive integer constants
c                  truncation required
c
 4971 continue
      ivtnum = 497
c
c      ****  test 497  ****
c
      if (iczero) 34970, 4970, 34970
 4970 continue
      ivcomp = 24/3/3
      go to 44970
34970 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44970, 4981, 44970
44970 if (ivcomp -2) 24970,14970,24970
14970 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4981
24970 ivfail = ivfail + 1
      ivcorr = 2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4981 continue
      ivtnum = 498
c
c      ****  test 498  ****
c
      if (iczero) 34980, 4980, 34980
 4980 continue
      ivcomp = 230/2/3
      go to 44980
34980 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44980, 4991, 44980
44980 if (ivcomp - 38) 24980,14980,24980
14980 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 4991
24980 ivfail = ivfail + 1
      ivcorr = 38
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 4991 continue
      ivtnum = 499
c
c      ****  test 499  ****
c
      if (iczero) 34990, 4990, 34990
 4990 continue
      ivcomp = 7151/3/10
      go to 44990
34990 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 44990, 5001, 44990
44990 if (ivcomp - 238) 24990,14990,24990
14990 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5001
24990 ivfail = ivfail + 1
      ivcorr = 238
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5001 continue
      ivtnum = 500
c
c      ****  test 500  ****
c
      if (iczero) 35000, 5000, 35000
 5000 continue
      ivcomp = 15248/51/13
      go to 45000
35000 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45000, 5011, 45000
45000 if (ivcomp - 22) 25000,15000,25000
15000 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5011
25000 ivfail = ivfail + 1
      ivcorr = 22
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5011 continue
      ivtnum = 501
c
c      ****  test 501  ****
c
      if (iczero) 35010, 5010, 35010
 5010 continue
      ivcomp = 27342/4/3
      go to 45010
35010 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45010, 5021, 45010
45010 if (ivcomp - 2278) 25010,15010,25010
15010 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5021
25010 ivfail = ivfail + 1
      ivcorr = 2278
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5021 continue
      ivtnum = 502
c
c      ****  test 502  ****
c
      if (iczero) 35020, 5020, 35020
 5020 continue
      ivcomp = 32767/2/1
      go to 45020
35020 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45020, 5031, 45020
45020 if (ivcomp - 16383) 25020,15020,25020
15020 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5031
25020 ivfail = ivfail + 1
      ivcorr = 16383
      write (i02,80004) ivtnum, ivcomp ,ivcorr
c
c     test 503 through test 507 - negative integer constants included
c                  no truncation required
c
 5031 continue
      ivtnum = 503
c
c      ****  test 503  ****
c
      if (iczero) 35030, 5030, 35030
 5030 continue
      ivcomp = -24/3/4
      go to 45030
35030 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45030, 5041, 45030
45030 if (ivcomp +2) 25030,15030,25030
15030 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5041
25030 ivfail = ivfail + 1
      ivcorr = -2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5041 continue
      ivtnum = 504
c
c      ****  test 504  ****
c
      if (iczero) 35040, 5040, 35040
 5040 continue
      ivcomp = 330/(-3)/2
      go to 45040
35040 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45040, 5051, 45040
45040 if (ivcomp + 55) 25040,15040,25040
15040 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5051
25040 ivfail = ivfail + 1
      ivcorr = -55
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5051 continue
      ivtnum = 505
c
c      ****  test 505  ****
c
      if (iczero) 35050, 5050, 35050
 5050 continue
      ivcomp = 15249/(-13)/(-51)
      go to 45050
35050 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45050, 5061, 45050
45050 if (ivcomp - 23) 25050,15050,25050
15050 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5061
25050 ivfail = ivfail + 1
      ivcorr = 23
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5061 continue
      ivtnum = 506
c
c      ****  test 506  ****
c
      if (iczero) 35060, 5060, 35060
 5060 continue
      ivcomp = -7150/(-2)/(-25)
      go to 45060
35060 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45060, 5071, 45060
45060 if (ivcomp + 143) 25060,15060,25060
15060 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5071
25060 ivfail = ivfail + 1
      ivcorr = -143
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5071 continue
      ivtnum = 507
c
c      ****  test 507  ****
c
      if (iczero) 35070, 5070, 35070
 5070 continue
      ivcomp = (-32766)/(-2)/(-3)
      go to 45070
35070 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45070, 5081, 45070
45070 if (ivcomp + 5461) 25070,15070,25070
15070 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5081
25070 ivfail = ivfail + 1
      ivcorr = -5461
      write (i02,80004) ivtnum, ivcomp ,ivcorr
c
c     test 508 through test 513 - negative integer constants included
c                       truncation required
c
 5081 continue
      ivtnum = 508
c
c      ****  test 508  ****
c
      if (iczero) 35080, 5080, 35080
 5080 continue
      ivcomp = -24/3/3
      go to 45080
35080 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45080, 5091, 45080
45080 if (ivcomp + 2) 25080,15080,25080
15080 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5091
25080 ivfail = ivfail + 1
      ivcorr = -2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5091 continue
      ivtnum = 509
c
c      ****  test 509  ****
c
      if (iczero) 35090, 5090, 35090
 5090 continue
      ivcomp = 230/(-2)/3
      go to 45090
35090 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45090, 5101, 45090
45090 if (ivcomp + 38) 25090,15090,25090
15090 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5101
25090 ivfail = ivfail + 1
      ivcorr = -38
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5101 continue
      ivtnum = 510
c
c      ****  test 510  ****
c
      if (iczero) 35100, 5100, 35100
 5100 continue
      ivcomp = 7151/(-3)/(-10)
      go to 45100
35100 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45100, 5111, 45100
45100 if (ivcomp - 238) 25100,15100,25100
15100 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5111
25100 ivfail = ivfail + 1
      ivcorr = 238
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5111 continue
      ivtnum = 511
c
c      ****  test 511  ****
c
      if (iczero) 35110, 5110, 35110
 5110 continue
      ivcomp = -15248/(-51)/(-13)
      go to 45110
35110 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45110, 5121, 45110
45110 if (ivcomp + 22) 25110,15110,25110
15110 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5121
25110 ivfail = ivfail + 1
      ivcorr = -22
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5121 continue
      ivtnum = 512
c
c      ****  test 512  ****
c
      if (iczero) 35120, 5120, 35120
 5120 continue
      ivcomp = (-27342)/(-4)/(-3)
      go to 45120
35120 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45120, 5131, 45120
45120 if (ivcomp + 2278) 25120,15120,25120
15120 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5131
25120 ivfail = ivfail + 1
      ivcorr = -2278
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5131 continue
      ivtnum = 513
c
c      ****  test 513  ****
c
      if (iczero) 35130, 5130, 35130
 5130 continue
      ivcomp = 32767/2/(-1)
      go to 45130
35130 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45130, 5141, 45130
45130 if (ivcomp + 16383) 25130,15130,25130
15130 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5141
25130 ivfail = ivfail + 1
      ivcorr = -16383
      write (i02,80004) ivtnum, ivcomp ,ivcorr
c
c     test 514 through test 519 - positive and negative signed integer
c           constants in arithmetic expression.
c
 5141 continue
      ivtnum = 514
c
c      ****  test 514  ****
c
      if (iczero) 35140, 5140, 35140
 5140 continue
      ivcomp = +24/(-3)/4
      go to 45140
35140 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45140, 5151, 45140
45140 if (ivcomp +2) 25140,15140,25140
15140 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5151
25140 ivfail = ivfail + 1
      ivcorr = -2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5151 continue
      ivtnum = 515
c
c      ****  test 515  ****
c
      if (iczero) 35150, 5150, 35150
 5150 continue
      ivcomp = 24/(+3)/(-4)
      go to 45150
35150 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45150, 5161, 45150
45150 if (ivcomp +2) 25150,15150,25150
15150 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5161
25150 ivfail = ivfail + 1
      ivcorr = -2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5161 continue
      ivtnum = 516
c
c      ****  test 516  ****
c
      if (iczero) 35160, 5160, 35160
 5160 continue
      ivcomp = -24/(-3)/(+4)
      go to 45160
35160 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45160, 5171, 45160
45160 if (ivcomp -2) 25160,15160,25160
15160 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5171
25160 ivfail = ivfail + 1
      ivcorr = 2
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5171 continue
      ivtnum = 517
c
c      ****  test 517  ****
c
      if (iczero) 35170, 5170, 35170
 5170 continue
      ivcomp = -16811/(-16812)/(+1)
      go to 45170
35170 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45170, 5181, 45170
45170 if (ivcomp - 0) 25170,15170,25170
15170 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5181
25170 ivfail = ivfail + 1
      ivcorr = 0
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5181 continue
      ivtnum = 518
c
c      ****  test 518  ****
c
      if (iczero) 35180, 5180, 35180
 5180 continue
      ivcomp = (-16811) / (+16811) / (+1)
      go to 45180
35180 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45180, 5191, 45180
45180 if (ivcomp +1) 25180,15180,25180
15180 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5191
25180 ivfail = ivfail + 1
      ivcorr = -1
      write (i02,80004) ivtnum, ivcomp ,ivcorr
 5191 continue
      ivtnum = 519
c
c      ****  test 519  ****
c
      if (iczero) 35190, 5190, 35190
 5190 continue
      ivcomp = (-335)/(+168)/(+1)
      go to 45190
35190 ivdele = ivdele + 1
      write (i02,80003) ivtnum
      if (iczero) 45190, 5201, 45190
45190 if (ivcomp + 1) 25190,15190,25190
15190 ivpass = ivpass + 1
      write (i02,80001) ivtnum
      go to 5201
25190 ivfail = ivfail + 1
      ivcorr = -1
      write (i02,80004) ivtnum, ivcomp ,ivcorr
c      ****    end of tests    ****
 5201 continue
c
c     write page footings and run summaries
99999 continue
      write (i02,90002)
      write (i02,90006)
      write (i02,90002)
      write (i02,90002)
      write (i02,90007)
      write (i02,90002)
      write (i02,90008)  ivfail
      write (i02,90009) ivpass
      write (i02,90010) ivdele
c
c
c     terminate routine execution
      stop
c
c     format statements for page headers
90000 format (1h1)
90002 format (1h )
90001 format (1h ,10x,34hfortran compiler validation system)
90003 format (1h ,21x,11hversion 1.0)
90004 format (1h ,10x,38hfor official use only - copyright 1978)
90005 format (1h ,5x,4htest,5x,9hpass/fail, 5x,8hcomputed,8x,7hcorrect)
90006 format (1h ,5x,46h----------------------------------------------)
90011 format (1h ,18x,17hsubset level test)
c
c     format statements for run summaries
90008 format (1h ,15x,i5,19h errors encountered)
90009 format (1h ,15x,i5,13h tests passed)
90010 format (1h ,15x,i5,14h tests deleted)
c
c     format statements for test results
80001 format (1h ,4x,i5,7x,4hpass)
80002 format (1h ,4x,i5,7x,4hfail)
80003 format (1h ,4x,i5,7x,7hdeleted)
80004 format (1h ,4x,i5,7x,4hfail,10x,i6,9x,i6)
80005 format (1h ,4x,i5,7x,4hfail,4x,e12.5,3x,e12.5)
c
90007 format (1h ,20x,20hend of program fm037)
      end

max_stars_repo_name: weiss/original-bsd

max_issues_repo_name: JamesLinus/original-bsd

max_forks_repo_name: JamesLinus/original-bsd


text:
module hello_interface
  character(len=13) :: string="Hello, world!"
  interface
    module function get() result(result_string)
      character(:), allocatable :: result_string
    end function
  end interface
end module

submodule(hello_interface) hello_implementation
contains
  module function get() result(result_string)
    character(:), allocatable :: result_string
    result_string = string
  end function
end submodule

  use hello_interface
  print *,len(get())
end

max_stars_repo_name: rouson/bug-reports

max_issues_repo_name: rouson/bug-reports

max_forks_repo_name: rouson/bug-reports


subroutine wind_vertical_interpolation(m,n,o,p,req,pres,&
	pressure_full_level,wind_on_model_level,wind_on_press_level)
implicit none
integer :: m,n,o,p,req,o_sub
integer :: x,y,s,t,plev
double precision :: pres(req),wind_on_model_level(m,n,o,p)
double precision :: pressure_full_level(m,n,o,p)
double precision :: delta_w,delta_p,grad_w_p,diff_p
double precision, intent(out) :: wind_on_press_level(m,n,req,p)
real :: arg = -1.0,NaN
NaN= sqrt(arg)
o_sub = o - 1
do plev=1,req
 do t=1,p
  do  x=1,m
   do  y=1,n
    do  s=1,o_sub
!above uppest level
     if(pres(plev) .LT. pressure_full_level(x,y,1,t)) then
      wind_on_press_level(x,y,plev,t) = NaN
     end if
! in between levels
     if(pres(plev) .GT. pressure_full_level(x,y,s,t) .AND. pres(plev) .LT. &
       pressure_full_level(x,y,s+1,t)) then
      delta_w = wind_on_model_level(x,y,s,t) - wind_on_model_level(x,y,s+1,t)
      delta_p = log(pressure_full_level(x,y,s,t)) &
       - log(pressure_full_level(x,y,s+1,t))
      grad_w_p = delta_w /delta_p
      diff_p = log(pres(plev)) - log(pressure_full_level(x,y,s,t))
      wind_on_press_level(x,y,plev,t) = wind_on_model_level(x,y,s,t) &
        +  grad_w_p * diff_p
     end if
! extrapolation below the ground
     if(pres(plev) .GT. pressure_full_level(x,y,o,t)) then
      wind_on_press_level(x,y,plev,t) = wind_on_model_level(x,y,o,t)
     end if
    end do
   end do
  end do
 end do
end do
end subroutine wind_vertical_interpolation


text:
c
c     file tcud34sp.f
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1999 by UCAR                 .
c  .                                                             .
c  .       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                      MUDPACK version 5.0                    .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... author and specialist
c
c          John C. Adams (National Center for Atmospheric Research)
c          email: johnad@ucar.edu, phone: 303-497-1213

c ... For MUDPACK 5.0 information, visit the website:
c     (http://www.scd.ucar.edu/css/software/mudpack)
c
c ... purpose
c
c     test program for the MUDPACK solver cud34sp
c
c ... required MUDPACK files
c
c     cud34sp.f, cud3sp.f, cudcom.f
c
c
c *********************************************************
c *********************************************************
c
c     sample program/test driver for cud34sp
c
c **********************************************************
************************************************************
c
c     a sample program/test driver for cud34sp is below. it can be
c     executed as an initial test.  the output is listed for the
c     test case described.
c
c     test cud34sp below by solving the complex separable elliptic pde
c
c        d(cmplx(x,-x)*dp/dx)/dx + d(cmplx(y,-y)*dp/dy)/dy +
c
c        d(cmplx(z,-z)*dp/dz)/dz + cmplx(1.,-1.)*(dp/dx + dp/dy + dp/dz) +
c
c        cmplx(-(x+y+z),x+y+z)*p(x,y,z) = r(x,y,z)
c
c     on the (x,y,z) region
c
c       1/4 < x < 3/4, 1/3 < y < 2/3, 1/5 < z < 4/5
c
c     with dirchlet boundary conditions at the upper x,y,z boundaries
c     and the mixed complex derirvative boundary conditions:
c
c          (1) dp/dx + cmplx(1.,1.)*p(0.5,y,z) = gbdxa(y,z) at x = 1/4
c
c          (2) dp/dy + cmplx(-1.,-1.)*p(x,0.5,z) = gbdyc(x,z) at y = 1/3
c
c          (3) dp/dz + cmplx(1.,-1.)*p(x,y,0.5) = gbdze(x,y) at z = 1/5
c
c     note the complex number multiplying "p" must be constant along
c     each surface.
c
c     use the exact solution
c
c          pe(x,y,z) = cmplx(exp(x+y),exp(y+z))
c
c     for testing and choose a 49 by 33 by 41 grid.  First cud3sp
c     is called with 5 cycles to ensure second-order discretization
c     level error is reached (a requirement for mudpack fourth-order
c     solvers).  Then cud34sp is called for a fourth-order estimate.
c
c ******************************************************
c     output (64 bit floating point arithmetic)
c *******************************************************
c
c     cud3sp test
c
c     input arguments
c     intl =  0 nxa =  2 nxb =  1 nyc =  2 nyd =  1
c     nze =  2 nzf =  1
c     ixp =  3 jyq =  2 kzr =  3
c     iex =  5 jey =  5 kez =  5
c     nx =  49 ny =  33 nz =  49 iguess =  0 maxcy =  5
c     method =  0 work space length input =  291605
c     xa =  0.25 xb =  0.75
c     yc =  0.33 yd =  0.67
c     ze =  0.20 zf =  0.80
c     tolmax =  0.000E+00
c
c     multigrid options
c     kcycle =  2
c     iprer =  2
c     ipost =  1
c     intpol =  3
c
c     discretization call to cud3sp intl =  0
c     ierror =  0 minimum work space =  291605
c
c     approximation call to cud3sp
c     intl =  1 method =  0 iguess =  0 maxcy =  5
c     ierror =  0
c     maximum error  =   0.287E-04
c
c     cud34sp test  ierror =  0
c     maximum error  =   0.592E-07
c
c ************************************************************
c     end of output
c ************************************************************
c
c
      program tcud34sp
      implicit none
c
c     set grid sizes with parameter statements
c
      integer iixp,jjyq,kkzr,iiex,jjey,kkez,llwork,nnx,nny,nnz
      parameter(iixp=3,jjyq=2,kkzr=3)
      parameter(iiex=5,jjey=5,kkez=5)
      parameter (nnx = iixp*2**(iiex-1)+1)
      parameter (nny = jjyq*2**(jjey-1)+1)
      parameter (nnz = kkzr*2**(kkez-1)+1)
c
c     set minimal work space (see tcud3sp.f)
c
      parameter (llwork = 291605)
c
c     dimension solution,right hand side, and work arrays
c
      complex phi(nnx,nny,nnz),rhs(nnx,nny,nnz),work(llwork)
      integer iprm(23),mgopt(4)
      real fprm(8)
c
c     put integer and floating point arguments names in contiguous
c     storeage for labelling purposes
c
      integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,
     +              kez,nx,ny,nz,iguess,maxcy,method,nwork,
     +              lwrkqd,itero
      common/itcd3sp/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,
     +              kez,nx,ny,nz,iguess,maxcy,method,nwork,
     +              lwrkqd,itero
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/ftcud3sp/xa,xb,yc,yd,ze,zf,tolmax,relmax
      real dlx,dly,dlz,x,y,z,errm
      complex cxx,cyy,czz,cx,cy,cz,ce,cex,cey,cez
      complex pxx,pyy,pzz,px,py,pz,pe
      integer i,j,k,ierror
      equivalence(intl,iprm)
      equivalence(xa,fprm)
c
c     declare coefficient and boundary condition input subroutines external
c
      external cfx,cfy,cfz,bndc
c
c     set for initial call
c
      intl = 0
c
c     set boundary condition flags
c
      nxa = 2
      nxb = 1
      nyc = 2
      nyd = 1
      nze = 2
      nzf = 1
c
c     set grid sizes from parameter statements
c
      ixp = iixp
      jyq = jjyq
      kzr = kkzr
      iex = iiex
      jey = jjey
      kez = kkez
      nx = nnx
      ny = nny
      nz = nnz
c
c     set  5 multigrid cycles
c
      maxcy = 5
c
c     set work space length approximation from parameter statement
c
      nwork = llwork
c
c     set point relaxation (only choice with cud3sp)
c
      method = 0
c
c     set full multigrid cycling beginning at coarsest grid
c
      iguess = 0
c
c     set end points of solution region in (x,y,z) space
c
      xa = 0.25
      xb = 0.75
      yc = 1.0/3.0
      yd = 2.0/3.0
      ze = 0.20
      zf = 0.80
c
c     set default multigrid options
c
      mgopt(1) = 2
      mgopt(2) = 2
      mgopt(3) = 1
      mgopt(4) = 3
c
c     set mesh increments
c
      dlx = (xb-xa)/float(nx-1)
      dly = (yd-yc)/float(ny-1)
      dlz = (zf-ze)/float(nz-1)
c
c     set for no error control
c
      tolmax = 0.0
c
c     set right hand side in rhs and phi to zero
c
      do k=1,nz
	z = ze+(k-1)*dlz
	call cfz(z,czz,cz,cez)
	do j=1,ny
	  y = yc+float(j-1)*dly
	  call cfy(y,cyy,cy,cey)
	  do i=1,nx
	    x = xa+float(i-1)*dlx
	    call cfx(x,cxx,cx,cex)
	    call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	    ce = cex+cey+cez
	    rhs(i,j,k) = cxx*pxx+cyy*pyy+czz*pzz+cx*px+cy*py+cz*pz+ce*pe
	    phi(i,j,k) = (0.,0.)
	  end do
	end do
      end do
c
c     set specified values at upper x,y,z boundaries in phi
c
      x = xb
      do k=1,nz
	z = ze+(k-1)*dlz
	do j=1,ny
	  y = yc+float(j-1)*dly
	  call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	  phi(nx,j,k) = pe
	end do
      end do
      y = yd
      do k=1,nz
	z = ze+(k-1)*dlz
	do i=1,nx
	  x = xa+float(i-1)*dlx
	  call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	  phi(i,ny,k) = pe
	end do
      end do
      z = zf
      do j=1,ny
	y = yc+(j-1)*dly
	do i=1,nx
	  x = xa+float(i-1)*dlx
	  call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	  phi(i,j,nz) = pe
	end do
      end do
      write(*,50)
   50 format(//' cud3sp test ')
c
c     print input arguments
c
      write(*,100)intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez,
     +            nx,ny,nz,iguess,maxcy,method,nwork,xa,xb,yc,yd,ze,zf,
     +            tolmax,mgopt(1),mgopt(2),mgopt(3),mgopt(4)
  100 format(/' input arguments ',
     +/'intl = ',i2,' nxa = ',i2,' nxb = ',i2,' nyc = ',i2,' nyd = ',i2,
     +/' nze = ',i2, ' nzf = ',i2,
     +/' ixp = ',i2,' jyq = ',i2,' kzr = ',i2,
     +/' iex = ',i2, ' jey = ',i2, ' kez = ',i2,
     +/' nx = ',i3,' ny = ',i3,' nz = ',i3, ' iguess = ',i2,
     +/' maxcy = ',i2,
     +/' method = ',i2, ' work space length input = ',i7,
     +/' xa = ',f5.2,' xb = ',f5.2,
     +/' yc = ',f5.2,' yd = ',f5.2,
     +/' ze = ',f5.2,' zf = ',f5.2,
     +/' tolmax = ' ,e10.3
     +//' multigrid options '
     +/' kcycle = ',i2
     +/' iprer = ',i2
     +/' ipost = ',i2
     +/' intpol = ',i2 )
c
c     discretize pde only
c
      write(*,104) intl
  104 format(/' discretization call to cud3sp', ' intl = ', i2)
      call cud3sp(iprm,fprm,work,cfx,cfy,cfz,bndc,rhs,phi,mgopt,ierror)
      write (*,105) ierror,iprm(21)
  105 format(' ierror = ',i2, ' minimum work space = ',i7)
      if (ierror.gt.0) call exit(0)
c
c     approximate pde
c
      intl = 1
      write(*,106) intl,method,iguess,maxcy
  106 format(/' approximation call to cud3sp ',
     +/' intl = ',i2, ' method = ',i2,' iguess = ',i2, ' maxcy = ',i2)
      call cud3sp(iprm,fprm,work,cfx,cfy,cfz,bndc,rhs,phi,mgopt,ierror)
      write (*,107) ierror
  107 format(' ierror = ',i2)
      if (ierror.gt.0) call exit(0)
c
c     compute and print maximum error
c
      call error(nx,ny,nz,phi,errm)
      write(*,108) errm
  108 format(' maximum error  =  ',e10.3)
c
c     attempt fourth-order estimate
c
      call cud34sp(work,phi,ierror)
      write(*,109) ierror
  109 format(/ ' cud34sp test ', ' ierror = ',i2)
      if (ierror.gt.0) call exit(0)
      call error(nx,ny,nz,phi,errm)
      write(*,108) errm
      end

      subroutine error(nx,ny,nz,phi,errm)
c
c     compute the error in the estimate in phi
c
      implicit none
      integer nx,ny,nz
      complex phi(nx,ny,nz)
      real xa,xb,yc,yd,ze,zf,tolmax,relmax,errm
      common/ftcud3sp/xa,xb,yc,yd,ze,zf,tolmax,relmax
      real dlx,dly,dlz,x,y,z
      complex pxx,pyy,pzz,px,py,pz,pe
      integer i,j,k
      dlx = (xb-xa)/(nx-1)
      dly = (yd-yc)/(ny-1)
      dlz = (zf-ze)/(nz-1)
      errm = 0.0
      do k=1,nz
	z = ze+(k-1)*dlz
	do j=1,ny
	  y = yc+float(j-1)*dly
	  do i=1,nx
	    x = xa+float(i-1)*dlx
	    call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	    errm = amax1(errm,abs(phi(i,j,k)-pe))
	  end do
	end do
      end do
      return
      end
c
c     complex coefficient subroutines
c
      subroutine cfx(x,cxx,cx,cex)
      implicit none
      real x
      complex cxx,cx,cex
      cxx = cmplx(x,-x)
      cx = (1.,-1.)
      cex = cmplx(-x,x)
      return
      end
      subroutine cfy(y,cyy,cy,cey)
      implicit none
      real y
      complex cyy,cy,cey
      cyy= cmplx(y,-y)
      cy = (1.,-1.)
      cey = cmplx(-y,y)
      return
      end
      subroutine cfz(z,czz,cz,cez)
      implicit none
      real z
      complex czz,cz,cez
      czz = cmplx(z,-z)
      cz = (1.,-1.)
      cez = cmplx(-z,z)
      return
      end

      subroutine bndc(kbdy,xory,yorz,cons,gbdy)
c
c     input complex derivative boundary conditions to cud3sp
c
      implicit none
      integer kbdy
      real xory,yorz,x,y,z
      complex cons,gbdy,pe,px,py,pz,pxx,pyy,pzz
      real xa,xb,yc,yd,ze,zf,tolmax,relmax
      common/ftcud3sp/xa,xb,yc,yd,ze,zf,tolmax,relmax
      if (kbdy.eq.1) then
c
c     x=xa surface
c
	x = xa
	y = xory
	z = yorz
	call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	cons = (1.,1.)
	gbdy = px+cons*pe
	return
      end if
      if (kbdy.eq.3) then
c
c     y=yc surface
c
	y = yc
	x = xory
	z = yorz
	call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	cons = (-1.,-1.)
	gbdy = py+cons*pe
	return
      end if
      if (kbdy.eq.5) then
c
c     z=ze surface
c
	z = ze
	x = xory
	y = yorz
	call exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
	cons = (1.,-1.)
	gbdy = pz + cons*pe
	return
      end if
      end

      subroutine exact(x,y,z,pxx,pyy,pzz,px,py,pz,pe)
c
c     this subroutine is used to set an exact solution for testing cud3sp
c     (i.e., setting the rhs, boundary conditions and computing the exact
c     error)
c
      implicit none
      real x,y,z,exy,eyz
      complex pe,px,py,pz,pxx,pyy,pzz
      exy = exp(x+y)
      eyz = exp(y+z)
      pe = cmplx(exy,eyz)
      px = cmplx(exy,0.)
      pxx = px
      py = pe
      pyy = pe
      pz = cmplx(0.,eyz)
      pzz = pz
      return
      end

max_stars_repo_name: utastudents/selalib

max_issues_repo_name: utastudents/selalib

max_forks_repo_name: utastudents/selalib


text:
: testkeys
	repeat
		key? dup .
		if ekey 27 - 0; drop then
		100 ms
	again ;

." bug 315 : key? holds its value also after releasing the key" cr
." This example will print a '0' if no key is pending; otherwise it will" cr
." print a '-1'.  Press ESC to quit" cr

testkeys cr bye

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
Measure C is a local ballot measure, the only one on the ballot in the March 2012 Election.  The text of the measure reads, Shall the Davis Joint Unified School District continue a local parcel tax to preserve existing classroom programs including math, science, English, and music; key schoolbased personnel including librarians, nurses, counselors, school technology and reading specialists; and other programs such as athletics, drama, student nutrition and class size reduction for a period of 5 years not to exceed the Base Annual Tax of $150.00 per unit for multidwelling parcels and $320.00 per parcel for all other parcels?

The Yes for Our Students campaign, which supports Measure C, hosted an Election Night party on March 6, 2012 beginning at 7 p.m. at the home of Kim Eichorn, 631 Rutgers Drive. All supporters of Measure C were welcome.

The No on Measure C campaign met for its party beginning at 7 p.m. Tuesday upstairs at Woodstocks Pizza, 219 G Street downtown.

Mar. 06, 2012, 11:08 PM.  Turnout: 17,210 of 43,138 Registered Voters (39.9%). Returns available at: http://www.yoloelections.org/returns/

See the The Vanguards http://davisvanguard.org/index.php?optioncom_content&viewarticle&id5155:measurecsurgestoresoundingwinatpolls&Itemid121 Measure C Surges to Resounding Win At Polls.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search