text:
MODULE forpy
!1-ExpandArrayR
!2-ExpandArrayI
!3-AppendR
!4-AppendI
!5-AppendInverseR
!6-AppendInverseI
!7-SetI
!8-SetF
!9-PrintArrayI
!10-PrintArrayR
!11-FindInArray
!12-FindInArrayF (SetF, should be the general one )
!13-FindFInArrayF (0)
!14-Sorting1Column
!15-SORTING
!16-int2str
!17-Interp
!18-sol2eq
!19-PolygonGeometry
!20-ReadASCII
!21-WriteASCII
IMPLICIT NONE

CONTAINS
    Function Round(x, precision)
		real :: x, defualt_precision, Round
		real, optional :: precision

		if (.not. present(precision)) then
			defualt_precision = 100.0
		else
			defualt_precision = precision
		end if
		Round = INT(x * defualt_precision) / defualt_precision
		return
	end function Round

	SUBROUTINE ExpandArrayR(OldArray, Dim)
	!------------------------------------------------------------------
	! ExpandArrayR(OldArray, Dim)
	!------------------------------------------------------------------
	! ExpandArrayR subroutine expand the dimension of a given 2D array
	! by adding one row at the bottom or one column at the right of the
	! given array
	!	Inpus:
	!		1- OldArray:[2D Real]
	!			1D array you want to append a value to it
	!		3-Dim:[Integer]
	!			1 to append the NewValue after the last row, 2 to append the the NewValue after
	!			the last column
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: rowsold, colsold, Dim
		Real, POINTER :: OldArray(:,:), NewArray(:,:)

		rowsold = size(OldArray,1)
		colsold = size(OldArray,2)
		if (Dim == 1) then
			allocate(NewArray(rowsold+1,colsold))
		else
			allocate(NewArray(rowsold,colsold+1))
		end if

		NewArray = 0

		if (Dim == 1) then
			if (rowsold == 1) then
				NewArray(1,:) = OldArray(1,:)
			else
				NewArray(1:rowsold,:) = OldArray(:,:)
			end if
		else
			if (colsold == 1) then
				NewArray(:,1) = OldArray(:,1)
			else
				NewArray(:,1:colsold) = OldArray(:,:)
			end if
		end if

		deallocate(OldArray)

		if (Dim == 1) then
			allocate(NewArray(rowsold+1,colsold))
		else
			allocate(NewArray(rowsold,colsold+1))
		end if

		OldArray = NewArray

		RETURN
	ENDSUBROUTINE ExpandArrayR

	SUBROUTINE Expand1DArrayR(OldArray)
	!------------------------------------------------------------------
	! ExpandArrayR(OldArray, Dim)
	!------------------------------------------------------------------
	! ExpandArrayR subroutine expand the dimension of a given 2D array
	! by adding one row at the bottom or one column at the right of the
	! given array
	!	Inpus:
	!		1- OldArray:[2D Real]
	!			1D array you want to append a value to it
	!		3-Dim:[Integer]
	!			1 to append the NewValue after the last row, 2 to append the the NewValue after
	!			the last column
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: rowsold
		Real, POINTER :: OldArray(:), NewArray(:)

		rowsold = size(OldArray,1)

		allocate(NewArray(rowsold+1))

		NewArray = -123456789

		if (rowsold == 1) then
			NewArray(1) = OldArray(1)
		else
			NewArray(1:rowsold) = OldArray(:)
		end if

		deallocate(OldArray)

		allocate(OldArray(rowsold+1))

		OldArray = NewArray

		RETURN
	ENDSUBROUTINE Expand1DArrayR

	SUBROUTINE Expand1DArrayI(OldArray)
	!------------------------------------------------------------------
	! ExpandArrayR(OldArray, Dim)
	!------------------------------------------------------------------
	! ExpandArrayR subroutine expand the dimension of a given 2D array
	! by adding one row at the bottom or one column at the right of the
	! given array
	!	Inpus:
	!		1- OldArray:[2D Real]
	!			1D array you want to append a value to it
	!		3-Dim:[Integer]
	!			1 to append the NewValue after the last row, 2 to append the the NewValue after
	!			the last column
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: rowsold
		integer, POINTER :: OldArray(:), NewArray(:)

		rowsold = size(OldArray,1)

		allocate(NewArray(rowsold+1))

		NewArray = -123456789

		if (rowsold == 1) then
			NewArray(1) = OldArray(1)
		else
			NewArray(1:rowsold) = OldArray(:)
		end if

		deallocate(OldArray)

		allocate(OldArray(rowsold+1))

		OldArray = NewArray

		RETURN
	ENDSUBROUTINE Expand1DArrayI

	SUBROUTINE ExpandArrayI(OldArray, Dim)
	!------------------------------------------------------------------
	! ExpandArrayI(OldArray, Dim)
	!------------------------------------------------------------------
	! ExpandArrayR subroutine expand the dimension of a given 2D array
	! by adding one row at the bottom or one column at the right of the
	! given array
	!	Inpus:
	!		1- OldArray:[2D Integer]
	!			1D array you want to append a value to it
	!		3-Dim:[Integer]
	!			1 to append the NewValue after the last row, 2 to append the the NewValue after
	!			the last column
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: rowsold, colsold, Dim
		INTEGER, POINTER :: OldArray(:,:), NewArray(:,:)

		rowsold = size(OldArray,1)
		colsold = size(OldArray,2)
		if (Dim == 1) then
			allocate(NewArray(rowsold+1,colsold))
		else
			allocate(NewArray(rowsold,colsold+1))
		end if

		NewArray = 0

		if (Dim == 1) then
			if (rowsold == 1) then
				NewArray(1,:) = OldArray(1,:)
			else
				NewArray(1:rowsold,:) = OldArray(:,:)
			end if
		else
			if (colsold == 1) then
				NewArray(:,1) = OldArray(:,1)
			else
				NewArray(:,1:colsold) = OldArray(:,:)
			end if
		end if

		deallocate(OldArray)

		if (Dim == 1) then
			allocate(OldArray(rowsold+1,colsold))
		else
			allocate(OldArray(rowsold,colsold+1))
		end if

		OldArray = NewArray

		RETURN
	ENDSUBROUTINE ExpandArrayI

	SUBROUTINE AppendR(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendR(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendR subroutine append a new value at the end of an array
	!
	!	Inpus:
	!		1- OldArray:[Real]
	!			1D array you want to append a value to it
	!		2-NewValue:[real]
	!			real value to append to the array, if you want to append a value of a different
	!			type(integer) conver the value using the intrinsic functions int, real
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: oldsize
		Real :: NewValue
		Real, POINTER :: OldArray(:), NewArray(:)

		oldsize = size(OldArray,1)
		allocate(NewArray(oldsize+1))

		NewArray = 0
		if (oldsize == 1) then
			NewArray(1) = OldArray(oldsize)
			NewArray(2) = NewValue
		else
			NewArray(1:oldsize) = OldArray
			NewArray(oldsize+1) = NewValue
		end if

		deallocate(OldArray)
		allocate(OldArray(oldsize+1))
		OldArray = NewArray

		RETURN
	ENDSUBROUTINE AppendR

	SUBROUTINE AppendI(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendI(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendI subroutine append a new value at the end of an array
	!
	!	Inpus:
	!		1- OldArray:[Integer]
	!			1D array you want to append a value to it
	!		2-NewValue:[real]
	!			real value to append to the array, if you want to append a value of a different
	!			type(integer) conver the value using the intrinsic functions int, real
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: oldsize
		Real :: NewValue
		Integer, POINTER :: OldArray(:), NewArray(:)

		oldsize = size(OldArray,1)
		allocate(NewArray(oldsize+1))

		if (oldsize == 1) then
			NewArray(1) = OldArray(oldsize)
			NewArray(2) = NewValue
		else
			NewArray(1:oldsize) = OldArray
			NewArray(oldsize+1) = NewValue
		end if

		deallocate(OldArray)
		allocate(OldArray(oldsize+1))
		OldArray = NewArray

		RETURN
	ENDSUBROUTINE AppendI

	SUBROUTINE AppendInverseR(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendInverseR(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendInverseR subroutine append a new value at the top of an array
	!
	!	Inpus:
	!		1- OldArray:[Real]
	!			1D array you want to append a value to it
	!		2-NewValue:[real]
	!			real value to append to the array, if you want to append a value of a different
	!			type(integer) conver the value using the intrinsic functions int, real
	!	Example:
	!
	!------------------------------------------------------------------
		INTEGER :: oldsize
		Real :: NewValue
		Real, POINTER :: OldArray(:), NewArray(:)

		! get the size of the array and assign the new array with +1 the old size
		oldsize = size(OldArray,1)
		allocate(NewArray(oldsize+1))

		NewArray = 0
		if (oldsize == 1) then
			NewArray(2) = OldArray(oldsize)
			NewArray(1) = NewValue
		else
			NewArray(2:oldsize+1) = OldArray
			NewArray(1) = NewValue
		end if
		! reallocate the old array and assign it with the new array
		deallocate(OldArray)
		allocate(OldArray(oldsize+1))
		OldArray = NewArray
		RETURN
	ENDSUBROUTINE AppendInverseR

	SUBROUTINE AppendInverseI(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendInverseI(OldArray, NewValue)
	!------------------------------------------------------------------
	! AppendInverseR subroutine append a new value at the top of an array
	!
	!	Inpus:
	!		1- OldArray:[Integer]
	!			1D array you want to append a value to it
	!		2-NewValue:[real]
	!			real value to append to the array, if you want to append a value of a different
	!			type(integer) conver the value using the intrinsic functions int, real
	!	Example:
	!
	!------------------------------------------------------------------

		INTEGER :: oldsize
		Real :: NewValue
		Integer, POINTER :: OldArray(:), NewArray(:)

		oldsize = size(OldArray,1)
		allocate(NewArray(oldsize+1))

		if (oldsize == 1) then
			NewArray(2) = OldArray(oldsize)
			NewArray(1) = NewValue
		else
			NewArray(2:oldsize+1) = OldArray
			NewArray(1) = NewValue
		end if
		deallocate(OldArray)
		allocate(OldArray(oldsize+1))
		OldArray = NewArray

		RETURN
	ENDSUBROUTINE AppendInverseI

	SUBROUTINE SetI(Array1, Array2)
	!------------------------------------------------------------------
	! SetI(Array1, Array2)
	!------------------------------------------------------------------
	! SetI subroutine takes an array full of values that might be repeated
	! and return an array of unique values (delete repeated values)
	!	Inpus:
	!		1-Array1:[Integer]
	!
	!	Example:
	!		INTEGER , POINTER :: Array1(:), Array2(:)
	!		allocate(Array1(7),Array2(7))
	!		Array1 = (/1,2,3,4,5,6,6/)
	!		call SetI(Array1, Array2)

		Integer, pointer :: Array1(:), Array2(:)
		Integer :: i, loc, CurrentValue, N

		N = 1
		allocate(Array2(1))
		Array2(1) = Array1(1)

		do i=2, size(Array1,1)
			CurrentValue = Array1(i)
			! search for the value i array2
			call FindInArrayF(real(CurrentValue),real(Array2),loc)
			! if you did not find it add it to the new array
			if (loc == -1) then
				call Expand1DArrayI(Array2)
				N = N + 1
				Array2(N) = CurrentValue
			end if
		end do

	END SUBROUTINE SetI

	SUBROUTINE SetF(Array1, Array2)
	!------------------------------------------------------------------
	! SetI(Array1, Array2)
	!------------------------------------------------------------------
	! SetI subroutine takes an array full of values that might be repeated
	! and return an array of unique values (delete repeated values)
	!	Inpus:
	!		1-Array1:[Integer]
	!
	!	Example:
	!		Real , POINTER :: Array1(:), Array2(:)
	!		allocate(Array1(7),Array2(7))
	!		Array1 = (/1,2,3,4,5,6,6/)
	!		call SetF(Array1, Array2)
	implicit none
		Real, pointer :: Array1(:), Array2(:)
		real :: CurrentValue
		Integer :: i, loc, N
		N = 1

		allocate(Array2(1))
		Array2(1) = Array1(1)

		do i=2, size(Array1,1)
			CurrentValue = Array1(i)
			! search for the value i array2
			call FindInArrayF(real(CurrentValue),Array2,loc)
			! if you did not find it add it to the new array
			if (loc == -1) then
				call Expand1DArrayR(Array2)
				N = N + 1
				Array2(N) = CurrentValue
			end if
		end do

	END SUBROUTINE SetF

	SUBROUTINE FindDuplicate()

		implicit none
		Real, pointer :: Array1(:), Array2(:)
		real :: CurrentValue
		Integer :: i, loc, N
		N = 1

		allocate(Array2(1))
		Array2(1) = Array1(1)

		do i=2, size(Array1,1)
			CurrentValue = Array1(i)
			! search for the value i array2
			call FindInArrayF(real(CurrentValue),Array2,loc)
			! if you did not find it add it to the new array
			if (loc == -1) then
				call Expand1DArrayR(Array2)
				N = N + 1
				Array2(N) = CurrentValue
			end if
		end do

	END SUBROUTINE FindDuplicate

	SUBROUTINE PrintArrayI(array, rows, columns, fname,decision, formatt)
	!------------------------------------------------------------------
	!! PrintArray(array, rows, columns, fileU,decision, formatt)
	!! PrintArray is used to print a 2D array on the screen or in a file
	!! if decision is 0 the array will be printed on screen if not array
	!! will be printed on file
	!!
	!! iNPUTS:
	!!		1-array:[integer] 2D
	!!				dimensions of the array  should be (rows, columns)
	!!		2-rows:[integer]
	!!				number of rows (size of the array)
	!!		3-columns:[integer]
	!!				number of columns (size of the array)
	!!		4-fileU:[integer]
	!!				file unit to open (not yet)
	!!		5-decision:[integer]
	!!				0 to write on screen or 1 to write on a file
	!!		6-formatt:[character]
	!!				format of how things are going to be written in the file
	!!Outputs:
	!!		1- content of the array are going to be eitherw ritten in the
	!!			screen or in a file
	!! Examples:
	!!		! free format
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,0, "")
	!!		! on the screen formatted
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,0,"(i3,2x,i6,2x,i6)")
	!!		! on file free formatted
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,1,'' )
	!!		! on file formatted
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,1,"(i3,2x,i6,2x,i6)" )
	!------------------------------------------------------------------
	! DECLARATION
	CHARACTER(*) :: formatt, fname
	! 1D
	INTEGER :: rows, columns,i, j, decision !,fileU
	! 2D
	INTEGER, DIMENSION(rows,columns) :: array

	! EXECUTION
	! check if decision 0 on screen and format length more than 1 write formatted
	if (decision == 0 .and. len(trim(formatt)) == 0 ) then
	! print on the screen free format
		do i = 1,rows
			write(*,*) (array(i,j), j = 1,columns)
		end do
	else if(decision == 0 .and. len(trim(formatt)) > 0 ) then
	! print on the screen formatted
		do i = 1,rows
			write(*,formatt) (array(i,j), j = 1,columns)
		end do
	else if (decision == 1 .and. len(trim(formatt)) == 0 ) then
	! print on file free formatted
		open(90000, file=fname, status='unknown')
		do i = 1,rows
			write(90000,*) (array(i,j), j = 1,columns)
		end do
	else
		! print on file free formatted
		open(90000, file=fname, status='unknown')
		do i = 1,rows
			write(90000,formatt) (array(i,j), j = 1,columns)
		end do
	end if
	close(90000)

	END SUBROUTINE

	SUBROUTINE PrintArrayR(array, rows, columns, fname,decision, formatt)
	!------------------------------------------------------------------
	!! PrintArray(array, rows, columns, fileU,decision, formatt)
	!! PrintArray is used to print a 2D array on the screen or in a file
	!! if decision is 0 the array will be printed on screen if not array
	!! will be printed on file
	!!
	!! iNPUTS:
	!! 		1-array:[Real] 2D
	!!				dimensions of the array  should be (rows, columns)
	!!		2-rows:[integer]
	!!				number of rows (size of the array)
	!!		3-columns:[integer]
	!!				number of columns (size of the array)
	!!		4-fileU:[integer]
	!!				file unit to open (not yet)
	!!		5-decision:[integer]
	!!				0 to write on screen or 1 to write on a file
	!!		6-formatt:[character]
	!!				format of how things are going to be written in the file
	!! Outputs:
	!!		1- content of the array are going to be eitherw ritten in the
	!!			screen or in a file
	!! Examples:
	!!		! free format
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,0, "")
	!!		! on the screen formatted
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,0,"(i3,2x,i6,2x,i6)")
	!!		! on file free formatted
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,1,'' )
	!!		! on file formatted
	!!		CALL PrintArrayI(info(1:3,1:3), 3, 3, 900,1,"(i3,2x,i6,2x,i6)" )
	!------------------------------------------------------------------

	! DECLARATION
	CHARACTER(*) :: formatt, fname
	! 1D
	INTEGER :: rows, columns,i, j, decision !fileU,
	! 2D
	REAL, DIMENSION(rows,columns) :: array

	! EXECUTION
	 ! check if decision 0 on screen and format length more than 1 write formatted
	if (decision == 0 .and. len(trim(formatt)) == 0 ) then
	! print on the screen free format
		do i = 1,rows
			write(*,*) (array(i,j), j = 1,columns)
		end do
	else if(decision == 0 .and. len(trim(formatt)) > 0 ) then
	! print on the screen formatted
		do i = 1,rows
			write(*,formatt) (array(i,j), j = 1,columns)
		end do
	else if (decision == 1 .and. len(trim(formatt)) == 0 ) then
		! print on file free formatted
		open(90000, file=fname, status='unknown')
		do i = 1,rows
			write(90000,*) (array(i,j), j = 1,columns)
		end do
	else
		! print on file free formatted
		open(90000, file=fname, status='unknown')
		do i = 1,rows
			write(90000,formatt) (array(i,j), j = 1,columns)
		end do
	end if
	close(90000)

	END SUBROUTINE PrintArrayR

	Subroutine FindInArray(TargetValue,Arraye,loc)
	!-------------------------------------------------------------------------------------------------
	! Subroutine FindInArray(TargetValue,Arraye,loc)
	! ----------------------------------------------------
	! This function find the location of an integer input value in an integer 1D array and if
	! the value does not exist in the array a value of -1 will be returned
	! as a location
	!
	!	1-Argument   : DESCRIPTION (UNIT)
	!		1-TargetValue:
	!			[Integer] the value you want to locate
	!
	!		2-Arraye:
	!			[Integer] 1D array that contains values

	!	2-Output:
	!		1- loc:
	!			[integer] the index of the row that the target value
	!			locates in the array
	!
	!	3- Examples:
	!			allocate(a(50))
	!			a=(/(i,i=1,50)/)
	!			call FindInArray(20,a,loc)
	!-------------------------------------------------------------------------------------------------
	! DECLARATION
		INTEGER :: TargetValue, i, loc
		INTEGER, INTENT(in) :: Arraye (:)

	! Execution
		loc = -1
		do i = 1, size(Arraye,1)
			if (Arraye(i) == TargetValue) then
				loc = i
				exit
			endif
		end do

	END SUBROUTINE FindInArray

	Subroutine FindInArrayF(TargetValue, Arr, loc, precision)
	!------------------------------------------------------------------------------------------------------
	! Subroutine FindInArray(TargetValue,Array,loc)
	! ----------------------------------------------------
	! This function find the location of an Rael input value in a Real 1D array and if
	! the value does not exist in the array a value of -1 will be returned
	! as a location, to use the same function to search for an integer value use the real
	! intrinsic function to convert the TargetValue into real
	! FindInArrayF(int(TargetValue),Arr,loc)
	!
	!	1-Argument   : DESCRIPTION (UNIT)
	!		1-TargetValue:
	!			[Integer] the value you want to locate
	!
	!		2-Array:
	!			[Real] 1D array that contains values
	!
	!	2-Output:
	!		1- loc:
	!			[integer] the index of the row that the target value
	!			locates in the array
	!
	!	3- Examples:
	!			allocate(a(50))
	!			a=(/(i,i=1,50)/)
	!			call FindInArray(20,a,loc)
	!------------------------------------------------------------------------------------------------------
	! DECLARATION
		!INTEGER :: TargetValue, i, loc
		INTEGER :: i, loc
		Real :: TargetValue, val1, val2
		Real, INTENT(in) :: Arr (:)
		real, optional :: precision

		loc = -1
		if (.not. present(precision)) then
			do i = 1, size(Arr,1)
				val1 = Round(Arr(i))
				val2 = Round(TargetValue)
				if (val1 == val2) then
					loc = i
					exit
				endif
			end do
		else
			do i = 1, size(Arr,1)
				val1 = Round(Arr(i),precision)
				val2 = Round(TargetValue,precision)
				if ( val1 == val2) then
					loc = i
					exit
				endif
			end do
		end if
	END SUBROUTINE FindInArrayF

	Subroutine FindFInArrayF(TargetValue,Arr,loc)
	!------------------------------------------------------------------------------------------------------
	! Subroutine FindInArray(TargetValue,Array,loc)
	! ----------------------------------------------------
	! This function find the location of a the closest number in real 1D array to Real input value
	! the value does not exist in the array a value of -1 will be returned
	! as a location
	!
	!	1-Argument   : DESCRIPTION (UNIT)
	!		1-TargetValue:
	!			[Real] the value you want to locate
	!
	!		2-Array:
	!			[Real] 1D array that contains values

	!	2-Output:
	!		1- loc:
	!			[integer] the index of the row that the target value
	!			locates in the array
	!
	!	3- Examples:
	!			allocate(a(50))
	!			a=(/(i,i=1,50)/)
	!			call FindInArray(20,a,loc)
	!------------------------------------------------------------------------------------------------------
	! DECLARATION
		REAL :: TargetValue
		INTEGER :: i, loc
		Real, INTENT(in) :: Arr (:)
		REAL,DIMENSION(:), ALLOCATABLE :: diff
	! Execution
		ALLOCATE(diff(size(Arr,1)))
		loc = -1
		do i = 1, size(Arr,1)
			diff(i) = abs( Arr(i) - TargetValue)
		end do

		loc = minloc(diff,1)

	END SUBROUTINE FindFInArrayF

	Subroutine FindIn2Darray(TargetValue, DirectUS, i)
		REAL :: TargetValue !, DS
		INTEGER :: i, loc
		!Real, INTENT(in) :: Arr (:)
		REAL, POINTER :: DirectUS(:,:)

		do i = 1, size(DirectUS,1)
			! search in each row if the TargetValue is found that is our row
			call FindInArrayF(real(TargetValue), DirectUS(i,:), loc)
			if (loc /= -1 ) then
				!DS = SubID(loc)
				EXIT
			end if
		end do

	END SUBROUTINE

	SUBROUTINE Sort1D(arr,n)
	!---------------------------------------------------------------------------------
	!Sorting1Column(a,m,n,column)
	!	Inputs:
	!		1-arr:
	!			[array] 2D array
	!		2-m:
	!			[INTEGER] number of rows in the array
	!		3-n:
	!			[INTEGER] number of columns in the array
	!		4-column:
	!			[INTEGER] the column you want to sort the array with respect to
	!
	!		Outputs:
	!			1- array:
	!				[array] sorted array
	!---------------------------------------------------------------------------------
		IMPLICIT NONE
		INTEGER :: row, swap_index, n !nrows,
		integer :: arr(n), tmp

		DO row = 1, n - 1
			swap_index = MINLOC( arr(row:n) , dim=1 ) + row - 1
			! Swap_index
			tmp = arr(row)
			! Swap_index
			arr(row) = arr(swap_index)
			arr(swap_index) = tmp
		END DO

	END SUBROUTINE Sort1D

	SUBROUTINE Sorting1Column(arr,m,n,column)
	!---------------------------------------------------------------------------------
	!Sorting1Column(a,m,n,column)
	!	Inputs:
	!		1-arr:
	!			[array] 2D array
	!		2-m:
	!			[INTEGER] number of rows in the array
	!		3-n:
	!			[INTEGER] number of columns in the array
	!		4-column:
	!			[INTEGER] the column you want to sort the array with respect to
	!
	!		Outputs:
	!			1- array:
	!				[array] sorted array
	!---------------------------------------------------------------------------------
		IMPLICIT NONE
		INTEGER :: column, row, nrows, swap_index, m, n
		REAL :: arr(m,n), tmp(n)
		nrows = size(arr,1 )

		! sort based on the first column
		DO row = 1, nrows-1
			! first column
			swap_index = MINLOC( arr(row:nrows,column) , dim=1 ) + row - 1
			! Swap_index
			tmp = arr(row,:)
			! Swap_index
			arr(row,:) = arr(swap_index,:)
			arr(swap_index,:) = tmp
		END DO

	END SUBROUTINE Sorting1Column

	SUBROUTINE SORTING(array,m,n,column1,column2)
	!---------------------------------------------------------------------------------
	!SORTING(a,m,n,column)
	!	Inputs:
	!		1-array:
	!			[REAL] 2D array
	!		2-m:
	!			[INTEGER] number of rows in the array
	!		3-n:
	!			[INTEGER] number of columns in the array
	!		4-column:
	!			[INTEGER] the column you want to sort the array with respect to
	!
	!		Outputs:
	!			1- array:
	!				[array] sorted array
	!---------------------------------------------------------------------------------

	! DECLARATION
		IMPLICIT NONE
		INTEGER :: column1, column2, row, nrows, m, n, ro, &
					no_unrepeated, start, endd
			REAL :: array(m,n)
		REAL, DIMENSION(:,:), ALLOCATABLE :: unrepeated, smaller_array!, second_sort
		REAL, DIMENSION(:), ALLOCATABLE :: dummy

	! EXECUTION
		! sort based on the first column
		nrows = SIZE( array,1 )
		!write(*,*) 'sorted column = ', column1
		CALL Sorting1Column(array,m,n,column1)
			!write(*,*) 'sorted 1D---------------------------------------------------------------------------------'

		!do row=1,m
		!	write(*,'(i4,f8.1,2x,f8.1)') row,array(row,:)
			!end do
		!write(*,*) 'done'
		! sort based on the second column
		allocate(dummy(m))
		ro = 2
		dummy = -1
		dummy(1)=array(1,column1)

			do row = 2, nrows
				if (array(row-1,column1) /= array(row,column1)) then
				dummy(ro) = array(row,column1)
					ro = ro + 1
				end if
		end do

			no_unrepeated = size(pack(dummy,dummy /= -1.0),1)

			!write(*,*) dummy
		!write(*,*) 'no_unrepeated',no_unrepeated

			allocate(unrepeated(no_unrepeated,2))

			unrepeated = 0
		! first column will be the number from the sorted column
		! the second column will be the number of rows that has the same
		! number in the first column
		unrepeated(:,1) = dummy(1:no_unrepeated)

		!write(*,*)unrepeated(:,1)

			deallocate (dummy)
		!write(*,*) unrepeated(:,1)
		ro = 1
		! initialize the first repeated number by 1
		unrepeated(ro,2) = 1
		!counter = 0

		do row = 2, nrows
		  ! if the number is different
			if (array(row-1,column1) /= array(row,column1)) then
				! go to the next row
			  ro = ro +1
				! initialize no of repeated number in the second column with 1
			  unrepeated(ro,2) = 1
			else
			  ! if the number is the same
			  ! increase the number in the second column with 1
				unrepeated(ro,2) = unrepeated(ro,2) +1
				end if
		end do
			start = 1

			do row = 1, no_unrepeated
			if (row == 1) then
				endd = unrepeated(row,2)
					allocate (smaller_array(start:endd,n))
					smaller_array(:,:) = array(start:endd,1:n)

				CALL Sorting1Column(smaller_array,size(smaller_array,1),size(smaller_array,2),column2)

					array(start:endd,1:n) = smaller_array
				!write(*,*) smaller_array!(:,1:n)
				!write(*,*)'-----------'
			else
				start = start + unrepeated(row-1,2)
				endd = start + unrepeated(row,2)-1
					allocate (smaller_array(start:endd,n))
					smaller_array(:,:) = array(start:endd,1:n)

				CALL Sorting1Column(smaller_array,size(smaller_array,1),size(smaller_array,2),column2)

					array(start:endd,1:n) = smaller_array
				!write(*,*) smaller_array!(:,1:n)
				!write(*,*)'-----------'
				end if

				deallocate(smaller_array)
		end do

	END SUBROUTINE SORTING

	FUNCTION int2str(num)
	!------------------------------------------------------------------------------------------------------
	! int2str(int)
	! ----------------------------------------------------
	! This function converts the integer number into a string type
	!
	!	1-Argument:
	! 		1-int:
	!			[integer] number of type integer
	!
	!	2-Output:
	!		1- int2str:
	!			[character] the string form of the inputed integer number
	!
	!	3- Examples:
	!				num_str=int2str(100)
	!------------------------------------------------------------------------------------------------------

	INTEGER, INTENT(IN) :: num
	CHARACTER*6 :: int2str,temp
	CHARACTER*6 :: fs
	! choosing the format of the number by checking the length
	if (num>=100000) then
		fs = '(I6)'
	elseif (num>=10000) then
		fs = '(I5)'
	elseif (num>=1000) then
		fs = '(I4)'
	elseif (num>=100) then
		fs = '(I3)'
	elseif (num>=10) then
		fs = '(I2)'
	else
		fs = '(I1)'
	end if

	write(temp,fs) num
	int2str = trim(temp)
	return
	END FUNCTION int2str

	FUNCTION str2int(string)
	!---------------------------------------------------------------------------
	! int2str(int)
	! ----------------------------------------------------
	! This function converts the integer number into a string type
	!
	!	1-Argument:
	! 		1-int:
	!			[integer] number of type integer
	!
	!	2-Output:
	!		1- int2str:
	!			[character] the string form of the inputed integer number
	!
	!	3- Examples:
	!				num_str=int2str(100)
	!---------------------------------------------------------------------------

	INTEGER :: str2int !intent(out)
	character(len=*), intent(in) :: string
	! choosing the format of the number by checking the length
	read(string,*) str2int

	return

	END FUNCTION str2int

	FUNCTION str2real(string)
		!---------------------------------------------------------------------------
		! int2str(int)
		! ----------------------------------------------------
		! This function converts the integer number into a string type
		!
		!	1-Argument:
		! 		1-int:
		!			[integer] number of type integer
		!
		!	2-Output:
		!		1- int2str:
		!			[character] the string form of the inputed integer number
		!
		!	3- Examples:
		!				num_str=int2str(100)
		!---------------------------------------------------------------------------

		real :: str2real !intent(out)
		character(len=*), intent(in) :: string
		! choosing the format of the number by checking the length
		read(string,*) str2real

		return

	END FUNCTION str2real

	SUBROUTINE Interp(InputVector,TargetTempRes,InputTempRes, InterpolatedVector)
	!------------------------------------------------------------------------------------------------------
	! Interp(InputVector,TargetTempRes,InputTempRes, InterpolatedVector)
	! ----------------------------------------------------
	! This function interpolate input time series from input temporal resolution (tsinput)
	! to another dt
	!
	!	1-Argument   : DESCRIPTION (UNIT)
	!		1-InputVector:
	!			[real] 1D array of values you want to interpolate
	!		2-dt:
	!			[integer] resulted time resolution (sec)
	!			 default is 60 sec
	!		3-tsInput:
	!			[integer] input time step (hrs)
	!			 default is 24 hrs
	!
	!	2-Output:
	!		1- hin:
	!			[real] 1D array of interpolated values(m)
	!
	!	3- Examples:
	!		HydrologicalTempRes = 24  ! (hourly)
	!		dt1D = 60		  !(min)
	!		! size(lflow(sSIM:eSIM,i))=508
	!		call Interp(lflow(sSIM:eSIM,i), dt1D, HydrologicalTempRes,  bnd_q)
	!------------------------------------------------------------------------------------------------------
	! DECLARATION
		REAL*4, DIMENSION(:), INTENT(IN) :: InputVector			! input InputVector
		INTEGER*4, INTENT(IN) :: TargetTempRes, InputTempRes	!time step of the input
		!REAL*4, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: hin !interpolated h
		REAL*4, POINTER :: InterpolatedVector(:) !interpolated h
		INTEGER :: i, j, correct

		correct = 24/InputTempRes							! 1
		! new size will be(508 day(size(InputVector)) * 24 *60)
		allocate(InterpolatedVector(86400/correct/TargetTempRes*(size(InputVector)-1)))		! (86400/1/60)*508= 731,520

		do i = 0, size(InputVector)-2 ! 0 to 507 days
			do j = 1, int(86400/correct/TargetTempRes) ! from 1 min to 1440 min
				InterpolatedVector(j+i*86400/correct/TargetTempRes) = InputVector(i+1) + &
                (InputVector(i+2) - InputVector(i+1))/(86400/correct)*j*TargetTempRes
			end do
		end do

		return
	END SUBROUTINE Interp

	FUNCTION sol2eq(a,b,c)
	!------------------------------------------------------------------------------------------------------
	! FUNCTION sol2eq(a,b,c)
	! ----------------------------------------------------
	! This function solves a quadratic equation at the form of aX^2+bX+c=0
	!
	!	1-Argument   : DESCRIPTION (UNIT)
	!		1-a:
	!			[real] coefficient of X^2
	!
	!		2-b:
	!			[real] coefficient of X
	!
	!		3-c:
	!			[real] Free term
	!
	!	2-Output:
	!		1- sol2eq:
	!			[real] only one solution of the equation (-b+sqrt(b^2-4ac))/2a
	!
	!	3- Examples:
	!		sol2eq(2,5,4)
	!------------------------------------------------------------------------------------------------------
	! DECLARATION
	!solution for the special case: a>0, delta>0, only the greater solution is needed
		REAL*4 :: sol2eq
		REAL*4, INTENT(IN) :: a, b, c
		REAL*4 :: delta
		delta = b**2 - 4*a*c
		sol2eq = (-b+sqrt(delta))/(2*a)
		!print*, sol2eq
		return
	END FUNCTION sol2eq

	FUNCTION PolygonGeometry(mp)
	!------------------------------------------------------------------------------------------------------
	! PolygonGeometry(mp)
	! ----------------------------------------------------
	! This function calculates the area and perimeter of the polygon using
	! the coordinates of the vortexes of polygon
	!
	!	1-Argument : DESCRIPTION (UNIT)
	! 		1-mp:
	!			[real] 2D array of containing the coordinates of the vortexes of
	!			the cross section with the x coord in the first column and the
	!			y coord in the second column
	!
	!	2-Output:
	!		1-PolygonGeometry:
	!			[real] 1D array contains the area and perimeter of the input polygon
	!
	!	3- Examples:
	!				mp(1,1:4)=(/0,1,5,6/)   ! x coordinates
	!				mp(2,1:4)=(/1,0,0,1/)   ! y coordinates
	!				prop=PolygonGeometry(mp)
	!				area=prop(1)
	!				peri=prop(2)
	!------------------------------------------------------------------------------------------------------

		!REAL*4, ALLOCATABLE, DIMENSION(:,:), INTENT(IN) :: mp
		REAL*4, POINTER :: mp(:,:)
		REAL*4, DIMENSION(2):: PolygonGeometry
		REAL*4 :: area, peri
		INTEGER :: i, j

		j = size(mp,1)
		area = 0.0
		peri = 0.0
		do i = 1, size(mp,1)-1
			area = area + mp(i,1)*mp(i+1,2) - mp(i+1,1)*mp(i,2)
			peri = peri + ((mp(i+1,1)-mp(i,1))**2 + (mp(i+1,2)-mp(i,2))**2)**0.5
		end do
		area = area + mp(j,1)*mp(1,2) - mp(1,1)*mp(j,2)
		area = area*0.5
		PolygonGeometry(:) = (/area, peri/)
		return
	END FUNCTION PolygonGeometry

	SUBROUTINE ReadASCII(Filename,Grid,Gridrows,Gridcols,LeftCornerX,Leftcornery,Nodata)
	!----------------------------------------------------------------------
	! populates the grid from the inputfile
	! clips the calculation grid from the base grid
	! based on the breach locations
	! and defines inflow cells
	!----------------------------------------------------------------
	!!   1-Argument   : DESCRIPTION (UNIT)
	!!
	!!		1-grd:
	!!			[Real] 2D array having the values of the grid you want to clip
	!!		2-grdrows:
	!!			[Integer] number of rows in the grid
	!!		3-grdcols:
	!!			[Integer] number of columns in the grid
	!!		4-filename:
	!!
	!!
	!!
	!!
	!----------------------------------------------------------------------
	CHARACTER*200		:: Filename
	DOUBLE PRECISION	:: Leftcornerx,Leftcornery,Nodata
	INTEGER				:: cellsize
	INTEGER				:: Gridrows,Gridcols
	REAL*4, Pointer		:: Grid(:,:)
	CHARACTER			:: dummy*15
	INTEGER				:: i!,j

	! open grid file
	open(9000,file=trim(Filename),action='read')
	read(9000,*) dummy,Gridcols
	read(9000,*) dummy,Gridrows
	read(9000,*) dummy,Leftcornerx
	read(9000,*) dummy,Leftcornery
	read(9000,*) dummy,cellsize
	read(9000,*) dummy,Nodata

	allocate(Grid(Gridrows,Gridcols))

	do i=1,Gridrows
		read(9000,*) Grid(i,1:Gridcols)
	end do

	close(9000)

	END SUBROUTINE ReadASCII

	SUBROUTINE WriteASCII(Filename,Grid,Gridrows,Gridcols,LeftCornerX,Leftcornery,cellsize,Nodata,formatt,&
				typeagree)
	!----------------------------------------------------------------------
	! populates the grid from the inputfile
	! clips the calculation grid from the base grid
	! based on the breach locations
	! and defines inflow cells
	!----------------------------------------------------------------
	!!	1-Argument   : DESCRIPTION (UNIT)
	!!		1-Filename:
	!!			[Real] 2D array having the values of the grid you want to clip
	!!		2-Grid:
	!!			[]
	!!		3-Gridrows:
	!!			[Integer] number of rows in the grid
	!!		4-Gridcols:
	!!			[Integer] number of columns in the grid
	!!		5-LeftCornerX:
	!!			[]
	!!		6-Leftcornery:
	!!			[]
	!!		7-cellsize:
	!!			[]
	!!		8-Nodata:
	!!			[]
	!!		9-formatt:
	!!			[String] string having the format of the number between two parentheses
	!!		10-typeagree:
	!!			[Integer] 1 if the array is float, 0 if the array is integer
	!----------------------------------------------------------------------
	CHARACTER*200		:: Filename,formatt
	DOUBLE PRECISION	:: Leftcornerx,Leftcornery,Nodata
	INTEGER				:: cellsize
	INTEGER				:: Gridrows,Gridcols
	REAL*4				:: Grid(Gridrows,Gridcols)
	!CHARACTER			:: dummy*15
	INTEGER				:: typeagree,i!,j

	! open grid file
	open(666,file=trim(Filename),action='write',STATUS="Unknown")
	write(666,'(A14,i6)') 'ncols         ', Gridcols
	write(666,'(A14,i6)') 'nrows         ', Gridrows
	write(666,'(A14,f16.4)') 'xllcorner     ', Leftcornerx
	write(666,'(A14,f16.4)') 'yllcorner     ', Leftcornery
	write(666,'(A14,i4)') 'cellsize      ', cellsize
	write(666,'(A14,i6)') 'NODATA_value  ', nint(Nodata)

	!allocate(Grid(Gridrows,Gridcols))
	if (typeagree ==1 ) then
	! if the type of the array is float as it is declared inside of the function
		do i=1,Gridrows
			!write(666,'('//trim(int2str(Gridcols))//'f6.2)') Grid(i,1:Gridcols)
			write(666,formatt) Grid(i,1:Gridcols)
		end do
	else
	! if the type of the array is integer in contrary to what is declared inside of the function
		do i=1,Gridrows
			write(666,formatt) int(Grid(i,1:Gridcols))
		end do
	end if

	close(666)

	END SUBROUTINE WriteASCII

	SUBROUTINE WriteData(wpath, step, SubID, BCWrite, QminWrite, HminWrite, bnd_q, bnd_h, RIMTempRes, &
						RRMTempRes, q_day, h_day)

		implicit none
		character*200 :: wpath, formattxt, dummy, fname
		integer ::BCWrite, QminWrite, HminWrite, step, SubID, RIMTempRes, RRMTempRes, m
		real*4, pointer :: bnd_h(:), bnd_q(:), h_day(:,:), q_day(:,:)

		if (BCWrite == 1) then
			open(10, file = trim(wpath)//'/results/USbnd/'//trim(int2str(SubID))//&
					'-'//trim(int2str(step))//'.txt',err=824)
			write(10,'(2F9.2)') (bnd_q(m), bnd_h(m), m= 1, size(bnd_q))
			close(10)
		end if
		! write the whole timeseries of Q and H----------------------------------------------------
		formattxt = "("//int2str(RIMTempRes*RRMTempRes)//'(F8.3,2x)'//")"
		dummy = trim(wpath)//'\results\'

		if (QminWrite == 1) then
			fname = trim(dummy)//'q\'//trim(int2str(int(SubID)))//'-q-'//trim(int2str(step))//'.txt'
			call PrintArrayR(q_day, size(q_day,1), size(q_day,2),fname,1, formattxt)
		end if

		if (HminWrite == 1) then
			fname = trim(dummy)//'h\'//trim(int2str(int(SubID)))//'-h-'//trim(int2str(step))//'.txt'
			call PrintArrayR(h_day, size(h_day,1), size(h_day,2),fname,1, formattxt)
		end if

		go to 1000
		824 write(*,*) 'file = '//'/USbnd/'//trim(int2str(SubID))//'-'//trim(int2str(step))//'.txt','  failed to open'

	1000 END SUBROUTINE WriteData

END MODULE forpy

max_stars_repo_name: MAfarrag/HydroFortran

max_issues_repo_name: MAfarrag/HydroFortran

max_forks_repo_name: MAfarrag/HydroFortran


text:
!
! $Id$
!
! Copyright 1989-2014 MINES ParisTech
!
! This file is part of PIPS.
!
! PIPS is free software: you can redistribute it and/or modify it
! under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! any later version.
!
! PIPS is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE.
!
! See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with PIPS.  If not, see <http://www.gnu.org/licenses/>.
!

!
! the main for the HPFC program (for both host and node)
!
      program MAIN
      include "hpfc_commons.h"
      call HPFC INIT MAIN
      if (MY TID.eq.HOST TID) then
       call HOST
      else
       call NODE
      end if
      end
!
! end of it
!

ext: f


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20100717 11:24:22 nbsp  I hope you get your bracelet back.  One of the most upsetting losses of an item Ive ever had is when my horned hat was stolen from a rather expensive restaurants coat room.  Like you, I had made the horns from sculpey, and it really upset me that I had put the time into something uniquely from and for me that somebody else had just grabbed and taken for their own.   Users/JabberWokky Evan JabberWokky Edwards

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!==============================================================================!
  subroutine Allocate_Memory(grid)
!------------------------------------------------------------------------------!
!   Allocates memory for the convertor.                                        !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Grid_Mod, only: Grid_Type,                &
                      Grid_Mod_Allocate_Nodes,  &
                      Grid_Mod_Allocate_Cells,  &
                      Grid_Mod_Allocate_Faces,  &
                      Grid_Mod_Allocate_New_Numbers
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!==============================================================================!

  ! Allocate memory
  ! =--> carefull: there is no checking!
  call Grid_Mod_Allocate_Nodes(grid, grid % n_nodes)
  call Grid_Mod_Allocate_Cells(grid, grid % n_cells,   grid % n_bnd_cells)
  call Grid_Mod_Allocate_Faces(grid, grid % n_cells*5)
  call Grid_Mod_Allocate_New_Numbers(grid,                &
                                     grid % n_nodes,      &
                                     grid % n_bnd_cells,  &
                                     grid % n_cells,      &
                                     grid % n_faces)

  allocate(grid % bnd_cond % color(-grid % n_bnd_cells:-1))
  grid % bnd_cond % color=0

  grid % n_copy = grid % n_faces  ! I believe it is n_cells * 5 at this point
  allocate(grid % bnd_cond % copy_c( -grid % n_bnd_cells:-1))
  grid % bnd_cond % copy_c = 0
  allocate(grid % bnd_cond % copy_s(2,grid % n_copy))
  grid % bnd_cond % copy_s=0

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
!----------------------------------------------------------------------------
! This file is part of UCLALES.
!
! UCLALES is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! UCLALES is distributsed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1999-2007, Bjorn B. Stevens, Dep't Atmos and Ocean Sci, UCLA
!----------------------------------------------------------------------------
!
module stat

  use ncio, only : open_nc, define_nc
  use grid, only : level
  use util, only : get_avg, get_cor, get_avg3, get_cor3, get_var3, get_csum
  implicit none
  private

  integer, parameter :: nvar1 = 27, nvar2 = 92

  integer, save      :: nrec1, nrec2, ncid1, ncid2, nv1=nvar1, nv2=nvar2
  real, save         :: fsttm, lsttm, nsmp = 0

  logical            :: sflg = .false.
  real               :: ssam_intvl = 30.   ! statistical sampling interval
  real               :: savg_intvl = 1800. ! statistical averaging interval

  character (len=7), save :: s1(nvar1)=(/                           &
       'time   ','cfl    ','maxdiv ','zi1_bar','zi2_bar','zi3_bar', & ! 1
       'vtke   ','sfcbflx','wmax   ','tsrf   ','ustar  ','shf_bar', & ! 7
       'lhf_bar','zi_bar ','lwp_bar','lwp_var','zc     ','zb     ', & !13
       'cfrac  ','lmax   ','albedo ','rwp_bar','prcp   ','pfrac  ', & !19
       'CCN    ','nrain  ','nrcnt  '/),                             &
       s2(nvar2)=(/                                                 &
       'time   ','zt     ','zm     ','dn0    ','u0     ','v0     ', & ! 1
       'fsttm  ','lsttm  ','nsmp   ','u      ','v      ','t      ', & ! 7
       'p      ','u_2    ','v_2    ','w_2    ','t_2    ','w_3    ', & !13
       't_3    ','tot_tw ','sfs_tw ','tot_uw ','sfs_uw ','tot_vw ', & !19
       'sfs_vw ','tot_ww ','sfs_ww ','km     ','kh     ','lmbd   ', & !25
       'lmbde  ','sfs_tke','sfs_boy','sfs_shr','boy_prd','shr_prd', & !31
       'trans  ','diss   ','dff_u  ','dff_v  ','dff_w  ','adv_u  ', & !37
       'adv_v  ','adv_w  ','prs_u  ','prs_v  ','prs_w  ','prd_uw ', & !43
       'storage','q      ','q_2    ','q_3    ','tot_qw ','sfs_qw ', & !49
       'rflx   ','rflx2  ','sflx   ','sflx2  ','l      ','l_2    ', & !55
       'l_3    ','tot_lw ','sed_lw ','cs1    ','cnt_cs1','w_cs1  ', & !61
       'tl_cs1 ','tv_cs1 ','rt_cs1 ','rl_cs1 ','wt_cs1 ','wv_cs1 ', & !67
       'wr_cs1 ','cs2    ','cnt_cs2','w_cs2  ','tl_cs2 ','tv_cs2 ', & !73
       'rt_cs2 ','rl_cs2 ','wt_cs2 ','wv_cs2 ','wr_cs2 ','Nc     ', & !79
       'Nr     ','rr     ','precip ','evap   ','frc_prc','prc_prc', & !85
       'frc_ran','hst_srf' /)

  real, save, allocatable   :: tke_sgs(:), tke_res(:), tke0(:), wtv_sgs(:),  &
       wtv_res(:), wrl_sgs(:), thvar(:), svctr(:,:), ssclr(:)

  public :: sflg, ssam_intvl, savg_intvl, statistics, init_stat, write_ps,   &
       acc_tend, updtst, sfc_stat, close_stat, fill_scalar, tke_sgs, sgsflxs,&
       sgs_vel, comp_tke, get_zi

contains
  !
  ! ---------------------------------------------------------------------
  ! INIT_STAT:  This routine initializes the statistical arrays which
  ! are user/problem defined.  Note that svctr is given 100 elements, and
  ! elements 90 and above are used for computing the TKE budget. Hence
  ! if (nvar2 >= 90 the program stops
  !
  subroutine init_stat(time, filprf, expnme, nzp)

    use grid, only : nxp, nyp, iradtyp
    use mpi_interface, only : myid

    character (len=80), intent (in) :: filprf, expnme
    integer, intent (in)            :: nzp
    real, intent (in)               :: time

    character (len=80) :: fname

    allocate (wtv_sgs(nzp),wtv_res(nzp),wrl_sgs(nzp))
    allocate (tke_res(nzp),tke_sgs(nzp),tke0(nzp),thvar(nzp))
    allocate (ssclr(nvar1),svctr(nzp,100))

    wtv_sgs(:) = 0.
    wtv_res(:) = 0.
    wrl_sgs(:) = 0.
    tke_res(:) = 0.
    tke_sgs(:) = 0.
    tke0(:)    = 0.

    svctr(:,:) = 0.
    ssclr(:)   = 0.                 ! changed from = -999. to = 0.

    select case(level)
    case (0)
       nv1 = 13
       nv2 = 58
    case (1)
       nv1 = 14
       nv2 = 58
    case (2)
       nv1 = 20
       nv2 = 83
       if (iradtyp == 3) nv1=21
    case default
       nv1 = nvar1
       nv2 = nvar2
    end select

    fname =  trim(filprf)//'.ts'
    if(myid == 0) print                                                  &
         "(//' ',49('-')/,' ',/,'  Initializing: ',A20)",trim(fname)
    call open_nc( fname, expnme, time, (nxp-4)*(nyp-4), ncid1, nrec1)
    call define_nc( ncid1, nrec1, nv1, s1)
    if (myid == 0) print *, '   ...starting record: ', nrec1

    fname =  trim(filprf)//'.ps'
    if(myid == 0) print                                                  &
         "(//' ',49('-')/,' ',/,'  Initiailzing: ',A20)",trim(fname)
    call open_nc( fname, expnme, time,(nxp-4)*(nyp-4), ncid2, nrec2)
    call define_nc( ncid2, nrec2, nv2, s2, n1=nzp)
    if (myid == 0) print *, '   ...starting record: ', nrec2

  end subroutine init_stat
  !
  ! ---------------------------------------------------------------------
  ! Subroutine Statistics:  This subroutine is the statistics driver
  ! it calls various other subroutines to compute and accumulate
  ! statistical quantities.  These are stored in two arrays:  SVCTR,
  ! and SSCLR (which accumulate scalar and vector statistics respectively
  !
  subroutine statistics(time)

    use grid, only : a_up, a_vp, a_wp, a_rc, a_theta, a_scr1, a_scr2   &
         , a_rp, a_tp, a_press, nxp, nyp, nzp, dzm, dzt, zm, zt, th00, umean &
         , vmean, dn0, precip, a_rpp, a_npp, albedo, CCN, iradtyp, a_rflx    &
         , a_sflx, albedo

    real, intent (in) :: time

    if (nsmp == 0.) fsttm = time
    nsmp=nsmp+1.
    ssclr(14:nvar1) = -999.
    !
    ! profile statistics
    !
    call accum_stat(nzp, nxp, nyp, a_up, a_vp, a_wp, a_tp, a_press, umean &
         ,vmean,th00)
    if (iradtyp == 3) then
       call accum_rad(nzp, nxp, nyp, a_rflx, sflx=a_sflx, alb=albedo)
    elseif (iradtyp > 0) then
       call accum_rad(nzp, nxp, nyp, a_rflx)
    end if
    if (level >=1) call accum_lvl1(nzp, nxp, nyp, a_rp)
    if (level >=2) call accum_lvl2(nzp, nxp, nyp, th00, dn0, zm, a_wp, &
         a_scr1, a_theta, a_tp, a_rc, a_scr2, a_rp)
    if (level >=3) call accum_lvl3(nzp, nxp, nyp, dn0, zm, a_rc, a_rpp,     &
         a_npp, precip, CCN)
    !
    ! scalar statistics
    !
    call set_ts(nzp, nxp, nyp, a_wp, a_theta, dn0, zt,zm,dzt,dzm,th00,time)
    if (level >=1) call ts_lvl1(nzp, nxp, nyp, dn0, zt, dzm, a_rp)
    if (level >=2) call ts_lvl2(nzp, nxp, nyp, a_rp, a_scr2, zt)

    call write_ts

  end subroutine statistics
  !
  ! -----------------------------------------------------------------------
  ! subroutine set_ts: computes and writes time sequence stats
  !
  subroutine set_ts(n1,n2,n3,w,th,dn0,zt,zm,dzt,dzm,th00,time)

    use defs, only : cp

    integer, intent(in) :: n1,n2,n3
    real, intent(in)    :: w(n1,n2,n3),th(n1,n2,n3)
    real, intent(in)    :: dn0(n1),zt(n1),zm(n1),dzt(n1),dzm(n1),th00,time

    integer :: k
    real    :: bf(n1)

    ssclr(1) = time
    ssclr(4) = get_zi(n1, n2, n3, 2, th, dzm, zt, 1.)   ! maximum gradient
    ssclr(5) = get_zi(n1, n2, n3, 3, th, thvar, zt, 1.) ! maximum variance
    !
    ! buoyancy flux statistics
    !
    ssclr(7) = 0.
    do k = 2,n1-2
       bf(k) = wtv_res(k) + wtv_sgs(k)
       ssclr(7) = ssclr(7) + (tke_res(k)+tke_sgs(k))*dn0(k)/dzt(k)
       svctr(k,33) = svctr(k,33) + wtv_sgs(k)*9.8/th00
    end do
    ssclr(6) = get_zi(n1, n2, n3, 4, th, bf, zm, 1.) ! minimum buoyancy flux

    ssclr(8) = bf(2)
    ssclr(9) = maxval(w)

    ssclr(12) = ssclr(12)*cp*(dn0(1)+dn0(2))*0.5

  end subroutine set_ts
  !
  ! -----------------------------------------------------------------------
  ! subroutine ts_lvl1: computes and writes time sequence stats; for the
  ! zi calculation setting itype=1 selects a concentration threshold
  !
  subroutine ts_lvl1(n1,n2,n3,dn0,zt,dzm,q)

    use defs, only : alvl

    integer, intent(in) :: n1,n2,n3
    real, intent(in)    :: q(n1,n2,n3)
    real, intent(in)    :: dn0(n1),zt(n1),dzm(n1)

    ssclr(13) = ssclr(13)*alvl*(dn0(1)+dn0(2))*0.5
    ssclr(14) = get_zi(n1, n2, n3, 1, q, dzm, zt, 0.5e-3)

  end subroutine ts_lvl1
  !
  ! -----------------------------------------------------------------------
  ! subroutine ts_lvl2: computes and writes time sequence stats
  !
  subroutine ts_lvl2(n1,n2,n3,rt,rs,zt)

    integer, intent(in) :: n1,n2,n3
    real, intent(in)    :: rt(n1,n2,n3),rs(n1,n2,n3), zt(n1)

    integer :: k,i,j
    real    :: cpnt, unit, xaqua

    ssclr(18)  = zt(n1)
    ssclr(19)  = 0.
    ssclr(20)  = 0.

    unit = 1./real((n2-4)*(n3-4))
    do j=3,n3-2
       do i=3,n2-2
          cpnt  = 0.
          do k=2,n1-2
             xaqua = rt(k,i,j) - rs(k,i,j)
             if (xaqua > 1.e-5) then
                ssclr(17) = max(ssclr(17),zt(k))
                ssclr(18) = min(ssclr(18),zt(k))
                cpnt = unit
                ssclr(20) = max(ssclr(20), xaqua)
             end if
          end do
          ssclr(19) = ssclr(19) + cpnt
       end do
    end do

    if (ssclr(18) == zt(n1)) ssclr(18) = -999.

  end subroutine ts_lvl2
  !
  !---------------------------------------------------------------------
  ! SUBROUTINE ACCUM_STAT: Accumulates various statistics over an
  ! averaging period for base (level 0) version of model
  !
  subroutine accum_stat(n1,n2,n3,u,v,w,t,p,um,vm,th00)

    integer, intent (in) :: n1,n2,n3
    real, dimension (n1,n2,n3), intent (in)    :: u, v, w, t, p
    real, intent (in)           :: um, vm, th00

    integer :: i,j,k
    real    :: a1(n1), b1(n1), c1(n1), d1(n1), a3(n1), b3(n1), x

    x = 1./real( (n3-4)*(n2-4))
    call get_avg3(n1,n2,n3, u,a1)
    call get_avg3(n1,n2,n3, v,b1)
    call get_avg3(n1,n2,n3, t,c1)
    call get_avg3(n1,n2,n3, p,d1)
    call get_var3(n1,n2,n3, t, c1, thvar)

    a3(:) = 0.
    b3(:) = 0.
    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             a3(k)=a3(k) + w(k,i,j)**3
             b3(k)=b3(k) + (t(k,i,j)-a1(k))**3
          end do
       end do
    end do

    do k=1,n1
       svctr(k,10)=svctr(k,10) + a1(k) + um
       svctr(k,11)=svctr(k,11) + b1(k) + vm
       svctr(k,12)=svctr(k,12) + c1(k) + th00
       svctr(k,13)=svctr(k,13) + d1(k)
       svctr(k,17)=svctr(k,17) + thvar(k)
       svctr(k,18)=svctr(k,18) + a3(k) * x
       svctr(k,19)=svctr(k,19) + b3(k) * x
    end do

  end subroutine accum_stat
  !
  !---------------------------------------------------------------------
  ! SUBROUTINE ACCUM_STAT: Accumulates various statistics over an
  ! averaging period for radiation variables
  !
  subroutine accum_rad(n1,n2,n3,rflx,sflx,alb)

    integer, intent (in) :: n1,n2,n3
    real, intent (in)    :: rflx(n1,n2,n3)
    real, optional, intent (in) :: sflx(n1,n2,n3), alb(n2,n3)

    integer :: k
    real    :: a1(n1),a2(n1)

    call get_avg3(n1,n2,n3,rflx,a1)
    call get_var3(n1,n2,n3,rflx,a1,a2)
    do k=1,n1
       svctr(k,55)=svctr(k,55) + a1(k)
       svctr(k,56)=svctr(k,56) + a2(k)
    end do

    if (present(sflx)) then
       call get_avg3(n1,n2,n3,sflx,a1)
       call get_var3(n1,n2,n3,sflx,a1,a2)
       do k=1,n1
          svctr(k,57)=svctr(k,57) + a1(k)
          svctr(k,58)=svctr(k,58) + a2(k)
       end do
       ssclr(21) = get_avg(1,n2,n3,1,alb)
    end if

  end subroutine accum_rad
  !
  !---------------------------------------------------------------------
  ! SUBROUTINE ACCUM_LVL1: Accumulates various statistics over an
  ! averaging period for moisture variable (smoke or total water)
  !
  subroutine accum_lvl1(n1,n2,n3,rt)

    integer, intent (in) :: n1,n2,n3
    real, intent (in)  :: rt(n1,n2,n3)

    integer :: i,j,k
    real    :: a1(n1),a2(n1),a3(n1)

    call get_avg3(n1,n2,n3,rt,a1)
    call get_var3(n1,n2,n3,rt,a1,a2)

    a3(:) = 0.
    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             a3(k) = a3(k) + (rt(k,i,j)-a1(k))**3
          end do
       end do
    end do

    do k=1,n1
       svctr(k,50)=svctr(k,50) + a1(k)*1000.
       svctr(k,51)=svctr(k,51) + a2(k)
       svctr(k,52)=svctr(k,52) + a3(k)/REAL((n2-4)*(n3-4))
    end do

  end subroutine accum_lvl1
  !
  !---------------------------------------------------------------------
  ! SUBROUTINE ACCUM_LVL2: Accumulates specialized statistics that depend
  ! on level 3 variables.
  !
  subroutine accum_lvl2(n1, n2, n3, th00, dn0, zm, w, tv, th, tl, &
       rl, rs, rt)

    use defs, only : ep2

    integer, intent (in) :: n1,n2,n3
    real, intent (in)                       :: th00
    real, intent (in), dimension(n1)        :: zm, dn0
    real, intent (in), dimension(n1,n2,n3)  :: w, th, tl, rl, rs, rt
    real, intent (out), dimension(n1,n2,n3) :: tv

    integer                   :: k, i, j, km1
    logical                   :: aflg
    real                      :: xy1mx
    real, dimension(n1)       :: a1, a2, a3, tvbar
    real, dimension(n2,n3)    :: scr, xy1, xy2

    !
    ! liquid water statistics
    !
    call get_avg3(n1,n2,n3,rl,a1)
    call get_var3(n1,n2,n3,rl,a1,a2)

    a3(:) = 0.
    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             a3(k) = a3(k) + (rl(k,i,j)-a1(k))**3
          end do
       end do
    end do

    do k=1,n1
       svctr(k,59)=svctr(k,59) + a1(k)*1000.
       svctr(k,60)=svctr(k,60) + a2(k)
       svctr(k,61)=svctr(k,61) + a3(k)/REAL((n2-4)*(n3-4))
    end do
    !
    ! do some conditional sampling statistics: cloud, cloud-core
    !
    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             tv(k,i,j) = th(k,i,j)*(1.+ep2*rt(k,i,j) - rl(k,i,j))
          end do
       end do
    end do
    call get_avg3(n1,n2,n3,tv,tvbar)

    xy1mx = 0.
    do k=1,n1
       aflg = .false.
       do j=3,n3-2
          do i=3,n2-2
             xy1(i,j) = 0.
             xy2(i,j) = 0.
             if (rt(k,i,j) > rs(k,i,j) + 0.01e-3) then
                aflg = .true.
                xy1(i,j)=1.
                if (tv(k,i,j) > tvbar(k)) xy2(i,j)=1.
             end if
          end do
       end do

       if (aflg) then
          svctr(k,64)=svctr(k,64)+get_avg(1,n2,n3,1,xy1)
          svctr(k,65)=svctr(k,65)+get_csum(1,n2,n3,1,xy1,xy1)
          svctr(k,66)=svctr(k,66)+get_csum(n1,n2,n3,k,w,xy1)
          svctr(k,67)=svctr(k,67)+get_csum(n1,n2,n3,k,tl,xy1)
          svctr(k,68)=svctr(k,68)+get_csum(n1,n2,n3,k,tv,xy1)
          svctr(k,69)=svctr(k,69)+get_csum(n1,n2,n3,k,rt,xy1)*1000.
          svctr(k,70)=svctr(k,70)+get_csum(n1,n2,n3,k,rl,xy1)*1000.

          svctr(k,74)=svctr(k,74)+get_avg(1,n2,n3,1,xy2)
          svctr(k,75)=svctr(k,75)+get_csum(1,n2,n3,1,xy2,xy2)
          svctr(k,76)=svctr(k,76)+get_csum(n1,n2,n3,k,w,xy2)
          svctr(k,77)=svctr(k,77)+get_csum(n1,n2,n3,k,tl,xy2)
          svctr(k,78)=svctr(k,78)+get_csum(n1,n2,n3,k,tv,xy2)
          svctr(k,79)=svctr(k,79)+get_csum(n1,n2,n3,k,rt,xy2)*1000.
          svctr(k,80)=svctr(k,80)+get_csum(n1,n2,n3,k,rl,xy2)*1000.

          do j=3,n3-2
             do i=3,n2-2
                scr(i,j)=(.5*(tl(k,i,j)+tl(k+1,i,j))+th00)*w(k,i,j)
             end do
          end do
          svctr(k,71)=svctr(k,71)+get_csum(1,n2,n3,1,scr,xy1)
          svctr(k,81)=svctr(k,81)+get_csum(1,n2,n3,1,scr,xy2)

          do j=3,n3-2
             do i=3,n2-2
                scr(i,j)=(.5*(tv(k,i,j)+tv(k+1,i,j)))*w(k,i,j)
             end do
          end do
          svctr(k,72)=svctr(k,72)+get_csum(1,n2,n3,1,scr,xy1)
          svctr(k,82)=svctr(k,82)+get_csum(1,n2,n3,1,scr,xy2)

          do j=3,n3-2
             do i=3,n2-2
                scr(i,j)=(.5*(rt(k,i,j)+rt(k+1,i,j)))*w(k,i,j)
             end do
          end do
          svctr(k,73)=svctr(k,73)+get_csum(1,n2,n3,1,scr,xy1)
          svctr(k,83)=svctr(k,83)+get_csum(1,n2,n3,1,scr,xy2)
       end if
    end do
    !
    ! water paths
    !
    do j=3,n3-2
       do i=3,n2-2
          scr(i,j) = 0.
          do k=1,n1
             km1=max(1,k-1)
             scr(i,j)=scr(i,j)+rl(k,i,j)*dn0(k)*(zm(k)-zm(km1))*1000.
          enddo
       end do
    end do
    ssclr(15) = get_avg(1,n2,n3,1,scr)
    ssclr(16) = get_cor(1,n2,n3,1,scr,scr)
  end subroutine accum_lvl2
  !
  !---------------------------------------------------------------------
  ! SUBROUTINE ACCUM_LVL3: Accumulates specialized statistics that depend
  ! on level 3 variables.
  !
  subroutine accum_lvl3(n1, n2, n3, dn0, zm, rc, rr, nr, rrate, CCN)

    use defs, only : alvl

    integer, intent (in) :: n1,n2,n3
    real, intent (in)                      :: CCN
    real, intent (in), dimension(n1)       :: zm, dn0
    real, intent (in), dimension(n1,n2,n3) :: rc, rr, nr, rrate

    integer                :: k, i, j, km1
    real                   :: nrsum, nrcnt, rrsum, rrcnt, xrain, xaqua
    real                   :: rmax, rmin
    real, dimension(n1)    :: a1
    real, dimension(n2,n3) :: scr1,scr2
    logical                :: aflg

    !
    ! conditionally average rain numbers, and droplet concentrations
    !

    call get_avg3(n1,n2,n3,rr,a1)
    nrsum = 0.
    nrcnt = 0.
    do k=1,n1-1
       aflg = .false.
       do j=3,n3-2
          do i=3,n2-2
             scr1(i,j) = 0.
             if (rr(k,i,j) > 0.001e-3) then
                aflg = .true.
                scr1(i,j) = 1.
                nrsum = nrsum + nr(k,i,j)
                nrcnt = nrcnt + 1.
             end if
          end do
       end do
       svctr(k,84)=svctr(k,84) + CCN*dn0(k)/1000. ! nc (1/litre)
       svctr(k,86)=svctr(k,86) + a1(k)*1000.
       if (aflg) then
          svctr(k,85)=svctr(k,85)+get_csum(n1,n2,n3,k,nr,scr1)
          svctr(k,91)=svctr(k,91)+get_avg(1,n2,n3,1,scr1)
       end if
    end do
    !
    ! conditionally average precip fluxes
    !
    do k=2,n1-2
       rrsum = 0.
       rrcnt = 0.
       aflg = .false.
       do j=3,n3-2
          do i=3,n2-2
             scr1(i,j) = 0.
             if (rrate(k,i,j) > 3.65e-5) then
                aflg = .true.
                scr1(i,j) = 1.
                rrsum = rrsum + rrate(k,i,j) * alvl * 0.5*(dn0(1)+dn0(2))
                rrcnt = rrcnt + 1.
             end if
          end do
       end do
       if (aflg) then
          if (k == 2 ) ssclr(24) = rrcnt
          svctr(k-1,90)=svctr(k-1,90)+get_csum(n1,n2,n3,k,rrate,scr1)
          svctr(k-1,89)=svctr(k-1,89)+get_avg(1,n2,n3,1,scr1)
       end if
    end do
    !
    ! histograms
    !
    do k=1,n1
       rrcnt = 0.
       do j=3,n3-2
          do i=3,n2-2
             rmin = max(6.2e-8,(k-1)*3.421e-5)
             rmax =  k * 3.421e-5
             if (rrate(2,i,j) > rmin .and. rrate(2,i,j) <= rmax) rrcnt=rrcnt+1.
          end do
       end do
       if (rrcnt > 0.) svctr(k,92)=svctr(k,92)+rrcnt
    end do
    !
    ! water paths
    !
    do j=3,n3-2
       do i=3,n2-2
          scr1(i,j) = 0.
          scr2(i,j) = 0.
          do k=1,n1
             km1=max(1,k-1)
             xrain = max(0.,rr(k,i,j))
             xaqua = max(xrain,rc(k,i,j))
             scr1(i,j)=scr1(i,j)+xaqua*dn0(k)*(zm(k)-zm(km1))*1000.
             scr2(i,j)=scr2(i,j)+xrain*dn0(k)*(zm(k)-zm(km1))*1000.
          enddo
       end do
    end do
    ssclr(15) = get_avg(1,n2,n3,1,scr1)
    ssclr(16) = get_cor(1,n2,n3,1,scr1,scr1)
    ssclr(22) = get_avg(1,n2,n3,1,scr2)
    ssclr(23) = get_avg(1,n2,n3,1,rrate(2,:,:))
    ssclr(25) = CCN*1.e-6 ! per cc
    ssclr(26) = nrsum
    ssclr(27) = nrcnt

  end subroutine accum_lvl3
  !
  ! ---------------------------------------------------------------------
  ! subroutine comp_tke: calculates some components of the turbulent
  ! kinetic energy budgets and velocity statistics
  !
  subroutine comp_tke(n1,n2,n3,dzm,th00,u,v,w,s,scr)

    integer, intent (in) :: n1,n2,n3
    real, intent (in)    :: dzm(n1),th00,u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3)
    real, intent (inout) :: s(n1,n2,n3)
    real, intent (out) :: scr(n1,n2,n3)

    integer :: k,kp1,i,j
    real    :: x1(n1), x2(n1), xx

    !
    ! ------
    ! Calculates buoyancy forcing
    !
    call get_buoyancy(n1,n2,n3,s,w,th00)
    !
    ! ------
    ! Estimates shear component of TKE budget
    !
    call get_shear(n1,n2,n3,u,v,w,dzm)
    !
    ! ------
    ! Calculates horizontal variances and resolved TKE
    !
    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             scr(k,i,j)=u(k,i,j)**2
          end do
       end do
    end do
    call get_avg3(n1,n2,n3,scr,x1)
    call get_avg3(n1,n2,n3,u,x2)
    do k=1,n1
       svctr(k,14) = svctr(k,14) + x1(k)
       xx = x1(k)-x2(k)**2
       tke_res(k)  = xx
    end do

    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             scr(k,i,j)=v(k,i,j)**2
          end do
       end do
    end do
    call get_avg3(n1,n2,n3,scr,x1)
    call get_avg3(n1,n2,n3,v,x2)
    do k=1,n1
       svctr(k,15) = svctr(k,15) + x1(k)
       xx = x1(k)-x2(k)**2
       tke_res(k)  = tke_res(k) + xx
    end do

    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             scr(k,i,j)=w(k,i,j)**2
          end do
       end do
    end do
    call get_avg3(n1,n2,n3,scr,x1)

    do k=1,n1
       svctr(k,16) = svctr(k,16)+x1(k)
       kp1 = min(k+1,n1)
       tke_res(k)  = 0.5*(0.5*(tke_res(k)+tke_res(kp1)) + x1(k))
       if (nsmp == 0) tke0(k) = tke_res(k)
    end do

  end subroutine comp_tke
  !
  ! ---------------------------------------------------------------------
  ! get_buoyancy:  estimates buoyancy production term in tke budget
  !
  subroutine get_buoyancy(n1,n2,n3,b,w,th00)

    use defs, only : g

    integer, intent(in) :: n1,n2,n3
    real, intent(in)    :: w(n1,n2,n3),th00
    real, intent(inout) :: b(n1,n2,n3)

    integer :: i,j,k,kp1

    do j=3,n3-2
       do i=3,n2-2
          do k=1,n1
             kp1 = min(k+1,n1)
             b(k,i,j) = (b(k,i,j) + b(kp1,i,j))
          end do
       end do
    end do
    call get_cor3(n1,n2,n3,b,w,wtv_res)
    do k=1,n1
       svctr(k,35) = svctr(k,35) + wtv_res(k)
       wtv_res(k) = wtv_res(k) * th00/g
    end do

  end subroutine get_buoyancy
  !
  ! ---------------------------------------------------------------------
  ! get_shear:  estimates shear production term in tke budget
  !
  subroutine get_shear(n1,n2,n3,u,v,w,dzm)

    integer, intent(in) :: n3,n2,n1
    real, intent(in)    :: w(n1,n2,n3),dzm(n1),u(n1,n2,n3),v(n1,n2,n3)

    real :: ub(n1), vb(n1)
    integer i,j,k
    real fact, uw_shear, vw_shear

    fact = 0.25/float((n2-4)*(n3-4))

    call get_avg3(n1,n2,n3,u,ub)
    call get_avg3(n1,n2,n3,v,vb)

    do j=3,n3-2
       do i=3,n2-2
          do k=2,n1-1
             uw_shear = -(u(k,i,j)-ub(k))*fact*(                              &
                  (w(k,i,j)  +w(k,i+1,j)  )*(ub(k+1)-ub(k)  )*dzm(k) +        &
                  (w(k-1,i,j)+w(k-1,i+1,j))*(ub(k)  -ub(k-1))*dzm(k-1))
             if (j > 1) vw_shear = -(v(k,i,j)-vb(k))*fact*(                   &
                  (w(k,i,j)  +w(k,i,j+1)  )*(vb(k+1)-vb(k)  )*dzm(k) +        &
                  (w(k-1,i,j)+w(k-1,i,j+1))*(vb(k)  -vb(k-1))*dzm(k-1))

             svctr(k,48) = svctr(k,48)+uw_shear
             svctr(k,36) = svctr(k,36)+uw_shear+vw_shear
          end do
       end do
    end do

  end subroutine get_shear
  !
  ! ----------------------------------------------------------------------
  ! Subroutine write_ts: writes the statistics file
  !
  subroutine write_ts

    use netcdf

    integer :: iret, n, VarID
    !
    ! define different dimensions
    !
    do n=1,nv1
       iret = nf90_inq_varid(ncid1, s1(n), VarID)
       iret = nf90_put_var(ncid1, VarID, ssclr(n), start=(/nrec1/))
       ssclr(n) = 0.
    end do
    iret = nf90_sync(ncid1)
    nrec1 = nrec1 + 1

  end subroutine write_ts
  !
  ! ----------------------------------------------------------------------
  ! Subroutine write_ps: writes the time averaged elements of the
  ! statistics file
  !
  subroutine  write_ps(n1,dn0,u0,v0,zm,zt,time)

    use netcdf
    use defs, only : alvl, cp

    integer, intent (in) :: n1
    real, intent (in)    :: time
    real, intent (in)    :: dn0(n1), u0(n1), v0(n1), zm(n1), zt(n1)

    integer :: iret, VarID, k, n, kp1

    lsttm = time
    do k=1,n1
       kp1 = min(n1,k+1)
       svctr(k,20) = (svctr(k,20)+svctr(k,21))*cp
       svctr(k,22) = svctr(k,22)+svctr(k,23)
       svctr(k,24) = svctr(k,24)+svctr(k,25)
       svctr(k,26) = svctr(k,26)+svctr(k,27)
       svctr(k,53) = (svctr(k,53)+svctr(k,54))*alvl
       svctr(k,21) = svctr(k,21)*cp
       svctr(k,54) = svctr(k,54)*alvl
       svctr(k,62) = svctr(k,62)*alvl
       svctr(k,37) = svctr(k,44) + svctr(k,47) +(                         &
            +svctr(k,45) + svctr(kp1,45) + svctr(k,46) + svctr(kp1,46)    &
            +svctr(k,42) + svctr(kp1,42) + svctr(k,43) + svctr(kp1,43)    &
            -svctr(k,36) - svctr(kp1,36)   )*0.5
       if (lsttm>fsttm) then
          svctr(k,49) = (tke_res(k) - tke0(k))/(lsttm-fsttm)
       else
          svctr(k,49) = 0.
       end if
       svctr(k,10:nv2) = svctr(k,10:nv2)/nsmp
    end do

    iret = nf90_inq_VarID(ncid2, s2(1), VarID)
    iret = nf90_put_var(ncid2, VarID, time, start=(/nrec2/))
    if (nrec2 == 1) then
       iret = nf90_inq_varid(ncid2, s2(2), VarID)
       iret = nf90_put_var(ncid2, VarID, zt, start = (/nrec2/))
       iret = nf90_inq_varid(ncid2, s2(3), VarID)
       iret = nf90_put_var(ncid2, VarID, zm, start = (/nrec2/))
       iret = nf90_inq_varid(ncid2, s2(4), VarID)
       iret = nf90_put_var(ncid2, VarID, dn0, start = (/nrec2/))
       iret = nf90_inq_varid(ncid2, s2(5), VarID)
       iret = nf90_put_var(ncid2, VarID, u0, start = (/nrec2/))
       iret = nf90_inq_varid(ncid2, s2(6), VarID)
       iret = nf90_put_var(ncid2, VarID, v0, start = (/nrec2/))
    end if

    iret = nf90_inq_VarID(ncid2, s2(7), VarID)
    iret = nf90_put_var(ncid2, VarID, fsttm, start=(/nrec2/))
    iret = nf90_inq_VarID(ncid2, s2(8), VarID)
    iret = nf90_put_var(ncid2, VarID, lsttm, start=(/nrec2/))
    iret = nf90_inq_VarID(ncid2, s2(9), VarID)
    iret = nf90_put_var(ncid2, VarID, nsmp,  start=(/nrec2/))

    do n=10,nv2
       iret = nf90_inq_varid(ncid2, s2(n), VarID)
       iret = nf90_put_var(ncid2,VarID,svctr(:,n), start=(/1,nrec2/),    &
            count=(/n1,1/))
    end do

    iret  = nf90_sync(ncid2)
    nrec2 = nrec2+1
    nsmp  = 0.

    do k=1,n1
       svctr(k,:) = 0.
    end do

  end subroutine write_ps
  !
  ! ----------------------------------------------------------------------
  ! subroutine: sfc_stat:  Updates statistical arrays with surface flux
  ! variables
  !
  subroutine sfc_stat(n2,n3,tflx,qflx,ustar,sst)

    integer, intent(in) :: n2,n3
    real, intent(in), dimension(n2,n3) :: tflx, qflx, ustar
    real, intent(in)    :: sst

    ssclr(10) = sst
    ssclr(11) = get_avg(1,n2,n3,1,ustar)

    ssclr(12) = get_avg(1,n2,n3,1,tflx)
    if (level >= 1) ssclr(13) = get_avg(1,n2,n3,1,qflx)

  end subroutine sfc_stat
  !
  ! ----------------------------------------------------------------------
  ! subroutine: fills scalar array based on index
  ! 1: cfl; 2 max divergence
  !
  subroutine fill_scalar(index,xval)

    integer, intent(in) :: index
    real, intent (in)   :: xval

    select case(index)
    case(1)
       ssclr(2) = xval
    case(2)
       ssclr(3) = xval
    end select

  end subroutine fill_scalar
  !
  ! ----------------------------------------------------------------------
  ! subroutine: calculates the dissipation for output diagnostics, if
  ! isgstyp equals 2 then le is passed in via diss
  !
  subroutine sgs_vel(n1,n2,n3,v1,v2,v3)

    integer, intent(in) :: n1,n2,n3
    real, intent(in)    :: v1(n1),v2(n1),v3(n1)

    svctr(:,23)=svctr(:,23)+v1(:)/float((n2-2)*(n3-2))
    svctr(:,25)=svctr(:,25)+v2(:)/float((n2-2)*(n3-2))
    svctr(:,27)=svctr(:,27)+v3(:)/float((n2-2)*(n3-2))

  end subroutine sgs_vel
  !
  ! --------------------------------------------------------------------------
  ! SGSFLXS: estimates the sgs rl and tv flux from the sgs theta_l and sgs r_t
  ! fluxes
  !
  subroutine sgsflxs(n1,n2,n3,level,rl,rv,th,flx,type)

    use defs, only : alvl, cp, rm, ep2

    integer, intent(in) :: n1,n2,n3,level
    real, intent(in)    :: rl(n1,n2,n3),rv(n1,n2,n3)
    real, intent(in)    :: th(n1,n2,n3),flx(n1,n2,n3)
    character (len=2)   :: type

    integer :: k,i,j
    real    :: rnpts      ! reciprical of number of points and
    real    :: fctl, fctt ! factors for liquid (l) and tv (t) fluxes

    if (type == 'tl') then
       wrl_sgs(:) = 0.
       wtv_sgs(:) = 0.
    end if
    rnpts = 1./real((n2-4)*(n3-4))
    !
    ! calculate fluxes assuming the possibility of liquid water.  if liquid
    ! water does not exist sgs_rl = 0.
    !
    if ( level >= 2 ) then
       do j = 3,n3-2
          do i = 3,n2-2
             do k = 1,n1-1
                if (rl(k+1,i,j) > 0.) then
                   fctt = rnpts*(1. + rv(k,i,j)*(1.+ep2 +ep2*rv(k,i,j)*alvl   &
                        /(rm*th(k,i,j))))                                     &
                        /(1.+(rv(k,i,j)*(alvl/th(k,i,j))**2)/(rm*cp))
                   select case (type)
                   case ('tl')
                      fctl =-rnpts/(rm*th(k,i,j)**2/(rv(k,i,j)*alvl)+alvl/cp)
                   case ('rt')
                      fctl =rnpts/(1.+(rv(k,i,j)*alvl**2)/(cp*rm*th(k,i,j)**2))
                      fctt = (alvl*fctt/cp - th(k,i,j)*rnpts)
                   end select
                   wrl_sgs(k) = wrl_sgs(k) + fctl*flx(k,i,j)
                   wtv_sgs(k) = wtv_sgs(k) + fctt*flx(k,i,j)
                else
                   select case (type)
                   case ('tl')
                      fctt = rnpts*(1. + ep2*rv(k,i,j))
                   case ('rt')
                      fctt = rnpts*(ep2*th(k,i,j))
                   end select
                   wtv_sgs(k) = wtv_sgs(k) + fctt*flx(k,i,j)
                end if
             end do
          end do
       end do
       !
       ! calculate fluxes for dry thermodynamics, i.e., wrl_sgs is by def
       ! zero
       !
    else
       do k = 1,n1
          wrl_sgs(k) = 0.
       end do
       do j = 3,n3-2
          do i = 3,n2-2
             do k = 1,n1-1
                if ( level >= 1) then
                   select case (type)
                   case ('tl')
                      fctt = rnpts * (1. + ep2*rv(k,i,j))
                   case ('rt')
                      fctt = rnpts * ep2*th(k,i,j)
                   end select
                else
                   fctt = rnpts
                end if
                wtv_sgs(k) = wtv_sgs(k) + fctt*flx(k,i,j)
             end do
          end do
       end do
    end if

  end subroutine sgsflxs
  !
  ! ----------------------------------------------------------------------
  ! subroutine fill_tend: fills arrays with current value of tendencies
  !
  subroutine acc_tend(n1,n2,n3,f1,f2,f3,t1,t2,t3,v1,v2,v3,ic,routine)

    integer, intent(in) :: n1,n2,n3,ic
    real, intent(in)    :: f1(n1,n2,n3),f2(n1,n2,n3),f3(n1,n2,n3)
    real, intent(in)    :: t1(n1,n2,n3),t2(n1,n2,n3),t3(n1,n2,n3)
    real, intent(inout) :: v1(n1),v2(n1),v3(n1)
    character (len=3)   :: routine

    integer :: k,ii
    real    :: x1(n1),x2(n1),x3(n1)

    call get_cor3(n1,n2,n3,f1,t1,x1)
    call get_cor3(n1,n2,n3,f2,t2,x2)
    call get_cor3(n1,n2,n3,f3,t3,x3)

    select case (routine)
    case ('sgs')
       ii = 39
    case ('adv')
       ii = 42
    end select

    select case (ic)
    case (1)
       do k=1,n1
          v1(k) = x1(k)
          v2(k) = x2(k)
          v3(k) = x3(k)
       end do
    case (2)
       do k=1,n1
          svctr(k,ii)   = svctr(k,ii)   + (x1(k)-v1(k))
          svctr(k,ii+1) = svctr(k,ii+1) + (x2(k)-v2(k))
          svctr(k,ii+2) = svctr(k,ii+2) + (x3(k)-v3(k))
       end do
    end select

  end subroutine acc_tend
  !
  !---------------------------------------------------------------------
  ! subroutine updtst: updates appropriate statistical arrays
  !
  subroutine updtst(n1,routine,nfld,values,ic)

    integer, intent(in)            :: n1,nfld,ic
    real, intent (in)              :: values(n1)
    character (len=3), intent (in) :: routine

    integer :: nn,k

    select case (routine)
    case("sgs")
       select case (nfld)
       case (-6)
          nn = 31 ! dissipation length-scale
       case (-5)
          nn = 30 ! mixing length
       case (-4)
          nn = 29 ! eddy diffusivity
       case (-3)
          nn = 28 ! eddy viscosity
       case (-2)
          nn = 38 ! dissipation
       case (-1)
          nn = 32 ! estimated sgs energy
       case (1)
          nn = 21 ! sgs tl flux
       case (2)
          nn = 54 ! sgs rt flux
       case default
          nn = 0
       end select
    case("adv")
       select case (nfld)
       case (-3)
          nn = 26 ! adv w flux
       case (-2)
          nn = 24 ! adv v flux
       case (-1)
          nn = 22 ! adv u flux
       case (0)
          nn = 62 ! adv rl flux
       case (1)
          nn = 20 ! adv tl flux
       case (2)
          nn = 53 ! adv rt flux
       case default
          nn = 0
       end select
    case("prs")
       select case (nfld)
       case (1)
          nn = 45 ! dpdx u corr
       case (2)
          nn = 46 ! dpdy v corr
       case (3)
          nn = 47 ! dpdz w corr
       case default
          nn = 0
       end select
    case("prc")
       select case (nfld)
       case (1)
          nn = 87
       case (2)
          nn = 88
       case (3)
          nn = 63
       case default
          nn = 0
       end select
    case default
       nn = 0
    end select

    if (nn > 0) then
       if (ic == 0) svctr(:,nn)=0.
       do k=1,n1
          svctr(k,nn)=svctr(k,nn)+values(k)
       enddo
    end if

  end subroutine updtst
  !
  ! -------------------------------------------------------------------------
  !
  integer function close_stat()

    use netcdf

    close_stat = nf90_close(ncid1) + nf90_close(ncid2)

  end function close_stat
  !
  ! -------------------------------------------------------------------------
  !
  real function get_zi (n1, n2, n3, itype, sx, xx, z, threshold)

    integer, intent (in) :: n1, n2, n3, itype
    real, intent (in)    :: xx(n1), z(n1), sx(n1,n2,n3), threshold

    integer :: i, j, k, kk
    real    :: zibar, sval, dmy, scr(n2,n3)

    get_zi = -999.
    select case(itype)
    case (1)
       !
       ! find level at which sx=threshold (xx is one over grid spacing)
       !
       zibar = 0.
       do j=3,n3-2
          do i=3,n2-2
             k = 2
             do while (k < n1-2 .and. sx(k,i,j) > threshold)
                k = k+1
             end do
             if (k == n1-2) zibar = -999.
             if (zibar /= -999.) zibar = zibar + z(k-1) +  &
                  (threshold - sx(k-1,i,j))/xx(k-1)     /  &
                  (sx(k,i,j) - sx(k-1,i,j) + epsilon(1.))
          end do
       end do
       if (zibar /= -999.) get_zi = zibar/real((n3-4)*(n2-4))

    case(2)
       !
       ! find level of maximum gradient (xx is one over grid spacing)
       !
       scr=0.
       do j=3,n3-2
          do i=3,n2-2
             sval = 0.
             do k=2,n1-5
                dmy = (sx(k+1,i,j)-sx(k,i,j))*xx(k)
                if (dmy > sval) then
                   sval = dmy
                   scr(i,j) = z(k)
                end if
             end do
          end do
       end do
       get_zi = get_avg(1,n2,n3,1,scr)

    case(3)
       !
       ! find level where xx is a maximum
       !
       sval = -huge(1.)
       kk = 1
       do k=2,n1
          if (xx(k) > sval) then
             kk = k
             sval = xx(k)
          end if
       end do
       get_zi = z(kk)

    case(4)
       !
       ! find level where xx is a maximum
       !
       sval = huge(1.)
       kk = 1
       do k=2,n1-2
          if (xx(k) < sval) then
             kk = k
             sval = xx(k)
          end if
       end do
       get_zi = z(kk)
    end select

  end function get_zi

end module stat

max_stars_repo_name: waqarnabi/tybec

max_issues_repo_name: waqarnabi/tybec

max_forks_repo_name: waqarnabi/tybec


text:
! This file was automatically generated by SWIG (http://www.swig.org).
! Version 4.0.2
!
! Do not make changes to this file unless you know what you are doing--modify
! the SWIG interface file instead.
module pointers_mod
 use, intrinsic :: ISO_C_BINDING
 implicit none
 private

 ! DECLARATION CONSTRUCTS
 type, bind(C) :: SwigArrayWrapper
  type(C_PTR), public :: data = C_NULL_PTR
  integer(C_SIZE_T), public :: size = 0
 end type
 public :: fill_with_zeros
 public :: accumulate

! WRAPPER DECLARATIONS
interface
subroutine swigc_fill_with_zeros(farg1) &
bind(C, name="_wrap_fill_with_zeros")
use, intrinsic :: ISO_C_BINDING
import :: swigarraywrapper
type(SwigArrayWrapper) :: farg1
end subroutine

function swigc_accumulate(farg1) &
bind(C, name="_wrap_accumulate") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
import :: swigarraywrapper
type(SwigArrayWrapper) :: farg1
integer(C_INT) :: fresult
end function

end interface

contains
 ! MODULE SUBPROGRAMS
subroutine SWIGTM_fin_double_Sb__SB_(finp, iminp)
  use, intrinsic :: ISO_C_BINDING
  real(C_DOUBLE), dimension(:), intent(in), target :: finp
  type(SwigArrayWrapper), intent(out) :: iminp
  integer(C_SIZE_T) :: sz
  real(C_DOUBLE), pointer :: imtemp

  sz = size(finp, kind=C_SIZE_T)
  if (sz > 0_c_size_t) then
    imtemp => finp(1)
    iminp%data = c_loc(imtemp)
  else
    iminp%data = c_null_ptr
  end if
  iminp%size = sz
end subroutine
subroutine fill_with_zeros(x)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), dimension(:), target :: x
type(SwigArrayWrapper) :: farg1

call SWIGTM_fin_double_Sb__SB_(x, farg1)
call swigc_fill_with_zeros(farg1)
end subroutine

subroutine SWIGTM_fin_int_Sb__SB_(finp, iminp)
  use, intrinsic :: ISO_C_BINDING
  integer(C_INT), dimension(:), intent(in), target :: finp
  type(SwigArrayWrapper), intent(out) :: iminp
  integer(C_SIZE_T) :: sz
  integer(C_INT), pointer :: imtemp

  sz = size(finp, kind=C_SIZE_T)
  if (sz > 0_c_size_t) then
    imtemp => finp(1)
    iminp%data = c_loc(imtemp)
  else
    iminp%data = c_null_ptr
  end if
  iminp%size = sz
end subroutine
function accumulate(arr) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), dimension(:), intent(in), target :: arr
integer(C_INT) :: fresult
type(SwigArrayWrapper) :: farg1

call SWIGTM_fin_int_Sb__SB_(arr, farg1)
fresult = swigc_accumulate(farg1)
swig_result = fresult
end function

end module

max_stars_repo_name: ExternalRepositories/shroud

max_issues_repo_name: ExternalRepositories/shroud

max_forks_repo_name: ExternalRepositories/shroud


text: Take a look around you. If you see a friend, you are among the luckiest of us. If not, take a deep breath. Reevaluate your life and return to Front Page this Wiki only when you have discovered your true inner spirituality.


text:
      MODULE vcoilpts
      USE stel_constants
      INTEGER, PARAMETER :: ncdim = 100
      INTEGER :: nwdim, nwdim1
      INTEGER, PARAMETER :: nfdim = 5
      END MODULE vcoilpts

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: ORNL-Fusion/LIBSTELL


text:
c               Copyright 1993 Colorado State University
c                       All Rights Reserved

C ... FSFUNC.F

      double precision function fsfunc(minrlp, sorpaf, sorpmx)

      implicit none

c ... Argument declarations
      double precision minrlp, sorpaf, sorpmx

c ... CALLED FROM:  SIMSOM
c ...               PSCHEM
c ...               GROWTH
c ...               SAVARP

c ... This function calculates the fraction of mineral P that is in
c ... solution.  FSFUNC is a functon of the maximum sorption potential
c ... of the soil and sorption affinity.  Both variables are site variables
c ... dependent on the soil mineralogy.

c ... New function by Alister Metherell   March 23, 1992

c ... MINERL  - Mineral phosphorous in the first layer.  (g P / m2)
c ... SORPMX  - The maximum P sorption potential for a soil.  (g P / m2)
c ... SORPAF  - The sorption affinity which controls the slope of the
c ...           sorption curve when mineral P = 0.  Mathematically defined
c ...           as the ratio of mineral P / sorption max when labile P =
c ...           sorbed P. (Range 1.0 - 2.0)

c ... Local variables
      double precision b, c, labile

      c = sorpmx * (2.0 - sorpaf)/2
      b = sorpmx - minrlp + c
      labile = (-b + sqrt(b*b + 4*c*minrlp))/2
      fsfunc = labile / minrlp

      if (fsfunc .lt. 0.0 .or. fsfunc .gt. 1.0) then
        write(*,*) 'Error in calculation of fsfunc, fsfunc = ', fsfunc
      endif

      return
      end


text:
      SUBROUTINE ztsint6 (IFLTAB, CA, NA, CB, NB, CC, NC, CF, NF, JULIS,
     1 ISITIM, JULE, IETIME, INTIN, INTOUT, IHEAD, NH, WRKSPC, VALUES,
     2 NVALS, CUNITS, CTYPE, ISTAT)
C
C
C     Given pathname components, a date/time window and time
C     interval, return interpolated data in users array.
C                      ------------
C
C     Time window may cross the data blocks actually stored.
C
C
C
C     Input:
C        IFLTAB:  Working DSS array used in zopen6 call.  Must be
C                 be dimensioned as INTEGER with 1200 words
C        CA-CF :  Pathname parts ins array form (for harris integer*6).
C        NA-NF :  Respective pathname part lenghts.
C        JULIS :  The julian startin date of the time window.
C        ISITIM:  The starting time of the time window, in minutes.
C        JULE  :  The julian ending date of the time window.
C        IETIME:  The ending time of the time window, in minutes
C        INTIN :  The interval of the data in the array  WRKSPC, in min
C                 5 minutes is the minimum interval.
C        IHEAD :  This is not currently being used.
C        NH    :  This is also not used.
C        WRKSPC:  Work space - area used to actually read data from dss.
C                 Dimension WRKSPC to be at least as large as the number
C                 of data in the time window specified.
C        NVALS :  Dimension limit of array  VALUES.
C
C
C     Output:
C        NVALS:   The number of interpolated data retrieved.  Note that
C                 this is also and input argument.
C        VALUES:  The interpolated data retrieved.
C        CUNITS:  Character string returning the units of the data.
C                 CUNITS must be declared CHARACTER*8
C        CTYPE:   Character string returning the type of the data
C                 (e.g., PER-AVER).  CTYPE must be declared CHARACTER*8
C        ISTAT:   Integer status parameter, indicating the
C                 successfullness of the retrieval.
C                 ISTAT = 0  All ok.
C                 ISTAT = 1  Some missing data (still ok)
C                 ISTAT = 2  Missing data blocks, but some data found
C                 ISTAT = 3  Combination of 1 and 2 (some data found)
C                 ISTAT = 4  No data found, although a pathname was read
C                 ISTAT = 5  No pathname(s) found
C                 ISTAT > 9  Illegal call to zrrts6
C
      INCLUDE 'zdssmz.h'
C
C
C
C        Parameters ----------------------------------------------------
      CHARACTER*(*) CA, CB, CC, CF
      CHARACTER*(*) CTYPE, CUNITS
      CHARACTER CE*16
      INTEGER IFLTAB(*), IHEAD(*)
      REAL WRKSPC(*), VALUES(*)
      CHARACTER PATH*400
C
C       Local Variables -----------------------------------------------
C
      CHARACTER CJTYPE(4)*8
      CHARACTER CDATE*20, CTIME*4
C
      DATA CJTYPE/'PER-AVER','PER-CUM ','INST-VAL','INST-CUM'/
C
C ----------------------------------------------------------------------
C
C
C
C     Copy starting date and time so they can be changed.
      JULS = JULIS
      ISTIME = ISITIM
C
      IF (MLEVEL.LT.9) GO TO 60
      WRITE(MUNIT,30)JULS,ISTIME,JULE,IETIME,INTIN,INTOUT
 30   FORMAT(T20,'----- ENTERING ztsint6-----------------',
     */,T10,'STARTING JULIAN DATE, TIME = ',3X,2I8,
     */,T10,'ENDING   JULIAN DATE, TIME = ',3X,2I8,
     */,T10,'TIME INTERVAL IN, OUT = ',2I7)
      WRITE(MUNIT,40)CA(1:4),CB(1:4),CC(1:4),CF(1:4),NA,NB,NC,NF
 40   FORMAT(T10,'FIRST 4 CHARS AND CHAR LENGTH OF CA, CB, CC, CF',
     */,T15,4(A4,4X),/,T15,4(I4,4X))
 60   CONTINUE
C
C     Check that time does not start or end on 00 mins (illegal time).
      IF (ISTIME.NE.0) GO TO 70
          ISTIME = 1440
          JULS = JULS - 1
 70   IF (IETIME.NE.0) GO TO 75
          IETIME = 1440
          JULE = JULE - 1
 75   CONTINUE
C
C     Initialize work space.
      DO 77 I=1,35
      WRKSPC(I) = -901.0
 77   CONTINUE
C
      NSD = 1
      ISTAT = 0
C
C     Calculate the number of periods asked for.
      NED = NOPERS(INTOUT,0,JULS,ISTIME,JULE,IETIME) + 1
      IF (NED.GT.NVALS) GO TO 910
C
C
C     Get alpha value of interval.
      I = 2
      CALL zgintl6 ( INTIN, CE, NUMDAT, I)
      IF (I.LT.0) GO TO 900
C
C     Get start date and time of first data to be read.
C
      INCBAK = MAX(INTOUT,INTIN)
      ILARGE = 1440
      INCBAK = MAX(INCBAK,ILARGE)
      IDUMMY = INCTIM(INCBAK,0,-2,JULS,ISTIME,JULSD,ISTIM)
C
      CALL zofset6(JULSD,ISTIM,INCBAK,1,IOFS)
      IF (MLEVEL.GE.8) WRITE (MUNIT, 79) JULSD, ISTIM, INCBAK, IOFS
 79   FORMAT (T8,'ztsint6: After zofset6, JULSD:',I8,',  ISTIM:',I8,/,
     * T10,'INCBAK:',I8,',  IOFS:',I8)
      INCFWD = MIN (INTOUT,INTIN)
      IDUMMY = INCTIM (INCFWD, 0, 1, JULSD, 0, JUL, ISTIM)
      IF ((ISTIM.EQ.1440).AND.(INCFWD.NE.1440)) THEN
      JUL = JUL + 1
      ISTIM = 5
      ENDIF
C
C
C     Determine the number of periods (a block may have an
C     incorrect amount of data.)
C     NDW = Number in WRKSPC, ND = number in VALUES array.
      NDW = NOPERS(INTIN,0,JULSD,0,JULE,IETIME) + 2
      ND = NOPERS(INTOUT,0,JULSD,0,JULE,IETIME)
C
      I = M2IHM(ISTIM,CTIME)
      CALL JULDAT (JUL,104,CDATE,N)
C     Form pathname.
      CALL zfpn(CA,NA,CB,NB,CC,NC,CDATE,9,CE,8,CF,NF,PATH,NPATH)
C     Read data from DSS file.
      CALL zrrts(IFLTAB,PATH(1:NPATH),CDATE,CTIME,NDW,WRKSPC(36),
     +CUNITS,CTYPE,IWOFS,ISTAT)
C
      IF (ISTAT.GE.9) GO TO 920
      IF (ISTAT.GE.4) GO TO 800
C
C     Determine data type.
C
      DO 80 I=1,4
      IFLAG = I
      IF (CTYPE .EQ. CJTYPE(IFLAG)) GO TO 85
 80   CONTINUE
C
C     Could not recognize type - Fatal Error.
      GO TO 930
C
 85   CONTINUE
C
C
      OFSETW = FLOAT(IWOFS)/FLOAT(INTIN)
C
      IF (MLEVEL.GE.8) WRITE (MUNIT, 86) CTYPE, IFLAG, IWOFS, INTIN,
     * OFSETW
 86   FORMAT ('  ztsint6:  Data Type -',A,'-   Flag:',I4,/,
     * T10,'IWOFS:',I8,',  INTIN:',I8,',   OFSETW:',F12.5)
C
      CALL zintbk6 (VALUES, NSD, NED, JULS, ISTIME, INTOUT, ND,
     * WRKSPC, NDW, JULSD, INTIN, OFSETW, ISTAT, IFLAG)
      IF (ISTAT.GE.9) GO TO 920
C
      NVALS = NED
C
C     Locate any missing data, and indicate in istat.
      DO 100 I=1,NVALS
      IF (VALUES(I).NE.-901.0.AND.VALUES(I).NE.-902.0)GO TO 110
 100  CONTINUE
      ISTAT=4
       GO TO 800
C
C     Locate first and last data.
 110  IFIRST = I
C
      DO 120 J=1,NVALS
      I = NVALS - J + 1
C
 120  IF(VALUES(I).NE.-901.0.AND.VALUES(I).NE.-902.0) GO TO 130
      GO TO 140
 130  CONTINUE
      K = IFIRST
C
C
 140  DO 150 I=K,NVALS
      IF (VALUES(I).EQ.-901.0)GO TO 160
 150  CONTINUE
       GO TO 800
C
 160  IF(ISTAT.EQ.2)ISTAT=3
      IF(ISTAT.EQ.0)ISTAT=1
C
 800  CONTINUE
      IF(MLEVEL.GE.9)WRITE(MUNIT,810)NVALS
 810  FORMAT(T20,'----- EXITING ztsint6, NO DATA READ =',I6)
C
      RETURN
C
C     --- Error Statements ---
C
 900   WRITE(MUNIT,901) INTIN
 901  FORMAT(/,'*****  ERROR-ztsint6-INTERVAL NON-STANDARD',I10/)
      ISTAT = 12
      GO TO 800
C
 910  WRITE(MUNIT,911)NED,NVALS
 911  FORMAT(//,'*****  ERROR - ztsint6 - NO PERIODS.GT.NVALS '  ,
     +/, 'NO OF PERIODS = ' ,I6,'   NVALS = ',I6)
      ISTAT = 11
      GO TO 800
C
 920  WRITE(MUNIT,921)ISTAT
 921  FORMAT('***** - ztsint6 - ERROR.  UNABLE TO',
     *' RETRIEVE DATA - STATUS = ',I3)
      GO TO 800
C
C
 930  CONTINUE
      CALL CHRLNB(PATH,N)
      WRITE(MUNIT,931)PATH(1:N),CUNITS,CTYPE
 931  FORMAT(//' ***** ERROR - ztsint6;  COULD NOT RECOGNIZE DATA TYPE',
     */,' PATHNAME = ',A,/,' UNITS = ',A,'  TYPE = ',A,/)
      ISTAT = 35
      GO TO 800
C
      END

max_stars_repo_name: sharkcooper/hec-dss

max_issues_repo_name: sharkcooper/hec-dss

max_forks_repo_name: sharkcooper/hec-dss


text:
      subroutine leshouche_inc_chooser()
c For a given nFKSprocess, it fills the c_leshouche_inc common block with the
c leshouche.inc information
      implicit none
      include 'nexternal.inc'
      include 'genps.inc'
      integer i,j,k
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
     &     icolup(2,nexternal,maxflow),niprocs
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
      logical firsttime
      data firsttime /.true./
      include 'leshouche_decl.inc'
      common/c_leshouche_idup_d/ idup_d
      save mothup_d, icolup_d, niprocs_d

c
      if (maxproc_used.gt.maxproc) then
         write (*,*) 'ERROR in leshouche_inc_chooser: increase maxproc',
     &        maxproc,maxproc_used
         stop
      endif
      if (maxflow_used.gt.maxflow) then
         write (*,*) 'ERROR in leshouche_inc_chooser: increase maxflow',
     &        maxflow,maxflow_used
         stop
      endif

      if (firsttime) then
        call read_leshouche_info(idup_d,mothup_d,icolup_d,niprocs_d)
        firsttime = .false.
      endif

      niprocs=niprocs_d(nFKSprocess)
      do j=1,niprocs
         do i=1,nexternal
            IDUP(i,j)=IDUP_D(nFKSprocess,i,j)
            MOTHUP(1,i,j)=MOTHUP_D(nFKSprocess,1,i,j)
            MOTHUP(2,i,j)=MOTHUP_D(nFKSprocess,2,i,j)
         enddo
      enddo
c
      do j=1,maxflow_used
         do i=1,nexternal
            ICOLUP(1,i,j)=ICOLUP_D(nFKSprocess,1,i,j)
            ICOLUP(2,i,j)=ICOLUP_D(nFKSprocess,2,i,j)
         enddo
      enddo
c
      return
      end

      subroutine read_leshouche_info(idup_d,mothup_d,icolup_d,niprocs_d)
C read the various information from the configs_and_props_info.dat file
      implicit none
      include "nexternal.inc"
      integer itmp_array(nexternal)
      integer i,j,k,l
      character *200 buff
      include 'leshouche_decl.inc'
      include 'nFKSconfigs.inc'
      include 'fks_info.inc'
      include 'born_maxamps.inc'
      integer idup(nexternal,maxproc)
      integer mothup(2,nexternal,maxproc)
      integer icolup(2,nexternal,maxflow)
      include 'born_leshouche.inc'
      if (fks_configs.eq.1) then
         if (pdg_type_d(1,fks_i_d(1)).eq.-21) then
c SPECIAL for [LOonly=QCD] process. Simply use the information from the
c born_leshouche.inc file.
            do j=1,maxproc
               do i=1,nexternal-1
                  idup_d(1,i,j)=idup(i,j)
                  mothup_d(1,1,i,j)=mothup(1,i,j)
                  mothup_d(1,2,i,j)=mothup(2,i,j)
               enddo
               idup_d(1,nexternal,j)=-21
               mothup_d(1,1,nexternal,j)=mothup(1,fks_j_d(1),j)
               mothup_d(1,2,nexternal,j)=mothup(2,fks_j_d(1),j)
            enddo
            do j=1,maxflow
               do i=1,nexternal-1
                  icolup_d(1,1,i,j)=icolup(1,i,j)
                  icolup_d(1,2,i,j)=icolup(2,i,j)
               enddo
               icolup_d(1,1,nexternal,j)=-99999 ! should not be used
               icolup_d(1,2,nexternal,j)=-99999
            enddo
            niprocs_d(1)=maxproc_used
            return
         endif
      endif

      open(unit=78, file='leshouche_info.dat', status='old')
      do while (.true.)
        read(78,'(a)',end=999) buff
        if (buff(:1).eq.'#') cycle
        if (buff(:1).eq.'I') then
        ! idup
        ! I  i   j   id1 ..idn -> IDUP_D(i,k,j)=idk
          read(buff(2:),*) i,j,(itmp_array(k),k=1,nexternal)
          do k=1,nexternal
            idup_d(i,k,j)=itmp_array(k)
          enddo
          niprocs_d(i)=j
        else if (buff(:1).eq.'M') then
        ! idup
        ! I  i   j   l   id1 ..idn -> MOTHUP_D(i,j,k,l)=idk
          read(buff(2:),*) i,j,l,(itmp_array(k),k=1,nexternal)
          do k=1,nexternal
            mothup_d(i,j,k,l)=itmp_array(k)
          enddo
        else if (buff(:1).eq.'C') then
        ! idup
        ! I  i   j   l   id1 ..idn -> ICOLUP_D(i,j,k,l)=idk
          read(buff(2:),*) i,j,l,(itmp_array(k),k=1,nexternal)
          do k=1,nexternal
            icolup_d(i,j,k,l)=itmp_array(k)
          enddo
        endif
      enddo
 999  continue
      close(78)

      return
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
C     =================
      SUBROUTINE TRKFIL
C     =================
C     This routine reconstructs the event tracking vertex,
C     the number of tracks, the number of vertex tracks and the
C     chi**2/ndf for the vertex fit. At present the vertex
C     parameters are stored for `CTD-only' and `regular mode',
C     all other track parameters (cal track matching etc.) are
C     determined from `regular tracking' mode.
C
C     A.Quadt, May  1997.
C              July 1997, CTD and global tracking modes added.
C     O.Ruske, January 1998.
C              added the three angles needed to determine the
C              vertex distribution: theta_perpf,b,r
C              (from CTD only)  see note book II, page 81
C     ----------------------------------------------------------

      IMPLICIT NONE

#include "partap.inc"
#include "vctvtx.inc"
#include "vctpar.inc"
#include "vctrhl.inc"
#include "zdskey.inc"
#include "zrevt.inc"
#include "vcevtctd.inc"
#include "vcevtreg.inc"
C
c#include "zvttup.inc"
#include "common.inc"
C
      INTEGER     NO_USE, Ierr, NrGoodCells, VzErr, I
      DATA        NO_USE / -999. /
      REAL        COMPF, COMPB, COMPR, PI
      REAL        VtxtFCAL
      REAL        VtxtFTD(3)

      PI =  ACOS(-1.0)
C ---------------------------------------------
C --- Unpack the tracking tables (CTD only) ---
C ---------------------------------------------

      IF (CouTab(VCEVTCTD).ne.0) CALL VCGETCTD(Ierr)

C -----------------------------------------
C --- reconstructed vertex from VCTRACK ---
C -----------------------------------------
      CHVCC   = NO_USE
      VCT_XVC = NO_USE
      VCT_YVC = NO_USE
      VCT_ZVC = NO_USE
      NVTRKC  = COUTAB(VCTPAR)
      NTRKC   = COUTAB(VCTRHL)

      IF (COUTAB(VCTVTX).GT.0) THEN
            CALL FETTAB(VCTVTX,ID,coutab(vctvtx))
            VCT_XVC = VCTVTX_V(1)
            VCT_YVC = VCTVTX_V(2)
            VCT_ZVC = VCTVTX_V(3)
            IF (VCTVTX_NDF .GT.0.AND.
     &          VCTVTX_CHI2.GT.0) THEN
              CHVCC = VCTVTX_CHI2/FLOAT(VCTVTX_NDF)
            ENDIF
        ENDIF
c      ENDIF

C ------------------------------------------------------------
C --- Unpack the tracking tables (regular mode; incl. FTD) ---
C ------------------------------------------------------------
c      IF (CouTab(VCEVTREG).ne.0) CALL VCGETREG(Ierr)

C -----------------------------------------
C --- reconstructed vertex from VCTRACK ---
C -----------------------------------------
c      CHVC   = NO_USE
c      VCT_XV = NO_USE
c      VCT_YV = NO_USE
c      VCT_ZV = NO_USE
c      IF (IErr.EQ.0) THEN
c        NVTRK  = COUTAB(VCTPAR)
c        NTRK   = COUTAB(VCTRHL)

c        IF (COUTAB(VCTVTX).GT.0) THEN
c            CALL FETTAB(VCTVTX,ID,1)
c            VCT_XV = VCTVTX_V(1)
c            VCT_YV = VCTVTX_V(2)
c            VCT_ZV = VCTVTX_V(3)
c            IF (VCTVTX_NDF .GT.0.AND.
c     &          VCTVTX_CHI2.GT.0) THEN
c              CHVC = VCTVTX_CHI2/FLOAT(VCTVTX_NDF)
c            ENDIF
c        ENDIF
c      ENDIF

C ----------------------------------------
C --- FCAL timing vertex from R.Pawlak ---
C ----------------------------------------
      FCAL_VTX  = VtxtFCAL(NrGoodCells,VzErr)
      FCAL_VTXE = VzErr

C ----------------------------------------
C --- FTD vertex from M.Eckert
C ----------------------------------------
      FTD_VTX  = -999.9
      FTD_VSEG = 0
      CALL FTDVTX(VTXTFTD,FTD_VSEG)
      FTD_VTX = VtxtFTD(3)

      RETURN
      END

ext: fpp