text:
!***************************************************************************************************
!>Subroutine initialMeshUniform():
!The subroutine creates meshes and sectors, and establishes the mapping relationship among meshes, sectors and processes.
!Each process is divided into eight sectors.
!NOTE: Each sector should have at least two meshes per dimension (x, y, z dimension).
!***************************************************************************************************
subroutine initialMesh()
	use mod_constants
	use mod_globalVariables
	use mod_structures

	implicit none
	include 'mpif.h'

	integer :: status(MPI_STATUS_SIZE), i, j, k, dir, localElem
	double precision :: length
	integer :: x, y, z, numXmin,numXmax,numYmin,numYmax,numZmin,numZmax	!used to determine numxLocal, numyLocal, numzLocal
	double precision :: tempCenter(6)

	length=meshLength
	numTotal=numx*numy*numz	!total cell in the system
	systemVol = dble(numTotal)*length**3d0

	!The  boundary coordinates of the system (xmin,xmax,ymin,ymax,zmin,zmax)
	do i=1,5,2
		myProc%globalCoord(i)=0d0
	end do
	myProc%globalCoord(2)=length*dble(numx)
	myProc%globalCoord(4)=length*dble(numy)
	myProc%globalCoord(6)=length*dble(numz)

	!The  boundary coordinates of the processor (xmin,xmax,ymin,ymax,zmin,zmax)
	myProc%localCoord(1)=dble(myProc%coords(1))*myProc%globalCoord(2)/dble(dims(1))    !<xmin
	myProc%localCoord(2)=myProc%localCoord(1)+myProc%globalCoord(2)/dble(dims(1))
	myProc%localCoord(3)=dble(myProc%coords(2))*myProc%globalCoord(4)/dble(dims(2))
	myProc%localCoord(4)=myProc%localCoord(3)+myProc%globalCoord(4)/dble(dims(2))
	myProc%localCoord(5)=dble(myProc%coords(3))*myProc%globalCoord(6)/dble(dims(3))
	myProc%localCoord(6)=myProc%localCoord(5)+myProc%globalCoord(6)/dble(dims(3))

	totalVolume=(myProc%localCoord(2)-myProc%localCoord(1))*(myProc%localCoord(4)-myProc%localCoord(3))*&
			(myProc%localCoord(6)-myProc%localCoord(5))

	!numxLocal, numyLocal, numzLocal are used to determine the size of the mesh inside the local processor.
	numxLocal=0
	numyLocal=0
	numzLocal=0

	!get numxLocal
	if(myProc%localCoord(1)==myProc%globalCoord(1) .AND. myProc%localCoord(2)==myProc%globalCoord(2)) then
		numXmin=0
		numXmax=numx
		numxLocal=numx
		tempCenter(1)=length/2d0
		tempCenter(2)=myProc%globalCoord(2)-length/2d0

	else if(myProc%localCoord(1)==myProc%globalCoord(1)) then	!at xmin

		numXmin=0
		numXmax = myProc%localCoord(2)/(length/2d0)
		if((dble(numXmax)*(length/2d0)) <=  myProc%localCoord(2) .AND. myProc%localCoord(2)< (dble(numXmax+1)*(length/2d0))) then
			if(mod(numXmax,2)==0) then
				numXmax = numXmax/2
			else
				numXmax = (numXmax+1)/2
			end if
		end if

		numxLocal=numXmax
		tempCenter(1)=length/2d0
		tempCenter(2)=dble(numXmax)*length-length/2d0

	else if(myProc%localCoord(2)==myProc%globalCoord(2)) then	!at xmax
		numXmin = myProc%localCoord(1)/(length/2d0)
		numXmax=numx
		if((dble(numXmin)*(length/2d0)) <=  myProc%localCoord(1) .AND. myProc%localCoord(1) < (dble(numXmin+1)*(length/2d0))) then
			if(mod(numXmin,2)==0) then
				numXmin = numXmin/2
			else
				numXmin = (numXmin+1)/2
			end if
		end if

		numxLocal=numx-numXmin
		tempCenter(1)=dble(numXmin)*length+length/2d0
		tempCenter(2)=myProc%globalCoord(2)-length/2d0

	else	!in the middle
		numXmin = myProc%localCoord(1)/(length/2d0)
		numXmax = myProc%localCoord(2)/(length/2d0)
		if((dble(numXmin)*(length/2d0)) <=  myProc%localCoord(1) .AND. myProc%localCoord(1) < &
				(dble(numXmin+1)*(length/2d0))) then
			if(mod(numXmin,2)==0) then
				numXmin = numXmin/2
			else
				numXmin = (numXmin+1)/2
			end if
		end if
		if((dble(numXmax)*(length/2d0)) <=  myProc%localCoord(2) .AND. myProc%localCoord(2) < &
				(dble(numXmax+1)*(length/2d0))) then
			if(mod(numXmax,2)==0) then
				numXmax = numXmax/2
			else
				numXmax = (numXmax+1)/2
			end if
		end if

		numxLocal = numXmax-numXmin
		tempCenter(1)=dble(numXmin)*length+length/2d0
		tempCenter(2)=dble(numXmax)*length-length/2d0

	end if

	!get numyLocal
	if(myProc%localCoord(3)==myProc%globalCoord(3) .AND. myProc%localCoord(4)==myProc%globalCoord(4)) then
		numYmin=0
		numYmax=numy
		numyLocal=numy
		tempCenter(3)=length/2d0
		tempCenter(4)=myProc%globalCoord(4)-length/2d0

	else if(myProc%localCoord(3)==myProc%globalCoord(3)) then	!at ymin
		numYmin=0
		numYmax = myProc%localCoord(4)/(length/2d0)
		if((dble(numYmax)*(length/2d0)) <=  myProc%localCoord(4) .AND. myProc%localCoord(4) < (dble(numYmax+1)*(length/2d0))) then
			if(mod(numYmax,2)==0) then
				numYmax = numYmax/2
			else
				numYmax = (numYmax+1)/2
			end if
		end if
		numyLocal=numYmax
		tempCenter(3)=length/2d0
		tempCenter(4)=dble(numYmax)*length-length/2d0

	else if(myProc%localCoord(4)==myProc%globalCoord(4)) then	!at ymax
		numYmin = myProc%localCoord(3)/(length/2d0)
		numYmax=numy
		if((dble(numYmin)*(length/2d0)) <=  myProc%localCoord(3) .AND. myProc%localCoord(3) < (dble(numYmin+1)*(length/2d0))) then
			if(mod(numYmin,2)==0) then
				numYmin = numYmin/2
			else
				numYmin = (numYmin+1)/2
			end if
		end if
		numyLocal=numy-numYmin
		tempCenter(3)=dble(numYmin)*length+length/2d0
		tempCenter(4)=myProc%globalCoord(4)-length/2d0

	else	!in the middle
		numYmin = myProc%localCoord(3)/(length/2d0)
		numYmax = myProc%localCoord(4)/(length/2d0)
		if((dble(numYmin)*(length/2d0)) <=  myProc%localCoord(3) .AND. myProc%localCoord(3) < &
				(dble(numYmin+1)*(length/2d0))) then
			if(mod(numYmin,2)==0) then
				numYmin = numYmin/2
			else
				numYmin = (numYmin+1)/2
			end if
		end if
		if((dble(numYmax)*(length/2d0)) <=  myProc%localCoord(4) .AND. myProc%localCoord(4) < &
				(dble(numYmax+1)*(length/2d0))) then
			if(mod(numYmax,2)==0) then
				numYmax = numYmax/2
			else
				numYmax = (numYmax+1)/2
			end if
		end if
		numyLocal = numYmax-numYmin
		tempCenter(3)=dble(numYmin)*length+length/2d0
		tempCenter(4)=dble(numYmax)*length-length/2d0
	end if

	!get numzLocal
	if(myProc%localCoord(5)==myProc%globalCoord(5) .AND. myProc%localCoord(6)==myProc%globalCoord(6)) then
		numZmin=0
		numZmax=numz
		numzLocal=numz
		tempCenter(5)=length/2d0
		tempCenter(6)=myProc%globalCoord(6)-length/2d0

	else if(myProc%localCoord(5)==myProc%globalCoord(5)) then	!at zmin
		numZmin=0
		numZmax = myProc%localCoord(6)/(length/2d0)
		if((dble(numZmax)*(length/2d0)) <=  myProc%localCoord(6) .AND. myProc%localCoord(6) < (dble(numZmax+1)*(length/2d0))) then
			if(mod(numZmax,2)==0) then
				numZmax = numZmax/2
			else
				numZmax = (numZmax+1)/2
			end if
		end if
		numzLocal=numZmax
		tempCenter(5)=length/2d0
		tempCenter(6)=dble(numZmax)*length-length/2d0

	else if(myProc%localCoord(6)==myProc%globalCoord(6)) then	!at zmax
		numZmin = myProc%localCoord(5)/(length/2d0)
		numZmax=numz
		if((dble(numZmin)*(length/2d0)) <=  myProc%localCoord(5) .AND. myProc%localCoord(5)< (dble(numZmin+1)*(length/2d0))) then
			if(mod(numZmin,2)==0) then
				numZmin = numZmin/2
			else
				numZmin = (numZmin+1)/2
			end if
		end if
		numzLocal=numz-numZmin
		tempCenter(5)=dble(numZmin)*length+length/2d0
		tempCenter(6)=myProc%globalCoord(6)-length/2d0

	else	!in the middle
		numZmin = myProc%localCoord(5)/(length/2d0)
		numZmax = myProc%localCoord(6)/(length/2d0)
		if((dble(numZmin)*(length/2d0)) <=  myProc%localCoord(5) .AND. myProc%localCoord(5) < &
				(dble(numZmin+1)*(length/2d0))) then
			if(mod(numZmin,2)==0) then
				numZmin = numZmin/2
			else
				numZmin = (numZmin+1)/2
			end if
		end if
		if((dble(numZmax)*(length/2d0)) <=  myProc%localCoord(6) .AND. myProc%localCoord(6) < &
				(dble(numZmax+1)*(length/2d0))) then
			if(mod(numZmax,2)==0) then
				numZmax = numZmax/2
			else
				numZmax = (numZmax+1)/2
			end if
		end if
		numzLocal = numZmax-numZmin
		tempCenter(5)=dble(numZmin)*length+length/2d0
		tempCenter(6)=dble(numZmax)*length-length/2d0
	end if

	!this variable is used at various points in SRSCD; lets us know the length fo myMesh(:)
	numCells=numxLocal*numyLocal*numzLocal	!total cells of this processor
	if(myProc%taskid==MASTER) then
		write(*,*) 'numx numy numz', numx, numy, numz
	end if
	!write(*,*) 'proc', myProc%taskid, 'numxLocal numyLocal numzLocal', numxLocal, numyLocal, numzLocal

	!If any processors don't have volume elements in them (too many procs), we create an error message
	if(numCells==0) then
		write(*,*) 'error processors with no volume elements'
		call MPI_ABORT(comm,ierr)
	end if

	!Step 5: Create meshes. And for each volume element in myMesh, assign coordinates and material number (and proc number)
	allocate(myMesh(numCells))
	localElem=0

	do k=1,numzLocal
		do j=1,numyLocal
			do i=1,numxLocal

				localElem=localElem+1

				x = numXmin+i
				y = numYmin+j
				z = numZmin+k

				myMesh(localElem)%coordinates(1)=tempCenter(1)+(i-1)*length
				myMesh(localElem)%coordinates(2)=tempCenter(3)+(j-1)*length
				myMesh(localElem)%coordinates(3)=tempCenter(5)+(k-1)*length
				myMesh(localElem)%material=1
				myMesh(localElem)%proc=myProc%taskid
				myMesh(localElem)%globalCell=(z-1)*numx*numy+(y-1)*numx+x

				myMesh(localElem)%length=length
				myMesh(localElem)%volume=length**3d0

				myMesh(localElem)%numNeighbors=1
				myMesh(localElem)%neighbors=0
				myMesh(localElem)%neighborProcs=-1
				!allocate(myMesh(localElem)%neighbors(6))
				!allocate(myMesh(localElem)%neighborProcs(6))
				!do dir=1,6
				!	myMesh(localElem)%numNeighbors(dir)=1
				!end do
			end do
		end do
	end do

	!Step 3d: assign neighbors and processor numbers for neighbors (connectivity in myMesh) - periodic or free surfaces
	call createMeshConnect(length)

end subroutine

!**************************************************************************************************
!>Subroutine createMeshConnect(length) (periodic boundary conditions or free surface)
!It identifies the mesh and processor of neighboring meshes for each mesh.
!The connectivity scheme is the same as in the global case, but neighboring processor numbers are used here.
!**************************************************************************************************
subroutine createMeshConnect(length)
	use mod_structures
	use mod_globalVariables
	implicit none
	include 'mpif.h'

	double precision, intent(in) :: length
	integer :: cell, localCell, maxElement, maxBndryCells
	integer :: i, dir, tag
	integer, allocatable :: send(:,:,:), recv(:,:)
	integer :: materialBuff(numCells,6)
	integer :: numSend(6)=0, numRecv(6)=0
	integer :: status(MPI_STATUS_SIZE)
	integer :: tempRecv
	double precision :: commTime1, commTime2

	maxBndryCells = numxLocal*numyLocal
	if((numyLocal*numzLocal) > maxBndryCells) then
		maxBndryCells = numyLocal*numzLocal
	end if
	if((numxLocal*numzLocal) > maxBndryCells) then
		maxBndryCells = numxLocal*numzLocal
	end if

	!allocate(send(2, numxLocal*numyLocal*numzLocal,6))
	allocate(send(2, maxBndryCells,6))

	do cell=1,numCells
		!Right (+x)
		if(mod(cell,numxLocal)==0) then !identify cell to the right
			myMesh(cell)%neighborProcs(1)=myProc%procNeighbor(1)
			if(myProc%procNeighbor(1)==-1) then	!+x free surface
				myMesh(cell)%neighbors(1)=0
			else if(myMesh(cell)%neighborProcs(1)==myProc%taskid) then	!+x periodic
				myMesh(cell)%neighbors(1)=cell-numxLocal+1
			else
				numSend(1)=numSend(1)+1			!count the number of items in buffer
				send(1,numSend(1),1)=cell		!localCellID in this processor
				send(2,numSend(1),1)=myMesh(cell)%material	!Material ID number that this element is composed of
			end if
		else
			myMesh(cell)%neighbors(1)=cell+1
			myMesh(cell)%neighborProcs(1)=myProc%taskid
		end if

		!Left (-x)
		if(mod(cell+numxLocal-1,numxLocal)==0) then !identify cell to the left
			myMesh(cell)%neighborProcs(2)=myProc%procNeighbor(2)
			if(myProc%procNeighbor(2)==-1) then	!-x free surface
				myMesh(cell)%neighbors(2)=0
			else if(myMesh(cell)%neighborProcs(2)==myProc%taskid) then
				myMesh(cell)%neighbors(2)=cell+numxLocal-1
			else
				numSend(2)=numSend(2)+1
				send(1,numSend(2),2)=cell
				send(2,numSend(2),2)=myMesh(cell)%material
			end if
		else
			myMesh(cell)%neighbors(2)=cell-1
			myMesh(cell)%neighborProcs(2)=myProc%taskid
		end if

		!Front (+y)
		if(mod(cell,numxLocal*numyLocal) > numxLocal*(numyLocal-1) .OR. mod(cell,numxLocal*numyLocal)==0) then
			myMesh(cell)%neighborProcs(3)=myProc%procNeighbor(3)
			if(myProc%procNeighbor(3)==-1) then
				myMesh(cell)%neighbors(3)=0
			else if(myMesh(cell)%neighborProcs(3)==myProc%taskid) then
				myMesh(cell)%neighbors(3)=cell-(numxLocal*(numyLocal-1))
			else
				numSend(3)=numSend(3)+1
				send(1,numSend(3),3)=cell
				send(2,numSend(3),3)=myMesh(cell)%material
			end if
		else
			myMesh(cell)%neighbors(3)=cell+numxLocal
			myMesh(cell)%neighborProcs(3)=myProc%taskid
		end if

		!Back (-y)
		if(mod(cell,numxLocal*numyLocal) <= numxLocal .AND. (mod(cell, numxLocal*numyLocal) /= 0 &
				.OR. numyLocal==1)) then
			myMesh(cell)%neighborProcs(4)=myProc%procNeighbor(4)
			if(myProc%procNeighbor(4)==-1) then
				myMesh(cell)%neighbors(4)=0
			else if(myMesh(cell)%neighborProcs(4)==myProc%taskid) then
				myMesh(cell)%neighbors(4)=cell+(numxLocal*(numyLocal-1))
			else
				numSend(4)=numSend(4)+1
				send(1,numSend(4),4)=cell
				send(2,numSend(4),4)=myMesh(cell)%material
			end if

		else
			myMesh(cell)%neighbors(4)=cell-numxLocal
			myMesh(cell)%neighborProcs(4)=myProc%taskid
		end if

		!Up (+z)
		if(mod(cell,numxLocal*numyLocal*numzLocal) > numxLocal*numyLocal*(numzLocal-1) &
				.OR. mod(cell, numxLocal*numyLocal*numzLocal)==0) then
			myMesh(cell)%neighborProcs(5)=myProc%procNeighbor(5)
			if(myProc%procNeighbor(5)==-1) then
				myMesh(cell)%neighbors(5)=0
			else if(myMesh(cell)%neighborProcs(5)==myProc%taskid) then
				myMesh(cell)%neighbors(5)=cell-(numxLocal*numyLocal*(numzLocal-1))
			else
				numSend(5)=numSend(5)+1
				send(1,numSend(5),5)=cell
				send(2,numSend(5),5)=myMesh(cell)%material
			end if

		else
			myMesh(cell)%neighbors(5)=cell+numxLocal*numyLocal
			myMesh(cell)%neighborProcs(5)=myProc%taskid
		end if

		!Down (-z)
		if(mod(cell,numxLocal*numyLocal*numzLocal) <= numxLocal*numyLocal &
				.AND. (mod(cell,numxLocal*numyLocal*numzLocal) /= 0 .OR. numzLocal==1)) then
			myMesh(cell)%neighborProcs(6)=myProc%procNeighbor(6)
			if(myProc%procNeighbor(6)==-1) then
				myMesh(cell)%neighbors(6)=0
			else if(myMesh(cell)%neighborProcs(6)==myProc%taskid) then
				myMesh(cell)%neighbors(6)=cell+(numxLocal*numyLocal*(numzLocal-1))
			else
				numSend(6)=numSend(6)+1
				send(1,numSend(6),6)=cell
				send(2,numSend(6),6)=myMesh(cell)%material
			end if

		else
			myMesh(cell)%neighbors(6)=cell-numxLocal*numyLocal
			myMesh(cell)%neighborProcs(6)=myProc%taskid
		end if
	end do

	maxElement=0
	!*******************************************************
	do dir = 1,6
		if(dir == 1 .OR. dir == 2) then			!<x
			numRecv(dir) = numyLocal*numzLocal
		else if(dir == 3 .OR. dir == 4) then	!<y
			numRecv(dir) = numxLocal*numzLocal
		else if(dir == 5 .OR. dir == 6) then	!<z
			numRecv(dir) = numxLocal*numyLocal
		end if
	end do

	do dir = 1,6

		if(mod(dir,2)==0) then
			tag = dir-1
		else
			tag = dir+1
		end if

		if(myProc%procNeighbor(dir) /= myProc%taskid .AND. myProc%procNeighbor(dir) /= -1) then
			commTime1=MPI_WTIME()
			call MPI_SEND(send(1,1,dir),numSend(dir)*2,MPI_INTEGER,myProc%procNeighbor(dir),dir,comm,ierr)
			commTime2=MPI_WTIME()
			commTimeSum=commTime2-commTime1
		end if

		if(myProc%procNeighbor(tag)/=myProc%taskid .AND. myProc%procNeighbor(tag) /= -1) then
			allocate(recv(2,numRecv(tag)))
			commTime1=MPI_WTIME()
			call MPI_RECV(recv,numRecv(tag)*2,MPI_INTEGER,myProc%procNeighbor(tag),dir,comm,status,ierr)
			commTime2=MPI_WTIME()
			commTimeSum=commTime2-commTime1

			do i=1, numRecv(tag)
				localCell=send(1,i,tag)
				myMesh(localCell)%neighbors(tag)=recv(1,i)
				materialBuff(localCell,tag)=recv(2,i)
			end do

			!if(myProc%taskid == 0) then
			!	do i = 1,numRecv(tag)
			!		write(*,*) 'recvDir', tag, 'numRecv', numRecv(tag), 'recvCell', recv(1,i)
			!	end do
			!end if

			deallocate(recv)
		end if
	end do
	deallocate(send)

	!***************************************************************************************************
	!Initializing myGhost with elements that are in neighboring processors that bound this one
	!***************************************************************************************************
	maxElement=0
	do i=1, numCells
		do dir=1, 6
			if(myMesh(i)%neighborProcs(dir) /= myProc%taskid) then	!we are pointing to a different proc
				if(myMesh(i)%neighbors(dir) > maxElement) then		!searching for the max element number in a neighbor
					maxElement=myMesh(i)%neighbors(dir)
				end if
			end if
		end do
	end do

	allocate(myGhost(maxElement,6))	!6 directions, maxElement elements in each direction (more than needed)
	!initialize myGhost with 0 in localNeighbor - signal that myGhost is not attached to anything
	do dir=1,6
		do i=1,maxElement
			myGhost(i,dir)%localNeighbor=0	!default, says that this is not a real element of myGhost.
			myGhost(i,dir)%proc=-10			!default, says that this is not a real element of myGhost.
			myGhost(i,dir)%material=0
			myGhost(i,dir)%length=0d0
			myGhost(i,dir)%volume=0d0
		end do
	end do

	do i=1,numCells
		do dir=1,6
			if(myMesh(i)%neighborProcs(dir) == -1) then										!this is a free surface
				!do nothing
			else if(myMesh(i)%neighborProcs(dir) /= myProc%taskid) then
				myGhost(myMesh(i)%neighbors(dir),dir)%proc=myMesh(i)%neighborProcs(dir)	!set proc # of elements in myGhost
				myGhost(myMesh(i)%neighbors(dir),dir)%length=length						!set length of elements in myGhost
				myGhost(myMesh(i)%neighbors(dir),dir)%volume=length**3d0					!set volume of elements in myGhost (changes with cascade addition)
				myGhost(myMesh(i)%neighbors(dir),dir)%material=materialBuff(i,dir)
				myGhost(myMesh(i)%neighbors(dir),dir)%localNeighbor=i
			end if
		end do
	end do

end subroutine

!***************************************************************************************
!>Function findgNeighborGID(globalID, dir)
!Inputs: globalID, direction
!Output: globalNeighborID of the globalID in the direaction
!***************************************************************************************
integer function findgNeighborGID(globalID, dir)
	use mod_structures
	use mod_globalVariables
	implicit none

	integer, intent(in) :: globalID, dir
	integer :: neighborID

	if(dir==1) then	!+x
		if(mod(globalID,numx)==0) then !identify cell to the right
			if(periods(1) .eqv. .true.) then	!periodic
				neighborID=globalID-numx+1
			else	!free surface
				neighborID=0
			end if
		else
			neighborID=globalID+1
		end if
	else if(dir==2) then
		if(mod(globalID+numx-1,numx)==0) then !identify cell to the left
			if(periods(1) .eqv. .true.) then
				neighborID=globalID+numx-1
			else
				neighborID=0
			end if
		else
			neighborID=globalID-1
		end if
	else if(dir==3) then
		if(mod(globalID,numx*numy) > numx*(numy-1) .OR. mod(globalID,numx*numy)==0) then	!identify cell to the front
			if(periods(2) .eqv. .true.) then
				neighborID=globalID-(numx*(numy-1))
			else
				neighborID=0
			end if
		else
			neighborID=globalID+numx
		end if
	else if(dir==4) then
		if(mod(globalID,numx*numy) <= numx .AND. (mod(globalID, numx*numy) /= 0 .OR. numy==1)) then	!identify cell to the back
			if(periods(2) .eqv. .true.) then
				neighborID=globalID+(numx*(numy-1))
			else
				neighborID=0
			end if
		else
			neighborID=globalID-numx
		end if
	else if(dir==5) then
		if(mod(globalID,numx*numy*numz) > numx*numy*(numz-1) .OR. mod(globalID, numx*numy*numz)==0) then	!identify cell to up
			if(periods(3) .eqv. .true.) then
				neighborID=globalID-(numx*numy*(numz-1))
			else
				neighborID=0
			end if
		else
			neighborID=globalID+numx*numy
		end if
	else if(dir==6) then
		if(mod(globalID,numx*numy*numz) <= numx*numy .AND. (mod(globalID,numx*numy*numz) /= 0 .OR. numz==1)) then	!identify cell to down
			if(periods(3) .eqv. .true.) then
				neighborID=globalID+(numx*numy*(numz-1))
			else
				neighborID=0
			end if
		else
			neighborID=globalID-numx*numy
		end if
	end if

	findgNeighborGID=neighborID

end function

max_stars_repo_name: MISA-SCD/misa-scd

max_issues_repo_name: MISA-SCD/misa-scd

max_forks_repo_name: MISA-SCD/misa-scd


text:
| Scene2 with some more utilities
| Reva OpenGL demo program
| Andrew Price, 2006

context: ~app

needs callbacks
needs util/auxstack
needs math/floats
~floats
needs ui/gl  needs ui/glu  needs ui/glut
~gl
~app
variable tmpf32
2variable tmpf64
: f>32 tmpf32 f!4  tmpf32 @ ;
: f>64 tmpf64 f!8  tmpf64 2@ ;
: 3f f>32 f>32 f>32 rot rot swap ;     		| do 3 f>32s

: 3floats  create f>32 f>32 f>32 , , , ; 	| Initialised vector of 3 32bit floats
: 4floats  f>32 3floats , ;  		 	| Initialised vector of 4 32bit floats
: -f1 f1 fnegate ;				| -1.0 on TOFS

0 variable, noargs
: blankargs "  " ;
: windowtitle " GL Test -- Scene" ;
blankargs drop variable, blankargs*
500 variable, width
500 variable, height
f0 f0 f0 f1  4floats light_ambient drop
f1 f1 f1 f1  4floats light_diffuse drop
f1 f1 f1 f1  4floats light_specular drop
f1 f1 f1 f0  4floats light_position drop

variable mywindow
: init	noargs blankargs* glutInit drop
	GLUT_RGB GLUT_DOUBLE GLUT_DEPTH or or  glutInitDisplayMode drop
	width @ height @ glutInitWindowSize drop
	100 100 glutInitWindowPosition drop
	windowtitle drop glutCreateWindow mywindow !
	GL_LIGHT0 GL_AMBIENT light_ambient glLightfv
	GL_LIGHT0 GL_DIFFUSE light_diffuse glLightfv
	GL_LIGHT0 GL_SPECULAR light_specular glLightfv
	GL_LIGHT0 GL_POSITION light_position glLightfv
	GL_LIGHTING glEnable
	GL_LIGHT0 glEnable
	GL_DEPTH_TEST glEnable
;

fvariable posx fvariable posy fvariable posz
-f1 0.75e0 f* posx f!  0.5e0 posy f!  f0 posz f!
variable key

: HandleKey ( key x y -- )
| 3dup rot ." Key: " dup $ff and . ." (" emit ." ) x: " swap . ." y: " .
rot $ff and key ! 2drop  key @
|	 ." Key " rot dup dup key !  emit 32 emit . . . cr
	key @
	case
		'w of posy f@ 0.05e0 f+ posy f! endof | w
		's of posy f@ 0.05e0 f- posy f! endof | s
		'd of posx f@ 0.05e0 f+ posx f! endof | d
		'a of posx f@ 0.05e0 f- posx f! endof | a
		'q of posz f@ 0.05e0 f+ posz f! endof | q
		'e of posz f@ 0.05e0 f- posz f! endof | e
		27 of bye endof
	endcase
	glutPostRedisplay drop ;	| ( key a b -- )

: HandleMouse  ." Mouse " . . . . cr ;

: DrawScene  | ." Draw" cr
	GL_COLOR_BUFFER_BIT GL_DEPTH_BUFFER_BIT or glClear
	GL_PROJECTION glMatrixMode
	glLoadIdentity
	-f1 2.5e0 f* f>64 2.5e0 f>64 -f1 2.5e0 f* f>64 2.5e0 f>64 -f1 10.0e0 f* f>64 10.0e0 f>64 glOrtho
	GL_MODELVIEW glMatrixMode
	glLoadIdentity
	glPushMatrix
	 20.0e0 f1 f0	3f f0 f>32 glRotatef
	 glPushMatrix
	  posx f@ posy f@ posz f@ 		3f glTranslatef
	  90.0e0 f>32 	f1 f0 f0	3f glRotatef
	  0.275e0 f>64 0.85e0 f>64 30 30 glutSolidTorus
	 glPopMatrix
	 glPushMatrix
	  -f1 0.75e0 f* f>32 -f1 0.5e0 f* f>32 0.0e0 f>32 glTranslatef
	  270.0e0 f>32 f1 f>32 f0 f>32 f0 f>32 glRotatef
	  1.0e0 f>64 2.0e0 f>64 15 15 glutSolidCone
	 glPopMatrix
	 glPushMatrix
	  0.75e0 f>32 0.0e0 f>32 -f1 f>32 glTranslatef
|	  1.0e0 f>64 30 30 glutSolidSphere
	  1.0e0 f>64 glutSolidTeapot
	 glPopMatrix
	glPopMatrix
	glFlush
	glutSwapBuffers
;

: Reshape ( width height -- )
   height ! width ! 					| ( width height -- )
	0 0 width @ height @ glViewport ." Reshape " width @ . height @ . cr
	GL_PROJECTION glMatrixMode
	glLoadIdentity
	-f1 2.5e0 f* f>64 2.5e0 f>64 -f1 2.5e0 f* f>64 2.5e0 f>64 -f1 10.0e0 f* f>64 10.0e0 f>64 glOrtho
	GL_MODELVIEW glMatrixMode
	glLoadIdentity

;

: SceneDrawer callback DrawScene ;
: KeyHandler callback 0 cb-param 1 cb-param 2 cb-param HandleKey ;
: MouseHandler callback  HandleMouse ;
: ReshapeHandler callback 0 cb-param 1 cb-param Reshape ;
: setup-callbacks 	['] SceneDrawer glutDisplayFunc drop 	['] KeyHandler glutKeyboardFunc drop
			['] MouseHandler glutMouseFunc drop 	['] ReshapeHandler glutReshapeFunc drop ;

init setup-callbacks ." press ESC to quit" cr glutMainLoop drop

ext: f

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
*dk,nbrsmooth
c -----------------------------------------------------------------

	subroutine nbrsmooth(cmo,mpary,mpno,rlxwt,ntimes,nwttry,useisn
     &  ,extnbr)

C -----------------------------------------------------------------
C
C  PURPOSE
C Smooth by moving nodes towards average position of neighbors.
C Limit definition of neighbor to those points connected by an edge
C  of an element to a node with a 'subset' of contraints/type
C  on the node of those on the neighbor and edge.
C Optionally check that the materials on the node are also a subset
C  of those on the neighbor and connecting edge.
C
C constraints:
c (1) above neighbor constraints
C (2) resultant element volumes must be positive
C
C caveats:
C (1) limited to tets for now
C (2) repeated iterations flatten curved surfaces, and shrink spheres
C (3) constaints do not know about analytic shapes
C
C
C  INPUT
C
C  cmo    name of mesh object
C  mpary  array of nodes numbers to smooth
C  mpno   number of nodes to smooth
C  rlxwt  weight for under-relaxed Laplacian smoothing
C         [recommended value: 0.5]
C  ntimes number of smoothing iterations
C         [recommended value: 1-5]
C  nwttry number of attempts to not tangle mesh by halving the
C         smoothing weights
C         [recommended value: 1-5]
C  useisn flag  0  ignore parent/child chains
C                  (interface node positions smoothed
C                   based on all interface neighbors)
C               1  use parent/child relationship
C                  (interface node positions smoothed
C                   based on only those neighbors
C                   who are connected along multimaterial edges
C                   with all the materials of the node)
C         [recommended value: 1; unlisted values same as 1]
C  extnbr flag  0  neighbors must be in mpary array
C               1  neighbors can be not in mpary array
C         [recommended value: 1; unlisted values same as 1]
C
C Rule for under-relaxed Laplacian smoothing is
C     x -> (1-rlxwt)*x + rlxwt*[ av_neighbor_position ]
C If a volume is inverted, I modify this to try nwttry times to
C reduce the smoothing weight of the nodes defining this volume by
C     rlxwt -> local_rlxwt*rlxwt,
C where local_rlxwt=rlxwt unless c rlxwt>1 or rlxwt<0,
C in which case local_rlxwt is set to 0.5 to avoid iterations blowing up.
C If locally reducing the weight fails to keep volumes positive
C after ntimes tries, I set the smoothing weight to zero on nodes
C defining inverted tets, and if this fails the smoothing is aborted.
C
C      CHANGE HISTORY -
C
C $Log: nbrsmooth.f,v $
C Revision 2.00  2007/11/05 19:46:02  spchu
C Import to CVS
C
CPVCS
CPVCS       Rev 1.5   14 Dec 2000 14:56:22   jtg
CPVCS    fixed error occuring when no neighbors
CPVCS
CPVCS       Rev 1.4   26 Jan 2000 16:57:24   jtg
CPVCS    now accepts mbndry=0
CPVCS
CPVCS       Rev 1.3   Mon Dec 21 18:51:22 1998   jtg
CPVCS    dfloat to changed to dble
CPVCS
CPVCS       Rev 1.2   Fri Oct 23 13:08:30 1998   dcg
CPVCS    replace lenptr with '*' - DEC compiler complaint
CPVCS
CPVCS       Rev 1.1   Wed Oct 21 01:09:52 1998   jtg
CPVCS    forgot to add PVCS log - corrected
c -----------------------------------------------------------------
c23456789012345678901234567890123456789012345678901234567890123456789012

      implicit none

      include 'local_element.h'

      include 'chydro.h'
      integer mpno
      integer mpary(mpno)

      integer nwttry,useisn,extnbr
      real*8 rlxwt

      pointer (ip_itet,itet),(ip_imt1,imt1),(ip_itp1,itp1)
     &  ,(ip_isn1,isn1),(ip_icr1,icr1),(ip_icontab,icontab)
     &  ,(ip_iparent,iparent),(ip_jtet,jtet),(ip_itettyp,itettyp)
      integer itet(4,*),imt1(*),itp1(*)
     &  ,isn1(*),icr1(*),icontab(50,*)
     &  ,iparent(*),jtet(4,*),itettyp(*)
c
      pointer (ip_xic,xic),(ip_yic,yic),(ip_zic,zic)
     &  ,(ip_dx,dx),(ip_dy,dy),(ip_dz,dz),(ip_weight,weight)
      real*8 xic(*),yic(*),zic(*)
     &  ,dx(*),dy(*),dz(*),weight(*)
c
      pointer (ip_nbr1,nbr1),(ip_nbor,nbor),(ip_link,link)
     &  ,(ip_neigh,neigh),(ip_neigh_off,neigh_off)
     &  ,(ip_neigh_n,neigh_n)
      integer nbr1(*),nbor(*),link(*)
     &  ,neigh(*),neigh_off(*)
     &  ,neigh_n(*)
c
      integer ntimes,itime,ii
     &  ,nnodes,nelements,nconbnd,mbndry
     &  ,length,itype,ierr,pos_vol,nwt,n_nbrs,length_add
     &  ,i,j,k,iel,i1(2),i2(2)
     &  ,ip1,ip2,ic1,ic2,jt1,jt2,is1,is2
c
      pointer (ip_child_nbr1,child_nbr1),(ip_child_nbor,child_nbor)
     &  ,(ip_child_link,child_link)
      integer child_nbr1(*),child_nbor(*),child_link(*)
      integer child_n_nbrs,child_length,n_dud
C
      real*8 xxx,rrr,x1(4),y1(4),z1(4),vvv,local_rlxwt
      character*32 cmo,isubname
      character*132 logmess
      real*8 zero,one,half,one_sixth
c
      zero=0.0d0
      one=1.0d0
      one_sixth=1.d0/6.0d0
      half=0.5d0

c -----------------------------------------------------------------
c -----------------------------------------------------------------

      isubname='nbrsmooth'

      call cmo_get_info('nnodes',cmo,nnodes,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info nnodes'
      call cmo_get_info('nelements',cmo,nelements,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info nelements'
      call cmo_get_info('itp1',cmo,ip_itp1,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info itp1'
      call cmo_get_info('imt1',cmo,ip_imt1,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info imt1'
      call cmo_get_info('isn1',cmo,ip_isn1,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info isn1'
      call cmo_get_info('icr1',cmo,ip_icr1,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info icr1'
      call cmo_get_info('itet',cmo,ip_itet,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info itet'
      call cmo_get_info('jtet',cmo,ip_jtet,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info jtet'
      call cmo_get_info('itettyp',cmo,ip_itettyp,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info itettyp'
      call cmo_get_info('xic',cmo,ip_xic,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info xic'
      call cmo_get_info('yic',cmo,ip_yic,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info yic'
      call cmo_get_info('zic',cmo,ip_zic,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info zic'
      call cmo_get_info('mbndry',cmo,mbndry,length,itype,ierr)
      if (ierr.ne.0) mbndry=0
      call cmo_get_info('nconbnd',cmo,nconbnd,length,itype,ierr)
      if (ierr.ne.0) then
	 nconbnd=0
      elseif (nconbnd.gt.0) then
	 call cmo_get_info('icontab',cmo,ip_icontab,length,itype,ierr)
	 if (ierr.ne.0) stop 'NBRSMOOTH: cmo_get_info icontab'
      endif
      itype=1
      length=nnodes
c   (use neigh_n to hold temp info...)
      call mmggetbk('neigh_n',isubname,ip_neigh_n,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk neigh_n'
      itype=2
      length=nnodes
      call mmggetbk('dx',isubname,ip_dx,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk dx'
      call mmggetbk('dy',isubname,ip_dy,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk dy'
      call mmggetbk('dz',isubname,ip_dz,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk dz'
      itype=2
      length=nnodes
      call mmggetbk('weight',isubname,ip_weight,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk weight'
      itype=1
      length=nnodes
      call mmggetbk("iparent",isubname,ip_iparent,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk iparent'
      call unpackpc(nnodes,itp1,isn1,iparent)
c
      itype=1
      length=nnodes
      call mmggetbk('nbr1',isubname,ip_nbr1,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk nbr1'
      if (useisn.ne.0) then
         call mmggetbk('child_nbr1',isubname,ip_child_nbr1
     &    ,length,itype,ierr)
         if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk child_nbr1'
      endif
      length=24*mpno+1000
      call mmggetbk('link',isubname,ip_link,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk link'
      call mmggetbk('nbor',isubname,ip_nbor,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk nbor'
      if (useisn.ne.0) then
         child_length=length
         call mmggetbk('child_link',isubname,ip_child_link
     &		,child_length,itype,ierr)
         if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk child_link'
         call mmggetbk('child_nbor',isubname,ip_child_nbor
     &		,child_length,itype,ierr)
         if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk child_nbor'
      else
          child_length=100
c   bigger than 24 so don't try to incr length..
      endif

c--------------------------------------------------------------------
c	 ........... (create neighbor table) .............
c
c initialize arrays for neighbor table
c nbr1 is "masked"
      do i=1,nnodes
	 nbr1(i)=-1
      enddo
      n_nbrs=0
      if (useisn.ne.0) then
         do i=1,nnodes
            child_nbr1(i)=0
         enddo
      endif
      child_n_nbrs=0
c
c unmask nbr1 for parents of points in mpary
      do ii=1,mpno
         i=iparent(mpary(ii))
         nbr1(i)=0
      enddo
c
      length_add=mpno+1000
c
c loop through elements to fill in neighbor table
      do iel=1,nelements
         if(itettyp(iel).ne.ifelmtet) then
            write(logmess,'(a)')
     *     'laplace smoothing only implemented on 3D tet meshes'
            call writloga('default',0,logmess,0,ierr)
            go to 9999
          endif
c
c  increment arrays if necessary
         if (n_nbrs+24.gt.length) then
            length=length+length_add
            call mmggetbk('link',isubname,ip_link,length,itype,ierr)
            if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk link'
            call mmggetbk('nbor',isubname,ip_nbor,length,itype,ierr)
            if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk nbor'
         endif
         if (child_n_nbrs+24.gt.child_length) then
            child_length=child_length+length_add
            call mmggetbk('child_link',isubname,ip_child_link
     &			,child_length,itype,ierr)
            if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk child_link'
            call mmggetbk('child_nbor',isubname,ip_child_nbor
     &			,child_length,itype,ierr)
            if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk child_nbor'
         endif
c
C  loop through edges of element iel
c   define endpoints (i1) and connected tets (i2)
c   at some point this should be generalized from "elmedge/elmface"...
c   as this is the only part of the code which is specific to tets
c   except dimension for itet,jtet assumed (4,nelements)
c   rather than using itetoff,etc
c   and formula for volume in volume check
         do k=1,6
            if (k.ge.1.and.k.le.4) then
               i1(1)=k
               i1(2)=1+mod(k,4)
               i2(1)=1+mod(k+1,4)
               i2(2)=1+mod(k+2,4)
            elseif (k.eq.5) then
               i1(1)=1
               i1(2)=3
               i2(1)=2
               i2(2)=4
c  if (k.eq.6) then
            else
               i1(1)=2
               i1(2)=4
               i2(1)=1
               i2(2)=3
            endif
            do i=1,2
               i1(i)=itet(i1(i),iel)
               i2(i)=jtet(i2(i),iel)
            enddo

c   loop through endpoints of edge k of element iel
            do i=1,2
c
c    check if OK re mask
c     (only unmasked nodes have neighbors)
               ip1=iparent(i1(i))
               ip2=iparent(i1(3-i))
               if (nbr1(ip1).lt.0) goto 250 ! not OK
               if (nbr1(ip2).lt.0 .and. extnbr.eq.0) goto 250 ! not OK
c
c    check if OK re itp
c     (point type of (child) node must be "lower" than neighbor)
c     Note: should switch these tests to proper "ifitp..." from chydro.h
               ip1=itp1(i1(i))
               ip2=itp1(i1(3-i))
               if (ip1.eq.0 .or. ip1.eq.ip2) then
                  continue
               elseif (ip1.eq.2
     &			.and.(ip2.eq.12.or.ip2.eq.13.or.ip2.eq.15)) then
                  continue
               elseif ((ip1.eq.10.or.ip1.eq.11)
     &			.and.(ip2.eq.12.or.ip2.eq.13)) then
                  continue
               elseif ((ip1.eq.12.or.ip1.eq.13).and.ip2.eq.15) then
                  continue
               elseif (ip1.eq.14.and.ip2.eq.15) then
                  continue
               else
c  not OK
                  goto 250
               endif

c    check if OK re jtet
c     (interface nodes only neighbors if edge is interface edge)
               if (ip1.ne.0) then ! ip1 is still point type
                  jt1=i2(1)
                  jt2=i2(2)
                  if ((jt1.lt.mbndry.and.jt2.lt.mbndry.and.mbndry.gt.0)
     &               .or.(jt1.gt.0.and.jt2.gt.0.and.mbndry.eq.0) ) then
                     goto 250
                  endif
               endif

c    check if OK re icr values
c     (constraints of node must be subset of constraints of neighbor)
c     Note: similar to material check, do I need to check that the
c     "bond constraints" also same or higher order? Or is the jtet check
c     that it's an interface bond enough? Answer is I probably should
c     check, but I don't know what the appropriate bond or face-based
c     variable to check the "bond constraint" is. -> DCG??
c     Probably in most cases the "useisn" check for "bond materials"
c     would insure the "bond constraints" are also consistent.
c     (Fails only if have interface bond connecting points on same
c      constraint surface with the bond not in the constraint surface,
c      and with the useisn flag on, the interface the bond lies on would
c      have to have all the same materials as the constraint surface.
c      A virtual interface on a curved surface between different
c      boundary conditions in the same material is an example requiring
c      one to check if the "bond constraints" are consistent
c      with the node constraints to properly assign neighbors and
c      avoid a "2nd neighbor" connection between points on the virtual
c      interface through either constraint surface.)
               ic1=icr1(i1(i))
               ic2=icr1(i1(3-i))
               if (ic1.eq.0 .or. ic1.eq.ic2) then
                  continue
               elseif (nconbnd.gt.0.and.ic1.gt.0.and.ic2.gt.0) then
                  do is1=1,icontab(1,ic1)
                     j=icontab(2+is1,ic1)
                     do is2=1,icontab(1,ic2)
                        if (icontab(2+is2,ic2).eq.j) goto 210
                     enddo
                     goto 250
210                  continue
                  enddo
               else
                  goto 250
               endif

c    if useisn flag on, check if OK re isn
c     (only neighbors if node materials subset of neighbor materials;
c     bond material check performed afterwards, using "child neighbors")
               if (useisn.ne.0) then
                  ip1=iparent(i1(i))
                  ip2=iparent(i1(3-i))
                  if (isn1(ip1).eq.0) then
                     continue
                  elseif (isn1(ip2).eq.0) then
                     goto 250
                  else
                     is1=isn1(ip1)
                     do while (is1.ne.ip1)
                        is2=isn1(ip2)
                        do while (is2.ne.ip2)
                           if (imt1(is1).eq.imt1(is2)) goto 220
                           is2=isn1(is2)
                        enddo
                        goto 250
220                     is1=isn1(is1)
                     enddo
                  endif
               endif

c    all checks OK: add i1(3-i) parent position to i1(i) parent neighbor table
               ip1=iparent(i1(i))
               ip2=iparent(i1(3-i))
               if (nbr1(ip1).eq.0) then
                  n_nbrs=n_nbrs+1
                  nbr1(ip1)=n_nbrs
                  nbor(n_nbrs)=ip2
                  link(n_nbrs)=0
c nbr1>0 as check would have failed if <0
               else
                  j=nbr1(ip1)
                  if (nbor(j).eq.ip2) goto 230
                  do while(link(j).ne.0)
                     j=link(j)
                     if (nbor(j).eq.ip2) goto 230
                  enddo
                  n_nbrs=n_nbrs+1
                  nbor(n_nbrs)=ip2
                  link(j)=n_nbrs
                  link(n_nbrs)=0
230               continue
               endif

c    if useisn flag on,
c    also add i1(3-i) child position to i1(i) child neighbor table
c    Note: ONLY for multimatl-multimatl connections
c     are child neighbors necessary to do multimatl edge check
              if (useisn.ne.0.and.isn1(ip1).gt.0) then
                  ip1=i1(i)
                  ip2=i1(3-i)
                  if (child_nbr1(ip1).eq.0) then
                     child_n_nbrs=child_n_nbrs+1
                     child_nbr1(ip1)=child_n_nbrs
                     child_nbor(child_n_nbrs)=ip2
                     child_link(child_n_nbrs)=0
                  else
                     j=child_nbr1(ip1)
                     if (child_nbor(j).eq.ip2) goto 240
                     do while(child_link(j).ne.0)
                        j=child_link(j)
                        if (child_nbor(j).eq.ip2) goto 240
                     enddo
                     child_n_nbrs=child_n_nbrs+1
                     child_nbor(child_n_nbrs)=ip2
                     child_link(j)=child_n_nbrs
                     child_link(child_n_nbrs)=0
240                  continue
                  endif
               endif

250            continue

C    end loop on endpoints
            enddo

c   end loop on edges
         enddo

c  end loop on elements
      enddo

c  if useisn flag on,
c  need to recheck multimatl-multimatl connections re edge materials
      if (useisn.ne.0) then
         n_dud=0
         do ip1=1,nnodes
            if (isn1(ip1).gt.0) then
               k=nbr1(ip1)
               is1=0
               do while (k.gt.0)
c                   isn1(ip2) must also be >0 if no error above
c                   => multimatl-multimatl connection
c                   => check that for every child of ip1,
c                   there is a child_nbor whose parent is ip2;
c                   and if not, remove from parent nbor table
                  ip2=nbor(k)
                  i=isn1(ip1)
                  do while (i.ne.ip1)
                     j=child_nbr1(i)
                     do while (j.ne.0)
                        if (iparent(child_nbor(j)).eq.ip2) goto 205
                        j=child_link(j)
                     enddo
                     n_dud=n_dud+1
                     k=link(k)
                     if (is1.eq.0) then
                        nbr1(ip1)=k
                     else
                        link(is1)=k
                     endif
                     goto 206
205                  continue
                     i=isn1(i)
c  end loop on isn1 chain of ip1
                  enddo
                  is1=k
                  k=link(k)
206               continue
               enddo
            endif
         enddo
         n_nbrs=n_nbrs-n_dud
      else
         n_dud=0
      endif
c

c redo linked list nbor into ordered list neigh:
c  actually could use nbor table and link, but create neigh table
c  as I wanted to look at for debugging: might eventually want to remove
c  to make more efficient?
      itype=1
      length=n_nbrs+1
      call mmggetbk('neigh',isubname,ip_neigh,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk neigh'
      itype=1
      length=nnodes
      call mmggetbk('neigh_off',isubname,ip_neigh_off
     &		,length,itype,ierr)
      if (ierr.ne.0) stop 'NBRSMOOTH: mmggetbk neigh_off'

      j=0
      do i=1,nnodes
         neigh_n(i)=0
         neigh_off(i)=j
         k=nbr1(i)
         do while (k.gt.0)
            j=j+1
            neigh_n(i)=neigh_n(i)+1
            neigh(j)=nbor(k)
            k=link(k)
        enddo
      enddo

c--------------------------------------------------------------------
ccc
cccc ........... (test nbr groups) .............
ccc
ccc	is1=0 ! so not used, but keep code here for debugging -> switch to idebug?
ccc
ccc	if (is1.gt.0) then
ccc
ccc	if (useisn.ne.0) write(*,*) 'removed ',n_dud,' nbrs as bad bond color'
ccc
ccc	is2=0
ccc	do i=1,nnodes
ccc	  if (neigh_n(i).gt.is2) is2=neigh_n(i)
ccc	enddo
ccc	do j=1,is2
ccc	 k=0
ccc	 do i=1,nnodes
ccc	  if (neigh_n(i).eq.j) k=k+1
ccc	 enddo
ccc	 if (k.gt.0) then
ccc	  write(*,'(a,i5,a,i8)') '# sites w ',j,' nbrs: ',k
ccc	 endif
ccc	enddo
ccc
cccc  just count number in isn1 chain vs "burn"??
cccc  also, should make attribute so can view,
cccc  and need to add neigh_type back to declarations
ccc	do i=1,nnodes
ccc	  neigh_type(i)=0
ccc	  if (itp1(i).ge.10.and.itp1(i).lt.20)
ccc     &		neigh_type(i)=neigh_type(i)+1
ccc	  if (icr1(i).ne.0)
ccc     &		neigh_type(i)=neigh_type(i)+icontab(1,icr1(i))
ccc	  if (isn1(i).ne.0) then
ccc	    j=isn1(i)
ccc	    if (itp1(i).eq.41.and.itp1(j).ge.10.and.itp1(j).lt.20)
ccc     &		neigh_type(i)=neigh_type(i)+1
ccc	    k=0
ccc	    do while (j.ne.i.and.j.ne.0)
ccc	     k=k+1
ccc	     j=isn1(j)
ccc	    enddo
ccc	    neigh_type(i)=neigh_type(i)+k
ccc	  endif
ccc	enddo
cccc   still not quite "spin" as icr....
ccc
ccc	endif
ccc
c--------------------------------------------------------------------

c ........... (local smooth - multi matl) .............
c   assuming that itp,icr checks on neigh table enough to fix that
c   can't arbitrarily move interior,exterior bdry defining nodes...

c check that the rlxwt won't lead to the iterative procedure blowing up
      local_rlxwt=rlxwt
      if (rlxwt.lt.zero .or. rlxwt.ge.one) local_rlxwt=half

c do the smooth the desired number of times
      do itime=1,ntimes

c   initialize weights and average neighbor positions for this itime of smoothing
c   Note: neigh_n already masked by mpary, and mpary could contain multiple children
         do i = 1,nnodes
            dx(i)=zero
            dy(i)=zero
            dz(i)=zero
            if (neigh_n(i).gt.0) then
               do j=1,neigh_n(i)
                  k=neigh(j+neigh_off(i))
                  dx(i)=dx(i)+xic(k)
                  dy(i)=dy(i)+yic(k)
                  dz(i)=dz(i)+zic(k)
               enddo
               dx(i)=dx(i)/dble(neigh_n(i))
               dy(i)=dy(i)/dble(neigh_n(i))
               dz(i)=dz(i)/dble(neigh_n(i))
               weight(i)=rlxwt
            else
               weight(i)=zero
            endif
         enddo
         nwt=nwttry

c   check volumes at new location and flag nodes bordering inverted elements
300      pos_vol=1
         do iel=1,nelements
            do j=1,4
              k=iparent(itet(j,iel))
	      xxx=one-weight(k)
	      rrr=weight(k)
	      x1(j)=xxx*xic(k)+rrr*dx(k)
	      y1(j)=xxx*yic(k)+rrr*dy(k)
	      z1(j)=xxx*zic(k)+rrr*dz(k)
	    enddo
	    vvv = ( (x1(4)-x1(1))*((y1(2)-y1(1))*(z1(3)-z1(1))
     &					-(z1(2)-z1(1))*(y1(3)-y1(1)))
     &		+(y1(4)-y1(1))*((z1(2)-z1(1))*(x1(3)-x1(1))
     &					-(x1(2)-x1(1))*(z1(3)-z1(1)))
     &		+(z1(4)-z1(1))*((x1(2)-x1(1))*(y1(3)-y1(1))
     &					-(y1(2)-y1(1))*(x1(3)-x1(1)))
     &		) * one_sixth
	    if (vvv.le.zero) then
               pos_vol=0
               do j=1,4
                  k=iparent(itet(j,iel))
                  neigh_n(k)=-iabs(neigh_n(k))
               enddo
            endif
         enddo

c   if inverting volume and if haven't tried to many times,
c    try to locally reduce weight for smoothing (and then go back to volume check);
c   otherwise punt
         if (pos_vol.eq.0) then
             nwt=nwt-1
             if (nwt.gt.0) then
                do i=1,nnodes
                   if (neigh_n(i).lt.0) then
                      neigh_n(i)=-neigh_n(i)
                      weight(i)=local_rlxwt*weight(i)
                   endif
                enddo
                goto 300
             elseif (nwt.eq.0) then
                do i=1,nnodes
                   if (neigh_n(i).lt.0) then
                      neigh_n(i)=-neigh_n(i)
                      weight(i)=zero
                   endif
                enddo
                goto 300
             else
                do i=1,nnodes
                   if (neigh_n(i).lt.0) neigh_n(i)=-neigh_n(i)
                   weight(i)=zero
                enddo
c  no need to reset or interpolate
                goto 400
             endif
          endif

c   reset parent positions (children done at end):
c   move parent node to weighted current location plus av of neighbors
c   Note: neigh_n already masked by mpary, and mpary could contain multiple children
          do i=1,nnodes
c  should be poisitive, but to be safe...
             if (neigh_n(i).ne.0) then
                neigh_n(i)=iabs(neigh_n(i))
                xxx=one-weight(i)
                rrr=weight(i)
                xic(i)=xxx*xic(i)+rrr*dx(i)
                yic(i)=xxx*yic(i)+rrr*dy(i)
                zic(i)=xxx*zic(i)+rrr*dz(i)
             endif
          enddo

c  ........... (interpolate) .............
c ADD INTERPOLATION IF NEEDED HERE...
c   could store "orig position" and interpolate at end, or interpolate here....

c   end loop on ntimes
      enddo
400   continue

c--------------------------------------------------------------------
c  ........... (fix child nodes) .............

c   spread new position of parent to child nodes
      do i=1,nnodes
         if (neigh_n(i).gt.0.and.itp1(i).eq.ifitpcup) then
            j=isn1(i)
            do while (j.ne.i.and.j.ne.0)
               xic(j)=xic(i)
               yic(j)=yic(i)
               zic(j)=zic(i)
ccc		 ! ! only if these are added as attributes re debugging...
ccc	         ! neigh_n(j)=neigh_n(i)
ccc	         ! dx(j)=dx(i)
ccc	         ! dy(j)=dy(i)
ccc	         ! dz(j)=dz(i)
ccc	         ! weight(j)=weight(i)
               j=isn1(j)
            enddo
         endif
      enddo

c--------------------------------------------------------------------

c  ........... (interpolate) .............
c ADD INTERPOLATION IF NEEDED HERE...

c--------------------------------------------------------------------

c  ........... (clean up/return) .............
9999  continue
c  release temp memory ............
      call mmrelprt(isubname,ierr)
      return

      end

c--------------------------------------------------------------------

max_forks_repo_name: daniellivingston/LaGriT


text:
Our mission:
To make a positive impact on the world, one person at a time, by sharing our gifts of healing, education, and genuine care.

We are very proud of our unique approach. We take the time to address the whole person, to educate, and to answer your questions. We understand the connection between the physical structure of your body, your nervous system, your biochemical environment, and your energetic/emotional health. We seek to achieve balance of these things through commonsense, natural means. We do this by combining gentle approaches such as specific spinal adjusting, Cranial Sacral therapy, passive motion, and energetic techniques.

We provide gentle adjusting for adults and children; including pregnant women by combining techniques such as:

Chiropractic Diversified spinal adjusting
Cranial Sacral Therapy
Passive motion (Leander flexiondistraction)
Energetic/emotional approaches

Offering holistic solutions to:

Headaches, Ear Infections, ADHD, allergies, asthma, colic, anxiety/depression, scoliosis, immune weakness, musculoskeletal injuries, Sciatica and more.

Gerner Chiropractic is one of many other Chiropractors in Davis.

	 Comments

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module CART_COE
! Description:
!    Routines to convert between Cartesian coordinates in the inertial frame
!    and classical orbital elements.
!
! Author:
!    Davide Amato
!    Space Dynamics Group - Technical University of Madrid
!    davideamato@email.arizona.edu
!
! ==============================================================================

use KINDS,      only: dk
use PHYS_CONST, only: pi,twopi,mzero
implicit none

contains

subroutine COE2CART(COE,R,V,GM,tolKE)
! Description:
!    Computes Cartesian state (R,V) from classical orbital elements COE.

! MODULES
use KEPLER, only: KESOLVE

! VARIABLES
implicit none
! Arguments
real(dk),intent(in)           ::  COE(1:6),GM
real(dk),intent(in),optional  ::  tolKE
real(dk),intent(out)          ::  R(1:3),V(1:3)
!logical,intent(out)           ::  KEflag
! Locals
real(dk)  ::  a,e,i,O,w,M,p
real(dk)  ::  EA
real(dk)  ::  sini,cosi,sinO,cosO,sinw,cosw,sinE,cosE,sinf,cosf
real(dk)  ::  R_pf(1:3),V_pf(1:3)
real(dk)  ::  ROT(1:3,1:3)

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

! CLASSICAL ORBITAL ELEMENTS
! N.B: Units of length have to be consistent between COE(1) and GM; COE(2:6) are
! in radians.
!
! COE(1): a, semi-major axis
! COE(2): e, eccentricity
! COE(3): i, inclination
! COE(4): O, right ascension of ascending node
! COE(5): w, argument of periapsis
! COE(6): M, mean anomaly

! Modification (15.01.2016): tolKE is the tolerance for the iterative solver of
! the Kepler equation. If not present, it is set to the machine epsilon

! ==============================================================================================
! 01. UNPACKING OF STATE VECTOR AND CHECK FOR SPECIAL CASES
! ==============================================================================================

a = COE(1); e = COE(2); i = COE(3)
O = COE(4); w = COE(5); M = COE(6)
p = a*(1._dk - e**2)

! Bring M to [-2pi,2pi]
M = mod(M,real(twopi,dk))

! Elliptical, equatorial (assume that 'w' is the true longitude of periapsis)
!if (abs(sin(i)) <= zero) then
!! Circular, inclined (assume that 'nu' is the argument of latitude)
!if (e <= zero) w = 0._dk
!! The case "Circular, equatorial" is covered now, assuming that 'nu' is the true longitude.

sini = sin(i); cosi = cos(i)
sinO = sin(O); cosO = cos(O)
sinw = sin(w); cosw = cos(w)
! Compute true anomaly
if (present(tolKE)) then
    EA = KESOLVE(e,M,tolKE)
else
    EA = KESOLVE(e,M,mzero)
end if
sinE = sin(EA); cosE = cos(EA)
sinf = (sqrt(1._dk - e**2)*sinE)/(1._dk - e*cosE)
cosf = (cosE - e)/(1._dk - e*cosE)

! ==============================================================================
! 02. COMPUTE POS, VEL IN PERifOCAL FRAME
! ==============================================================================

R_pf = p*[cosf/(1._dk+e*cosf), sinf/(1._dk+e*cosf), 0._dk]
V_pf = sqrt(GM/p)*[-sinf, e + cosf, 0._dk]

! ==============================================================================
! 03. COMPUTE ROTATION MATRIX
! ==============================================================================

ROT(1,1) =  cosO*cosw - sinO*sinw*cosi
ROT(2,1) =  sinO*cosw + cosO*sinw*cosi
ROT(3,1) =  sinw*sini

ROT(1,2) = -cosO*sinw - sinO*cosw*cosi
ROT(2,2) = -sinO*sinw + cosO*cosw*cosi
ROT(3,2) =  cosw*sini

ROT(1,3) =  sinO*sini
ROT(2,3) = -cosO*sini
ROT(3,3) =  cosi

! ==============================================================================================
! 04. ROTATE VECTORS TO INERTIAL FRAME
! ==============================================================================================

R = matmul(ROT,R_pf)
V = matmul(ROT,V_pf)

end subroutine COE2CART

subroutine CART2COE(R,V,COE,mu)
! Description:
!    Computes classical orbital elements COE from Cartesian coordinates R_in,V_in.
!    Units have to be consistent.
!
! ==============================================================================

! VARIABLES
implicit none
! Arguments
real(dk),intent(in)   ::  R(1:3),V(1:3)
real(dk),intent(in)   ::  mu
real(dk),intent(out)  ::  COE(1:6)
! Locals
real(dk)              ::  p,a,e             ! Orbital elements
real(dk)              ::  inc,RAAN,AoP,nu,M ! Orbital elements
real(dk)              ::  En                ! Total energy
real(dk)              ::  h(1:3),n(1:3)     ! Angular momentum and nodal vector
real(dk)              ::  ecc(1:3)          ! Eccentricity vector
real(dk)              ::  sinn,cosn
real(dk)              ::  EA,sinE,cosE
real(dk)              ::  GA                ! Gudermannian anomaly
real(dk)              ::  rmag,nmag,hmag
real(dk)              ::  zero

! ==============================================================================
!
!                                            EXECUTION
!
! ==============================================================================

! CLASSICAL ORBITAL ELEMENTS
! COE(1): a,  semi-major axis
! COE(2): e,  eccentricity
! COE(3): i,  inclination
! COE(4): Om, right ascension of ascending node
! COE(5): w,  argument of periapsis OR longitude of periapsis
! COE(6): M,  mean anomaly

! =============================================================================
! 02. COMPUTE p,a,ecc,e
! =============================================================================
zero = 2._dk*mzero

rmag   = sqrt(dot_product(R,R))

h(:)   = CROSS_PRODUCT(R,V)
n(:)   = CROSS_PRODUCT([0._dk,0._dk,1._dk],h)
ecc(:) = ((DOT_PRODUCT(V,V) - mu/rmag)*R - DOT_PRODUCT(R,V)*V)/mu
e      = sqrt(dot_product(ecc,ecc))

nmag   = sqrt(dot_product(n,n))

En     = 0.5_dk*DOT_PRODUCT(V,V) - mu/rmag

if (ABS(En) > zero) then
    a  = -mu/(2._dk*En)
    p  = a*(1._dk-e**2)

else ! Exactly parabolic orbit, should be extremely rare
    a  = HUGE(0._dk)
    p  = DOT_PRODUCT(h,h)/mu

end if

! =============================================================================
! 03. COMPUTE inc,RAAN,AoP,nu
! =============================================================================

hmag = sqrt(dot_product(h,h))
inc  = acos(h(3)/hmag)
RAAN = ACOS(n(1)/nmag); if (n(2) < 0._dk) RAAN = twopi - RAAN
AoP  = ACOS(DOT_PRODUCT(n,ecc)/(nmag*e)); if (ecc(3) < 0._dk) AoP  = twopi - AoP

nu   = ACOS(DOT_PRODUCT(ecc,R)/e/rmag)
if (DOT_PRODUCT(R,V) < 0._dk) then
    nu   = twopi - nu

end if

! SPECIAL CASES:
! Elliptical, equatorial. AoP = true longitude of periapsis
if (ABS(ABS(h(3)/hmag)-1._dk) < zero) then
    RAAN = 0._dk
    AoP  = ACOS(ecc(1)/e); if (ecc(2) < 0._dk) AoP = twopi - AoP
end if
! Circular, inclined.   nu = argument of latitude
if (e <= zero) then
    AoP = 0._dk
    nu  = ACOS(DOT_PRODUCT(n,R)/nmag/rmag); if (r(3) < 0._dk) nu = twopi - nu
    M = nu
end if
! Circular, equatorial. nu = true longitude.
if (ABS(ABS(h(3)/hmag)-1._dk) < zero .AND. e <= zero) then
    RAAN = 0._dk
    AoP  = 0._dk
    nu   = ACOS(R(1)/rmag); if (r(2) < 0._dk) nu = twopi - nu
end if

! Mean anomaly (see Battin, Ch. 4)
sinn = sin(nu); cosn = cos(nu)
if (e < 1._dk) then
  sinE = (sinn*sqrt(1._dk - e**2))/(1._dk + e*cosn)
  cosE = (e + cosn)/(1._dk + e*cosn)
  EA = atan2(sinE,cosE)
  M = EA - e*sinE

else if (e > 1._dk) then
  GA = 2._dk * atan( sqrt( (e - 1._dk) / (e + 1._dk) ) * tan(0.5_dk * nu) )
  M = e * tan(GA) - log( tan( 0.5_dk * GA + 0.25_dk * pi ))

end if

! Wrap M to [0,360]
M = mod(M + twopi,twopi)

COE = [a,e,inc,RAAN,AoP,M]

contains

function CROSS_PRODUCT(a,b)
! Description:
!    Computes cross product of a*b.
! ==============================================================================
!                                            PREAMBLE
! ==============================================================================

! VARIABLES
implicit none
! Function definition
real(dk)            :: CROSS_PRODUCT(1:3)
! Arguments
real(dk)            :: a(1:3),b(1:3)

! ==============================================================================
!                                            EXECUTION
! ==============================================================================
CROSS_PRODUCT(1) = a(2)*b(3) - a(3)*b(2)
CROSS_PRODUCT(2) = a(3)*b(1) - a(1)*b(3)
CROSS_PRODUCT(3) = a(1)*b(2) - a(2)*b(1)
end function CROSS_PRODUCT

end subroutine CART2COE

end module CART_COE

max_stars_repo_name: isaiasAaraya/TLEanalysis

max_issues_repo_name: isaiasAaraya/TLEanalysis

max_forks_repo_name: isaiasAaraya/TLEanalysis


text:
subroutine djpi2(dj, lmax, exitstatus)
!------------------------------------------------------------------------------
!
!   This subroutine computes
!     j
!   d     (pi/2)
!     m N
!
!   for all posible values of m (>=0) and N for 0<l<lmax.
!   The output array corresponds to dj(N,m,j). (I'm not positive
!   about the N and m ordering).
!
!   Note that this algorithm is easy to modify to work on a
!   single l (just need a loop to compute f1) since there are
!   no l recursions involved.
!
!   Calling Parameters
!
!       IN
!           lmax    Maximum spherical harmonic degree to be computed.
!
!       OUT
!           dj  Rotation matrix with dimensions (lmax+1, lmax+1, lmax+1).
!
!       OPTIONAL (OUT)
!           exitstatus  If present, instead of executing a STOP when an error
!                       is encountered, the variable exitstatus will be
!                       returned describing the error.
!                       0 = No errors;
!                       1 = Improper dimensions of input array;
!                       2 = Improper bounds for input variable;
!                       3 = Error allocating memory;
!                       4 = File IO error.
!
!   History
!
!       1. Based on routine from Guy Masters (July 16, 1993)
!       2. Modified by Mark Simons (July 19, 1993)
!       3. Turned into readable f95 code by Mark Wieczorek (August, 2003).
!
!   Copyright (c) 2005-2019, SHTOOLS
!   All rights reserved.
!
!------------------------------------------------------------------------------
    use ftypes

    implicit none

    integer, intent(in) :: lmax
    real(dp), intent(out) :: dj(:,:,:)
    integer, intent(out), optional :: exitstatus
    integer :: i, l, n, lp1, j, m, isn, np
    real(dp) :: f((lmax+1)*8), f1, f2, g1, g2, en2, fl2p1

    if (present(exitstatus)) exitstatus = 0

    if (size(dj(:,1,1)) < lmax+1 .or. size(dj(1,:,1)) < lmax + 1 .or. &
            size(dj(1,1,:)) < lmax +1) then
        print*, "Error --- djpi2"
        print*, "DJ must be dimensioned (LMAX+1, LMAX+1, LMAX+1) " // &
                "where LMAX is ", lmax
        print*, "Input array is dimensioned as ", size(dj(:,1,1)), &
                size(dj(1,:,1)), size(dj(1,1,:))
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    end if

    dj = 0.0_dp

    dj(1,1,1) = 1.0_dp
    dj(1,1,2) = 0.0_dp
    dj(1,2,2) = -1.0_dp / sqrt(2.0_dp)
    dj(2,1,2) = -dj(1,2,2)
    dj(2,2,2) = 0.5_dp
    f1        = 0.5_dp

    do l = 2, lmax

        lp1 = l + 1
        fl2p1 = l + lp1

        do i = 1, l
            f(i) = dsqrt(i * (fl2p1-i))
        end do

            f1 = f1 * (l+l-1.0_dp) / (l+l)

        ! Do N = 0 terms
        dj(lp1,1,lp1) = -dsqrt(f1)
        dj(l,1,lp1)   = 0.0_dp

        do i = 2, l
            j = lp1 - i
            dj(j,1,lp1) = -f(i-1) * dj(j+2,1,lp1) / f(i)
        end do

        ! Do positive N terms (bottom triangle)
        f2 = f1
        g1 = l
        g2 = lp1

        do n = 1, l
            np = n + 1
            en2 = n + n
            g1 = g1 + 1.0_dp
            g2 = g2 - 1.0_dp
            f2 = f2 * g2 / g1
            dj(lp1,np,lp1) = -dsqrt(f2)
            dj(l,np,lp1) = dj(lp1,np,lp1) * en2 / f(1)

            do i = 2, l - n
                j = lp1 - i
                dj(j,np,lp1) = ( en2 * dj(j+1,np,lp1) &
                               - f(i-1) * dj(j+2,np,lp1) ) / f(i)
            end do

        end do

        ! Fill upper triangle and fix signs
        do j = 1, l
            do m = j, l
                dj(j,m+1,lp1) = dj(m+1,j,lp1)
            end do
        end do

        isn = 1 + mod(l,2)

        do np = 1, lp1
            do i = isn, lp1, 2
                dj(i,np,lp1) = -dj(i,np,lp1)
            end do
        end do

    end do

end subroutine djpi2

max_stars_repo_name: dilkins/SHTOOLS

max_issues_repo_name: dilkins/SHTOOLS

max_forks_repo_name: dilkins/SHTOOLS


text:
! Fortran code for computing response spectra.
! Originial version is writting by: Leonardo Ramirez-Guzman
! The code base is tailored for F2py.
!

subroutine max_osc_response(acc, npts, dt, period, csi, max_disp, &
    max_vel, max_acc)

    implicit none
    integer, intent(in) :: npts
    real, intent(in) :: acc(npts)
    real :: a(npts), v(npts), d(npts)
    real, intent(in) :: dt, period, csi
    real, intent(out) :: max_disp, max_vel, max_acc
    integer :: i=1
    real, parameter :: PI = 3.14159265359
    real :: w, ww, csicsi, dcsiw, rcsi, csircs, wd, ueskdt, dcsiew, um2csi, e, &
     s, c0, ca, cb, cc, cd, cap, cbp, ccp, cdp

    d(1) = 0.
    v(1) = 0.

    w = (2*PI)/period
    ww = w ** 2
    csicsi = csi ** 2
    dcsiw = 2. * csi * w

    rcsi = sqrt(1-csicsi)
    csircs = csi/rcsi
    wd = w * rcsi
    ueskdt = -1/(ww*dt)
    dcsiew = 2 * csi / w
    um2csi = (1.-2*csicsi)/wd
    e=exp(-w*dt*csi)
    s=sin(wd*dt)
    c0=cos(wd*dt)

    a(1) = -ww*d(1)-dcsiw*v(1)

    ca = e * (csircs*s+c0)
    cb = e * s/wd
    cc = (e*((um2csi-csircs*dt)*s -(dcsiew+dt)*c0)+dcsiew)*ueskdt
    cd = (e*(-um2csi*s+dcsiew*c0)+dt-dcsiew)*ueskdt
    cap = -cb * ww
    cbp = e * (c0-csircs*s)
    ccp = (e*((w*dt/rcsi+csircs)*s+c0)-1.)*ueskdt
    cdp = (1.-ca)*ueskdt

    do i=2,npts
        d(i) = ca*d(i-1)+cb*v(i-1)+cc*acc(i-1)+cd*acc(i)
        v(i) = cap*d(i-1)+cbp*v(i-1)+ccp*acc(i-1)+cdp*acc(i)
        a(i) = -ww*d(i)-dcsiw*v(i)
    end do

    max_disp = 0
    max_vel = 0
    max_acc = 0

    do i=1,npts
        if (abs(a(i)) .gt. max_acc) then
            max_acc = abs(a(i))
        end if
    end do

    do i=1,npts
        if (abs(v(i)) .gt. max_vel) then
            max_vel = abs(v(i))
        end if
    end do

    do i=1,npts
        if (abs(d(i)) .gt. max_disp) then
            max_disp = abs(d(i))
        end if
    end do

end subroutine max_osc_response

subroutine calc_response(acc, npts, dt, periods, n_prds,csi, max_disps, &
     max_vels, max_accs)
    implicit none
    integer, intent(in) :: npts, n_prds
    real, intent(in) :: acc(npts)
    real, intent(in) :: dt, csi
    real, intent(in) :: periods(n_prds)
    real, intent(out) :: max_disps(n_prds), max_vels(n_prds), max_accs(n_prds)
    integer :: i=1

    do i=1,n_prds
        call max_osc_response(acc, npts, dt, periods(i), csi, max_disps(i),&
         max_vels(i), max_accs(i))
    end do

end subroutine calc_response

max_stars_repo_name: Naeemkh/ts-process

max_issues_repo_name: Naeemkh/ts-process

max_forks_repo_name: Naeemkh/tsprocess


text:
        subroutine gensav (ientmp, mattmp, ien,    mater)
c
c----------------------------------------------------------------------
c
c  This routine saves the element block data.
c
c input:
c  ientmp (npro,nshl)   : nodal connectivity
c  mattmp (npro)        : material type flag
c
c output:
c  ien    (npro,nshl)   : nodal connectivity
c  mater  (npro)        : material type flag
c
c
c Zdenek Johan, Winter 1992.
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension   ientmp(npro,nshl),
     &              mattmp(npro),           ien(npro,nshl),
     &              mater(npro)
c
c.... save the element data
c
        do i = 1, nshl
          ien(:,i) = ientmp(:,i)
        enddo
c
        mater = mattmp
c
c.... end
c
        return
        end


text:
c     This version of plot2 takes a UNIX SIOSEIS seismic plot file from a
c  SIOSEIS disc file and writes directly to the Versatec plotter.  PLOT2 opens
c  all files for the user and assumes the Versatec is /dev/lpv0
c     This program is a beast because UNIX doesn't allow the use of the VERSATEC
c  SSP feature (simultaneous print/plot), therefore this program ors in a
c  rasterized version of the annotation.
c     Another problem is that CR/LF are need.
c
c  Copyright (C) 1985
c  Paul Henkart, Scripps Institution of Oceanography, La Jolla, Ca. 92093
c

      PARAMETER (nchars=80)
      INTEGER*2 ibuf(5000)
      CHARACTER*80 token
      CHARACTER*90 prtbuf
      DATA prtbuf/' '/

    1 PRINT *,' Enter the value used for SIOSEIS plot parameter nibs.'
      READ *, itype
      IF( itype .NE. 160 .AND. itype .NE. 200 .AND.
     *    itype .NE. 8122. AND. itype .NE. 8222 .AND.
     &    itype .NE. 3436 .AND.
     *    itype .NE. 201 .AND. itype .NE. 5732 .AND.
     *    itype .NE. 7222 .AND. itype .NE. 7224 .AND.
     *    itype .NE. 7236 .AND. itype .NE. 7244 .AND.
     *    itype .NE. 7422 .AND. itype .NE. 7424 .AND.
     *    itype .NE. 7436 .AND. itype .NE. 7444 .AND.
     *    itype .NE. 8625 .AND. itype .NE. 3444 ) THEN
          PRINT *,' Incorrect answer, try again.'
          GOTO 1
      ENDIF
      IF( itype .EQ. 160 ) mwrds = 90                                   ! the number of 32 bit words per raster line
      IF( itype .EQ. 200 ) mwrds = 66
      IF( itype .EQ. 201 ) mwrds = 146
      IF( itype .EQ. 3444 ) mwrds = 538
      IF( itype .EQ. 5732 ) mwrds = 147
      IF( itype .EQ. 8122 ) mwrds = 66
      IF( itype .EQ. 8222 ) mwrds = 132
      IF( itype .EQ. 7222 ) mwrds = 132
      IF( itype .EQ. 7224 ) mwrds = 144
      IF( itype .EQ. 7225 ) mwrds = 147
      IF( itype .EQ. 7236 ) mwrds = 220
      IF( itype .EQ. 7244 ) mwrds = 269
      IF( itype .EQ. 7422 ) mwrds = 264
      IF( itype .EQ. 7424 ) mwrds = 288
      IF( itype .EQ. 7425 ) mwrds = 294
      IF( itype .EQ. 7436 ) mwrds = 440
      IF( itype .EQ. 7444 ) mwrds = 538
      IF( itype .EQ. 8122 ) mwrds = 132*2
      IF( itype .EQ. 8222 ) mwrds = 264*2
      IF( itype .EQ. 8242 ) mwrds = 512*2
      IF( itype .EQ. 8625 ) mwrds = 294
      IF( itype .EQ. 9242 ) mwrds = 500*2
      IF( itype .EQ. 3436 ) mwrds = 856 / 2
      nbytes = mwrds * 4                                                ! the number of bytes per raster line
      PRINT *,' Enter the name of the input SIOSEIS plot file:'
      CALL rdline
      CALL getoke( token, nchars )
      token( nchars+1:nchars+1) = CHAR(0)
      CALL getfil( 4, ilun, token, istat )                               ! open the plot file
      CALL attach
c****
c****   Get the Versatec toner flowing and get the plot out of the dried toner
c****  by printing 20 blank lines.  FF (Form Feed) doesn't get the toner flowing
C****
      CALL vprint                                                       ! set UNIX to put the Versatec in print mode
      prtbuf = ' '
      prtbuf(10:10) = CHAR(10)                                          ! put in a CR
      DO 20 i = 1, 20
         CALL vwrite ( prtbuf, 10 )                                     ! write 10 characters
   20 CONTINUE
c****
c****    PRINT 44 LINES.  Line 1 is the date of the plot.  lines 2-41 are the
c****  SEGY tape header/SIOSEIS comments.  42-44 contain some plotting parameters.
c****
      DO 50 i = 1,44
         CALL rddisc( ilun, prtbuf, 20, istat )                          ! read 20 words (80 characters)
         prtbuf(81:81) = CHAR(10)                                       ! put in a CR
         CALL vwrite( prtbuf, 81 )
   50 CONTINUE
      prtbuf = ' '                                                      ! make a gap between the header and the data
      prtbuf(10:10) = CHAR(10)                                          ! put in a CR
      DO 70 i = 1, 10
         CALL vwrite( prtbuf, 10 )                                      ! write 10 characters
   70 CONTINUE
      CALL vplot
c****
  100 CONTINUE
      CALL rddisc( ilun, ibuf, mwrds, istat )
      IF( istat .NE. mwrds ) GOTO 9999
      CALL vwrite( ibuf, nbytes )
      GOTO 100
c****
c****    END OF PLOT FILE
c****
 9999 CONTINUE
      CALL vprint
      prtbuf = ' '
      prtbuf(10:10) = CHAR(10)                                          ! put in a CR
      DO 10000 i = 1, 60
         CALL vwrite( prtbuf, 10 )                                      ! write 10 characters
10000 CONTINUE
      CALL detach
      PRINT *,' finished the seismic plot.'

      END


program n_body_test
    use n_bodies
    implicit none
    integer(8), parameter :: n = 2, n_times = 100000
    real :: q = 1e-3
    real(8), dimension(n) :: masses
    real(8), dimension(n, 3) :: i_positions, i_velocities, poses, vels
    real(8), dimension(n_times, n, 3) :: velocities, positions
    real(8), dimension(n_times) :: times
    real(8) :: dt, t_tot, g, tolerance
    logical :: update_com
    character(256) :: status
    integer(8) :: i, j, k, q_dt
    masses = (/1. - q, q/)
    i_positions(1, :) = (/-q, 0., 0./)
    i_positions(2, :) = (/1. - q, 0., 0./)
    i_velocities(1, :) = (/0., -q, 0./)
    i_velocities(2, :) = (/0., 1. - q, 0./)
    dt = 1e-3
    t_tot = n_times * dt
    tolerance = 1e-10
    q_dt = 4
    g = 1.
    update_com = .true.

    print *, "Euler"
    call n_body_model("euler", n, n_times, masses, i_positions, i_velocities, &
                      t_tot, dt, g, update_com, times, positions, velocities, &
                      status)
    print *, status

    print *, "RK5"
    call n_body_model("rk5", n, n_times, masses, i_positions, i_velocities, &
                  t_tot, dt, g, update_com, times, positions, velocities, &
                  status)
    print *, status
!
    print *, "RKF45"
    status = "Ok"
    call adaptive_n_body_model("rkf45", n, n_times, masses, i_positions, &
            i_velocities, t_tot, tolerance, dt, q_dt, g, update_com, times, &
            positions, velocities, status)
    print *, status
!    do i=1, n_times
!        print *, positions(i, 2, :)
!    end do

end program n_body_test


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis
         subroutine readcsc(i4time_cur,imax,jmax,cscnew)
C
         integer  imax,jmax

         real csc(imax,jmax),lcv(imax,jmax),csctot(imax,jmax)
         real csc2nd(imax,jmax)
         real snow_total(imax,jmax)
         real snow_accum(imax,jmax)
         real missval,cscnew(imax,jmax)
         integer i4time,i4time_cur
         integer i4time_csc(imax,jmax),i4time_snow(imax,jmax)
         CHARACTER*24 TIME
c ************************************************************
c
         do j=1,jmax
            do i=1,imax
               snow_total(i,j)=0.0
            enddo
         enddo
        call get_laps_cycle_time(laps_cycle_time,istatus)
        if(istatus.ne.1)then
           write(6,*)'Error - get_laps_cycle_time'
           return
        endif
        call get_r_missing_data(missval,istatus)
        if(istatus.ne.1)then
           write(6,*)'Error - get_r_missing_data'
           return
        endif

        i4time=i4time_cur-48*3600 !force 48 hours

c ************************************************************
c initialize arrays
c
        csctot=missval
        csc2nd=missval
        i4time_csc = -999
        i4time_snow = -999
c
c now loop through past 48 hours, looking only at csc files
c
        rmxlcv=0.0
        rmnlcv=100.
        icsc=0
        ncycle_times=48*int(3600./float(laps_cycle_time))
        do itime=1,ncycle_times
c
           call cv_i4tim_asc_lp(i4time,time,istatus)
c           if(time(13:14).eq.'13')goto 52
           write(6,101)time
 101       format(1x,'First csc loop data search time: ',a17)
c
	   CALL GETLAPSLCV(I4TIME,LCV,CSC,IMAX,JMAX,ISTATUS)
	   IF (ISTATUS .NE. 1)go to 522

           icsc=icsc+1
           nspts=0
           nrpts=0
           do i=1,imax
           do j=1,jmax
             if(csc(i,j).le.1.1)then
                csctot(i,j)=csc(i,j)
                nspts=nspts+1
             endif
             if(csc(i,j).le.0.1)then
                nrpts=nrpts+1
                csctot(i,j)=0.0
             endif
             rmxlcv=max(lcv(i,j),rmxlcv)
             rmnlcv=min(lcv(i,j),rmnlcv)
           enddo
           enddo
           nmpts=(imax*jmax)-nspts-nrpts
           nspts=nspts-nrpts
           write(6,501)itime,nmpts,nspts,nrpts
 501       format(' Cycle ',i4,' nmiss/nsnow/nreject ',3i9)
 522       continue
c
           i4time=i4time+laps_cycle_time
        enddo
c
        if(icsc.eq.0)then
        write(6,*) '**** No lcv data available over 48 hours ****'
          do i=1,imax
          do j=1,jmax
             csctot(i,j)=missval
          enddo
          enddo
        endif
c
        call analyze(csctot,cscnew,imax,jmax,missval)
c
c now loop through past 48 hours, looking for latest snow_total and
c  latest csc obs. snow_total only applies to a point - it is not
c  spread horizontally. If there are then any later csc obs at that
c  point, then the csctot point is set to the csc ob.
c
        i4time=i4time_cur-48*3600
        ncycle_times=48*int(3600./float(laps_cycle_time))
        write(6,*)'---------------------------------------'
        write(6,*)

        do itime=1,ncycle_times
c
           call cv_i4tim_asc_lp(i4time,time,istatus)
c          if(time(13:14).eq.'13')goto 52
           write(6,102)time
 102       format(1x,'Second loop data search time: ',a17)
c
c First get csc field again. Will not analyze - just use this to
c   write over any previous time's snow_total obs.
	   CALL GETLAPSLCV(I4TIME,LCV,CSC,IMAX,JMAX,ISTATUS)
	   IF (ISTATUS .NE. 1)goto 52

           do i=1,imax
           do j=1,jmax
             if(csc(i,j).le.1.1)then
                 cscnew(i,j)=csc(i,j)
                 csc2nd(i,j)=csc(i,j)
                 i4time_csc(i,j) = i4time
             endif
	     if(csc(i,j).le.0.1)then
                 cscnew(i,j)=0.0
                 csc2nd(i,j)=0.0
                 i4time_csc(i,j) = i4time
             endif
           enddo
           enddo
 52        continue
           i4time=i4time+laps_cycle_time
        enddo
c
c now get snow_total obs
c J.Smart (11-3-97). Compute 48 hour snow accumulation instead of
c                    snow total to avoid false snow cover from "old"
c                    snow total accumulation.
c Recent snow is given more weight
c
        i4time=i4time_cur-48*3600
        ncycle_times=48*int(3600./float(laps_cycle_time))
        write(6,*)
        write(6,*)'Accumulate hrly Snow for 48 hour total '
        write(6,*)'---------------------------------------'

        do itime=1,ncycle_times

           call cv_i4tim_asc_lp(i4time,time,istatus)
           write(6,103)time
 103       format(1x,'Snow Accum data search time: ',a17)

           CALL GETLAPSL1S(I4TIME,snow_accum,imax,jmax,1,ISTATUS)
           IF (ISTATUS .EQ. 1)then
             do j=1,jmax
             do i=1,imax
               snow_total(i,j)=snow_total(i,j)
     &       +(snow_accum(i,j)*(float(itime)/float(ncycle_times)))
               if(snow_accum(i,j) .gt. 0.)then
                 i4time_snow(i,j) = i4time
               endif
             enddo
             enddo
             call array_range(snow_accum,imax,jmax,rmin1,rmax1,missval)
             call array_range(snow_total,imax,jmax,rmin2,rmax2,missval)
             write(6,111)rmin1,rmax1,rmin2,rmax2
111          format(' L1S/S01, snow total ranges ',4f9.5)
           else
             write(6,*)'WARNING - L1S/S01 unavailable'
           endif

           i4time=i4time+laps_cycle_time
        enddo
c
c Adjust snow cover field based upon snow total. ------
c   if snowfall total >= 2cm, then set snow cover to 1.
c   if snowfall total  = 1cm, then set snow cover to 0.5
c   if snowfall total between 1 and 2cm, then ramp snow cover 0.5 to 1.
c
c We added an extra condition that snowfall is more recent than satellite
c data.

        nrecent_snow = 0
        nrecent_csc2nd = 0
        n_neither = 0
        do i=1,imax
        do j=1,jmax
c         if(csc2nd(i,j) .eq. missval)then
c         if(i4time_csc(i,j) .lt. i4time_cur-86400)then
          if(i4time_snow(i,j) .gt. i4time_csc(i,j))then
            nrecent_snow = nrecent_snow + 1

            if(snow_total(i,j) .ge. 0.02)cscnew(i,j)=1.0
            if((snow_total(i,j).ge.0.01).and.(snow_total(i,j).lt.0.02))
     1         cscnew(i,j)=snow_total(i,j)/0.02

c           if(snow_total(i,j).ge.0.02.and.cscnew(i,j).gt.0.25)
c    1         cscnew(i,j)=1.0
c           if((snow_total(i,j).ge.0.01).and.(snow_total(i,j).lt.0.02))
c    1         cscnew(i,j)=snow_total(i,j)/0.02

          elseif(csc2nd(i,j) .ne. missval)then ! test csc2nd is recent
              nrecent_csc2nd = nrecent_csc2nd + 1
          else
              n_neither = n_neither + 1
          endif

        enddo
        enddo

        write(6,*)' recent snow(S01)/csc(sat)/none = ',nrecent_snow
     1                                                ,nrecent_csc2nd
     1                                                ,n_neither
c
c       do j=jmax,1,-1
c          write(6,76)j,(cscnew(i,j),i=1,18)
c 76       format(1x,i2,1x,18f4.1)
c       enddo
c
        return
        end
C-------------------------------------------------------------------------------
C
      SUBROUTINE GETLAPSLCV(I4TIME,LCV,CSC,IMAX,JMAX,ISTATUS)
C
      Integer*4 imax,jmax
      Integer*4 kmax
      parameter(kmax = 2)
c
      INTEGER*4 I4TIME,LVL(KMAX),I,J,ERROR(2),ISTATUS
C
      REAL*4 lcv(imax,jmax),csc(imax,jmax)
      Real*4 readv(imax,jmax,kmax)
C
      CHARACTER*150 LDIR
      CHARACTER*31 EXT
      CHARACTER*3 VAR(KMAX)
      CHARACTER*4 LVL_COORD(KMAX)
      CHARACTER*10 UNITS(KMAX)
      CHARACTER*125 COMMENT(KMAX)
C
C-------------------------------------------------------------------------------
C
	ERROR(1)=1
	ERROR(2)=0
C
        call get_directory('lcv',ldir,len)

	EXT='lcv'
	VAR(1)='LCV'
	VAR(2)='CSC'
	LVL(1)=0
	LVL(2)=0
C
C ****  Read LAPS lcv and csc.
C
	CALL READ_LAPS_DATA(I4TIME,LDIR,EXT,IMAX,JMAX,KMAX,KMAX,VAR,LVL,
     1      LVL_COORD,UNITS,COMMENT,readv,ISTATUS)
C
	IF (ISTATUS .NE. 1) THEN
		ISTATUS=ERROR(2)
		RETURN
	ENDIF
C
        do j=1,jmax
        do i=1,imax
          csc(i,j)=readv(i,j,2) !Cloud analysis implied snow cover
          lcv(i,j)=readv(i,j,1) !LAPS cloud cover
        enddo
        enddo
c
c        print *,' '
c        print *,'lcv field'
c        do j=1,jmax
c          write(6,20)j,(lcv(i,j),i=21,40)
c        enddo
c        print *,' '
c        print *,'csc field'
c        do j=1,jmax
c          write(6,20)j,(csc(i,j),i=21,40)
c        enddo
c        print *,'csc(21,55)=',csc(21,55)
c 20     format(1x,i2,1x,20f4.1)
c
	ISTATUS=ERROR(1)
	RETURN
C
	END
C
C-------------------------------------------------------------------------------
C
      SUBROUTINE GETLAPSL1S(I4TIME,snow_accum,imax,jmax,kmax,
     &ISTATUS)
C
      integer*4 imax,jmax,kmax
c
      INTEGER*4 I4TIME,LVL(KMAX),I,J,ERROR(2),ISTATUS
C
      REAL*4 snow_accum(imax,jmax),readv(imax,jmax,kmax)
C
      CHARACTER*150 LDIR
      CHARACTER*31 EXT
      CHARACTER*3 VAR(KMAX)
      CHARACTER*4 LVL_COORD(KMAX)
      CHARACTER*10 UNITS(KMAX)
      CHARACTER*125 COMMENT(KMAX)
C
C-------------------------------------------------------------------------------
C
	ERROR(1)=1
	ERROR(2)=0
C
        call get_directory('l1s',ldir,len)

	EXT='L1S'
	VAR(1)='S01'
	LVL(1)=0
C
C ****  Read LAPS l1s
C
	CALL READ_LAPS_DATA(I4TIME,LDIR,EXT,IMAX,JMAX,KMAX,KMAX,VAR,LVL,
     1      LVL_COORD,UNITS,COMMENT,readv,ISTATUS)
C
	IF (ISTATUS .NE. 1) THEN
		ISTATUS=ERROR(2)
		RETURN
	ENDIF
C
        do j=1,jmax
        do i=1,imax
          snow_accum(i,j)=readv(i,j,1)
        enddo
        enddo
c
c
	ISTATUS=ERROR(1)
	RETURN
C
	END
C
C
C-------------------------------------------------------------------------------
         subroutine analyze(csctot,cscnew,imax,jmax,missval)
C
         integer imax,jmax
         integer nboxes, nruns
         parameter(nboxes=100,nruns=5)
c
         real csctot(imax,jmax),cscnew(imax,jmax)
         dimension iimin(nruns,nboxes),
     1       iimax(nruns,nboxes),jjmin(nruns,nboxes),
     2       jjmax(nruns,nboxes),nbtot(nruns)
         real*4 missval
c
c set up boundaries for boxes to do averaging over
c
         do i=1,imax
         do j=1,jmax
            cscnew(i,j)=csctot(i,j)
         end do
         end do
         nbtot(1)=1
         nbtot(2)=4
         nbtot(3)=9
         nbtot(4)=36
         nbtot(5)=100
         call setboxes(imax,jmax,iimin,iimax,jjmin,jjmax,
     1    nbtot,nruns,nboxes)
c
c **************************************************************
c
c Don't know ahead of time how much missing data, so:
c   1) Take average for whole grid, apply to missing points (create
c        cscnew grid).
c   2) Divide grid into 4 boxes. Take average from original grid
c        in each box. If there are some non-missing points in the
c        box, apply average to missing points (modify cscnew grid).
c   3) Divide grid into 9 boxes. Take average from original grid
c        in each box. If there are some non-missing points in the
c        box, apply average to missing points (modify cscnew grid).
c     -> Actually do this again for 36 and 100 box subdivisions.
c        Thus the cscnew grid will contain the average values from
c        the smallest boxes (from the 100 box subdivision). The
c        missing grid points receive this value while the non-missing
c        grid points get overwritten with the original csc values.
c   4) Apply original grid non-missing values to cscnew grid.
c
       do nrun=1,nruns
c
       do nn=1,nbtot(nrun)
         total=0.
         sum=0.
         do i=iimin(nrun,nn),iimax(nrun,nn)
         do j=jjmin(nrun,nn),jjmax(nrun,nn)
           if(csctot(i,j).le.1.1)then
             sum=sum+csctot(i,j)
             total=total+1.
           endif
         enddo
         enddo
         if(total.eq.0.)then
c          print *,'all values missing for this box'
c          print *,'nrun,nbox=',nrun,nn
c          if(nrun.eq.1)then
c             print *,'ALL VALUES MISSING - CANNOT CONTINUE'
c             stop
c          endif
          go to 50
       else
          avrg=sum/total
c          print *,'nrun,nbox=',nrun,nn
c          print *,'box average=',avrg,' npts=',total
          do i=iimin(nrun,nn),iimax(nrun,nn)
          do j=jjmin(nrun,nn),jjmax(nrun,nn)
             cscnew(i,j)=avrg
          enddo
          enddo
       endif
 50    continue
c
       enddo
c
c       print *,'nrun,nbox=',nrun,nn
c       do j=jmax,1,-1
c       write(6,75)j,(cscnew(i,j),i=1,imax)
c       enddo
c 75    format(1x,i2,61f2.1)
c
       enddo
c
c now overwrite non-missing points
c
       do i=1,imax
       do j=1,jmax
         if(csctot(i,j).le.1.1)cscnew(i,j)=csctot(i,j)
       enddo
       enddo
c
        return
        end
c
c ***********************************************************
c
         subroutine setboxes(imax,jmax,iimin,iimax,jjmin,jjmax,
     1    nbtot,nruns,nboxes)
c
c set up boundaries for boxes to do averaging over
c
         dimension iimin(nruns,nboxes),
     1       iimax(nruns,nboxes),jjmin(nruns,nboxes),
     2       jjmax(nruns,nboxes),nbtot(nruns)
c
         do nn=1,nruns
           idir=sqrt(real(nbtot(nn)))
           inc=imax/idir+0.99
c           print *,'nn,idir,inc=',nn,idir,inc
           nbox=1
           do i=1,idir
             do j=1,idir
               iimin(nn,nbox)=1+(i-1)*inc
               iimax(nn,nbox)=min(imax,1+i*inc)
               jjmin(nn,nbox)=1+(j-1)*inc
               jjmax(nn,nbox)=min(jmax,1+j*inc)
               nbox=nbox+1
             enddo
           enddo
c           do nbox=1,nbtot(nn)
c             write(6,10)nn,nbox,inc,iimin(nn,nbox),iimax(nn,nbox),
c    1        jjmin(nn,nbox),jjmax(nn,nbox)
c10          format(1x,'nn,nbox,inc=',3i3,'imin,max,jmin,max=',4i4)
c           enddo
         enddo
c
         return
         end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      subroutine nitvol

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine estimates daily mineralization (NH3 to NO3)
!!    and volatilization of NH3

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    curyr         |none          |current year of simulation
!!    hru_dafr(:)   |km**2/km**2   |fraction of watershed area in HRU
!!    ihru          |none          |HRU number
!!    nyskip        |none          |number of years to skip output
!!                                 |summarization and printing
!!    sol_fc(:,:)   |mm H2O        |amount of water available to plants in soil
!!                                 |layer at field capacity (fc - wp)
!!    sol_nh3(:,:)  |kg N/ha       |amount of nitrogen stored in the ammonium
!!                                 |pool in soil layer
!!    sol_nly(:)    |none          |number of layers in soil profile
!!    sol_no3(:,:)  |kg N/ha       |amount of nitrogen stored in the
!!                                 |nitrate pool in soil layer
!!    sol_st(:,:)   |mm H2O        |amount of water stored in the soil layer
!!                                 |on any given day (less wp water)
!!    sol_tmp(:,:)  |deg C         |daily average temperature of soil layer
!!    sol_wpmm(:,:) |mm H20        |water content of soil at -1.5 MPa (wilting
!!                                 |point)
!!    sol_z(:,:)    |mm            |depth to bottom of soil layer
!!    wshd_nitn     |kg N/ha       |average annual amount of nitrogen moving
!!                                 |from the NH3 to the NO3 pool by
!!                                 |nitrification in the watershed
!!    wshd_voln     |kg N/ha       |average annual amount if nitrogen lost by
!!                                 |ammonia volatilization in watershed
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    sol_nh3(:,:)  |kg N/ha       |amount of nitrogen stored in the ammonium
!!                                 |pool in soil layer
!!    sol_no3(:,:)  |kg N/ha       |amount of nitrogen stored in the
!!                                 |nitrate pool in soil layer
!!    wshd_nitn     |kg N/ha       |average annual amount of nitrogen moving
!!                                 |from the NH3 to the NO3 pool by
!!                                 |nitrification in the watershed
!!    wshd_voln     |kg N/ha       |average annual amount if nitrogen lost by
!!                                 |ammonia volatilization in watershed
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    akn         |
!!    akv         |
!!    cecf        |none          |volatilization CEC factor
!!    dmidl       |
!!    dpf         |
!!    j           |none          |HRU number
!!    k           |none          |counter (soil layer)
!!    rnit        |kg N/ha       |amount of nitrogen moving from the NH3 to the
!!                               |NO3 pool (nitrification) in the layer
!!    rnv         |
!!    rvol        |kg N/ha       |amount of nitrogen lost from the NH3 pool due
!!                               |to volatilization
!!    sw25        |
!!    swf         |
!!    swwp        |
!!    tf          |
!!    xx          |
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Exp, Max

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
      use parm

      integer :: j, k
      real :: sw25, swwp, swf, xx, dmidl, dpf, akn, akv, rnv, rnit, rvol
      real :: tf
      real :: cecf = 0.15

      j = 0
      j = ihru

      do k = 1, sol_nly(j)
        tf = 0.
        tf = .41 * (sol_tmp(k,j) - 5.) / 10.

        if (sol_nh3(k,j) > 0. .and. tf >= 0.001) then
          sw25 = 0.
          swwp = 0.
          sw25 = sol_wpmm(k,j) + 0.25 * sol_fc(k,j)
          swwp = sol_wpmm(k,j) + sol_st(k,j)
          if (swwp < sw25) then
            swf = 0.
            swf = (swwp - sol_wpmm(k,j)) /(sw25 - sol_wpmm(k,j))
          else
            swf = 1.
          endif

          if (k == 1) then
            xx = 0.
          else
            xx = 0.
            xx = sol_z(k-1,j)
          endif

          dmidl = 0.
          dpf = 0.
          akn = 0.
          akv = 0.
          rnv = 0.
          rnit = 0.
          rvol = 0.
          dmidl = (sol_z(k,j) + xx) / 2.
          dpf = 1. - dmidl / (dmidl + Exp(4.706 - .0305 * dmidl))
          akn = tf * swf
          akv = tf * dpf * cecf
          rnv = sol_nh3(k,j) * (1. - Exp(-akn - akv))
          rnit = 1. - Exp(-akn)
          rvol = 1. - Exp(-akv)

          !! calculate nitrification (NH3 => NO3)
	    !! apply septic algorithm only to active septic systems
          if(k/=i_sep(j).or.isep_opt(j)/= 1) then  ! J.Jeong for septic, biozone layer
             if (rvol + rnit > 1.e-6) then
               rvol = rnv * rvol / (rvol + rnit)
               rnit = rnv - rvol
               if (rnit < 0.) rnit = 0.
               sol_nh3(k,j) = Max(1.e-6, sol_nh3(k,j) - rnit)
             endif
             if (sol_nh3(k,j) < 0.) then
               rnit = rnit + sol_nh3(k,j)
               sol_nh3(k,j) = 0.
             endif
             sol_no3(k,j) = sol_no3(k,j) + rnit

             !! calculate ammonia volatilization
             sol_nh3(k,j) = Max(1.e-6, sol_nh3(k,j) - rvol)
             if (sol_nh3(k,j) < 0.) then
               rvol = rvol + sol_nh3(k,j)
               sol_nh3(k,j) = 0.
             endif

             !! summary calculations
             if (curyr > nyskip) then
               wshd_voln = wshd_voln + rvol * hru_dafr(j)
               wshd_nitn = wshd_nitn + rnit * hru_dafr(j)
             end if
          end if
        end if

      end do

      return
      end