text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                                !!
!!  This file is part of SciFT project                                               !!
!!  Copyright (c) 2015-2015 Nestor F. Aguirre (nfaguirrec@gmail.com)                 !!
!!                                                                                !!
!!  Redistribution and use in source and binary forms, with or without               !!
!!  modification, are permitted provided that the following conditions are met:      !!
!!                                                                                !!
!!  1. Redistributions of source code must retain the above copyright notice, this   !!
!!     list of conditions and the following disclaimer.                              !!
!!  2. Redistributions in binary form must reproduce the above copyright notice,     !!
!!     this list of conditions and the following disclaimer in the documentation     !!
!!     and/or other materials provided with the distribution.                        !!
!!  3. Neither the name of the copyright holders nor the names of its contributors   !!
!!     may be used to endorse or promote products derived from this software         !!
!!     without specific prior written permission.                                    !!
!!                                                                                !!
!!  The copyright holders provide no reassurances that the source code provided      !!
!!  does not infringe any patent, copyright, or any other intellectual property      !!
!!  rights of third parties.  The copyright holders disclaim any liability to any    !!
!!  recipient for claims brought against recipient by any third party for            !!
!!  infringement of that parties intellectual property rights.                       !!
!!                                                                                !!
!!  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND  !!
!!  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED    !!
!!  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE           !!
!!  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR  !!
!!  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES   !!
!!  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;     !!
!!  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND      !!
!!  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT       !!
!!  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS    !!
!!  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                     !!
!!                                                                                !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module FourierTransform3D_
	use FFTW3_
	use GOptions_
	use Math_
	use Grid3D_
	use CNFunction3D_
	use RNFunction3D_
	use RandomUtils_
	use FourierTransform_
	use FourierTransform2D_
	implicit none
	private

	public :: &
		FourierTransform3D_omegaGrid, &
		FourierTransform3D_xyzGrid, &
		FourierTransform3D_phase, &
		FourierTransform3D_iphase, &
		FourierTransform3D_shift, &
		FourierTransform3D_ishift, &
		FourierTransform3D_plan, &
		FourierTransform3D_execute, &
		FourierTransform3D_destroyPlan, &
		FourierTransform3D_dft, &
		FourierTransform3D_idft, &
		FourierTransform3D_fft, &
		FourierTransform3D_ifft, &
		FourierTransform3D_nft, &
		FourierTransform3D_inft, &
		FourierTransform3D_test

	type, public :: FourierTransform3D
		class(CNFunction3D), pointer, private :: iFunc
		integer(8), private :: planF, planB

		integer :: nPoints(3)
		type(Grid3D) :: x
		type(Grid3D) :: omega

		contains
			procedure :: init
			final :: destroy
			procedure :: str
			procedure :: show
			procedure :: execute
! 			procedure :: filter
	end type FourierTransform3D

    interface FourierTransform3D_omegaGrid
		module procedure FourierTransform3D_omegaGridFromData
		module procedure FourierTransform3D_omegaGridFromXGrid
    end interface FourierTransform3D_omegaGrid

    interface FourierTransform3D_xyzGrid
		module procedure FourierTransform3D_xyzGridFromData
		module procedure FourierTransform3D_xyzGridFromOmegaGrid
    end interface FourierTransform3D_xyzGrid

    interface FourierTransform3D_phase
        module procedure FourierTransform3D_phase_realArray
        module procedure FourierTransform3D_phase_complexArray
        module procedure FourierTransform3D_phase_RNFunction3D
        module procedure FourierTransform3D_phase_CNFunction3D
    end interface FourierTransform3D_phase

    interface FourierTransform3D_iphase
        module procedure FourierTransform3D_iphase_realArray
        module procedure FourierTransform3D_iphase_complexArray
        module procedure FourierTransform3D_iphase_CNFunction3D
    end interface FourierTransform3D_iphase

    interface FourierTransform3D_shift
        module procedure FourierTransform3D_shift_Grid
        module procedure FourierTransform3D_shift_RNFunction3D
        module procedure FourierTransform3D_shift_CNFunction3D
    end interface FourierTransform3D_shift

    interface FourierTransform3D_ishift
        module procedure FourierTransform3D_ishift_Grid
        module procedure FourierTransform3D_ishift_CNFunction3D
    end interface FourierTransform3D_ishift

    interface FourierTransform3D_plan
        module procedure FourierTransform3D_plan_Array
        module procedure FourierTransform3D_plan_CNFunction3D
    end interface FourierTransform3D_plan

    interface FourierTransform3D_dft
        module procedure FourierTransform3D_dft_CArray
    end interface FourierTransform3D_dft

    interface FourierTransform3D_idft
        module procedure FourierTransform3D_idft_Array
    end interface FourierTransform3D_idft

    interface FourierTransform3D_fft
        module procedure FourierTransform3D_fft_CNFunction3D
    end interface FourierTransform3D_fft

    interface FourierTransform3D_ifft
        module procedure FourierTransform3D_ifft_CNFunction3D
    end interface FourierTransform3D_ifft

    interface FourierTransform3D_nft
        module procedure FourierTransform3D_nft_CNFunction3D
    end interface FourierTransform3D_nft

    interface FourierTransform3D_inft
        module procedure FourierTransform3D_inft_CNFunction3D
    end interface FourierTransform3D_inft

	contains

	!>
	!! @brief Constructor
	!! @todo En algunas ocasiones es importante hacer la transfomada de no todos los puntos
	!!
	subroutine init( this, iFunc, domain, oFunc )
		class(FourierTransform3D) :: this
		class(CNFunction3D), target, intent(in) :: iFunc
		integer(8), intent(in), optional :: domain
		class(CNFunction3D), target, optional, intent(in) :: oFunc

		integer(8) :: effDomain

		integer :: nx, ny, nz

		if( associated(this.iFunc) ) nullify(this.iFunc)
		this.iFunc => iFunc
		this.nPoints = iFunc.nPoints() ! << Hay que cambiar para hacer la transformada de algunos puntos

		nx = this.nPoints(1)
		ny = this.nPoints(2)
		nz = this.nPoints(3)

		effDomain = FourierTransform_SPATIAL_DOMAIN
		if( present(domain) ) effDomain = domain

		if( present(oFunc) ) call GOptions_error( "oFunc parameter is not implamented yet", "FourierTransform3D.init()" )

		call dfftw_plan_dft_3d( this.planF, nx, ny, nz, this.iFunc.fArray, this.iFunc.fArray, FFTW_FORWARD, FFTW_ESTIMATE )
		call dfftw_plan_dft_3d( this.planB, nx, ny, nz, this.iFunc.fArray, this.iFunc.fArray, FFTW_BACKWARD, FFTW_ESTIMATE )

		if( effDomain == FourierTransform_SPATIAL_DOMAIN ) then

			this.x = this.iFunc.xyzGrid ! << Hay que cambiar para hacer la transformada de algunos puntos
			this.omega = FourierTransform3D_omegaGrid( this.x )

		else if( effDomain == FourierTransform_FREQUENCY_DOMAIN ) then

			call GOptions_warning( "FourierTransform_FREQUENCY_DOMAIN have not been tested yet", "FourierTransform.init()" )
			this.omega = this.iFunc.xyzGrid ! << Hay que cambiar para hacer la transformada de algunos puntos
			this.x = FourierTransform3D_xyzGrid( this.omega )

		end if
	end subroutine init

	!>
	!! @brief Destructor
	!!
	subroutine destroy( this )
		type(FourierTransform3D) :: this

		! Si los activo el programa produce:
		! 	forrtl: error (76): Abort trap signal
! 		call dfftw_destroy_plan( this.planF )
! 		call dfftw_destroy_plan( this.planB )

		nullify(this.iFunc)

		this.nPoints = -1
	end subroutine destroy

	!>
	!! @brief
	!!
	function str( this ) result( output )
		class(FourierTransform3D) :: this
		character(len=200) :: output

		integer :: fmt
		character(len=200) :: strBuffer

		output = ""

		output = trim(output)//"<FourierTransform3D:"

! 		output = trim(output)//"min="
! 		fmt = int(log10(this.min+1.0))+1
! 		write(strBuffer, "(f<fmt+7>.6)") this.min
! 		output = trim(output)//trim(strBuffer)
!
! 		output = trim(output)//",max="
! 		fmt = int(log10(this.max+1.0))+1
! 		write(strBuffer, "(f<fmt+7>.6)") this.max
! 		output = trim(output)//trim(strBuffer)
!
! 		output = trim(output)//",h="
! 		fmt = int(log10(this.h+1.0))+1
! 		write(strBuffer, "(f<fmt+7>.6)") this.h
! 		output = trim(output)//trim(strBuffer)
!
! 		output = trim(output)//",size="
! 		fmt = int(log10(float(this.size+1)))+1
! 		write(strBuffer, "(i<fmt>)") this.size
! 		output = trim(output)//trim(strBuffer)

		output = trim(output)//">"
	end function str

	!>
	!! @brief
	!!
	subroutine show( this, unit )
		class(FourierTransform3D) :: this
		integer, optional, intent(in) :: unit

		integer :: effunit

		if( present(unit) ) then
			effunit = unit
		else
			effunit = 6
		end if

		write(effunit,"(a)") trim(this.str())
	end subroutine show

	!>
	!! @brief
	!!
	subroutine execute( this, sgn, sync, shift )
		class(FourierTransform3D) :: this
		integer, intent(in) :: sgn
		logical, optional, intent(in) :: sync
		logical, optional, intent(in) :: shift

		logical :: effSync
		logical :: effShift

		effSync = .false.
		if( present(sync) ) effSync = sync

		effShift = .false.
		if( present(shift) ) effShift = shift

		if( sgn == FourierTransform_FORWARD ) then

			call dfftw_execute( this.planF )

			if( effSync ) then
				this.iFunc.xyzGrid = this.omega

				if( effShift ) then
					call FourierTransform3D_phase( this.iFunc )
					call FourierTransform3D_shift( this.iFunc )
				end if
			end if

		else if ( sgn == FourierTransform_BACKWARD ) then

			if( effSync ) then
				if( effShift ) then
					call FourierTransform3D_ishift( this.iFunc )
					call FourierTransform3D_phase( this.iFunc )
				end if
			end if

			call dfftw_execute( this.planB )

			if( effSync ) then
				this.iFunc.xyzGrid = this.x
			end if

			this.iFunc = this.iFunc/real(this.nPoints(1)*this.nPoints(2)*this.nPoints(3),8)

		else
			call GOptions_error( "Bad value for sgn", "FourierTransform3D.execute()" )
		end if
	end subroutine execute

	!>
	!! @brief
	!!
	function FourierTransform3D_omegaGridFromData( n, h, order ) result( output )
		integer, intent(in) :: n(3)
		real(8), optional, intent(in) :: h(3)
		integer, optional, intent(in) :: order
		type(Grid3D) :: output

		integer :: i

		do i=1,3
			output.component(i) = FourierTransform_omegaGrid( n(i), h(i), order=order )
		end do
	end function FourierTransform3D_omegaGridFromData

	!>
	!! @brief
	!!
	function FourierTransform3D_omegaGridFromXGrid( xyzGrid, order ) result( output )
		type(Grid3D), intent(in) :: xyzGrid
		integer, optional, intent(in) :: order
		type(Grid3D) :: output

		output = FourierTransform3D_omegaGridFromData( xyzGrid.nPointsVec(), xyzGrid.stepSize(), order=order )
	end function FourierTransform3D_omegaGridFromXGrid

	!>
	!! @brief
	!!
	function FourierTransform3D_xyzGridFromData( n, h, order ) result( output )
		integer, intent(in) :: n(3)
		real(8), intent(in) :: h(3)
		integer, optional, intent(in) :: order
		type(Grid3D) :: output

		integer :: i

		do i=1,3
			output.component(i) = FourierTransform_xGrid( n(i), h(i), order=order )
		end do
	end function FourierTransform3D_xyzGridFromData

	!>
	!! @brief
	!!
	function FourierTransform3D_xyzGridFromOmegaGrid( omegaGrid, order ) result( output )
		type(Grid3D), intent(in) :: omegaGrid
		integer, optional, intent(in) :: order
		type(Grid3D) :: output

		real(8) :: dx(3)

		! dx = 2*pi/n/dp
		dx = 2.0_8*MATH_PI/real(omegaGrid.nPoints(),8)/omegaGrid.stepSize()

		output = FourierTransform3D_xyzGridFromData( omegaGrid.nPoints(), dx, order=order )
	end function FourierTransform3D_xyzGridFromOmegaGrid

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_phase_realArray( iArray, oArray )
		real(8) :: iArray(:)
		real(8), optional :: oArray(:)

! 		integer :: i, j, n, p0
!
! 		n = size( iArray )
! 		p0 = floor( n/2.0_8 )+2
!
! 		i = 1
! 		do j=p0,n
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**j
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**j
! 			end if
!
! 			i = i+1
! 		end do
!
! 		do j=1,p0-1
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**j
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**j
! 			end if
!
! 			i = i+1
! 		end do
	end subroutine FourierTransform3D_phase_realArray

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_phase_complexArray( iArray, oArray )
		complex(8) :: iArray(:)
		complex(8), optional :: oArray(:)

! 		integer :: i, j, n, p0
!
! 		n = size( iArray )
! 		p0 = floor( n/2.0_8 )+2
!
! 		i = 1
! 		do j=p0,n
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**j
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**j
! 			end if
!
! 			i = i+1
! 		end do
!
! 		do j=1,p0-1
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**j
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**j
! 			end if
!
! 			i = i+1
! 		end do
	end subroutine FourierTransform3D_phase_complexArray

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_iphase_realArray( iArray, oArray )
		real(8) :: iArray(:)
		real(8), optional :: oArray(:)

! 		integer :: i, j, n, p2
!
! 		n = size( iArray )
! 		p2 = n-floor(n/2.0_8)
!
! 		i = 1
! 		do j=p2,n
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			end if
!
! 			i = i+1
! 		end do
!
! 		do j=1,p2-1
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			end if
!
! 			i = i+1
! 		end do
	end subroutine FourierTransform3D_iphase_realArray

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_iphase_complexArray( iArray, oArray )
		complex(8) :: iArray(:)
		complex(8), optional :: oArray(:)

! 		integer :: i, j, n, p2
!
! 		n = size( iArray )
! 		p2 = n-floor(n/2.0_8)
!
! 		i = 1
! 		do j=p2,n
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			end if
!
! 			i = i+1
! 		end do
!
! 		do j=1,p2-1
! 			if( present(oArray) ) then
! 				oArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			else
! 				iArray(i) = iArray(i)*(-1.0_8)**(j-1)
! 			end if
!
! 			i = i+1
! 		end do
	end subroutine FourierTransform3D_iphase_complexArray

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_phase_RNFunction3D( iFunc, oFunc )
		class(RNFunction3D) :: iFunc
		class(RNFunction3D), optional :: oFunc

		integer :: i

		if( iFunc.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should not be equally spaced", &
				"FourierTransform3D_phase_RNFunction3D()" &
			)
		end if

		if( present(oFunc) ) then
			oFunc = iFunc

			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_phase( oFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_phase( oFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_phase( oFunc.fArray(:,:,i) )
			end do
		else
			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_phase( iFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_phase( iFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_phase( iFunc.fArray(:,:,i) )
			end do
		end if
	end subroutine FourierTransform3D_phase_RNFunction3D

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_phase_CNFunction3D( iFunc, oFunc )
		class(CNFunction3D) :: iFunc
		class(CNFunction3D), optional :: oFunc

		integer :: i

		if( iFunc.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should not be equally spaced", &
				"FourierTransform3D_phase_CNFunction3D()" &
			)
		end if

		if( present(oFunc) ) then
			oFunc = iFunc

			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_phase( oFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_phase( oFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_phase( oFunc.fArray(:,:,i) )
			end do
		else
			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_phase( iFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_phase( iFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_phase( iFunc.fArray(:,:,i) )
			end do
		end if
	end subroutine FourierTransform3D_phase_CNFunction3D

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_iphase_CNFunction3D( iFunc, oFunc )
		class(CNFunction3D) :: iFunc
		class(CNFunction3D), optional :: oFunc

! 		if( .not. iFunc.isEquallyspaced() ) then
! 			call GOptions_warning( &
! 				"Grid should be equally spaced", &
! 				"FourierTransform3D_ishift_Grid()" &
! 			)
! 		end if
!
! 		if( present(oFunc) ) then
! 			oFunc = iFunc
!
! 			call FourierTransform3D_iphase_complexArray( oFunc.yArray )
! 		else
! 			call FourierTransform3D_iphase_complexArray( iFunc.yArray )
! 		end if
	end subroutine FourierTransform3D_iphase_CNFunction3D

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_shift_Grid( iGrid, oGrid )
		class(Grid3D) :: iGrid
		class(Grid3D), optional :: oGrid

		integer :: i

		if( iGrid.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should not be equally spaced", &
				"FourierTransform3D_shift_Grid()" &
			)
		end if

		if( present(oGrid) ) then
			oGrid = iGrid

			do i=1,3
				call FourierTransform_shift( oGrid.component(i) )
			end do
		else
			do i=1,3
				call FourierTransform_shift( iGrid.component(i) )
			end do
		end if

	end subroutine FourierTransform3D_shift_Grid

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_ishift_Grid( iGrid, oGrid )
		class(Grid3D) :: iGrid
		class(Grid3D), optional :: oGrid

		integer :: i

		if( .not. iGrid.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should be equally spaced", &
				"FourierTransform3D_ishift_Grid()" &
			)
		end if

		if( present(oGrid) ) then
			oGrid = iGrid

			do i=1,3
				call FourierTransform_ishift( oGrid.component(i) )
			end do
		else
			do i=1,3
				call FourierTransform_ishift( iGrid.component(i) )
			end do
		end if

	end subroutine FourierTransform3D_ishift_Grid

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_shift_RNFunction3D( iFunc, oFunc )
		class(RNFunction3D) :: iFunc
		class(RNFunction3D), optional :: oFunc

		integer :: i

		if( iFunc.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should not be equally spaced", &
				"FourierTransform3D_shift_Grid()" &
			)
		end if

		if( present(oFunc) ) then
			oFunc = iFunc

			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_shift( oFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_shift( oFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_shift( oFunc.fArray(:,:,i) )
			end do

			call FourierTransform3D_shift( oFunc.xyzGrid )
		else
			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_shift( iFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_shift( iFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_shift( iFunc.fArray(:,:,i) )
			end do

			call FourierTransform3D_shift( iFunc.xyzGrid )
		end if

	end subroutine FourierTransform3D_shift_RNFunction3D

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_shift_CNFunction3D( iFunc, oFunc )
		class(CNFunction3D) :: iFunc
		class(CNFunction3D), optional :: oFunc

		integer :: i

		if( iFunc.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should not be equally spaced", &
				"FourierTransform3D_shift_Grid()" &
			)
		end if

		if( present(oFunc) ) then
			oFunc = iFunc

			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_shift( oFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_shift( oFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_shift( oFunc.fArray(:,:,i) )
			end do

			call FourierTransform3D_shift( oFunc.xyzGrid )
		else
			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_shift( iFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_shift( iFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_shift( iFunc.fArray(:,:,i) )
			end do

			call FourierTransform3D_shift( iFunc.xyzGrid )
		end if

	end subroutine FourierTransform3D_shift_CNFunction3D

	!>
	!! @brief
	!!
	subroutine FourierTransform3D_ishift_CNFunction3D( iFunc, oFunc )
		class(CNFunction3D) :: iFunc
		class(CNFunction3D), optional :: oFunc

		integer :: i

		if( .not. iFunc.isEquallyspaced() ) then
			call GOptions_warning( &
				"Grid should be equally spaced", &
				"FourierTransform3D_ishift_Grid()" &
			)
		end if

		if( present(oFunc) ) then
			oFunc = iFunc

			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_ishift( oFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_ishift( oFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_ishift( oFunc.fArray(:,:,i) )
			end do

			call FourierTransform3D_ishift( oFunc.xyzGrid )
		else
			do i=1,iFunc.nPoints(1)
				call FourierTransform2D_ishift( iFunc.fArray(i,:,:) )
			end do

			do i=1,iFunc.nPoints(2)
				call FourierTransform2D_ishift( iFunc.fArray(:,i,:) )
			end do

			do i=1,iFunc.nPoints(3)
				call FourierTransform2D_ishift( iFunc.fArray(:,:,i) )
			end do

			call FourierTransform3D_ishift( iFunc.xyzGrid )
		end if

	end subroutine FourierTransform3D_ishift_CNFunction3D

	!>
	!! @brief
	!!
	function FourierTransform3D_plan_Array( iArray, sgn, oArray ) result( plan )
		complex(8) :: iArray(:,:,:)
		integer, intent(in) :: sgn
		complex(8), optional :: oArray(:,:,:)
		integer(8) :: plan

		if( present(oArray) ) then
			call dfftw_plan_dft_3d( plan, size(iArray,dim=1), size(iArray,dim=2), size(iArray,dim=3), iArray, oArray, sgn, FFTW_ESTIMATE )
		else
			call dfftw_plan_dft_3d( plan, size(iArray,dim=1), size(iArray,dim=2), size(iArray,dim=3), iArray, iArray, sgn, FFTW_ESTIMATE )
		end if
	end function FourierTransform3D_plan_Array

	!>
	!! @brief
	!!
	function FourierTransform3D_plan_CNFunction3D( iFunc, sgn, oFunc ) result( plan )
		type(CNFunction3D) :: iFunc
		integer, intent(in) :: sgn
		type(CNFunction3D), optional :: oFunc
		integer(8) :: plan

		if( present(oFunc) ) then
			call dfftw_plan_dft_3d( plan, iFunc.nPoints(1), iFunc.nPoints(2), iFunc.nPoints(3), iFunc.fArray, oFunc.fArray, sgn, FFTW_ESTIMATE )
		else
			call dfftw_plan_dft_3d( plan, iFunc.nPoints(1), iFunc.nPoints(2), iFunc.nPoints(3), iFunc.fArray, iFunc.fArray, sgn, FFTW_ESTIMATE )
		end if
	end function FourierTransform3D_plan_CNFunction3D

	!>
	!!
	!!
	subroutine FourierTransform3D_execute( plan )
		integer(8), intent(in) :: plan

		call dfftw_execute( plan )
	end subroutine FourierTransform3D_execute

	!>
	!!
	!!
	subroutine FourierTransform3D_destroyPlan( plan )
		integer(8), intent(in) :: plan

		call dfftw_destroy_plan( plan )
	end subroutine FourierTransform3D_destroyPlan

	!>
	!!
	!!
	subroutine FourierTransform3D_dft_CArray( iArray, oArray, sgn )
		complex(8) :: iArray(:,:,:)
		complex(8), optional :: oArray(:,:,:)
		integer, optional :: sgn

		integer :: effSgn

		integer(8) :: plan

		effSgn = FourierTransform_FORWARD
		if( present(sgn) ) effSgn = sgn

		if( present(oArray) ) then
			if( size(iArray) /= size(oArray) ) &
				call GOptions_error( "iArray and oArray have not the same size", "FourierTransform3D_dft_CArray" )

			call dfftw_plan_dft_3d( plan, size(iArray,dim=1), size(iArray,dim=2), size(iArray,dim=3), iArray, oArray, effSgn, FFTW_ESTIMATE )
			call dfftw_execute( plan )
		else
			call dfftw_plan_dft_3d( plan, size(iArray,dim=1), size(iArray,dim=2), size(iArray,dim=3), iArray, iArray, effSgn, FFTW_ESTIMATE )
			call dfftw_execute( plan )
		end if

		call dfftw_destroy_plan( plan )
	end subroutine FourierTransform3D_dft_CArray

	!>
	!!
	!!
	subroutine FourierTransform3D_idft_Array( iArray, oArray, sgn )
		complex(8) :: iArray(:)
		complex(8), optional :: oArray(:)
		integer, optional :: sgn

		stop "FourierTransform3D_idft_Array is not implemented"
! 		integer :: effSgn
!
! 		integer(8) :: plan
!
! 		effSgn = FourierTransform3D_BACKWARD
! 		if( present(sgn) ) effSgn = sgn
!
! 		if( present(oArray) ) then
! 			if( size(iArray) /= size(oArray) ) &
! 				call GOptions_error( "iArray and oArray have not the same size", "FourierTransform3D_idft_Array" )
!
! 			call dfftw_plan_dft_1d( plan, size(iArray), iArray, oArray, effSgn, FourierTransform3DW_ESTIMATE )
! 			call dfftw_execute( plan )
!
! 			oArray = oArray/real( size(iArray), 8 )
! 		else
! 			call dfftw_plan_dft_1d( plan, size(iArray), iArray, iArray, effSgn, FourierTransform3DW_ESTIMATE )
! 			call dfftw_execute( plan )
!
! 			iArray = iArray/real( size(iArray), 8 )
! 		end if
!
! 		call dfftw_destroy_plan( plan )
	end subroutine FourierTransform3D_idft_Array

	!>
	!!
	!!
	function FourierTransform3D_fft_CNFunction3D( iFunc, sgn ) result( oFunc )
		type(CNFunction3D) :: iFunc
		integer, optional :: sgn
		type(CNFunction3D) :: oFunc

		integer :: n(3)
		real(8) :: dx(3), dp(3)
		type(Grid3D) :: xyzGrid

		! @todo Check for checkEquallyspaced
		n = iFunc.nPoints()
		dx = iFunc.xyzGrid.stepSize()
		dp = 2.0_8*Math_PI/dx/real(n,8)

		oFunc = iFunc
		call FourierTransform3D_dft( iFunc.fArray, oFunc.fArray, sgn=FourierTransform_FORWARD )
		oFunc.xyzGrid = FourierTransform3D_omegaGrid( n, dx, order=FourierTransform_SORDER )

		call FourierTransform3D_phase( oFunc )
		call FourierTransform3D_shift( oFunc )

		oFunc = oFunc*dx(1)*dx(2)*dx(3)/sqrt(2.0_8*Math_PI)**3
	end function FourierTransform3D_fft_CNFunction3D

	!>
	!!
	!!
	function FourierTransform3D_ifft_CNFunction3D( iFunc, sgn ) result( oFunc )
		type(CNFunction3D) :: iFunc
		integer, optional :: sgn
		type(CNFunction3D) :: oFunc

		integer :: n(3)
		real(8) :: dx(3), dp(3)
		type(Grid3D) :: xyzGrid

		integer :: i

		! @todo Check for checkEquallyspaced
		n = iFunc.nPoints()
		dp = iFunc.xyzGrid.stepSize()
		dx = 2.0_8*Math_PI/dp/real(n,8)

		oFunc = iFunc

		call FourierTransform3D_ishift( oFunc )
		call FourierTransform3D_phase( oFunc )

		call FourierTransform3D_dft( oFunc.fArray, oFunc.fArray, sgn=FourierTransform_BACKWARD )
		oFunc.xyzGrid = FourierTransform3D_xyzGrid( n, dx, order=FourierTransform_NORDER )

		oFunc = oFunc*dp(1)*dp(2)*dp(3)/sqrt(2.0_8*Math_PI)**3
	end function FourierTransform3D_ifft_CNFunction3D

	!>
	!! @brief
	!!
	function FourierTransform3D_nft_CNFunction3D( iFunc, sgn ) result( oFunc )
		type(CNFunction3D), intent(in) :: iFunc
		integer, optional, intent(in) :: sgn
		type(CNFunction3D) :: oFunc

		integer :: effSgn

		integer :: iw, jw, kw, ix, jx, kx, n(3)
		real(8) :: dx(3), dp(3)
		type(Grid3D) :: xyzGrid

		effSgn = FourierTransform_FORWARD
		if( present(sgn) ) effSgn = sgn

		! @todo Check for checkEquallyspaced
		n = iFunc.nPoints()
		dx = iFunc.xyzGrid.stepSize()
		dp = 2.0_8*Math_PI/dx/real(n,8)

		oFunc = iFunc
		oFunc.xyzGrid = FourierTransform3D_omegaGrid( n, dx, order=FourierTransform_NORDER )
		call oFunc.xyzGrid.show()

		do kw=1,n(3); do jw=1,n(2); do iw=1,n(1)
			oFunc.fArray(iw,jw,kw) = 0.0_8
			do kx=1,n(3); do jx=1,n(2); do ix=1,n(1)
				oFunc.fArray(iw,jw,kw) = oFunc.fArray(iw,jw,kw) &
					+ iFunc.fArray(ix,jx,kx)*exp( real(effSgn,8)*Math_I*sum( iFunc.xyzGrid.at(ix,jx,kx)*oFunc.xyzGrid.at(iw,jw,kw) ) )
			end do; end do; end do
		end do; end do; end do

		oFunc = oFunc*dx(1)*dx(2)*dx(3)/(2.0_8*Math_PI)
	end function FourierTransform3D_nft_CNFunction3D

	!>
	!! @brief
	!!
	function FourierTransform3D_inft_CNFunction3D( iFunc, sgn ) result( oFunc )
		type(CNFunction3D), intent(in) :: iFunc
		integer, optional, intent(in) :: sgn
		type(CNFunction3D) :: oFunc

		integer :: effSgn

		integer :: iw, jw, kw, ix, jx, kx, n(3)
		real(8) :: dx(3), dp(3)
		type(Grid3D) :: xyzGrid

		effSgn = FourierTransform_BACKWARD
		if( present(sgn) ) effSgn = sgn

		! @todo Check for checkEquallyspaced
		n = iFunc.nPoints()
		dp = iFunc.xyzGrid.stepSize()
		dx = 2.0_8*Math_PI/dp/real(n,8)

		oFunc = iFunc
		oFunc.xyzGrid = FourierTransform3D_xyzGrid( n, dx, order=FourierTransform_NORDER )

		do kw=1,n(3); do jw=1,n(2); do iw=1,n(1)
			oFunc.fArray(iw,jw,kw) = 0.0_8
			do kx=1,n(3); do jx=1,n(2); do ix=1,n(1)
				oFunc.fArray(iw,jw,kw) = oFunc.fArray(iw,jw,kw) &
					+ iFunc.fArray(ix,jx,kx)*exp( real(effSgn,8)*Math_I*sum( iFunc.xyzGrid.at(ix,jx,kx)*oFunc.xyzGrid.at(iw,jw,kw) ) )
			end do; end do; end do
		end do; end do; end do

		oFunc = oFunc*dp(1)*dp(2)*dp(3)/(2.0_8*Math_PI)
	end function FourierTransform3D_inft_CNFunction3D

	!>
	!! This is neccesary only for FourierTransform3D_test()
	!! Maxima:
	!!   f(x) := exp(-0.5*x**2);
	!!
	function funcGaussian( x, y, z ) result( output )
		real(8), intent(in) :: x, y, z
		complex(8) :: output

		real(8) :: alpha

		alpha = 5.0_8
		output = exp(-alpha*( x**2 + y**2 + z**2 ))
	end function funcGaussian

	!>
	!! This is neccesary only for FourierTransform3D_test()
	!! Maxima:
	!!   f(x) := exp(-0.5*x**2);
	!!
	!! La convencion de normalizacion de acuerdo con wikipedia es:
    !!   Fourier transform unitary, angular frequency
	!!
	function funcFGaussian( x ) result( output )
		real(8), intent(in) :: x
		complex(8) :: output

		real(8) :: alpha

		alpha = 0.5_8
! 		output = exp(-x**2/4.0_8/alpha)/sqrt(2.0_8*alpha)
		output = exp(-(x-5.0_8)**2/4.0_8/alpha)/sqrt(2.0_8*alpha) + 0.4_8*exp(-(x+2.0_8)**2/10.0_8/alpha)/sqrt(2.0_8*alpha)*sin(2.0*x) &
			+ 0.4_8*Math_I*exp(-(x-2.0_8)**2/10.0_8/alpha)/sqrt(2.0_8*alpha)*cos(2.0*x)
	end function funcFGaussian

	!>
	!! This is neccesary only for FourierTransform3D_test()
	!!
	function funcRectangular( x, y, z ) result( output )
		real(8), intent(in) :: x, y, z
		complex(8) :: output

		real(8) :: a1, a2, a3

		a1 = 0.5_8
		a2 = 1.0_8
		a3 = 1.5_8

		output = Math_ubox( a1*x )*Math_ubox( a2*y )*Math_ubox( a3*z )
	end function funcRectangular

	!>
	!! This is neccesary only for FourierTransform3D_test()
	!!
	function FfuncRectangular( w1, w2, w3 ) result( output )
		real(8), intent(in) :: w1, w2, w3
		complex(8) :: output

		real(8) :: a1, a2, a3

		a1 = 0.5_8
		a2 = 1.0_8
		a3 = 1.5_8

		output = Math_nsinc( w1/(2.0_8*Math_PI*a1) ) &
					*Math_nsinc( w2/(2.0_8*Math_PI*a2) ) &
					*Math_nsinc( w3/(2.0_8*Math_PI*a3) )&
					/sqrt(2.0_8*Math_PI)**3/a1/a2/a3
	end function FfuncRectangular

	!>
	!! @bief test
	!!
	subroutine FourierTransform3D_test()
		type(Grid3D) :: xyzGrid, omegaGrid
		type(CNFunction3D) :: funcA, dFuncA, fftFuncA, fftDFuncA
		type(CNFunction3D) :: funcB, dFuncB, fftFuncB, fftDFuncB
		type(CNFunction3D) :: funcAB, dFuncAB, fftFuncAB, fftDFuncAB
		type(RNFunction3D) :: aSpectrum, pSpectrum
		type(FourierTransform3D) :: fft
		real(8) :: exactValue
		real(8) :: value
		integer :: i, j
		real(8), allocatable :: array(:,:,:)
		complex(8), allocatable :: cArray(:,:,:)
		integer(8) :: planF, planB
		character :: cBuffer

		write(*,*) "----------------------------------"
		write(*,*) " Verifing omega Grid constructors "
		write(*,*) "----------------------------------"

		xyzGrid = FourierTransform3D_omegaGridFromData( [10,11,12], [0.1_8,0.2_8,0.3_8], FourierTransform_NORDER )
		call xyzGrid.show()
		call xyzGrid.save()

		xyzGrid = FourierTransform3D_omegaGridFromData( [10,11,12], [0.1_8,0.2_8,0.3_8], FourierTransform_SORDER )
		call xyzGrid.show()
		call xyzGrid.save()

		write(*,*) "-------------------------------"
		write(*,*) " Verifing xy Grid constructors "
		write(*,*) "-------------------------------"

		xyzGrid = FourierTransform3D_xyzGridFromData( [10,11,12], [0.1_8,0.2_8,0.3_8], FourierTransform_NORDER )
		call xyzGrid.show()
		call xyzGrid.save()

		xyzGrid = FourierTransform3D_xyzGridFromData( [10,11,12], [0.1_8,0.2_8,0.3_8], FourierTransform_SORDER )
		call xyzGrid.show()
		call xyzGrid.save()

		return
		call GOptions_doYouWantToContinue()

! 		write(*,*) "-------------------------------------------"
! 		write(*,*) " FourierTransform3D of a gaussian function "
! 		write(*,*) "-------------------------------------------"
!
! 		call xyzGrid.init( [-2.5_8,-2.5_8], [2.5_8,2.5_8], size=[90,110] )
! 		call funcA.init( xyzGrid, funcRectangular )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save( ".func", format=BLKS_FORMAT )
!
! 		funcB = FourierTransform3D_fft( funcA )
! 		write(*,"(A)") "exact (.Ffunc) = "
! 		call funcB.show()
! 		call funcB.save( ".Ffunc", format=BLKS_FORMAT )
!
! 		call xyzGrid.init( [-54.677675_8,-67.241507_8], [55.920349_8,68.486720_8], size=[90,110] )
! 		call funcB.init( xyzGrid, FfuncRectangular )
! ! 		funcB = FourierTransform3D_nft( funcA )
! 		write(*,"(A)") "exact (.Fexact) = "
! 		call funcB.show()
! 		call funcB.save( ".Fexact", format=BLKS_FORMAT )
!
! ! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! ! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! ! 		call system( "echo splot \"//achar(34)//".Fexact\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".Fexact\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .Ffunc .Fexact" )

! 		write(*,*) "-----------------------------"
! 		write(*,*) " iFourierTransform3D of a gaussian function"
! 		write(*,*) "-----------------------------"
!
! 		call xyzGrid.init( [-2.5_8,-2.5_8], [2.5_8,2.5_8], size=[90,110] )
! 		call funcA.init( xyzGrid, funcRectangular )
! 		funcA = FourierTransform3D_nft( funcA )
! 		write(*,"(A)") "input (.Ffunc) = "
! 		call funcA.show()
! 		call funcA.save( ".Ffunc", format=BLKS_FORMAT )
!
! 		write(*,"(A)") "iFourierTransform3D (.func) = "
! 		funcB = funcA
! 		funcB = FourierTransform3D_ifft( funcB )
! 		call funcB.show()
! 		call funcB.save( ".func", format=BLKS_FORMAT )
!
! 		funcB = FourierTransform3D_inft( funcA )
! 		write(*,"(A)") "exact (.exact) = "
! 		call funcB.show()
! 		call funcB.save( ".exact", format=BLKS_FORMAT )
!
! ! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! ! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! ! 		call system( "echo splot \"//achar(34)//".exact\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".exact\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .Ffunc .func .exact" )
!
! 		write(*,*) "-------------------------------"
! 		write(*,*) " iFourierTransform3D( FourierTransform3D(func) ). CNFunction3D "
! 		write(*,*) "-------------------------------"
!
! 		call xyzGrid.init( [-2.5_8,-2.5_8], [2.5_8,2.5_8], size=[90,110] )
! 		call funcA.init( xyzGrid, funcRectangular )
!
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save( ".func", format=BLKS_FORMAT )
!
! 		write(*,"(A)") "FourierTransform3D (.Ffunc) = "
! 		funcB = FourierTransform3D_fft( funcA )
! 		call funcB.show()
! 		call funcB.save( ".Ffunc", format=BLKS_FORMAT )
!
! 		write(*,"(A)") "iFourierTransform3D (.iFFfunc) = "
! 		funcA = FourierTransform3D_ifft( funcB )
! 		call funcA.show()
! 		call funcA.save( ".iFFfunc", format=BLKS_FORMAT )
!
! ! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! ! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! ! 		call system( "echo splot \"//achar(34)//".iFFfunc\"//achar(34)//" u 1:2:3 w l, \"//achar(34)//"\"//achar(34)//" u 1:2:4 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".iFFfunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .Ffunc .iFFfunc" )

! 		write(*,*) "--------------------------------"
! 		write(*,*) " FourierTransform3D( iFourierTransform3D(Ffunc) ). CNFunction3D "
! 		write(*,*) "--------------------------------"
!
! 		call xyzGrid.init( [-2.5_8,-2.5_8], [2.5_8,2.5_8], size=[90,110] )
! 		call funcA.init( xyzGrid, funcRectangular )
! 		funcA = FourierTransform3D_nft( funcA )
!
! 		write(*,"(A)") "input (.Ffunc) = "
! 		call funcA.show()
! 		call funcA.save( ".Ffunc", format=BLKS_FORMAT )
!
! 		write(*,"(A)") "iFourierTransform3D (.iFfunc) = "
! 		funcA = FourierTransform3D_ifft( funcA )
! 		call funcA.show()
! 		call funcA.save( ".iFfunc", format=BLKS_FORMAT )
!
! 		write(*,"(A)") "FourierTransform3D (.FiFfunc) = "
! 		funcA = FourierTransform3D_fft( funcA )
! 		call funcA.show()
! 		call funcA.save( ".FiFfunc", format=BLKS_FORMAT )
!
! 		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".iFfunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".FiFfunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .Ffunc .iFfunc .FiFfunc" )

		write(*,*) "-----------------------------------------------------"
		write(*,*) " Verifing FourierTransform3D CNFunction3D with plans "
		write(*,*) "-----------------------------------------------------"

		call xyzGrid.init( [-2.5_8,-2.5_8,-2.5_8], [2.5_8,2.5_8,2.5_8], size=[90,110,200] )
		call funcA.init( xyzGrid, funcRectangular )
		write(*,"(A)") "input (.func) = "
		call funcA.show()
		call funcA.save( "func.cube", format=AUTO_FORMAT )

		planF = FourierTransform3D_plan( funcA, FourierTransform_FORWARD )
		planB = FourierTransform3D_plan( funcA, FourierTransform_BACKWARD )

		write(*,"(A)") "FourierTransform3D (.Ffunc) = "
		call FourierTransform3D_execute( planF )
		funcA.xyzGrid = FourierTransform3D_omegaGrid( funcA.xyzGrid, order=FourierTransform_SORDER )
! 		call FourierTransform3D_phase( funcA )
		call FourierTransform3D_shift( funcA )
		call funcA.show()
! 		call funcA.save( ".Ffunc", format=BLKS_FORMAT )
		call funcA.save( "Ffunc.cube", format=AUTO_FORMAT )
		call funcA.save( "Ffunc.n3df", format=AUTO_FORMAT )
		stop
		write(*,"(A)") "iFourierTransform3D (.iFFfunc) = "
		call FourierTransform3D_ishift( funcA )
		call FourierTransform3D_phase( funcA )
		call FourierTransform3D_execute( planB )
		funcA.xyzGrid = FourierTransform3D_xyzGrid( funcA.xyzGrid, order=FourierTransform_NORDER )
		funcA = funcA/real( funcA.nPoints(1)*funcA.nPoints(2), 8 )
		call funcA.show()
		call funcA.save( ".iFFfunc", format=BLKS_FORMAT )

		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
		call system( "echo splot \"//achar(34)//".iFFfunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
		call GOptions_doYouWantToContinue()
		call system( "rm .func .Ffunc .iFFfunc" )

		call FourierTransform3D_destroyPlan( planF )
		call FourierTransform3D_destroyPlan( planB )

! 		write(*,*) "------------------------------------------"
! 		write(*,*) " iFourierTransform3D( FourierTransform3D(func) ). CNFunction3D with plans"
! 		write(*,*) "------------------------------------------"
!
! 		call xyzGrid.init( [-2.5_8,-2.5_8], [2.5_8,2.5_8], size=[90,110] )
! 		call funcA.init( xyzGrid, funcRectangular )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save( ".func", format=BLKS_FORMAT )
!
! 		planF = FourierTransform3D_plan( funcA, FourierTransform_FORWARD )
! 		planB = FourierTransform3D_plan( funcA, FourierTransform_BACKWARD )
!
! 		call FourierTransform3D_execute( planF )
! 		call FourierTransform3D_execute( planB )
! 		funcA = funcA/real( funcA.nPoints(1)*funcA.nPoints(2), 8 )
!
! 		call funcA.show()
! 		call funcA.save( ".iFFfunc", format=BLKS_FORMAT )
!
! 		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call system( "echo splot \"//achar(34)//".iFFfunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .iFFfunc" )
!
! 		call FourierTransform3D_destroyPlan( planF )
! 		call FourierTransform3D_destroyPlan( planB )

		write(*,*) "-----------------------------------------"
		write(*,*) " Verifing FourierTransform3D CNFunction3D oriented object"
		write(*,*) "-----------------------------------------"

		call xyzGrid.init( [-2.5_8,-2.5_8,-2.5_8], [2.5_8,2.5_8,2.5_8], size=[90,110,130] )
		call funcA.init( xyzGrid, funcRectangular )
		write(*,"(A)") "input (.func) = "
		call funcA.show()
		call funcA.save( ".func", format=BLKS_FORMAT )

		call fft.init( funcA, FourierTransform_SPATIAL_DOMAIN )

		call fft.execute( FourierTransform_FORWARD, sync=.true., shift=.true. )
		write(*,"(A)") "FourierTransform3D (.Ffunc) = "
		call funcA.show()
		call funcA.save( ".Ffunc", format=BLKS_FORMAT )

		call fft.execute( FourierTransform_BACKWARD, sync=.true., shift=.true. )
		write(*,"(A)") "iFourierTransform3D (.iFFfunc) = "
		call funcA.show()
		call funcA.save( ".iFFfunc", format=BLKS_FORMAT )

		call system( "echo splot \"//achar(34)//".func\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
		call system( "echo splot \"//achar(34)//".Ffunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
		call system( "echo splot \"//achar(34)//".iFFfunc\"//achar(34)//" u 1:2:3 w l | gnuplot -p" )
		call GOptions_doYouWantToContinue()
		call system( "rm .func .Ffunc .iFFfunc" )

! 		write(*,*) "--------------------------------------"
! 		write(*,*) " Second derivative via FourierTransform3D with plans"
! 		write(*,*) "--------------------------------------"
!
! 		call xyzGrid.init( -3.0_8*Math_PI, 3.0_8*Math_PI, 1001 )
! 		omegaGrid = FourierTransform3D_omegaGrid( xyzGrid )
! 		call funcA.init( xyzGrid, funcTest )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save(".func")
!
! 		planF = FourierTransform3D_plan( funcA, FourierTransform3D_FORWARD )
! 		planB = FourierTransform3D_plan( funcA, FourierTransform3D_BACKWARD )
!
! 		call FourierTransform3D_execute( planF )
!
! 		funcA.yArray = ( Math_I*fft.omega )**2*funcA.yArray
!
! 		write(*,"(A)") "iFourierTransform3D (.dfunc) = "
! 		call FourierTransform3D_execute( planB )
! 		funcA = funcA/real( funcA.nPoints(), 8 )
! 		call funcA.show()
! 		call funcA.save(".dfunc")
!
! 		call funcB.init( xyzGrid, d2funcTest )
! 		write(*,"(A)") "exact (.exact) = "
! 		call funcB.show()
! 		call funcB.save(".exact")
!
! ! 		call system( "echo plot \"//achar(34)//".func\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call system( "echo plot \"//achar(34)//".dfunc\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call system( "echo plot \"//achar(34)//".exact\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .dfunc .exact" )
!
! 		call FourierTransform3D_destroyPlan( planF )
! 		call FourierTransform3D_destroyPlan( planB )
!
! 		write(*,*) "------------------------------------------------------"
! 		write(*,*) " Second derivative via FourierTransform3D CNFunction3D oriented object"
! 		write(*,*) "------------------------------------------------------"
!
! 		call xyzGrid.init( -3.0_8*Math_PI, 3.0_8*Math_PI, 1001 )
! 		call funcA.init( xyzGrid, funcTest )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save(".func")
!
! 		call fft.init( funcA, FourierTransform_SPATIAL_DOMAIN )
!
! 		call fft.execute( FourierTransform3D_FORWARD )
!
! 		funcA.yArray = ( Math_I*fft.omega )**2*funcA.yArray
!
! 		call fft.execute( FourierTransform3D_BACKWARD )
!
! 		write(*,"(A)") "iFourierTransform3D (.dfunc) = "
! 		call funcA.show()
! 		call funcA.save(".dfunc")
!
! 		call funcB.init( xyzGrid, d2funcTest )
! 		write(*,"(A)") "exact (.exact) = "
! 		call funcB.show()
! 		call funcB.save(".exact")
!
! ! 		call system( "echo plot \"//achar(34)//".func\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call system( "echo plot \"//achar(34)//".dfunc\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call system( "echo plot \"//achar(34)//".exact\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .dfunc .exact" )
!
! 		write(*,*) "----------"
! 		write(*,*) " Spectrum "
! 		write(*,*) "----------"
!
! ! 		call xyzGrid.init( -3.0_8*Math_PI, 3.0_8*Math_PI, 1001 )
! 		call xyzGrid.init( -15.0_8, 15.0_8, 1001 )
! 		call funcA.init( xyzGrid, funcTest )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save(".func")
!
! 		aSpectrum = FourierTransform3D_powerSpectrum( funcA )
!
! 		write(*,"(A)") "aspec (.asfunc) = "
! 		call aSpectrum.show()
! 		call aSpectrum.save(".asfunc")
!
! 		pSpectrum = FourierTransform3D_phaseSpectrum( funcA )
!
! 		write(*,"(A)") "pspec (.psfunc) = "
! 		call pSpectrum.show()
! 		call pSpectrum.save(".psfunc")
!
! 		call system( "echo plot \"//achar(34)//".func\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! 		call system( "echo plot \"//achar(34)//".asfunc\"//achar(34)//" u 1:2 w l lw 1.5 | gnuplot -p" )
! 		call system( "echo plot \"//achar(34)//".psfunc\"//achar(34)//" u 1:2 w l lw 1.5 | gnuplot -p" )
! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .asfunc .psfunc" )
!
! 		write(*,*) "--------------------------------------------------------"
! 		write(*,*) " Filtering by amplitude. FourierTransform3D CNFunction3D oriented object"
! 		write(*,*) "--------------------------------------------------------"
!
! 		call xyzGrid.init( -3.0_8*Math_PI, 3.0_8*Math_PI, 1001 )
! 		call funcA.init( xyzGrid, funcTestWithNoise )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save(".func")
!
! 		call fft.init( funcA, FourierTransform_SPATIAL_DOMAIN )
!
! 		call fft.execute( FourierTransform3D_FORWARD )
!
! 		where( abs( funcA.yArray ) <= 20.0_8 ) funcA.yArray = 0.0_8
!
! 		call fft.execute( FourierTransform3D_BACKWARD )
!
! 		write(*,"(A)") "iFourierTransform3D (.ffunc) = "
! 		call funcA.show()
! 		call funcA.save(".ffunc")
!
! ! 		call system( "echo plot \"//achar(34)//".func\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call system( "echo plot \"//achar(34)//".ffunc\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .ffunc" )
!
! 		write(*,*) "--------------------------------------------------------"
! 		write(*,*) " Filtering by frequency. FourierTransform3D CNFunction3D oriented object"
! 		write(*,*) "--------------------------------------------------------"
!
! 		call xyzGrid.init( -3.0_8*Math_PI, 3.0_8*Math_PI, 1001 )
! 		call funcA.init( xyzGrid, funcTestWithNoise )
! 		write(*,"(A)") "input (.func) = "
! 		call funcA.show()
! 		call funcA.save(".func")
!
! 		call fft.init( funcA, FourierTransform_SPATIAL_DOMAIN )
!
! 		call fft.execute( FourierTransform3D_FORWARD )
!
! 		do i=1,funcA.nPoints()
! 			if( abs( fft.omega(i) ) > 10.0_8 ) then
! 				funcA.yArray(i) = 0.0_8
! 			end if
! 		end do
!
! 		call fft.execute( FourierTransform3D_BACKWARD )
!
! 		write(*,"(A)") "iFourierTransform3D (.ffunc) = "
! 		call funcA.show()
! 		call funcA.save(".ffunc")
!
! ! 		call system( "echo plot \"//achar(34)//".func\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call system( "echo plot \"//achar(34)//".ffunc\"//achar(34)//" u 1:2 w l lw 1.5, \"//achar(34)//"\"//achar(34)//" u 1:3 w l lw 1.5 | gnuplot -p" )
! ! 		call GOptions_doYouWantToContinue()
! 		call system( "rm .func .ffunc" )

	end subroutine FourierTransform3D_test

end module FourierTransform3D_

max_stars_repo_name: nfaguirrec/scift


text:
	SUBROUTINE ND_GGDS ( gdsarr, cprj, kx, ky, grdout, rnvblk,
     +			    iret )
C************************************************************************
C* ND_GGDS								*
C*									*
C* This subroutine takes the GDS information and constructs a 		*
C* navigation block for a GEMPAK grid file.				*
C*									*
C* ND_GGDS ( GDSARR, CPRJ, KX, KY, GRDOUT, RNVBLK, IRET )		*
C*									*
C* Input parameters:							*
C*	GDSARR (10)	REAL		GRIB GDS projection info	*
C*									*
C* Output parameters:							*
C*	CPRJ		CHAR*		Grid projection			*
C*	KX		INTEGER		Number of points in x dir	*
C*	KY		INTEGER		Number of points in y dir	*
C*	GRDOUT  (4)	REAL		Grid corners			*
C*	RNVBLK  (13)	REAL		Grid navigation block		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = invalid navigation	*
C**									*
C* Log:									*
C* T. Piper/SAIC	10/02	Modified from NAGGDS for GRIB2		*
C************************************************************************
	REAL		gdsarr (*), grdout (*), rnvblk (*)
	CHARACTER*(*)	cprj
C*
	LOGICAL		angflg
C------------------------------------------------------------------------
	iret = 0
C
C*	Set the projection.
C
	IF  ( NINT ( gdsarr (1) ) .eq. 0 )  THEN
	    cprj = 'CED'
	ELSE IF  ( NINT ( gdsarr (1) ) .eq. 10 )  THEN
	    cprj = 'MER'
	ELSE IF  ( NINT ( gdsarr (1) ) .eq. 20 )  THEN
	    cprj = 'STR'
	ELSE IF  ( NINT ( gdsarr (1) ) .eq. 30 )  THEN
	    cprj = 'LCC'
	ELSE
	    iret = -4
	    RETURN
	END IF
C
C*	Set the number of columns and rows.
	kx = NINT ( gdsarr (2) )
	ky = NINT ( gdsarr (3) )
C
C*	Set the lat/lon values of the corners.
C
	grdout (1) = gdsarr (4)
	grdout (2) = gdsarr (5)
	grdout (3) = gdsarr (6)
	grdout (4) = gdsarr (7)
C
C*	Set the projection angles.
C
	angle1 = gdsarr ( 8)
	angle2 = gdsarr ( 9)
	angle3 = gdsarr (10)
	angflg = .true.
C
C*	Fill navigation block.
C
	CALL GR_MNAV  ( cprj, kx, ky, grdout (1), grdout (2),
     +			grdout (3), grdout (4), angle1, angle2,
     +			angle3, angflg, rnvblk, ier )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine makbar ( procesdef, notot , syname, nocons, constants,
     +                    nopa     , paname, nofun , funame, nosfun   ,
     +                    sfname   , nodisp, diname, novelo, vename   ,
     +                    noq3     , laswi , no_act, actlst, noinfo   ,
     +                    nowarn   , nerror)

      ! Checks which processes can be activated

      use dlwq_data
      use processet
      use timers       !   performance timers

      implicit none

      ! declaration of arguments

      type(procespropcoll)      :: procesdef       ! all processes
      integer                   :: notot           ! number of substances
      character(len=*)          :: syname(*)       ! substance name
      integer                   :: nocons          ! number of constants
      type(t_dlwq_item)   , intent(inout) :: constants       !< delwaq constants list
      integer                   :: nopa            ! number of parameters
      character(len=*)          :: paname(*)       ! parameter names
      integer                   :: nofun           ! number of functions
      character(len=*)          :: funame(*)       ! function names
      integer                   :: nosfun          ! number of segment functions
      character(len=*)          :: sfname(*)       ! segment function names
      integer                   :: nodisp          ! number of dispersions
      character(len=*)          :: diname(*)       ! dispersion names
      integer                   :: novelo          ! number of velocities
      character(len=*)          :: vename(*)       ! velocity names
      integer                   :: noq3            ! number of exhcanges in third direction
      logical                   :: laswi           ! active only switch
      integer                   :: no_act          ! number of active processes
      character(len=*)          :: actlst(*)       ! active processes names
      integer                   :: noinfo          ! number of informative messages
      integer                   :: nowarn          ! number of warnings
      integer                   :: nerror          ! number of errors

      ! local decalarations

      integer                   :: nproc           ! number of processes
      integer                   :: iproc           ! loop counter processes
      integer                   :: iproc2          ! second loop counter processes
      type(procesprop), pointer :: proc1           ! process description
      type(procesprop), pointer :: proc2           ! description second process
      integer                   :: ivalip          ! index variable in pmsa
      character(len=20)         :: valnam          ! variable name
      character(len=50)         :: valtxt          ! variable description
      integer                   :: iflux           ! index flux
      integer                   :: i_input         ! index input item
      integer                   :: ioutput         ! index output item
      integer                   :: iact            ! index in active list
      integer                   :: imolev          ! monitoring level
      character(len=100)        :: line            ! line buffer for output
      integer, parameter        :: mismax = 50     ! maximum number of missing variables per process
      integer                   :: nmis            ! actual number of missing variables
      integer                   :: imis            ! index number of missing variables
      character(len=20)         :: misnam(mismax)  ! name missing variables
      character(len=50)         :: mistxt(mismax)  ! description missing variables
      logical                   :: dhrmis          ! function checks on missing value
      logical                   :: iok             ! indicates if its ok
      integer                   :: i_star                ! index of * in name
      integer(4) :: ithndl = 0
      if (timon) call timstrt( "makbar", ithndl )

      write ( line , '(a)' ) '# Determining which processes can be switched on'
      call monsys( line , 2 )
      line = ' '
      call monsys( line , 2 )
      call getmmo(imolev)

      nproc = procesdef%cursize
      do iproc = 1 , nproc

         proc1 => procesdef%procesprops(iproc)
         if (proc1%sfrac_type .eq. SFRAC_DUPLICATED_ORIGINAL) then
            !
            ! prevent the original version of duplicated processes from showing up in
            ! warning list
            !
            proc1%linvok = .false.
            cycle
         endif

         nmis = 0
         call dhswtr( proc1%swtransp, noq3  , iok )
         if (.not. iok ) then
            write (line,'(4a)') ' Input for [',proc1%name,'] ', proc1%text(1:50)
            call monsys( line , 4 )
            write (line,'(a)') ' process for different model dimensions'
            call monsys( line , 4 )
            proc1%linvok = .false.
            goto 550
         endif

         write (line,'(4a)') ' Input for [',proc1%name,'] ', proc1%text(1:50)
         call monsys( line , 4 )
         proc1%linvok = .true.

         ! check input items

         do i_input = 1 , proc1%no_input

            if ( proc1%input_item(i_input)%type .eq. IOTYPE_SEGMENT_INPUT ) then

               if ( .not.dhrmis(proc1%input_item(i_input)%actdef) .and. imolev .lt. 7 ) cycle
               valnam = proc1%input_item(i_input)%name
               valtxt = proc1%input_item(i_input)%item%text
               write(line,'(4a)') '       [',valnam,'] ',valtxt
               call monsys( line , 7 )

   10          continue

               ! specified in input?

               call valpoi ( notot , nopa     , nosfun, syname, nocons,
     +                       nofun , constants, paname, funame, sfname,
     +                       valnam, ivalip   , line  )

               ! output of previous proces ? , is this switched on , switch it on

               if ( ivalip .eq. -1 ) then

                  ! output of previous proces ? , including own pre-workspace

                  do iproc2 = 1 , iproc

                     proc2 => procesdef%procesprops(iproc2)

                     ! flux of process

                     call zoekio ( valnam, proc2%no_fluxoutput, proc2%fluxoutput, 20, iflux)
                     if ( iflux .gt. 0 ) then
                        if ( proc2%linvok ) then
                           write ( line , '(3a)' ) '       Using flux from proces [',proc2%name,']'
                           ivalip = -2
                           exit
                        else
                           write ( line , '(3a)' ) '  NOT  Using flux from proces [',proc2%name,']'
                        endif
                     endif

                     ! output variable of process

                     call zoekio ( valnam, proc2%no_output, proc2%output_item, 20, ioutput, IOTYPE_SEGMENT_OUTPUT)
                     if ( ioutput .gt. 0 ) then
                        if ( proc2%linvok ) then
                           write ( line , '(3a)' ) '       Using output from proces [',proc2%name,']'
                           ivalip = -3
                           exit
                        else
                           write(line,'(3a)') '  NOT  Using output from proces [',proc2%name,']'
                        endif
                     endif

                  enddo
               endif
               if ( ivalip .eq. -1 ) then

                  ! if this is a fraction input then first look for the generic name

                  i_star = index(valnam,'*')
                  if ( i_star .gt. 1 ) then
                     valnam(i_star:) = ' '
                     write(line,'(a)') '       fraction specific input not found, trying generic name'
                     call monsys( line , 7 )
                     write(line,'(4a)') '       [',valnam,'] ',valtxt
                     call monsys( line , 7 )
                     goto 10
                  endif

                  if ( dhrmis(proc1%input_item(i_input)%actdef) )then
                     write(line,'(a)') '       not found'
                     proc1%linvok = .false.
                     nmis = nmis + 1
                     if ( nmis .le. mismax ) then
                        misnam(nmis) = valnam
                        mistxt(nmis) = valtxt
                     endif
                  else
                     write(line,'(a)') '       can use default value'
                  endif
               endif
               call monsys( line , 7 )

            endif
         enddo

         ! on exchange level, seperate loop, then report file stays compatible

         do i_input = 1 , proc1%no_input

            if ( proc1%input_item(i_input)%type .eq. IOTYPE_EXCHANG_INPUT ) then

               if ( .not.dhrmis(proc1%input_item(i_input)%actdef) .and. imolev .lt. 7 ) cycle
               valnam = proc1%input_item(i_input)%name
               valtxt = proc1%input_item(i_input)%item%text
               write(line,'(4a)') '       [',valnam,'] ',valtxt
               call monsys( line , 7 )

   20          continue

               ! specified in input?

               call vxlpoi ( nocons, nofun , nodisp, novelo, constants,
     +                       funame, diname, vename, valnam, ivalip   ,
     +                       line  )

               ! output of previous proces ? , is this switched on , switch it on

               if ( ivalip .eq. -1 ) then

                  ! output of previous proces ? , including own pre-workspace

                  do iproc2 = 1 , iproc

                     proc2 => procesdef%procesprops(iproc2)

                     ! xoutput variable of process

                     call zoekio ( valnam, proc2%no_output, proc2%output_item, 20, ioutput, IOTYPE_EXCHANG_OUTPUT)
                     if ( ioutput .gt. 0 ) then
                        if ( proc2%linvok ) then
                           write ( line , '(3a)' ) '       Using output from proces [',proc2%name,']'
                           ivalip = -3
                           exit
                        else
                           write(line,'(3a)') '  NOT  Using output from proces [',proc2%name,']'
                        endif
                     endif

                  enddo
               endif
               if ( ivalip .eq. -1 ) then

                  ! if this is a fraction input then first look for the generic name

                  i_star = index(valnam,'*')
                  if ( i_star .gt. 1 ) then
                     valnam(i_star:) = ' '
                     write(line,'(a)') '       fraction specific input not found, trying generic name'
                     call monsys( line , 7 )
                     write(line,'(4a)') '       [',valnam,'] ',valtxt
                     call monsys( line , 7 )
                     goto 20
                  endif

                  if ( dhrmis(proc1%input_item(i_input)%actdef) )then
                     write(line,'(a)') '       not found'
                     proc1%linvok = .false.
                     nmis = nmis + 1
                     if ( nmis .le. mismax ) then
                        misnam(nmis) = valnam
                        mistxt(nmis) = valtxt
                     endif
                  else
                     write(line,'(a)') '       can use default value'
                  endif
               endif
               call monsys( line , 7 )
            endif
         enddo

  550    continue
         if ( laswi ) then
            call zoek ( proc1%name, no_act, actlst, 20    , iact  )
            if ( iact .gt. 0 ) then
               if ( proc1%linvok ) then
                  proc1%active = .true.
                  write(line,'(a)') '   Process is activated'
                  call monsys( line , 4 )
               else
                  nowarn = nowarn + 1
                  write(line,'(a)') '   WARNING : activated process can NOT be switched on'
                  call monsys( line , 4 )
                  do imis = 1 , min(nmis,mismax)
                     write(line,'(4a)') '   Not found:[',misnam(imis),'] ', mistxt(imis)
                     call monsys( line , 4 )
                  enddo
                  if ( nmis .gt. mismax ) then
                     write(line,'(a)') '   and more ...'
                     call monsys( line , 4 )
                  endif
               endif
            else
               proc1%linvok = .false.
               write(line,'(a)') '   Process is not activated'
               call monsys( line , 4 )
            endif
         else
            if ( proc1%linvok ) then
               write(line,'(a)') '   Proces can be switched on'
            else
               write(line,'(a)') '   Proces can NOT be switched on'
            endif
            call monsys( line , 4 )
         endif
         line = ' '
         call monsys( line , 4 )
      enddo
      call monsys( line , 2 )

      if (timon) call timstop( ithndl )
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      subroutine mkherm3(fun,x,nx,y,ny,z,nz,fherm)
C
C  create a data set for Hermite interpolation, from evaluation of
C  function and derivatives.  Function of 3 indep. coordinates.
C
C  input:
C
      external fun                      ! passed real function(x,y,z)
      real x(nx)                        ! x coordinate array  (1st dim)
      real y(ny)                        ! y coordinate array  (2nd dim)
      real z(nz)                        ! z coordinate array  (3rd dim)
C
C  the passed function fun must have the interface:
C
C        real function <name>(x,y,z,dfdx,dfdy,dfdz,
C               d2fdxdy,d2fdxdz,d2fdydz,d3fdxdydz)
C  where x,y,z are input, the function returns the function value,
C  and the arguments dfdx, dfdy and dfdz return as output the function
C  derivative at the point (x,y,z).
C
C  output:
C
      real fherm(0:7,nx,ny,nz)          ! function data & derivatives
C
C  fherm(0,i,j,k) = function value f at x(i),y(j),z(k)
C  fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k)
C  fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k)
C  fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k)
C  fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k)
C  fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k)
C  fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k)
C  fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k)
C
C----------------------------
C
      do iz=1,nz
         do iy=1,ny
            do ix=1,nx
               fherm(0,ix,iy,iz)=
     >            fun(x(ix),y(iy),z(iz),dfdx,dfdy,dfdz,
     >                d2fdxdy,d2fdxdz,d2fdydz,d3fdxyz)
               fherm(1,ix,iy,iz)=dfdx
               fherm(2,ix,iy,iz)=dfdy
               fherm(3,ix,iy,iz)=dfdz
               fherm(4,ix,iy,iz)=d2fdxdy
               fherm(5,ix,iy,iz)=d2fdxdz
               fherm(6,ix,iy,iz)=d2fdydz
               fherm(7,ix,iy,iz)=d3fdxyz
            enddo
         enddo
      enddo
C
      return
      end

max_stars_repo_name: ORNL-Fusion/RFSciDAC-testing

max_issues_repo_name: ORNL-Fusion/RFSciDAC-testing

max_forks_repo_name: ORNL-Fusion/RFSciDAC-testing


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Sparse Stoichiometric Data Structures File
!
! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
!
! File                 : kpp_achem_gas_StoichiomSP.f90
! Time                 : Thu Jun 11 11:33:46 2015
! Working directory    : /gpfsm/dnb32/adarmeno/models/geos-5/heracles-UNSTABLE-DEV_MICROPHYSICS/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphV
! Equation file        : kpp_achem_gas.kpp
! Output root filename : kpp_achem_gas
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MODULE kpp_achem_gas_StoichiomSP

  USE kpp_achem_gas_Precision
  PUBLIC
  SAVE

! Row-compressed sparse data for the Jacobian of reaction products JVRP

  INTEGER, PARAMETER, DIMENSION(6) :: CROW_JVRP = (/ &
       1,  2,  3,  4,  5,  6 /)

  INTEGER, PARAMETER, DIMENSION(5) :: ICOL_JVRP = (/ &
       1,  1,  1,  2,  3 /)

  INTEGER, PARAMETER, DIMENSION(5) :: IROW_JVRP = (/ &
       1,  2,  3,  4,  5 /)

!  Stoichiometric Matrix in Compressed Column Sparse Format

  INTEGER, PARAMETER, DIMENSION(6) :: CCOL_STOICM = (/ &
       1,  3,  6,  8, 10, 11 /)

  INTEGER, PARAMETER, DIMENSION(10) :: IROW_STOICM = (/ &
       1,  2,  1,  2,  4,  1,  2,  2,  5,  3 /)

  INTEGER, PARAMETER, DIMENSION(10) :: ICOL_STOICM = (/ &
       1,  1,  2,  2,  2,  3,  3,  4,  4,  5 /)

  REAL(kind=dp), PARAMETER, DIMENSION(10) :: STOICM = (/ &
       -1.000000e+00_dp,  1.000000e+00_dp,  -1.000000e+00_dp,  7.500000e-01_dp,  2.500000e-01_dp, &
       -1.000000e+00_dp,  1.000000e+00_dp,  -1.000000e+00_dp,  1.000000e+00_dp,  -1.000000e+00_dp /)

END MODULE kpp_achem_gas_StoichiomSP


! { dg-do compile }
! { dg-options "-O3 -fwhole-file -fdump-tree-optimized" }
! { dg-add-options bind_pic_locally }
!
! Check that inlining of functions declared BEFORE usage works.
! If yes, then the dump does not contain a call to F().
!

INTEGER FUNCTION f()
  f = 42
END FUNCTION

PROGRAM main
  INTEGER :: a, f
  a = f()
  print *, a, f()
END PROGRAM

! { dg-final { scan-tree-dump-times "= f \\(\\)" 0 "optimized" } }


text:
      PROGRAM PGKEX03
C
C  Use inquiry functions to implement a subroutine (DRWTXT below)
C  that puts out text using Normalized Device Coordinates (NDC)
C  for positioning.  Invoke Autograph to draw a linear/log plot
C  and then call DRWTXT two times to label the Autograph plot.
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)

      DIMENSION X(100),Y(100)
C
C  Open GKS, open and activate a workstation.
C
      CALL GOPKS (IERRF, ISZDUM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C  Define a small color table.
C
      CALL GSCR(IWKID,0,1.,1.,1.)
      CALL GSCR(IWKID,1,0.,0.,1.)
      CALL GSCR(IWKID,2,1.,0.,0.)
C
C  Generate a straight line with 100 points.
C
      DO 1 I=1,100
        X(I) = I
        Y(I) = 10.*I
  1   CONTINUE
C
C  Use SET to define normalization transformation 1 with linear
C  scaling in the X direction and log scaling in the Y direction.
C
      CALL SET(.15,.85,.15,.85,1.,100.,10.,1000.,2)
C
C  Set line color to red.
C
      CALL GSPLCI(2)
C
C  Initialize the AUTOGRAPH entry EZXY so that the frame is not
C  advanced and the Y axis is logarithmic.  Turn off axis labels.
C
      CALL DISPLA(2,0,2)
      CALL ANOTAT(' ',' ',0,0,0,0)
C
C  Output the polyline (X,Y) using EZXY.
C
      CALL EZXY(X,Y,100,' ')
C
C  Put out a couple of labels (DRWTXT uses NDC space).
C
      CALL DRWTXT(.50,.07,'The X axis is linear',-7,.025,0.)
      CALL DRWTXT(.07,.50,'The Y axis is log',-7,.025,90.)
C
C  Terminate the picture, deactivate and close the CGM workstation,
C  and close GKS.
C
      CALL FRAME
      CALL GDAWK(IWKID)
      CALL GCLWK(IWKID)
      CALL GCLKS
C
      STOP
      END
      SUBROUTINE DRWTXT(X,Y,TXT,IFNT,CHGT,ANG)
C
C  This subroutine draws the text string in TXT at position (X,Y) using
C  font IFNT with character height CHGT (specified in NDC) and text
C  angle ANG degrees.  The position (X,Y) is in NDC. This subroutine
C  is isolated from any GKS attribute settings in the calling program
C  by using inquiry functions to save all settings on entry and restore
C  all settings on exit.  The text is aligned as (center, half) and is
C  drawn in the foreground color.
C
      INTEGER   ERRIND
      DIMENSION CLRECT(4)
      CHARACTER*(*) TXT
      DATA PI  /3.1415927/
C
C  Inquire and save the state of all attributes that will be used in
C  this subroutine.  These will be restored on exit.
C
C   Clipping
      CALL GQCLIP(ERRIND,ICLIPO,CLRECT)
C
C   Character up vector.
      CALL GQCHUP(ERRIND,CHUPXO,CHUPYO)
C
C   Text alignment.
      CALL GQTXAL(ERRIND,ILNHO,ILNVO)
C
C   Text font.
      CALL GQTXFP(ERRIND,ITXFO,ITXPO)
C
C   Character height.
      CALL GQCHH(ERRIND,CHHO)
C
C  Get and save the existing normalization transformation information,
C  including the log scaling parameter..
C
      CALL GETSET(XV1,XV2,YV1,YV2,XW1,XW2,YW1,YW2,LS)
C
C  Use NDC space for drawing TXT.
C
      CALL GSELNT(0)
C
C  Define the text font.
C
      CALL GSTXFP(IFNT,2)
C
C  Set the character height.
C
      CALL GSCHH(CHGT)
C
C  Set the text alignment to (center, half).
C
      CALL GSTXAL(2,3)
C
C  Select the foreground color.
C
      CALL GSTXCI(1)
C
C  Define the character up vector in accordance with ANG (recall that
C  the up vector is perpendicular to the text path).
C
      RANG = (ANG+90.)*(2.*PI/360.)
      CALL GSCHUP(COS(RANG),SIN(RANG))
C
C  Draw the text string in TXT.
C
      CALL GTX(X,Y,TXT)
C
C  Restore the original normalization transformation.
C
      CALL SET(XV1,XV2,YV1,YV2,XW1,XW2,YW1,YW2,LS)
C
C  Restore all other attributes.
C
      CALL GSCLIP(ICLIPO)
      CALL GSCHUP(CHUPXO,CHUPYO)
      CALL GSTXAL(ILNHO,ILNVO)
      CALL GSTXFP(ITXFO,ITXPO)
      CALL GSCHH (CHHO)
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
Lower Freeborn Hall is a fancy way of saying the basement under Freeborn Hall.  Editors at The California Aggie and KDVS employees practically live down here.

Lower Freeborn is home to various campus student organizations organizations and services:
  AS Papers
  KDVS
 The Pantry
  KDVS Recordings
  The California Aggie
  El Rodeo
  SPACs business office  if youre looking for SPAC, you probably want to be on the 4th floor of the MU.
  The Memorial Union MU business office (MUBO)
  DPI Graphics & Design Services
  A Learn to Dance dance studio called The Space
  Lockers for student organizations
  Project Compost
  CoHo and Freeborn Hall storage rooms and industrial dishwashing room
  The cleanest bathrooms on campus. Seriously, I go in there barefoot.

Lower Freeborn can be accessed via:
  The south side of the building
  The north side of the building
  The MU Games Area  head down past the bowling lockers
  The Memorial Union MU  elevators and stairs

The ceiling of Lower Freeborn Hall, in between the basement and the first level, buildering can be accessed if:
  You have a flathead screwdriver and a short ladder
  No one else is around
  Youre not afraid of asbestos

Live in the KDVS Hallway
  The hallways of Lower Freeborn have also been used as a recording space by KDVS Engineers on occasion because of the natural reverb and interesting acoustics.
  Notable Performaces Include: Christina Carter and Gown (Officially released under the title Weve by Digitalis), Peter Walker, Jack Rose, Pumice, Tom Carter and Joel Pickard, Mi and Lau, and Dead Western

Departed:
  University Haircutters

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
enum E { }
constant a = E.X

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
      program hello
        print "(a)", "Hello, World!"
      end program hello

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
** Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
** See https://llvm.org/LICENSE.txt for license information.
** SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

*   Miscellaneous intrinsics which may appear in parameter statements
*   used in non-parameter context

	program ff01
	complex c1, c2
	double complex cd1, cd2
	double precision d1, d2, d3, d4, d5, d6, d7, d8, d9
	character ch
	logical l1, l2, l3, l4

	parameter (N=40)

	integer result(N)
	integer expect(N)
	common result, expect

	i1 = abs(2.5)	! tests 1 - 2
	i2 = abs(-2.5)

	x1 = abs(2.5)	! tests 3 - 6
	x2 = abs(-4.5e+3)
	d1 = abs(2.5)
	d2 = abs(-5.4999d+4)
	i3 = nint(x1)
	i4 = nint(-x2)
	i5 = idnint(d1)
	i6 = idnint(-d2)

	c1 = (2.54, -4.89)	! tests 7-9
	c2 = conjg(c1)
	x3 = aimag(c2)

	cd1 = (4.49, -10.2)	! tests 10-12
	cd2 = conjg(c2)
	d3 = dimag(cd2)

	d4 = dprod(x1, x2)	! test 13

	i7 = max(i1, i2, i3, i4, i5, i6)	! tests 14-15
	i8 = min(i1, i2, i3, i4, i5, i6, i7)

	x4 = max(x1, x2, x3)			! tests 16-17
	x5 = min(x1, x2, x3, x4)

	d5 = max(d1, d2, d3, d4)		! tests 18-19
	d6 = min(d1, d2, d3, d4, d5)

	i9 = iand('aaaaaaaa'x, '84848484'x)	! tests 20-23
	i10 = ior('78787878'x, '87878787'x)
	i11 = not('f0f0f0f0'x)
	i12 = ieor('bbbbbbbb'x, '8af8afcd'x)

	i13 = mod(13, -3)			! tests 24-26
	x6 = mod(-13.3, 5.0)
	d7 = mod(15.3d+2, 1.53d+2)

	i14 = ichar('a')			! test 27

	i15 = dim(20,15)			! tests 28-33
	i16 = dim(1,1)
	x7 = dim(20.0,15.0)
	x8 = dim(1.0,1.0)
	d8 = dim(20.0d+1, 15.0d+1)
	d9 = dim(1.0d+1, 1.0d+1)

	i17 = ishft('a0000000'x, 1)		! tests 34-35
	i18 = ishft('80000001'x, -1)

	ch = char(32)			! test 36

	l1 = lge('abcd', 'abcd')		! tests 37-40
	l2 = lgt('abcd', 'abc')
	l3 = lle('abc', 'abc')
	l4 = llt('abc', 'abcd')

	result(1) = i1
	result(2) = i2

	result(3) = i3
	result(4) = i4
	result(5) = i5
	result(6) = i6

	result(7) = c1
	result(8) = x3
	result(9) = cd1

	result(10) = cd1
	result(11) = d3
	result(12) = cd2

	result(13) = d4

	result(14) = i7
	result(15) = i8

	result(16) = x4
	result(17) = x5

	result(18) = d5
	result(19) = d6

	result(20) = i9
	result(21) = i10
	result(22) = i11
	result(23) = i12

	result(24) = i13
	result(25) = x6
	result(26) = d7

	result(27) = i14

	result(28) = i15
	result(29) = i16
	result(30) = x7
	result(31) = x8
	result(32) = d8
	result(33) = d9

	result(34) = i17
	result(35) = i18

	result(36) = ichar(ch)

	result(37) = iand(l1, 1)
	result(38) = iand(l2, 1)
	result(39) = iand(l3, 1)
	result(40) = iand(l4, 1)

	call check(result, expect, N)
	data result/N*-1/
	data expect/
     + 2, 2,					! tests 1 - 2
     + 3, -4500, 3, -54999,			! tests 3 - 6
     + 2, 4, 4,					! tests 7-9
     + 4, -4, 2,				! tests 10-12
     + 11250,					! test 13
     + 3, -54999,				! tests 14-15
     + 4500, 2,					! tests 16-17
     + 54999, -4,				! tests 18-19
     + -2139062144, -1, 252645135, 826479734,	! tests 20-23
     + 1, -3, 0,				! tests 24-26
     + 97,					! test 27
     + 5, 0, 5, 0, 50, 0,			! tests 28-33
     + 1073741824, 1073741824,			! tests 34-35
     + 32,					! test 36
     + 1, 1, 1, 1 /				! tests 37-40
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
*----------------------------------------------------------------------*
      subroutine init_contr(contr)
*----------------------------------------------------------------------*
*     initialize all pointers and set all mxvtx, mxarc, mxfac variables
*     to zero (describing the current size of allocated space for
*     vertices, arcs, and factorization info)
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'def_contraction.h'

      type(contraction), intent(inout) ::
     &     contr

      contr%mxvtx = 0
      contr%mxarc = 0
      contr%mxxarc = 0
      contr%mxfac = 0
      contr%nvtx = 0
      contr%narc = 0
      contr%nxarc = 0
      contr%nfac = 0
      contr%dagger = .false.
      contr%unique_set = .false.
      nullify(contr%vertex)
      nullify(contr%arc)
      nullify(contr%xarc)
      nullify(contr%joined)
      nullify(contr%svertex)
      nullify(contr%inffac)
      nullify(contr%vtx)
      nullify(contr%topo)
      nullify(contr%xlines)
      contr%index_info = .false.
      nullify(contr%contr_string)
      nullify(contr%result_string)
      contr%nidx = 0
      contr%nxidx = 0
      contr%total_sign = 0
      contr%eqvl_fact = 1d0

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
C Copyright(c) 1998, Space Science and Engineering Center, UW-Madison
C Refer to "McIDAS Software Acquisition and Distribution Policies"
C in the file  mcidas/data/license.txt

C *** $Id: nvxrect.dlm,v 1.13 1999/03/31 21:36:59 rodr Exp $ ***
C       ************************************************************************
C       * RECTINI(IFUNC,IPARMS) -- Initializes this modules Nav routine data
C       *  IFUNC -- ???
C       *  IPARMS -- array of Navblock entries as follows
C       * Projection : T
C       *    WORD 1 -- 4 character projection type RECT MERC LAMB etc..
C 	*         2 -- Image Row Number	of the Center of the image
C 	*         3 -- Latitude(degrees)*10000 of Image Row (2)
C 	*         4 -- Image Column Number of the Center of the image
C   	*         5 -- Longitude *10000 of 4
C 	*         6 -- Lat(degrees per line)*10000
C 	*         7 -- Lon(degrees per element)*10000
C 	*         8 -- radius of the planet in meters
C 	*         9 -- eccentricity of the planet * 1000000
C       *        10 -- coordinate type >=0 panetodetic <0 panetocentric
C 	*        11 -- longitude convention >=0 west positive <0 west neg.
C  	*        12 -- LAT power of 10 -- if 0 leave at 4 else this power
C  	*        13 -- LON power of 10 -- if 0 leave at 4 else this power
C  	*        14 -- Degress_per_line power
C 	*        15 -- Degrees_per_element power
C 	*        16 -- Degrees_radian power
C 	*        17 -- Degrees_eccentricity power
C 	*	 ** -- next for added to IMGREMAP but not currently used anyway
C 	*        21 -- ecor -- element number of upper left hand corner
C 	*        22 -- lcor -- line number of upper left hand corner
C 	*        23 -- esize -- size in elements of coverage (width)
C 	*        24 -- lsize -- size in lines of coverage (height)
C	*
C       *************************************************************************

      INTEGER FUNCTION RECTINI(IFUNC,IPARMS)
      DIMENSION IPARMS(*)
      REAL*8 DRAD,DECC

      INTEGER IPOWLAT
      INTEGER IPOWLON
      INTEGER IPOWDLIN		! Degrees Per Line Power
      INTEGER IPOWDELE		! Degrees Per Element Power
      INTEGER IPOWRAD		! Degrees Radian power
      INTEGER IPOWECC		! Eccentricity Power

      CHARACTER*4 CLIT
      COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST ! ,ULLON

      IF (IFUNC.EQ.1) THEN
	IF (IPARMS(1).NE.LIT('RECT')) GO TO 900
	  ITYPE=1
	  XROW=IPARMS(2)
	  IPOWLAT=IPARMS(12)
	  IF (IPOWLAT .EQ. 0) IPOWLAT=4 		! default is 10000 (10^4)
	  ZSLAT=IPARMS(3)/10.**IPOWLAT  		! REAL Latitude
	  XCOL=IPARMS(4)
	  IPOWLON=IPARMS(13)
	  IF (IPOWLON .EQ. 0) IPOWLON=4
	  ZSLON=IPARMS(5)/10.**IPOWLON  		! REAL Longitude
	  IPOWDLIN=IPARMS(14)
	  IF (IPOWDLIN .EQ. 0) IPOWDLIN=4
	  ZDLAT=IPARMS(6)/10.**IPOWDLIN  		! REAL Degrees_per_line_latitude
	  IPOWDELE=IPARMS(15)
	  IF (IPOWDELE .EQ. 0) IPOWDELE=4
	  ZDLON=IPARMS(7)/10.**IPOWDELE  		! REAL Degrees_per_line_longitude
	  IPOWRAD=IPARMS(16)
	  IF (IPOWRAD .EQ. 0) IPOWRAD=3
	  DRAD=IPARMS(8)/10.D0**IPOWRAD             ! REAL Radius of the planet in meters
	  R=DRAD
	  IPOWECC=IPARMS(17)
	  IF (IPOWECC .EQ. 0) IPOWECC=6
	  DECC=IPARMS(9)/10.D0**IPOWECC             ! REAL Eccentricity
	  IWEST=IPARMS(11)                          ! West positive vs. West negative
	  IF(IWEST.GE.0) IWEST=1
          CALL LLOPT(DRAD,DECC,IWEST,IPARMS(10))    ! Initialze LLCART code

          IF (XCOL.EQ.1) THEN		! special case of XCOL not located at image center
	     ZSLON=ZSLON-180.0*IWEST    ! -- so assume it's the left edge(duh)
	  ENDIF

	ELSE IF (IFUNC.EQ.2) THEN
	   IF(INDEX(CLIT(IPARMS(1)),'XY').NE.0) ITYPE=1
	   IF(INDEX(CLIT(IPARMS(1)),'LL').NE.0) ITYPE=2
	ENDIF

	RECTINI=0
	RETURN
900 	CONTINUE
	RECTINI=-1
	RETURN
	END

C     ***************************************************************************
C     * RECTSAE -- Line/Element to Lat/Lon					*
C     ***************************************************************************
      INTEGER FUNCTION RECTSAE(XLIN,XELE,XDUM,XLAT,XLON,Z)
      COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST !,ULLON

      XLDIF=XROW-XLIN
      if(xcol.eq.1) then
        XEDIF=IWEST*(XELE-XCOL)
        XLON=ZSLON+180*IWEST-XEDIF*ZDLON
      else
        XEDIF=IWEST*(XCOL-XELE)
        XLON=ZSLON+XEDIF*ZDLON
      endif
      XLAT=ZSLAT+XLDIF*ZDLAT
      IF(XLAT.GT.90. .OR. XLAT.LT.-90.) GO TO 900
C --  All we want is to keep the XLON within 180 degrees of the center LON on both sides
      IF(XLON.GT.(ZSLON+180.0)) GO TO 900
      IF(XLON.LT.(ZSLON-180.0)) GO TO 900
      IF(XLON.LT.-180.0) THEN
        XLON=XLON+360.0
      ELSE IF (XLON.GT.180) THEN
        XLON=XLON-360.0
      ENDIF

      IF(ITYPE.EQ.1) THEN
         YLAT=XLAT
         YLON=XLON
         CALL LLCART(YLAT,YLON,XLAT,XLON,Z)
      ENDIF

      RECTSAE=0
      RETURN
900   CONTINUE
      RECTSAE=-1
      RETURN
      END

C     ***************************************************************************
C     * RECTSAE -- Lat/Lon to Line/Element 					*
C     ***************************************************************************
      INTEGER FUNCTION RECTEAS(ZLAT,ZLON,Z,XLIN,XELE,XDUM)
      CHARACTER*12 cfe
      COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST ! ,ULLON

      XLAT=ZLAT
      XLON=ZLON
      IF(ITYPE.EQ.1) THEN
         X=XLAT
         Y=XLON
         CALL CARTLL(X,Y,Z,XLAT,XLON)
      ENDIF
C --  Keep XLON within 180.0 degrees of ZSLON
      IF(XLON.GT.(ZSLON+180.0)) THEN
        XLON=XLON-360.0
      ELSEIF(XLON.LT.ZSLON-180.0) THEN
        XLON=XLON+360.0
      ENDIF

      XLIN=XROW-(XLAT-ZSLAT)/ZDLAT
      IF (XCOL.EQ.1) THEN
C --    Need to adjust for the fact that XCOL is not really at the image center
        XELE=XCOL-(XLON-ZSLON-180.0*IWEST)/(ZDLON*IWEST)
      ELSE
        XELE=XCOL-(XLON-ZSLON)/(ZDLON*IWEST)
      ENDIF

      RECTEAS=0
      RETURN
      END
      INTEGER FUNCTION RECTOPT(IFUNC,XIN,XOUT)
      COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST
      REAL*4 XIN(*),XOUT(*)
      CHARACTER*4 CLIT,CFUNC
C
C IFUNC= 'SPOS'    SUBSATELLITE LAT/LON
C
C        XIN - NOT USED
C        XOUT - 1. STANDARD LATITUDE
C             - 2. NORMAL LONGITUDE
C
C
C IFUNC= 'ORAD'  OBLATE RADIUS
C
C        XIN - LATITUDE
C        XOUT - RADIUS IN KM
C
      CFUNC=CLIT(IFUNC)
      RECTOPT=0
      IF(CFUNC.EQ.'SPOS') THEN
ccc         XOUT(1)=ZSLAT
ccc         XOUT(2)=ZSLON
C
C--- NO SAT. POSITION FOR THIS PROJECTION
C
         RECTOPT = -1
C
      ELSE IF(CFUNC.EQ.'ORAD') THEN
         CALL LLOBL(XIN,XOUT)
      ELSE
         RECTOPT=1
      ENDIF
      RETURN
      END


text:
      SUBROUTINE COMPF(NUMPAR,NUMVAR,VALVAR,SUMERR,DIFFS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES/NOLIST'
      DIMENSION DIFFS(*), VALVAR(*)
      COMMON /GEOVAR/ NVAR, LOC(2,MAXPAR), XPARAM(MAXPAR)
***************************************************************************
*
*  THESE COMMON BLOCKS HOLD THE ADDRESSES OF THE PARAMETERS TO
*        BE OPTIMISED.
*
      COMMON /LOCVAR/ LOCVAR(2,100)
      DIMENSION GRAD(100), COORD(3,NUMATM)
      TIME1=SECOND()
      SUMERR=0.D0
          CALL COMPFG(VALVAR,.TRUE.,ESCF,.TRUE.,DIFFS,.TRUE.)
                  SUMERR=DOT(DIFFS,DIFFS,NUMVAR)
C#      WRITE(6,'('' SSQ=  '',F9.3)')SUMERR
C#      WRITE(6,'(6F12.4)')(DIFFS(I),I=1,NVAR)
      RETURN
      END

ext: for

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive