text:
	SUBROUTINE GH_UHST ( flpart, nstrm, strmid, advno, iret )
C************************************************************************
C* GH_UHST								*
C*									*
C* This subroutine updates the history file using the storm identifiers *
C* and advisory numbers processed in the current run.                   *
C*									*
C* GH_UHST ( FLPART, NSTRM, STRMID, ADVNO, IRET )                       *
C*									*
C* Input parameters:                                                    *
C*      FLPART		CHAR*		History file name (w/out year)  *
C*	NSTRM		INTEGER		Number of new storm advisories  *
C*	STRMID (*)	CHAR*		Storm ids for new advisories    *
C* 	ADVNO (*)	CHAR*		New advisory numbers            *
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return             *
C*					 -1 = error opening history file*
C*									*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 5/01  						*
C* D. Kidwell/NCEP	10/01	Added input parm flpart                 *
C* D. Kidwell/NCEP	 4/02   Used format for lunhst integer i/o      *
C* S. Gilbert/NCEP	 1/06   Allow for four character advisory num   *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	strmid (*), advno (*), flpart
C*
	CHARACTER*(MXFLSZ)      filnam
C*
	CHARACTER	strmhs (500)*13, strnew (500)*13
	LOGICAL		done
C-----------------------------------------------------------------------
	iret  = 0
C
C*	Open the history file.
C
	CALL ST_LSTR ( flpart, lenf, ier )
	filnam = flpart ( :lenf) // '.' // strmid ( 1 ) ( 5:8 )
	CALL FL_SWOP ( filnam, lunhst, ier )
	IF ( ier .ne. 0 ) THEN
	    iret = -1
	    RETURN
	END IF
C
C*	Read the history file.
C
	READ (lunhst,200,IOSTAT=iostat) numhst
200	FORMAT ( I4 )
	IF ( iostat .ne. 0 ) numhst = 0
	IF ( numhst .gt. 0 ) THEN
	    READ (lunhst,100) ( strmhs (ii), ii = 1, numhst )
100         FORMAT (A)
	  ELSE
C
C*	    There are no entries in the history file.  The current
C*	    file(s) will be added.
C
	    strmhs ( 1 ) = ' '
	END IF
C
C*	Loop to update current storm advisories and add new storms.
C
	CALL FL_REWD ( lunhst, ier )
	nnew = 0
	DO nn = 1, nstrm
	    done = .false.
	    ii   = 1
	    DO WHILE ( .not. done )
		IF ( strmid ( nn ) .eq. strmhs ( ii ) ( :8 ) ) THEN
C
C*		    Update the advisory number for an existing storm.
C
		    strmhs ( ii ) ( 10:13 ) = advno ( nn ) ( :4 )
		    done = .true.
		  ELSE
		    ii = ii + 1
		    IF ( ii .gt. numhst ) THEN
C
C*			Add a new storm.
C
			nnew = nnew + 1
			strnew ( nnew ) = strmid (nn) ( :8 ) // '_' //
     +					  advno (nn) ( :4 )
			done = .true.
		    END IF
		END IF
	    END DO
	END DO
C
C*	Write out the updated history file, with new storms first.
C
	numtot = numhst + nnew
	WRITE (lunhst,200) numtot
	DO ii = 1, nnew
	    WRITE (lunhst,100) strnew ( ii )
	END DO
	DO ii = 1, numhst
	    WRITE (lunhst,100) strmhs ( ii )
	END DO
	CALL FL_CLOS ( lunhst, 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:
!************************************************************

SUBROUTINE zludcmp(A,det,n)
IMPLICIT NONE

INTEGER*4, INTENT(IN) :: n
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(IN) :: A
COMPLEX*16, INTENT(OUT) :: det

! ******
! For input variables with zero-based indexing:
! Calculation of det(A) using LU decomposition. Note that in general,
! numerical overflow may result unless LOG(det(A)) is formed instead.

COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: A_aux
INTEGER, DIMENSION(:), ALLOCATABLE :: ipiv
INTEGER :: m,i,info

ALLOCATE(A_aux(1:n,1:n))
ALLOCATE(ipiv(1:n))

A_aux(1:n,1:n) = A(0:n-1,0:n-1)

m = n
CALL ZGETRF(m,m,A_aux,m,ipiv,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRF failed in zludcmp, code: ',info
	STOP

END IF

det = 1.d0

DO i = 1,n

	IF (ipiv(i) /= i) det = -det * A_aux(i,i)
	IF (ipiv(i) == i) det = det * A_aux(i,i)

END DO

DEALLOCATE(A_aux,ipiv)

END SUBROUTINE zludcmp

!************************************************************

SUBROUTINE zludcmp_ratio(A,detratio,ldetref,n)
IMPLICIT NONE

INTEGER*4, INTENT(IN) :: n
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(IN) :: A
DOUBLE PRECISION, INTENT(IN) :: ldetref
COMPLEX*16 :: ldet
COMPLEX*16, INTENT(OUT) :: detratio

! ******
! For input variables with zero-based indexing:
! Calculation of det(A)/detref = det(A)/exp(ldetref) using LU decomposition,
! where detref is positive

COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: A_aux
INTEGER, DIMENSION(:), ALLOCATABLE :: ipiv
INTEGER :: m,i,info

ALLOCATE(A_aux(1:n,1:n))
ALLOCATE(ipiv(1:n))

A_aux(1:n,1:n) = A(0:n-1,0:n-1)

m = n
CALL ZGETRF(m,m,A_aux,m,ipiv,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRF failed in zludcmp, code: ',info
	STOP

END IF

ldet = 0.d0

DO i = 1,n

	IF (ipiv(i) /= i) ldet = ldet + LOG(-A_aux(i,i))
	IF (ipiv(i) == i) ldet = ldet + LOG(A_aux(i,i))

END DO

detratio = exp(ldet-ldetref)

!det = 1.d0
!
!DO i = 1,n
!
!	IF (ipiv(i) /= i) det = -det * A_aux(i,i)
!	IF (ipiv(i) == i) det = det * A_aux(i,i)
!
!END DO

DEALLOCATE(A_aux,ipiv)

END SUBROUTINE zludcmp_ratio

!************************************************************

SUBROUTINE zludcmp_log(A,ldet,n)
IMPLICIT NONE

INTEGER*4, INTENT(IN) :: n
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(IN) :: A
COMPLEX*16, INTENT(OUT) :: ldet

! ******
! For input variables with zero-based indexing:
! Calculation of log(det(A)) using LU decomposition.

COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: A_aux
INTEGER, DIMENSION(:), ALLOCATABLE :: ipiv
INTEGER :: m,i,info

ALLOCATE(A_aux(1:n,1:n))
ALLOCATE(ipiv(1:n))

A_aux(1:n,1:n) = A(0:n-1,0:n-1)

m = n
CALL ZGETRF(m,m,A_aux,m,ipiv,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRF failed in zludcmp, code: ',info
	STOP

END IF

ldet = 0.d0

DO i = 1,n

	IF (ipiv(i) /= i) ldet = ldet + LOG(-A_aux(i,i))
	IF (ipiv(i) == i) ldet = ldet + LOG(A_aux(i,i))

END DO

!det = 1.d0
!
!DO i = 1,n
!
!	IF (ipiv(i) /= i) det = -det * A_aux(i,i)
!	IF (ipiv(i) == i) det = det * A_aux(i,i)
!
!END DO

DEALLOCATE(A_aux,ipiv)

END SUBROUTINE zludcmp_log

!***********************************************************

SUBROUTINE zdoinv(A,x,det,n)
IMPLICIT NONE

INTEGER*4, INTENT(IN) :: n
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(IN) :: A
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(OUT) :: x
COMPLEX*16, INTENT(OUT) :: det

! ******
! For input variables with zero-based indexing:
! Calculation of A^-1 from A*x = I including det(A).

COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: A_aux,B_aux
INTEGER, DIMENSION(:), ALLOCATABLE :: ipiv
INTEGER :: m,nRhs,i,info

ALLOCATE(A_aux(1:n,1:n))
ALLOCATE(ipiv(1:n))

A_aux(1:n,1:n) = A(0:n-1,0:n-1)

m = n
CALL ZGETRF(m,m,A_aux,m,ipiv,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRF failed in zdoinv, code: ',info
	STOP

END IF

det = 1.d0

DO i = 1,n

	IF (ipiv(i) /= i) det = -det * A_aux(i,i)
	IF (ipiv(i) == i) det = det * A_aux(i,i)

END DO

nRhs = n
ALLOCATE(B_aux(1:n,1:nRhs))
B_aux = 0.d0

DO i = 1,nRhs

	B_aux(i,i) = 1.d0

END DO

CALL ZGETRS('N',m,nRhs,A_aux,m,ipiv,B_aux,m,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRS failed in zdoinv, code: ',info
	STOP

END IF

x(0:n-1,0:n-1) = B_aux(1:n,1:n)

DEALLOCATE(A_aux,B_aux,ipiv)

END SUBROUTINE zdoinv

!***********************************************************

SUBROUTINE zdoinv_log(A,x,ldet,n)
IMPLICIT NONE

INTEGER*4, INTENT(IN) :: n
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(IN) :: A
COMPLEX*16, DIMENSION(0:n-1,0:n-1), INTENT(OUT) :: x
COMPLEX*16, INTENT(OUT) :: ldet

! ******
! For input variables with zero-based indexing:
! Calculation of A^-1 from A*x = I including det(A).

COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: A_aux,B_aux
INTEGER, DIMENSION(:), ALLOCATABLE :: ipiv
INTEGER :: m,nRhs,i,info

ALLOCATE(A_aux(1:n,1:n))
ALLOCATE(ipiv(1:n))

A_aux(1:n,1:n) = A(0:n-1,0:n-1)

m = n
CALL ZGETRF(m,m,A_aux,m,ipiv,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRF failed in zdoinv, code: ',info
	STOP

END IF

ldet = 0.d0

DO i = 1,n

	IF (ipiv(i) /= i) ldet = ldet + LOG(-A_aux(i,i))
	IF (ipiv(i) == i) ldet = ldet + LOG(A_aux(i,i))

END DO

! det = 1.d0
!
! DO i = 1,n
!
!	IF (ipiv(i) /= i) det = -det * A_aux(i,i)
!	IF (ipiv(i) == i) det = det * A_aux(i,i)
!
! END DO

nRhs = n
ALLOCATE(B_aux(1:n,1:nRhs))
B_aux = 0.d0

DO i = 1,nRhs

	B_aux(i,i) = 1.d0

END DO

CALL ZGETRS('N',m,nRhs,A_aux,m,ipiv,B_aux,m,info)

IF (info /= 0) THEN

	WRITE (*,*) 'ZGETRS failed in zdoinv, code: ',info
	STOP

END IF

x(0:n-1,0:n-1) = B_aux(1:n,1:n)

DEALLOCATE(A_aux,B_aux,ipiv)

END SUBROUTINE zdoinv_log

!*************************************************************

max_stars_repo_name: cpmoca/LectureNotesPhysics

max_issues_repo_name: cpmoca/LectureNotesPhysics

max_forks_repo_name: cpmoca/LectureNotesPhysics


text:
      program vacf

c*********************************************************************
c
c     dl_poly routine to calculate velocity autocorrelation function
c     for selected atoms
c
c     copyright daresbury laboratory 1995
c     author  w.smith jan 1995
c
c*********************************************************************

      implicit real*8(a-h,o-z)

      parameter (mxatms=1080,npts=128,nstr=250)
      character*8 atname
      character*8 atmnam
      character*40 fname
      character*80 title
      dimension vxx(mxatms),vyy(mxatms),vzz(mxatms)
      dimension vacf(npts,nstr),vcr0(npts,nstr,3),vcf(npts)
      dimension idv(npts),nvm(npts),key(npts)
      complex*16 aaa(npts),www(npts)

      write(*,'(a)')'Velocity Autocorrelation Program'

c     read the control variables

c     number of atoms
      read(*,*)natms
c     working length of vacf array
      read(*,*)nvacf
c     interval between origins
      read(*,*)iovacf
c     time interval between successive configurations
      read(*,*)tstep
c     name of species for autocorrelation
      read(*,*)atname
c     name of selected HISTORY file
      read(*,*)fname
c     number of configurations in HISTORY file
      read(*,*)nconf

c     check on specified control variables

      if(nvacf.gt.npts)then

         write(*,'(a,2i5)')'error - vacf array too large',nvacf,npts
         stop

      endif
      if(mod(nvacf,iovacf).ne.0)then

         nvacf=iovacf*(nvacf/iovacf)
         write(*,'(a,i5)')'warning - vacf array reset to ',nvacf

      endif

      novacf=nvacf/iovacf

c     initialise velocity autocorrelation arrays

      lor=0
      mor=0
      nsvacf=0
      do j=1,npts
         nvm(j)=0
         do i=1,nstr
            vacf(j,i)=0.d0
            vcr0(j,i,1)=0.d0
            vcr0(j,i,2)=0.d0
            vcr0(j,i,3)=0.d0
         enddo
      enddo

c     open HISTORY file and read headers

      open (10,file=fname,form='formatted')

      read(10,'(a80)')title
      write(*,'(a80)')title
      read(10,*)levcfg,imcon

      if(levcfg.eq.0)then

         write(*,'(a)')'error - HISTORY file contains no velocities'
         stop

      endif

      do iconf=1,nconf

         read(10,*,end=100)

         if(imcon.gt.0)then

            read(10,*,end=100)
            read(10,*,end=100)
            read(10,*,end=100)

         endif

         nsampl=0
         do i=1,natms

            read(10,*,end=100) atmnam
            read(10,*)
            read(10,*) xxx,yyy,zzz
            if(levcfg.gt.1)read(10,*)
            if(atmnam.eq.atname)then

               nsampl=nsampl+1

               if(nsampl.le.nstr)then

                  vxx(nsampl)=xxx
                  vyy(nsampl)=yyy
                  vzz(nsampl)=zzz

               endif

            endif

         enddo

         if(nsampl.gt.nstr)then

            write(*,'(a)')'error - too many atoms of required type'
            stop

         endif

c     calculate velocity autocorrelation function

         if(mod(nsvacf,iovacf).eq.0)then

            lor=min(lor+1,novacf)
            mor=mod(mor,novacf)+1
            idv(mor)=1

            do i=1,nsampl

               vcr0(mor,i,1)=vxx(i)
               vcr0(mor,i,2)=vyy(i)
               vcr0(mor,i,3)=vzz(i)

            enddo

         endif

         nsvacf=nsvacf+1

         do l=1,lor

            m=idv(l)
            idv(l)=m+1
            nvm(m)=nvm(m)+1

            do i=1,nsampl

               vacf(m,i)=vacf(m,i)+vxx(i)*vcr0(l,i,1)+
     x              vyy(i)*vcr0(l,i,2)+vzz(i)*vcr0(l,i,3)

            enddo

         enddo

      enddo

  100 continue

      last=min(nsvacf,nvacf)

c     normalise velocity autocorrelation functions

      do i=1,nsampl

         rnorm=dble(nvm(1))/vacf(1,i)

         do j=1,last

            vacf(j,i)=rnorm*vacf(j,i)/dble(nvm(j))

         enddo

      enddo

c     calculate normalised vacf

      vint=0.0d0
      do j=1,last

         vsum=0.0d0
         do i=1,nsampl

            vsum=vsum+vacf(j,i)

         enddo

         vcf(j)=vsum/dble(nsampl)
         vint=vint+vcf(j)
         time=tstep*dble(j-1)
         write(*,'(1p,3e15.7)')time,vcf(j),vint

      enddo
      write(*,'(a)')'&'

c     initialise FFT arrays

      do i=1,npts

         if(i.le.nvacf)then

            aaa(i)=cmplx(vcf(i),0.d0)

         else

            aaa(i)=(0.d0,0.d0)

         endif

      enddo
      aaa(1)=0.5d0*aaa(1)

c     perform power spectrum

      call fft(1,-1,npts,key,aaa,www,aaa)
      call fft(0,-1,npts,key,aaa,www,aaa)

      do i=1,npts

         freq=dble(i-1)/(dble(npts)*tstep)
         write(*,'(1p,2e15.7)')freq,real(aaa(i))

      enddo
      write(*,'(a)')'&'

      stop
      end
      subroutine fft(ind,isw,ndiv,key,aaa,wfft,bbb)
c***********************************************************************
c
c     fast fourier transform routine
c
c     copyright daresbury laboratory 1994
c     author w smith
c
c***********************************************************************

      implicit real*8(a-h,o-z)

      logical check
      complex*16 aaa(ndiv),bbb(ndiv),wfft(ndiv),ttt
      dimension key(ndiv)
      data tpi/6.2831853072d0/
   10 format(1h0,'error - number of points not a power of two')

c
c     check that array is of suitable length
      nt=1
      check=.true.
      do i=1,20
         nt=2*nt
         if(nt.eq.ndiv)then
            check=.false.
            nu=i
         endif
      enddo
      if(check)then
         write(*,10)
         stop
      endif

      if(ind.gt.0)then
c
c     set reverse bit address array

         do kkk=1,ndiv
            iii=0
            jjj=kkk-1
            do j=1,nu
               jj2=jjj/2
               iii=2*(iii-jj2)+jjj
               jjj=jj2
            enddo
            key(kkk)=iii+1
         enddo
c
c     initialise complex exponential factors

         tpn=tpi/dble(ndiv)
         arg=0.d0
         np1=ndiv+1
         np2=ndiv/2
         wfft(1)=(1.d0,0.d0)
         do i=1,np2
            arg=tpn*dble(i)
            wfft(i+1)=cmplx(cos(arg),sin(arg))
            wfft(np1-i)=conjg(wfft(i+1))
         enddo

         return
      endif

c
c     take conjugate of exponentials if required

      if(isw.lt.0)then

         do i=1,ndiv
            wfft(i)=conjg(wfft(i))
         enddo

      endif

c
c     take copy input array

      do i=1,ndiv
         bbb(i)=aaa(i)
      enddo
c
c     perform fourier transform

      kkk=0
      nu1=nu-1
      np2=ndiv/2
      do l=1,nu

  100    do i=1,np2
            iii=key(kkk/2**nu1+1)
            kk1=kkk+1
            k12=kk1+np2
            ttt=bbb(k12)*wfft(iii)
            bbb(k12)=bbb(kk1)-ttt
            bbb(kk1)=bbb(kk1)+ttt
            kkk=kkk+1
         enddo
         kkk=kkk+np2
         if(kkk.lt.ndiv)go to 100
         kkk=0
         nu1=nu1-1
         np2=np2/2

      enddo
c
c     unscramble the fft using bit address array

      do kkk=1,ndiv
         iii=key(kkk)
         if(iii.gt.kkk)then
            ttt=bbb(kkk)
            bbb(kkk)=bbb(iii)
            bbb(iii)=ttt
         endif
      enddo
c
c     restore exponentials to unconjugated values if necessary

      if(isw.lt.0)then

         do i=1,ndiv
            wfft(i)=conjg(wfft(i))
         enddo

      endif

      return
      end


text:
! NAME
!
!     pmca [ pmca ] - Probability and Mass Current Analysis Module
!
! SYNOPSIS
!
!     USE :: PMCA
!
! DESCRIPTION
!
!     PMCA is a custom Fortran module written to compute the probability
!     and mass currents from a single-particle and/or single-component
!     Bose-Einstein condensate wave function.
!
! OPTIONS
!
! SEE ALSO
!
! BUGS
!
! HISTORY
!
!     Created to study the quantum Compton generator data, where, in
!     particular, the more simple angular momentum expectation values
!     used initially led to large disagreements with the theory due to
!     the precense of many vortex-like structures and other complicating
!     factors in the raw simulation data.
!
! NOTES
!
! AUTHOR
!
!     Marty Kandes, Ph.D.
!     Computational & Data Science Research Specialist
!     High-Performance Computing User Services Group
!     San Diego Supercomputer Center
!     University of California, San Diego
!
! COPYRIGHT
!
!     Copyright (c) 2014, 2015, 2016, 2017, 2018, 2019, 2020 Martin Charles Kandes
!
! LAST UPDATED
!
!     Saturday, January 4th, 2019
!
! ----------------------------------------------------------------------

      MODULE pmca

! --- MODULE DECLARATIONS ----------------------------------------------

      USE, INTRINSIC :: ISO_FORTRAN_ENV
      USE :: MPI
      USE :: GRID

! --- MODULE DEFINITIONS -----------------------------------------------
!
!     ISO_FORTRAN_ENV is the intrinsic Fortran module that provides
!        information about the run-time environment.
!
!     MPI is the standard Message Passing Interface module.
!
!     GRID is a custom Fortran module written to build the regular
!        Cartesian grid of points that defines the spatial computational
!        domain on which external potentials, wave functions and other
!        physical quantities may be computed.
!
! ----------------------------------------------------------------------

      IMPLICIT NONE

!      INCLUDE 'mpif.h' ! if MPI module is not supported, use header file.

      PRIVATE

! --- VARIABLE DECLARATIONS --------------------------------------------

      REAL, PRIVATE :: pmcaI12 = 0.0
      REAL, PRIVATE :: pmcaI34 = 0.0
      REAL, PRIVATE :: pmcaI41 = 0.0
      REAL, PRIVATE :: pmcaI23 = 0.0

      REAL, PRIVATE :: pmcaI1 = 0.0
      REAL, PRIVATE :: pmcaI2 = 0.0
      REAL, PRIVATE :: pmcaI3 = 0.0
      REAL, PRIVATE :: pmcaI4 = 0.0

      REAL, PRIVATE :: pmcaV12 = 0.0
      REAL, PRIVATE :: pmcaV34 = 0.0
      REAL, PRIVATE :: pmcaV41 = 0.0
      REAL, PRIVATE :: pmcaV23 = 0.0

      REAL, PRIVATE :: pmcaV1 = 0.0
      REAL, PRIVATE :: pmcaV2 = 0.0
      REAL, PRIVATE :: pmcaV3 = 0.0
      REAL, PRIVATE :: pmcaV4 = 0.0

! --- VARIABLE DEFINITIONS ---------------------------------------------
!
! --- SUBROUTINE DECLARATIONS ------------------------------------------

      PUBLIC :: pmca_compute_currents
      PUBLIC :: pmca_write_currents

      PUBLIC :: pmca_compute_velocities
      PUBLIC :: pmca_write_velocities

      PUBLIC :: pmca_density
      PUBLIC :: pmca_phase
      PUBLIC :: pmca_velocity
      PUBLIC :: pmca_current_density
      PUBLIC :: pmca_quantum_potential
      PUBLIC :: pmca_vorticity

      PRIVATE :: pmca_grad_f_real_3d_rect_cd2
      PRIVATE :: pmca_grad_f_cmplx_3d_rect_cd2
      PRIVATE :: pmca_div_f_cmplx_3d_rect_cd2
      PRIVATE :: pmca_curl_f_cmplx_3d_rect_cd2
      PRIVATE :: pmca_laplacian_f_cmplx_3d_rect_cd2

! --- SUBROUTINE DEFINITIONS -------------------------------------------
!
! --- FUNCTION DECLARATIONS --------------------------------------------

      PRIVATE :: pmca_current_3d_rect
      PRIVATE :: pmca_velocity_3d_rect

! --- FUNCTION DEFINITIONS ---------------------------------------------
!
! ----------------------------------------------------------------------

      CONTAINS

! ----------------------------------------------------------------------

      SUBROUTINE pmca_compute_currents(mpiRank, mpiMaster, mpiReal, &
         & mpiError, pmcaQuadRule, nXa, nXb, nXbc, nYa, nYb, nYbc, nZa,&
         & nZb, nZbc, nZ, dX, dY, dZ, Xa, Ya, Za, Zc, J3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: mpiRank, mpiMaster, mpiReal
      INTEGER, INTENT(INOUT) :: mpiError
      INTEGER, INTENT(IN) :: pmcaQuadRule
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc, nZ

      REAL, INTENT (IN) :: dX, dY, dZ

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: Xa
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Ya
      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(IN) :: Za
      REAL, DIMENSION(  1 - nZbc : nZ  + nZbc), INTENT(IN) :: Zc

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc) :: J3

      REAL :: pmcaI = 0.0 ! temporary, local slab current variable

      IF (pmcaQuadRule == 1) THEN

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc,&
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, dX, dY, Ya(nYb), &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI12, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, Xa(nXa), 0.0, dY, dY, &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI23, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, 0.0, 0.0, Ya(nYa), 0.0, &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI34, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, Xa(nXb), 0.0, 0.0, &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI41, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, Xa(nXb), dY, &
            & Ya(nYb), 0.0, 0.0, Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI1, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, Xa(nXa), 0.0, dY, &
            & Ya(nYb), 0.0, 0.0, Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI2, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, Xa(nXa), 0.0, Ya(nYa), &
            & 0.0, 0.0, 0.0, Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI3, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaI = pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, Xa(nXb), Ya(nYa), &
            & 0.0, 0.0, 0.0, Xa, Ya, Za, Zc, J3)
         CALL MPI_REDUCE(pmcaI, pmcaI4, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

      ELSE

         WRITE(UNIT=ERROR_UNIT, FMT = *) 'gpse: pmca_compute_currents: &
            & ERROR - pmcaQuadRule is not supported.'
         STOP

      END IF

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_write_currents(mpiRank, mpiMaster, tN)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: mpiRank, mpiMaster

      REAL, INTENT (IN) :: tN

!     Write computed currents to file from mpiMaster
      IF ( mpiRank == mpiMaster ) THEN

         OPEN(UNIT=990, FILE='pmca-currents.output', ACCESS='APPEND', &
            & ACTION='WRITE', FORM='FORMATTED', STATUS='UNKNOWN')
         WRITE(UNIT=990, FMT='(9(F23.15))') tN, pmcaI1, &
            & pmcaI12, pmcaI2, pmcaI23, pmcaI3, pmcaI34, pmcaI4, &
            & pmcaI41
         CLOSE(UNIT=990, STATUS='KEEP')

      END IF

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_compute_velocities(mpiRank, mpiMaster, mpiReal, &
         & mpiError, pmcaQuadRule, nXa, nXb, nXbc, nYa, nYb, nYbc, nZa,&
         & nZb, nZbc, nZ, dX, dY, dZ, Xa, Ya, Za, Zc, V3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: mpiRank, mpiMaster, mpiReal
      INTEGER, INTENT(INOUT) :: mpiError
      INTEGER, INTENT(IN) :: pmcaQuadRule
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc, nZ

      REAL, INTENT (IN) :: dX, dY, dZ

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: Xa
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Ya
      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(IN) :: Za
      REAL, DIMENSION(  1 - nZbc : nZ  + nZbc), INTENT(IN) :: Zc

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc) :: V3

      REAL :: pmcaV = 0.0 ! temporary, local slab velocity variable

      IF (pmcaQuadRule == 1) THEN

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc,&
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, dX, dY, Ya(nYb), &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV12, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, Xa(nXa), 0.0, dY, dY, &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV23, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, 0.0, 0.0, Ya(nYa), 0.0, &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV34, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, Xa(nXb), 0.0, 0.0, &
            & Zc(1), Zc(nZ), Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV41, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, Xa(nXb), dY, &
            & Ya(nYb), 0.0, 0.0, Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV1, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, Xa(nXa), 0.0, dY, &
            & Ya(nYb), 0.0, 0.0, Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV2, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, Xa(nXa), 0.0, Ya(nYa), &
            & 0.0, 0.0, 0.0, Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV3, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

         pmcaV = pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, nYbc, &
            & nZa, nZb, nZbc, nZ, dX, dY, dZ, dX, Xa(nXb), Ya(nYa), &
            & 0.0, 0.0, 0.0, Xa, Ya, Za, Zc, V3)
         CALL MPI_REDUCE(pmcaV, pmcaV4, 1, mpiReal, MPI_SUM, &
            & mpiMaster, MPI_COMM_WORLD, mpiError)

      ELSE

         WRITE(UNIT=ERROR_UNIT, FMT = *) 'gpse: pmca_compute_velocities: &
            & ERROR - pmcaQuadRule is not supported.'
         STOP

      END IF

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_write_velocities(mpiRank, mpiMaster, tN)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: mpiRank, mpiMaster

      REAL, INTENT (IN) :: tN

!     Write computed currents to file from mpiMaster
      IF ( mpiRank == mpiMaster ) THEN

         OPEN(UNIT=991, FILE='pmca-velocities.output', ACCESS='APPEND',&
            & ACTION='WRITE', FORM='FORMATTED', STATUS='UNKNOWN')
         WRITE(UNIT=991, FMT='(9(F23.15))') tN, pmcaV1, pmcaV12, &
            & pmcaV2, pmcaV23, pmcaV3, pmcaV34, pmcaV4, pmcaV41
         CLOSE(UNIT=991, STATUS='KEEP')

      END IF

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_density(nXa, nXb, nXbc, nYa, nYb, nYbc, nZa, nZb,&
         & nZbc, Psi3, Rho3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa
      INTEGER, INTENT(IN) :: nXb
      INTEGER, INTENT(IN) :: nXbc
      INTEGER, INTENT(IN) :: nYa
      INTEGER, INTENT(IN) :: nYb
      INTEGER, INTENT(IN) :: nYbc
      INTEGER, INTENT(IN) :: nZa
      INTEGER, INTENT(IN) :: nZb
      INTEGER, INTENT(IN) :: nZbc

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: Rho3

      INTEGER :: j, k, l

      Rho3 = 0.0

!$OMP PARALLEL DO IF (nZa /= nZb) DEFAULT(SHARED) SCHEDULE(STATIC)
      DO l = nZa, nZb
!$OMP    PARALLEL DO IF (nZa == nZb) DEFAULT(SHARED) SCHEDULE(STATIC)
         DO k = nYa, nYb
            DO j = nXa, nXb
               Rho3(j,k,l) = ABS(Psi3(j,k,l))**2
            END DO
         END DO
!$OMP    END PARALLEL DO
      END DO
!$OMP END PARALLEL DO

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_phase(nXa, nXb, nXbc, nYa, nYb, nYbc, nZa, nZb,&
         & nZbc, Psi3, Phi3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa
      INTEGER, INTENT(IN) :: nXb
      INTEGER, INTENT(IN) :: nXbc
      INTEGER, INTENT(IN) :: nYa
      INTEGER, INTENT(IN) :: nYb
      INTEGER, INTENT(IN) :: nYbc
      INTEGER, INTENT(IN) :: nZa
      INTEGER, INTENT(IN) :: nZb
      INTEGER, INTENT(IN) :: nZbc

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: Phi3

      INTEGER :: j, k, l

      Phi3 = 0.0

!$OMP PARALLEL DO IF (nZa /= nZb) DEFAULT(SHARED) SCHEDULE(STATIC)
      DO l = nZa, nZb
!$OMP    PARALLEL DO IF (nZa == nZb) DEFAULT(SHARED) SCHEDULE(STATIC)
         DO k = nYa, nYb
            DO j = nXa, nXb
               Phi3(j,k,l) = ATAN2(AIMAG(Psi3(j,k,l)), REAL(Psi3(j,k,l)))
            END DO
         END DO
!$OMP    END PARALLEL DO
      END DO
!$OMP END PARALLEL DO

      RETURN
      END

! ----------------------------------------------------------------------

      SUBROUTINE pmca_velocity(nXa, nXb, nXbc, nYa, nYb, nYbc, nZa, &
         & nZb, nZbc, dX, dY, dZ, Rho3, J3, V3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa
      INTEGER, INTENT(IN) :: nXb
      INTEGER, INTENT(IN) :: nXbc
      INTEGER, INTENT(IN) :: nYa
      INTEGER, INTENT(IN) :: nYb
      INTEGER, INTENT(IN) :: nYbc
      INTEGER, INTENT(IN) :: nZa
      INTEGER, INTENT(IN) :: nZb
      INTEGER, INTENT(IN) :: nZbc

      REAL, INTENT(IN) :: dX
      REAL, INTENT(IN) :: dY
      REAL, INTENT(IN) :: dZ

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Rho3

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: J3

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: V3

      INTEGER :: i, j, k, l

      V3 = 0.0

      !$OMP PARALLEL DO IF (nZa /= nZb) DEFAULT(SHARED) SCHEDULE(STATIC)
      DO l = nZa, nZb
!$OMP    PARALLEL DO IF (nZa == nZb) DEFAULT(SHARED) SCHEDULE(STATIC)
         DO k = nYa, nYb
            DO j = nXa, nXb
               DO i = 1, 3
                  V3(i,j,k,l) = J3(i,j,k,l) / Rho3(j,k,l)
               END DO
            END DO
         END DO
!$OMP    END PARALLEL DO
      END DO
!$OMP END PARALLEL DO

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_current_density(nXa, nXb, nXbc, nYa, nYb, nYbc, &
         & nZa, nZb, nZbc, dX, dY, dZ, Psi3, GradPsi3, J3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      REAL, INTENT(IN) :: dX, dY, dZ

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      COMPLEX, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: GradPsi3

      COMPLEX, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: J3

      INTEGER :: j, k, l

      J3 = CMPLX(0.0,0.0)

      CALL pmca_grad_f_cmplx_3d_rect_cd2(nXa, nXb, nXbc, nYa, nYb, &
         & nYbc, nZa, nZb, nZbc, dX, dY, dZ, Psi3, GradPsi3)
      J3(1,:,:,:) = CONJG(Psi3) * GradPsi3(1,:,:,:)
      J3(2,:,:,:) = CONJG(Psi3) * GradPsi3(2,:,:,:)
      J3(3,:,:,:) = CONJG(Psi3) * GradPsi3(3,:,:,:)
      CALL pmca_grad_f_cmplx_3d_rect_cd2(nXa, nXb, nXbc, nYa, nYb, &
         & nYbc, nZa, nZb, nZbc, dX, dY, dZ, CONJG(Psi3), GradPsi3)
      J3(1,:,:,:) = J3(1,:,:,:) - Psi3 * GradPsi3(1,:,:,:)
      J3(2,:,:,:) = J3(2,:,:,:) - Psi3 * GradPsi3(2,:,:,:)
      J3(3,:,:,:) = J3(3,:,:,:) - Psi3 * GradPsi3(3,:,:,:)
      J3 = CMPLX(0.0,-0.5) * J3

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_quantum_potential()
      IMPLICIT NONE

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_vorticity()
      IMPLICIT NONE

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_grad_f_real_3d_rect_cd2(nXa, nXb, nXbc, nYa, &
         & nYb, nYbc, nZa, nZb, nZbc, dX, dY, dZ, F3, GradF3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      REAL, INTENT(IN) :: dX, dY, dZ

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: F3

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: GradF3

      INTEGER :: j, k, l

      GradF3 = 0.0

!$OMP PARALLEL DO IF ( nZa /= nZb ) DEFAULT ( SHARED ) SCHEDULE ( STATIC )
      DO l = nZa , nZb
!$OMP    PARALLEL DO IF ( nZa == nZb ) DEFAULT ( SHARED ) SCHEDULE ( STATIC )
         DO k = nYa , nYb
            DO j = nXa , nXb
               GradF3(1,j,k,l) = (0.5/dX) * (F3(j+1,k,l) - F3(j-1,k,l))
               GradF3(2,j,k,l) = (0.5/dY) * (F3(j,k+1,l) - F3(j,k-1,l))
               GradF3(3,j,k,l) = (0.5/dZ) * (F3(j,k,l+1) - F3(j,k,l-1))
            END DO
         END DO
!$OMP    END PARALLEL DO
      END DO
!$OMP END PARALLEL DO

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_grad_f_cmplx_3d_rect_cd2(nXa, nXb, nXbc, nYa, &
         & nYb, nYbc, nZa, nZb, nZbc, dX, dY, dZ, F3, GradF3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      REAL, INTENT(IN) :: dX, dY, dZ

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: F3

      COMPLEX, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: GradF3

      INTEGER :: j, k, l

      GradF3 = CMPLX(0.0,0.0)

!$OMP PARALLEL DO IF ( nZa /= nZb ) DEFAULT ( SHARED ) SCHEDULE ( STATIC )
      DO l = nZa , nZb
!$OMP    PARALLEL DO IF ( nZa == nZb ) DEFAULT ( SHARED ) SCHEDULE ( STATIC )
         DO k = nYa , nYb
            DO j = nXa , nXb
               GradF3(1,j,k,l) = CMPLX(0.5/dX,0.0) * (F3(j+1,k,l) - F3(j-1,k,l))
               GradF3(2,j,k,l) = CMPLX(0.5/dY,0.0) * (F3(j,k+1,l) - F3(j,k-1,l))
               GradF3(3,j,k,l) = CMPLX(0.5/dZ,0.0) * (F3(j,k,l+1) - F3(j,k,l-1))
            END DO
         END DO
!$OMP    END PARALLEL DO
      END DO
!$OMP END PARALLEL DO

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_div_f_cmplx_3d_rect_cd2()
      IMPLICIT NONE

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_curl_f_cmplx_3d_rect_cd2()
      IMPLICIT NONE

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE pmca_laplacian_f_cmplx_3d_rect_cd2()
      IMPLICIT NONE

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      REAL FUNCTION pmca_current_3d_rect(nXa, nXb, nXbc, nYa, nYb, &
         & nYbc, nZa, nZb, nZbc, nZ, dX, dY, dZ, xI, xF, yI, yF, zI, &
         & zF, Xa, Ya, Za, Zc, J3) RESULT(pmcaI)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc, nZ

      REAL, INTENT(IN) :: dX, dY, dZ

      REAL, INTENT(IN) :: xI, xF
      REAL, INTENT(IN) :: yI, yF
      REAL, INTENT(IN) :: zI, zF

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: Xa
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Ya
      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(IN) :: Za
      REAL, DIMENSION(  1 - nZbc : nZ  + nZbc), INTENT(IN) :: Zc

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: J3

      INTEGER :: j, k, l

      INTEGER :: nXiG = -1, nXfG = -1
      INTEGER :: nYiG = -1, nYfG = -1
      INTEGER :: nZiG = -1, nZfG = -1
      INTEGER :: nZiL = -1, nZfL = -1

      pmcaI = 0.0

!     Search for global location of integration limits in x, y, and z
      nXiG = grid_linear_search(nXa, nXb, nXbc, xI, Xa)
      nXfG = grid_linear_search(nXa, nXb, nXbc, xF, Xa)

      nYiG = grid_linear_search(nYa, nYb, nYbc, yI, Ya)
      nYfG = grid_linear_search(nYa, nYb, nYbc, yF, Ya)

      nZiG = grid_linear_search(1, nZ, nZbc, zI, Zc)
      nZfG = grid_linear_search(1, nZ, nZbc, zF, Zc)

!     1                        nZiG           nZfG                     nZ
!     |---------------|--------*-------|------*-------|----------------|
!     nZa1        nZb1 nZa2        nZb2 nZa3      nZb3 nZa4         nZa4
!     |------*-----*--|----------------|--------------|----------------|
!            nZiG  nZfG
!     |------*--------|----------------|--------------|-----*----------|
!            nZiG                                           nZfG

!     Search for slab-local location of integration limits in z
      IF ((nZa <= nZiG).AND.(nZiG <= nZb)) THEN
         nZiL = nZiG
      ELSE IF ((nZa > nZiG).AND.(nZa < nZfG)) THEN
         nZiL = nZa
      END IF

      IF ((nZa <= nZfG).AND.(nZfG <= nZb)) THEN
         nZfL = nZfG
      ELSE IF ((nZb > nZiG).AND.(nZb < nZfG)) then
         nZfL = nZb
      END IF

!     Compute surface integrals
      IF (nXiG == nXfG) THEN

         DO l = nZiL, nZfL
            DO k = nYiG, nYfG
               pmcaI = pmcaI + J3(1,nXiG,k,l)
            END DO
         END DO
         pmcaI = pmcaI * dY * dZ

      ELSE IF (nYiG == nYfG) THEN

         DO l = nZiL, nZfL
            DO j = nXiG, nXfG
               pmcaI = pmcaI + J3(2,j,nYiG,l)
            END DO
         END DO
         pmcaI = pmcaI * dX * dZ

      ELSE IF (nZiG == nZfG) THEN

         IF ((nZiG == nZiL).AND.(nZfG == nZfL)) THEN

            DO k = nYiG, nYfG
               DO j = nXiG, nXfG
                  pmcaI = pmcaI + J3(3,j,k,nZiG)
               END DO
            END DO
            pmcaI = pmcaI * dX * dY

         END IF

      ELSE

         WRITE(UNIT=ERROR_UNIT, FMT = *) 'gpse: pmca_current_3d_rect: &
            & ERROR - Surface integral must currently lie in one of &
            & the Cartesian planes.'
         STOP

      END IF

      RETURN
      END FUNCTION

! ----------------------------------------------------------------------

      REAL FUNCTION pmca_velocity_3d_rect(nXa, nXb, nXbc, nYa, nYb, &
         & nYbc, nZa, nZb, nZbc, nZ, dX, dY, dZ, xI, xF, yI, yF, zI, &
         & zF, Xa, Ya, Za, Zc, V3) RESULT(pmcaV)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc, nZ

      REAL, INTENT(IN) :: dX, dY, dZ

      REAL, INTENT(IN) :: xI, xF
      REAL, INTENT(IN) :: yI, yF
      REAL, INTENT(IN) :: zI, zF

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: Xa
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Ya
      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(IN) :: Za
      REAL, DIMENSION(  1 - nZbc : nZ  + nZbc), INTENT(IN) :: Zc

      REAL, DIMENSION(3, &
                    & nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: V3

      INTEGER :: j, k, l

      INTEGER :: nXiG = -1, nXfG = -1
      INTEGER :: nYiG = -1, nYfG = -1
      INTEGER :: nZiG = -1, nZfG = -1
      INTEGER :: nZiL = -1, nZfL = -1

      pmcaV = 0.0

!     Search for global location of integration limits in x, y, and z
      nXiG = grid_linear_search(nXa, nXb, nXbc, xI, Xa)
      nXfG = grid_linear_search(nXa, nXb, nXbc, xF, Xa)

      nYiG = grid_linear_search(nYa, nYb, nYbc, yI, Ya)
      nYfG = grid_linear_search(nYa, nYb, nYbc, yF, Ya)

      nZiG = grid_linear_search(1, nZ, nZbc, zI, Zc)
      nZfG = grid_linear_search(1, nZ, nZbc, zF, Zc)

!     Search for slab-local location of integration limits in z
      IF ((nZa <= nZiG).AND.(nZiG <= nZb)) THEN
         nZiL = nZiG
      ELSE IF ((nZa > nZiG).AND.(nZa < nZfG)) THEN
         nZiL = nZa
      END IF

      IF ((nZa <= nZfG).AND.(nZfG <= nZb)) THEN
         nZfL = nZfG
      ELSE IF ((nZb > nZiG).AND.(nZb < nZfG)) then
         nZfL = nZb
      END IF

!     Compute surface integrals
      IF (nXiG == nXfG) THEN

         DO l = nZiL, nZfL
            DO k = nYiG, nYfG
               pmcaV = pmcaV + V3(1,nXiG,k,l)
            END DO
         END DO
         pmcaV = pmcaV * dY * dZ
         pmcaV = pmcaV / ((Ya(nYfG) - Ya(nYiG)) * (Zc(nZfG) - Zc(nZiG)))

      ELSE IF (nYiG == nYfG) THEN

         DO l = nZiL, nZfL
            DO j = nXiG, nXfG
               pmcaV = pmcaV + V3(2,j,nYiG,l)
            END DO
         END DO
         pmcaV = pmcaV * dX * dZ
         pmcaV = pmcaV / ((Xa(nXfG) - Xa(nXiG)) * (Zc(nZfG) - Zc(nZiG)))

      ELSE IF (nZiG == nZfG) THEN

         IF ((nZiG == nZiL).AND.(nZfG == nZfL)) THEN

            DO k = nYiG, nYfG
               DO j = nXiG, nXfG
                  pmcaV = pmcaV + V3(3,j,k,nZiG)
               END DO
            END DO
            pmcaV = pmcaV * dX * dY
            pmcaV = pmcaV / ((Xa(nXfG) - Xa(nXiG)) * (Ya(nYfG) - Ya(nYiG)))

         END IF

      ELSE

         WRITE(UNIT=ERROR_UNIT, FMT = *) 'gpse: pmca_velocity_3d_rect: &
            & ERROR - Surface integral must currently lie in one of &
            & the Cartesian planes.'
         STOP

      END IF

      RETURN
      END FUNCTION

! ----------------------------------------------------------------------

      END MODULE

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

max_stars_repo_name: mkandes/gpse

max_issues_repo_name: mkandes/gpse

max_forks_repo_name: mkandes/gpse


text:
c\BeginDoc
c
c\Name: pdnorm2
c
c Message Passing Layer: BLACS
c
c\Description:
c
c\Usage:
c  call pdnorm2 ( COMM, N, X, INC )
c
c\Arguments
c  COMM    BLACS Communicator for the processor grid.  (INPUT)
c
c\SCCS Information:
c FILE: norm2.F   SID: 1.2   DATE OF SID: 2/22/96
c
c-----------------------------------------------------------------------
c
      Double precision function pdnorm2 ( comm, n, x, inc )
c
c     %------------------------------%
c     | BLACS Variables and Routines |
c     %------------------------------%
c
      integer    comm
      external   dgsum2d, dgamx2d
c
c     %------------------%
c     | Scalar Arguments |
c     %------------------%
c
      integer      n, inc
c
c     %-----------------%
c     | Array Arguments |
c     %-----------------%
c
      Double precision
     &             x(n)
c
c     %---------------%
c     | Local Scalars |
c     %---------------%
c
      Double precision
     &             max, buf, zero
      parameter    ( zero = 0.0 )
c
c     %---------------------%
c     | Intrinsic Functions |
c     %---------------------%
c
      intrinsic    abs, sqrt
c
c     %--------------------%
c     | External Functions |
c     %--------------------%
c
      Double precision
     &             dnrm2
      External     dnrm2
c
c     %-----------------------%
c     | Executable Statements |
c     %-----------------------%
c
      pdnorm2 = dnrm2( n, x, inc)
c
      max = pdnorm2
      call dgamx2d( comm, 'All', ' ', 1, 1, max, 1, ra, ca,
     &              -1, -1, -1 )
      if ( max .eq. zero ) then
         pdnorm2 = zero
      else
         pdnorm2 = (pdnorm2/max)**2.0
         call dgsum2d( comm, 'All', ' ', 1, 1, pdnorm2, 1, -1, -1 )
         pdnorm2 = max * sqrt(abs(pdnorm2))
      endif
c
c     %----------------%
c     | End of pdnorm2 |
c     %----------------%
c
      return
      end

max_stars_repo_name: JamiePringle/ROMScode-PringleBBLpaper

max_issues_repo_name: JamiePringle/ROMScode-PringleBBLpaper

max_forks_repo_name: JamiePringle/ROMScode-PringleBBLpaper


module linalg
  use iso_fortran_env, only: wp => real64
  implicit none

  public  :: linsolve_quick, linsolve

contains
  subroutine linsolve_quick(n, a, nrhs, b, x)

    ! Quick wrapper around linsolve

    integer,  intent(in)                          :: n, nrhs
    real(wp), intent(in),     dimension(n, n)     :: a
    real(wp), intent(in),     dimension(n, nrhs)  :: b
    real(wp), intent(out),    dimension(n, nrhs)  :: x

    integer,                  dimension(n)        :: P
    real(wp),                 dimension(n, n)     :: LU

    call linsolve (n, a, nrhs, b, x, LU, P, .False.)
    return
  end subroutine linsolve_quick

  subroutine linsolve (n, a, nrhs, b, x, LU, P, toggle)
    ! This routine is a wrapper dgesv, splitting it into its two primary
    ! components:
    !             dgetrf - Decomposes A into P*L*U
    !             dgetrs - Uses P*L*U to solve for x (Ax=b => (P*L*U)x=b)
    !
    ! Splitting these two up like this allows you to save the decomposed
    ! version of 'a' so that you don't have to do it again. If 'toggle' is
    ! equal to true, then the decomposition has already occured and LU can be
    ! trusted - OK to skip dgetrf

    ! Dummy variables
    integer,  intent(in)                          :: n, nrhs
    integer,  intent(inout),  dimension(n)        :: P
    real(wp), intent(in),     dimension(n, nrhs)  :: b
    real(wp), intent(in),     dimension(n, n)     :: a
    real(wp), intent(out),    dimension(n, nrhs)  :: x
    real(wp), intent(inout),  dimension(n, n)     :: LU
    logical,  intent(in)                          :: toggle

    ! Local variables
    integer                                       :: info
    integer,                  dimension(n)        :: my_P
    real(wp),                 dimension(n, n)     :: my_a
    real(wp),                 dimension(n, nrhs)  :: my_b

    my_a = a
    my_b = b

    if ( .not. toggle ) then
      call dgetrf (n, n, my_a, n, my_P, info)

      if ( info /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'DGETRF'
        write ( *, '(a,i8)' ) '  Factorization failed, INFO = ', info
        stop
      end if

      LU  = my_a
      P   = my_P

    end if

    call dgetrs ('No transpose', n, nrhs, LU, n, P, my_b, n, info)

    if ( info /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DGETRS'
      write ( *, '(a,i8)' ) '  Back substitution failed, INFO = ', info
      stop
    end if

    x = my_b

    return
  end subroutine linsolve
end module linalg


text:
The 1991 Davis Bike Plan set the stage for many improvements, including:
  The Dave Pelz Bike Overpass Dave Pelz overpass over Interstate 80
  The Putah Parkway underpass under I80 and the RR tracks
  Bike Lanes on Covell west of 113
and many others.

Each capital improvement (big project that would cost a lot of money) was planned years in advance, and funding sources were identified.

Compare to 2001 http://www.city.davis.ca.us/pw/pdfs/01bikeplan.pdf Bike Plan.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      INTEGER FUNCTION IDT_IEFUND(Pl,Ire)

      IMPLICIT NONE
      INTEGER i , ipla , iple , Ire , j
      DOUBLE PRECISION Pl
      SAVE

C*****IEFUN CALCULATES A MOMENTUM INDEX

      INCLUDE 'inc/dtflka'

      INCLUDE 'inc/hndrun'
      INCLUDE 'inc/hnredv'
      INCLUDE 'inc/hnreac'

      ipla = IEIi(Ire) + 1 + 1
      iple = IEIi(Ire+1)
      IF ( Pl.LT.0. ) THEN
         DO i = ipla , iple
            j = i - ipla + 1
            IF ( -Pl.LE.UMO(i) ) GOTO 100
         END DO
         i = iple
         IF ( EFTes.LE.40.D0 ) THEN
            EFTes = EFTes + 1.0D0

            IF ( LPRi.GT.4 ) WRITE (LOUt,99010) Pl , i
         END IF
      ELSE
         DO i = ipla , iple
            j = i - ipla + 1
            IF ( Pl.LE.PLAbf(i) ) GOTO 100
         END DO
         i = iple
         IF ( EFTes.LE.40.D0 ) THEN
            EFTes = EFTes + 1.0D0

            IF ( LPRi.GT.4 ) WRITE (LOUt,99010) Pl , j
         END IF
      END IF
 100  IDT_IEFUND = i
      RETURN
99010 FORMAT (' PLAB OR -ECM=',E12.4,' IS OUT OF CONSIDERED RANGE',
     &        ' IEFUN=',I5)
      END FUNCTION

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
C$Procedure  ZZEKKEY  ( EK, determine key column )

      SUBROUTINE ZZEKKEY ( HANDLE,  SEGDSC,  NROWS,   NCNSTR,  CLIDXS,
     .                     DSCLST,  OPS,     DTYPES,  CHRBUF,  CBEGS,
     .                     CENDS,   DVALS,   IVALS,   ACTIVE,  KEY,
     .                     KEYDSC,  BEGIDX,  ENDIDX,  FOUND          )

C$ Abstract
C
C     Determine the key column to use when searching an EK segment
C     for rows matching query constraints.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     EK
C
C$ Keywords
C
C     EK
C     PRIVATE
C
C$ Declarations

      INCLUDE 'ekbool.inc'
      INCLUDE 'ekcoldsc.inc'
      INCLUDE 'ekopcd.inc'
      INCLUDE 'ekqlimit.inc'
      INCLUDE 'eksegdsc.inc'
      INCLUDE 'ektype.inc'

      INTEGER               HANDLE
      INTEGER               SEGDSC ( SDSCSZ )
      INTEGER               NROWS
      INTEGER               NCNSTR
      INTEGER               CLIDXS ( * )
      INTEGER               DSCLST ( CDSCSZ, * )
      INTEGER               OPS    ( * )
      INTEGER               DTYPES ( * )
      CHARACTER*(*)         CHRBUF
      INTEGER               CBEGS  ( * )
      INTEGER               CENDS  ( * )
      DOUBLE PRECISION      DVALS  ( * )
      INTEGER               IVALS  ( * )
      LOGICAL               ACTIVE ( * )
      INTEGER               KEY
      INTEGER               KEYDSC ( * )
      INTEGER               BEGIDX
      INTEGER               ENDIDX
      LOGICAL               FOUND

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     HANDLE     I   Handle of EK file containing segment.
C     SEGDSC     I   Segment descriptor.
C     NROWS      I   Number of rows in segment.
C     NCNSTR     I   Number of relational constraints in query.
C     CLIDXS     I   Column attribute table indices for columns.
C     DSCLST     I   Array of column descriptors for constraints.
C     OPS        I   Operations used in query constraints.
C     DTYPES     I   Data types of scalar values used in constraints.
C     CHRBUF     I   Buffer containting query tokens.
C     CBEGS      I   Begin indices of character query tokens.
C     CENDS      I   End indices of character query tokens.
C     DVALS      I   D.p. values used in query constraints.
C     IVALS      I   Integer values used in query constraints.
C     ACTIVE    I-O  Array of flags indicating applicable constraints.
C     KEY        O   Index of key column.
C     KEYDSC     O   Descriptor of key column.
C     BEGIDX     O   Begin index of candidate row set.
C     ENDIDX     O   End index of candidate row set.
C     FOUND      O   Flag indicating whether a key column was found.
C     MAXCON     P   Maximum number of constraints allowed in query.
C     CDSCSZ     P   Column descriptor size.
C
C$ Detailed_Input
C
C     HANDLE         is the handle of the EK file containing the
C                    segment currently being searched for matching rows.
C
C     SEGDSC         the descriptor of the segment.
C
C     NROWS          is the number of rows in the segment designated by
C                    HANDLE and SEGDSC.
C
C     NCNSTR         is the number of input relational constraints.
C                    The input arrays CLIDXS, DSCLST, OPS, CHRBUF,
C                    CBEGS, CENDS, DVALS, and IVALS define these
C                    constraints.  Not all of the constraints may
C                    be applicable; the applicable constraints are
C                    identified by the input argument ACTIVE, which
C                    is described below.  Each *applicable* constraint
C                    has the form
C
C                       <table1.column1> <op> <value>
C
C     CLIDXS         is an array of column indices; the Nth index
C                    identifies the column on the left hand side of the
C                    Nth constraint.  Each index indicates the ordinal
C                    position of the attribute information for the
C                    corresponding column within the column attribute
C                    list for the column's parent table.  See the
C                    local variable declarations in EKBSR for further
C                    information on the column attribute list.
C
C     DSCLST         is an array of column descriptors for the columns
C                    referenced in the input constraints.  The Ith
C                    descriptor corresponds to the Ith constraint.
C
C
C     OPS            are relational operators used in the input
C                    constraints.  The elements of OPS are any of the
C                    integer parameters
C
C                       EQ, GE, GT, LE, LT, NE, LIKE, UNLIKE
C
C                    The Nth element of OPS corresponds to the Nth
C                    constraint.
C
C     DTYPES         is an array of data type codes for the values on
C                    the right hand sides of the input constraints.
C                    The Ith element of DTYPES applies to the Ith
C                    constraint.
C
C     CHRBUF,
C     CBEGS,
C     CENDS          are, respectively, a string containing character
C                    tokens representing values on the right hand sides
C                    of query constraints, and arrays of begin and end
C                    indices of these tokens within CHRBUF.  If the Nth
C                    constraint has a character value on the right hand
C                    side, that value is CHRBUF( CBEGS(N) : CENDS(N) ).
C                    For constraints whose right hand sides do not
C                    specify character values, the corresponding
C                    elements of CBEGS and CENDS are not used.
C
C     DVALS,
C     IVALS          are, respectively, arrays of double precision and
C                    integer values appearing on the right hand sides of
C                    input constraints.  The contents of DVALS and IVALS
C                    are meaningful only for those constraints whose
C                    right hand sides specify values having these data
C                    types.
C
C     ACTIVE         is an array of logical flags indicating which
C                    constraints are currently applicable.  The Nth
C                    element of ACTIVE indicates whether or not to apply
C                    the Nth constraint:  if ACTIVE(N) is .TRUE., the
C                    constraint is applicable, otherwise it isn't.
C
C                    The elements of the other input arguments that
C                    define constraints are defined when the
C                    corresponding element of ACTIVE is .TRUE.  For
C                    example, when the second constraint is not active,
C                    the second column descriptor in DSCLST may not be
C                    defined.
C
C                    Only constraints relating column entries to literal
C                    values may be active.
C
C$ Detailed_Output
C
C     ACTIVE         indicates, on output, which constraints are still
C                    active.  All constraints satisfied by the candidate
C                    row set are turned off on output.
C
C     KEY            is the index of the key column.  This index is
C                    taken from the input argument CLIDXS.
C
C     KEYDSC         is the column descriptor for the key column.  Note
C                    that this descriptor indicates whether the key
C                    column is indexed.
C
C     BEGIDX,
C     ENDIDX         are, respectively, begin and end indices for the
C                    candidate matching rows in the segment being
C                    searched.  These indices refer to positions in the
C                    key column's index:  the candidate rows are pointed
C                    to index elements having indices ranging from
C                    BEGIDX to ENDIDX, inclusively.  The actual
C                    candidate rows are referred to with one level of
C                    indirection.
C
C                    If the constraints on the key column entirely
C                    eliminate all rows in the segment, the returned
C                    values of BEGIDX and ENDIDX are, respectively, 1
C                    and 0.
C
C     FOUND          is a logical flag indicating whether a key column
C                    was determined.  The other outputs of this routine
C                    are valid only if a key column was found.  This
C                    routine will fail to find a key column if there are
C                    no active constraints on indexed columns.
C
C$ Parameters
C
C     MAXCON         is the maximum number of constraints that may
C                    be used in a query.
C
C     CDSCSZ         is the size of a column descriptor.
C
C$ Exceptions
C
C     1)  If the segment contains no indexed columns on which there are
C         active constraints, the output argument FOUND is set to
C         .FALSE.  The other output arguments are undefined in this
C         case.
C
C     2)  If the constraints on the key column entirely eliminate all
C         rows in the segment, the returned values of BEGIDX and ENDIDX
C         are, respectively, 1 and 0.
C
C     3)  If the number of input constraints is out of range, the error
C         SPICE(INVALIDCOUNT) is signalled.
C
C$ Files
C
C     See the description of the input argument HANDLE.
C
C$ Particulars
C
C     The EKSRCH algorithm for finding rows matching a given set
C     of constraints attempts to use constraints on indexed columns
C     to enable the matching process to be performed efficiently.
C     The idea is to find the indexed column whose constraints limit
C     the possible set of matching rows to the smallest number; then
C     to linearly search through this set of candidate rows to see
C     which ones satisfy the remaining applicable constraints.  The
C     column used to initially limit the set of candidate rows is
C     called the `key column'.  The constraints on the key column that
C     are of interest are ones involving order relations or equality:
C     these constraints use the operators
C
C        EQ  GE  GT  LE  LT
C
C     Note that the NE operator is not of much use here.
C
C     The set of candidate rows simultaneously satisifies all such
C     constraints on the key column, and therefore is the intersection
C     of the set of rows satisfying each such constraint.  This method
C     of selecting candidate rows can rapidly eliminate large numbers of
C     rows from consideration, because the index on the key column can
C     be employed in finding rows that match constraints involving order
C     relations:  the start and end indices of such rows can be found
C     by a binary, rather than linear, search.
C
C     A segment may have multiple indexed columns on which there are
C     constraints involving order or equality relations; in this case
C     the column whose constraints are most restrictive is selected as
C     the key column.
C
C     It may also happen that a segment contains no indexed columns.
C     In such a case, the key column is not useful for narrowing the
C     set of candidate rows.  The first column of the segment is
C     arbitrarily selected as the key column in this case.
C
C$ Examples
C
C     See EKSRCH.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.2.0, 26-JUL-1996 (NJB)
C
C        Added check of FAILED after calls to the EK search functions.
C
C-    SPICELIB Version 1.1.0, 17-APR-1996 (WLT)
C
C        Removed spurious periods that appeared at the
C        end of lines 524 and 577 in previous edition.
C
C-    Beta Version 1.0.0, 23-OCT-1995 (NJB)
C
C
C-&

C$ Revisions
C
C-    SPICELIB Version 1.2.0, 26-JUL-1996 (NJB)
C
C        Added check of FAILED after calls to the EK search functions.
C-&

C
C     SPICELIB functions
C
      INTEGER               CARDI
      INTEGER               ORDI

      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Other functions
C
      INTEGER               ZZEKILLE
      INTEGER               ZZEKILLT

C
C     Local parameters
C
      INTEGER               LBCELL
      PARAMETER           ( LBCELL = -5 )

C
C     Local variables
C
      INTEGER               B
      INTEGER               BEST
      INTEGER               COL
      INTEGER               CONMAP ( MAXCON )
      INTEGER               DTYPE
      INTEGER               E
      INTEGER               ELTIDX
      INTEGER               I
      INTEGER               IDXSET ( LBCELL : MAXCON )
      INTEGER               J
      INTEGER               LASTLE
      INTEGER               LASTLT
      INTEGER               MAXPTR
      INTEGER               MINPTR
      INTEGER               NMATCH

      LOGICAL               ELIM
      LOGICAL               FND
      LOGICAL               INDEXD

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'ZZEKKEY' )
      END IF

C
C     There's no key column to begin with.
C
      FOUND = .FALSE.

      IF (  ( NCNSTR .LT. 0 ) .OR. ( NCNSTR .GT. MAXCON )  ) THEN

         CALL SETMSG ( 'The number of constraints was #; valid '  //
     .                 'range is 0:#'                              )
         CALL ERRINT ( '#',   NCNSTR                               )
         CALL ERRINT ( '#',   MAXCON                               )
         CALL SIGERR ( 'SPICE(INVALIDCOUNT)'                       )
         CALL CHKOUT ( 'ZZEKKEY'                                   )
         RETURN

      END IF

C
C     Make a set out of the indices of indexed columns referenced
C     in active constraints.  Maintain a mapping from each column
C     to the index of some constraint that references that column.
C
      CALL SSIZEI ( MAXCON, IDXSET )

      DO I = 1, NCNSTR

         IF ( ACTIVE(I) ) THEN

            INDEXD  =  DSCLST(IXTIDX,I) .NE. IFALSE

            IF ( INDEXD ) THEN
               CALL INSRTI ( CLIDXS(I), IDXSET )
            END IF

         END IF

      END DO

      DO I = 1, CARDI(IDXSET)

         FND = .FALSE.
         J   =  1

         DO WHILE (  ( J .LE. NCNSTR ) .AND. ( .NOT. FND )  )

            IF (   ACTIVE(J) .AND. ( CLIDXS(J) .EQ. IDXSET(I) )  ) THEN
               FND       =   .TRUE.
               CONMAP(I) =   J
            ELSE
               J         =   J + 1
            END IF

         END DO

      END DO

C
C     We finish up now if there are no indexed columns
C     on which there are active constraints.
C
      IF ( CARDI(IDXSET) .EQ. 0 ) THEN

         CALL CHKOUT ( 'ZZEKKEY' )
         RETURN

      END IF

C
C     For each column in the `indexed' set, find out how many
C     candidate rows we'd have if we picked that column as the key
C     column.  If we find that the constraints on some column eliminate
C     all matching rows, we can stop.
C
      BEGIDX  =  1
      ENDIDX  =  NROWS
      BEST    =  IDXSET(1)
      NMATCH  =  NROWS

      ELIM    = .FALSE.
      ELTIDX  =  1

      DO WHILE (  ( ELTIDX .LE. CARDI(IDXSET) )  .AND. ( .NOT. ELIM )  )
C
C        Get the attribute list pointer for the current column.
C
         COL   = IDXSET(ELTIDX)

C
C        Set the initial values of MINPTR, MAXPTR, and NMATCH
C
         MINPTR       =  1
         MAXPTR       =  NROWS
         I            =  1

         DO WHILE (  ( I .LE. NCNSTR ) .AND. ( .NOT. ELIM )  )
C
C           For each constraint, increase MINPTR or decrease MAXPTR
C           if the constraint allows us to do so.
C

            IF (  ( CLIDXS(I) .EQ. COL ) .AND. ( ACTIVE(I) )  )  THEN
C
C              The Ith constraint is active and applies to this column.
C
C              If the column has character type, set the bounds of the
C              token on the right hand side of the constraint.
C              Otherwise, set the bounds to default valid values to
C              avoid subscript bounds errors.
C
               DTYPE  =  DSCLST( TYPIDX, I )

               IF ( DTYPE .EQ. CHR ) THEN
                  B = CBEGS(I)
                  E = CENDS(I)
               ELSE
                  B  =  1
                  E  =  1
               END IF

C
C              At this point, MINPTR and MAXPTR are in the range
C              1:NROWS, and MINPTR is less than or equal to MAXPTR.
C

               IF ( OPS(I) .EQ. LT ) THEN
C
C                 Find the index of the pointer to the last row
C                 whose value in this column is less than the
C                 value cited in the Ith constraint.
C
                  LASTLT  =  ZZEKILLT ( HANDLE,         SEGDSC,
     .                                  DSCLST(1,I),    NROWS,
     .                                  DTYPES(I),      CHRBUF(B:E),
     .                                  DVALS(I),       IVALS(I)  )
C
C                 If all column elements were greater than or equal
C                 to the specified value, MAXPTR will be set to zero.
C
                  MAXPTR  =  MIN ( LASTLT, MAXPTR )
                  ELIM    =  MAXPTR .EQ. 0

               ELSE IF ( OPS(I) .EQ. LE ) THEN
C
C                 Find the index of the pointer to the last row
C                 whose value in this column is less or equal to
C                 the value cited in the Ith constraint.
C
                  LASTLE  =  ZZEKILLE ( HANDLE,         SEGDSC,
     .                                  DSCLST(1,I),    NROWS,
     .                                  DTYPES(I),      CHRBUF(B:E),
     .                                  DVALS(I),       IVALS(I)    )

                  MAXPTR  =  MIN ( LASTLE, MAXPTR )
                  ELIM    =  MAXPTR .EQ. 0

               ELSE IF ( OPS(I) .EQ. EQ ) THEN
C
C                 Find both the pointer to the last row whose
C                 value in this column is less than the value cited in
C                 the Ith constraint, and the pointer to the last row
C                 whose value in this column is less than or equal to
C                 the value cited in the Ith constraint.  The
C                 successor of the former pointer, together with
C                 the latter pointer, bound the range of pointers
C                 to possible matching rows.
C
                  LASTLT  =  ZZEKILLT ( HANDLE,         SEGDSC,
     .                                  DSCLST(1,I),    NROWS,
     .                                  DTYPES(I),      CHRBUF(B:E),
     .                                  DVALS(I),       IVALS(I)    )

                  LASTLE  =  ZZEKILLE ( HANDLE,         SEGDSC,
     .                                  DSCLST(1,I),    NROWS,
     .                                  DTYPES(I),      CHRBUF(B:E),
     .                                  DVALS(I),       IVALS(I)    )

                  IF ( LASTLT .LT. LASTLE ) THEN
C
C                    There is at least one row whose value in the
C                    current column matches the value cited in the Ith
C                    constraint, and LASTLE is the index of the pointer
C                    to the last such row.  The successor of LASTLT is
C                    the first pointer to such a row (even if LASTLT is
C                    zero).
C
                     MINPTR = MAX ( LASTLT+1, MINPTR )
                     MAXPTR = MIN ( LASTLE,   MAXPTR )

                  ELSE
C
C                    No rows match this constraint.
C
                     ELIM  =  .TRUE.

                  END IF

               ELSE IF ( OPS(I) .EQ. GT ) THEN
C
C                 Find the index of the pointer to the last row
C                 whose value in this column is less or equal to
C                 the value cited in the Ith constraint.  The index of
C                 the pointer to the first row satisfying all of the
C                 constraints on this column is the successor of
C                 this pointer or a greater pointer.
C
                  LASTLE  =  ZZEKILLE ( HANDLE,         SEGDSC,
     .                                  DSCLST(1,I),    NROWS,
     .                                  DTYPES(I),      CHRBUF(B:E),
     .                                  DVALS(I),       IVALS(I)    )

                  MINPTR  =  MAX ( LASTLE + 1,  MINPTR )
                  ELIM    =  LASTLE .EQ. NROWS

               ELSE IF ( OPS(I) .EQ. GE ) THEN
C
C                 Find the index of the pointer to the last row
C                 whose value in this column is less than the
C                 value cited in the Ith constraint.  The index of the
C                 pointer to the first row satisfying all of the
C                 constraints on this column is the successor of
C                 this pointer or a greater pointer.
C
                  LASTLT  =  ZZEKILLT ( HANDLE,         SEGDSC,
     .                                  DSCLST(1,I),    NROWS,
     .                                  DTYPES(I),      CHRBUF(B:E),
     .                                  DVALS(I),       IVALS(I)    )

                  MINPTR  =  MAX ( LASTLT + 1, MINPTR )
                  ELIM    =  LASTLT .EQ. NROWS

               END IF
C
C              We've checked the Ith constraint to see whether
C              it applied to the current column, and if it did,
C              we adjusted MINPTR and MAXPTR to reflect the
C              constraint.
C
            END IF
C
C           We've applied the Ith constraint, if it was active.
C
            IF ( MINPTR .GT. MAXPTR ) THEN
               ELIM  =  .TRUE.
            END IF

            IF ( .NOT. ELIM ) THEN
               I  =  I  +  1
            END IF

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'ZZEKKEY' )
               RETURN
            END IF

         END DO

C
C        We've applied all of active, applicable constraints to column
C        COL.  If these constraints did not eliminate the current
C        segment entirely, save the number of candidate rows we'd have
C        if we kept this column as the key column.
C
         IF ( .NOT. ELIM ) THEN

            NMATCH  =  MAXPTR - MINPTR + 1

            IF (  NMATCH .LT. ( ENDIDX - BEGIDX + 1 )  ) THEN
C
C              This is our new key column, until a better one comes
C              along.
C
               BEST    =  COL
               BEGIDX  =  MINPTR
               ENDIDX  =  MAXPTR

            END IF

            ELTIDX  =  ELTIDX + 1

         END IF

      END DO

      IF ( ELIM ) THEN
C
C        If the segment was eliminated by constraints on the last column
C        we looked at, set BEGIDX and ENDIDX to indicate the absence of
C        matching rows.
C
         KEY     =  COL
         BEGIDX  =  1
         ENDIDX  =  0

      ELSE
C
C        BEST, BEGIDX, and ENDIDX are set to reflect the key column.
C        Set KEY and grab the descriptor of the key column.
C
         KEY     =  BEST

      END IF

      I   =  CONMAP(  ORDI( KEY, IDXSET )  )

      CALL MOVEI ( DSCLST(1,I), CDSCSZ, KEYDSC )

C
C     De-activate constraints on the key column that we've already
C     applied.
C
      DO I = 1, NCNSTR

         IF (  ACTIVE(I)  .AND.  ( CLIDXS(I) .EQ. KEY )  )  THEN

            IF (       ( OPS(I) .EQ. LT )
     .           .OR.  ( OPS(I) .EQ. LE )
     .           .OR.  ( OPS(I) .EQ. EQ )
     .           .OR.  ( OPS(I) .EQ. GE )
     .           .OR.  ( OPS(I) .EQ. GT )  ) THEN

C
C              This constraint is met by the candidate rows; we can
C              turn it off.
C
               ACTIVE(I)  =  .FALSE.

            END IF

         END IF

      END DO

C
C     At this point, we've found a key column.
C
      FOUND = .TRUE.

      CALL CHKOUT ( 'ZZEKKEY' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework