text:
! ###################################################################
! Copyright (c) 2013-2021, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - 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.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! 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 HOLDER 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.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:io.f90
!--------------------------------------------------------------------------
!
! MODULE: io
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief IO routines
!
!> @details  Several IO routines; the thinking is that this might be the only file that needs
!> to be rewritten if the package IO is modified.  In particular, in this version of the package,
!> all program input, with the exception of the mkxtal.f90 program and some other simple
!> programs, is done by means of namelist files.  This allows the entire package to be
!> controlled from inside the IDL environment or from the standard command line.
!
!> @date 01/05/99 MDG 1.0 original
!> @date 05/19/01 MDG 2.0 f90 version
!> @date 03/19/13 MDG 3.0 major changes in IO routines (use of interface)
!> @date 05/16/13 MDG 3.1 added stdout as an option to run from IDL
!> @date 06/05/14 MDG 4.0 changed stdout to regular argument; removed global "mess" declaration
!> @date 03/29/18 MDG 4.0 changed stdout handling to instrinsic io definitions in fortran 2003
!> @date 10/02/19 MDG 4.1 changed $ descriptor to advance="no" f03 standard
!--------------------------------------------------------------------------

module io

use local
use, intrinsic :: iso_fortran_env, only : stdin=>input_unit, &
                                          stdout=>output_unit, &
                                          stderr=>error_unit

public

interface ReadValue
        module procedure ReadValueIntShort
        module procedure ReadValueIntLong
        module procedure ReadValueRealSingle
        module procedure ReadValueRealDouble
        module procedure ReadValueString
        module procedure ReadValueStringArray
end interface ReadValue

interface WriteValue
        module procedure WriteValueIntShort
        module procedure WriteValueIntLong
        module procedure WriteValueIntLongLong
        module procedure WriteValueRealSingle
        module procedure WriteValueRealDouble
        module procedure WriteValueRealComplex
        module procedure WriteValueString
end interface WriteValue

contains

!--------------------------------------------------------------------------
!
! SUBROUTINE: Message
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief dump a message to standard output
!
!> @details Simple routine to print a string on the standard output, with optional formatting
!> instructions, for instance if one wants an empty line before (frm='(/A)') or after (frm='(A/)')
!> the string.  Note that one can include the name of the optional variable in the subroutine
!> call, as in:
!> call Message('this is a string', frm='(//A//)' , stdout = 22)
!> this makes it clear that frm and stdout are optional variables.
!
!> @param mess message string
!> @param frm optional string formatting command
!> @param stdout optional output unit identifier
!
!> @date 01/05/99 MDG 1.0 original
!> @date 05/19/01 MDG 2.0 f90 version
!> @date 03/19/13 MDG 3.0 made argument optional and introduced default format
!> @date 06/05/14 MDG 4.0 added stdout and mess as mandatory arguments
!> @date 03/29/18 MDG 4.1 removed stdout argument
!--------------------------------------------------------------------------
recursive subroutine Message(mess,frm,advance)
!DEC$ ATTRIBUTES DLLEXPORT :: Message

character(*),INTENT(IN)                 :: mess         !< message string
character(*),OPTIONAL,INTENT(IN)        :: frm          !< optional formatting string
character(*),OPTIONAL,INTENT(IN)        :: advance      !< optional formatting string

! default format or not ?
if (PRESENT(frm)) then
 if (present(advance)) then
   write (stdout,fmt=frm,advance="no") trim(mess)
 else
   write (stdout,fmt=frm) trim(mess)
 end if
else    ! default output format: a simple string
 write (stdout,fmt="(A)") trim(mess)
end if

end subroutine Message

! ###################################################################
! reading routines
! ###################################################################

! ###################################################################
!
!  subroutine ReadValueString
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief read a string from standard input (unit = 5)
!
!> @param Qstring question string
!> @param rd_string string to be read
!> @param frm optional format string

!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine ReadValueString( Qstring, rd_string, frm)
!DEC$ ATTRIBUTES DLLEXPORT :: ReadValueString

character(*),INTENT(IN)                         :: Qstring
character(*),INTENT(OUT)                        :: rd_string
character(*),INTENT(IN),OPTIONAL                :: frm

call Message(Qstring, frm = "(' ',A,' ')",advance="no")

if (PRESENT(frm)) then
  read (stdin, fmt=frm) rd_string
else
  read (stdin,*) rd_string
end if

end subroutine ReadValueString

! ###################################################################
!
!  subroutine ReadValueStringArray
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief read an array of strings from standard input (unit = 5)
!
!> @param Qstring question string
!> @param rd_string string to be read
!> @param num number of strings in array
!> @param frm optional format string

!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine ReadValueStringArray(Qstring, rd_string, num, frm)
!DEC$ ATTRIBUTES DLLEXPORT :: ReadValueStringArray

character(*),INTENT(IN)                         :: Qstring
character(1),INTENT(OUT)                        :: rd_string(num)
integer(kind=irg),INTENT(IN)                    :: num
character(*),INTENT(IN),OPTIONAL                :: frm

integer(kind=irg)                               :: i

call Message(Qstring, frm = "(' ',A,' ')",advance="no")

if (PRESENT(frm)) then
  do i=1,num
    read (stdin, fmt=frm) rd_string(i)
  end do
else
  do i=1,num
    read (stdin,*) rd_string(i)
  end do
end if

end subroutine ReadValueStringArray

! ###################################################################
!
!  subroutine ReadValueIntShort
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief read one or more short integers
!
!> @param Qstring question string
!> @param rd_int integer to be read
!> @param num optional number of integers to be read

!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine ReadValueIntShort(Qstring, rd_int, num)
!DEC$ ATTRIBUTES DLLEXPORT :: ReadValueIntShort

character(*), INTENT(IN)                        :: Qstring
integer(kind=ish),INTENT(OUT)                   :: rd_int(*)
integer(kind=irg),INTENT(IN),OPTIONAL           :: num

integer(kind=irg)                               :: i

call Message(Qstring, frm = "(' ',A,' ')",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
  read (stdin,*) (rd_int(i),i=1,num)
else
  read (stdin,*) rd_int(1)
end if

end subroutine ReadValueIntShort

! ###################################################################
!
!  subroutine ReadValueIntLong
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief read one or more regular (4-byte) integers
!
!> @param Qstring question string
!> @param rd_int integer to be read
!> @param num optional number of integers to be read

!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine ReadValueIntLong(Qstring, rd_int, num)
!DEC$ ATTRIBUTES DLLEXPORT :: ReadValueIntLong

character(*), INTENT(IN)                        :: Qstring
integer(kind=irg),INTENT(OUT)                   :: rd_int(*)
integer(kind=irg),INTENT(IN),OPTIONAL           :: num

integer(kind=irg)                               :: i

call Message(Qstring, frm = "(' ',A,' ')",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
  read (stdin,*) (rd_int(i),i=1,num)
else
  read (stdin,*) rd_int(1)
end if

end subroutine ReadValueIntLong

! ###################################################################
!
!  subroutine ReadValueRealSingle
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief read one or more regular (4-byte) reals
!
!> @param Qstring question string
!> @param rd_real integer to be read
!> @param num optional number of integers to be read

!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine ReadValueRealSingle(Qstring, rd_real, num)
!DEC$ ATTRIBUTES DLLEXPORT :: ReadValueRealSingle

character(*), INTENT(IN)                        :: Qstring
real(kind=sgl),INTENT(OUT)                      :: rd_real(*)
integer(kind=irg),INTENT(IN),OPTIONAL           :: num

integer(kind=irg)                               :: i

call Message(Qstring, frm = "(' ',A,' ')",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
  read (stdin,*) (rd_real(i),i=1,num)
else
  read (stdin,*) rd_real(1)
end if

end subroutine ReadValueRealSingle

! ###################################################################
!
!  subroutine ReadValueRealDouble
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief read one or more regular (4-byte) reals
!
!> @param Qstring question string
!> @param rd_real integer to be read
!> @param num optional number of integers to be read
!> @param stdout optional output unit identifier

!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine ReadValueRealDouble(Qstring, rd_real, num)
!DEC$ ATTRIBUTES DLLEXPORT :: ReadValueRealDouble

character(*), INTENT(IN)                        :: Qstring
real(kind=dbl),INTENT(OUT)                      :: rd_real(*)
integer(kind=irg),INTENT(IN),OPTIONAL           :: num

integer(kind=irg)                               :: i

call Message(Qstring, frm = "(' ',A,' ')",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
  read (stdin,*) (rd_real(i),i=1,num)
else
  read (stdin,*) rd_real(1)
end if

end subroutine ReadValueRealDouble

! ###################################################################
! writing routines
! ###################################################################

! ###################################################################
!
!  subroutine WriteValueString
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief write a string
!
!> @param Qstring question string
!> @param out_string output string
!> @param frm optional formatting argument
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueString(Qstring, out_string, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueString

character(*),INTENT(IN)                         :: Qstring
character(*),INTENT(IN)                         :: out_string
character(*),INTENT(IN),OPTIONAL                :: frm
character(*),INTENT(IN),OPTIONAL                :: advance

! send Qstring to the output only if it is non-zero length
if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

if (PRESENT(frm)) then
  if (present(advance)) then
    call Message(out_string, frm = frm, advance="no")
  else
    call Message(out_string, frm = frm)
  end if
else
 call Message(out_string, frm = "(A)")
end if

end subroutine WriteValueString

! ###################################################################
!
!  subroutine WriteValueIntShort
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief write one or more short integers
!
!> @param Qstring question string
!> @param out_int output string
!> @param num optional number of integers
!> @param frm optional formatting argument
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueIntShort(Qstring, out_int, num, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueIntShort

character(*), INTENT(IN)                        :: Qstring
integer(kind=ish),INTENT(IN)                    :: out_int(*)
character(*),INTENT(IN),OPTIONAL                :: frm
integer(kind=irg),INTENT(IN),OPTIONAL           :: num
character(*),INTENT(IN),OPTIONAL                :: advance

if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
 if (PRESENT(frm)) then
   if (present(advance)) then
    write (stdout, fmt=frm, advance="no") (out_int(i),i=1,num)
  else
    write (stdout, fmt=frm) (out_int(i),i=1,num)
  end if
 else
  write (stdout,*) (out_int(i),i=1,num)
 end if
else
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") out_int(1)
  else
    write (stdout, fmt=frm) out_int(1)
  end if
 else
  write (stdout,*) out_int(1)
 end if
end if

end subroutine WriteValueIntShort

! ###################################################################
!
!  subroutine WriteValueIntLong
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief write one or more 4-byte integers
!
!> @param Qstring question string
!> @param out_int output string
!> @param num optional number of integers
!> @param frm optional formatting argument
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueIntLong(Qstring, out_int, num, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueIntLong

character(*), INTENT(IN)                        :: Qstring
integer(kind=irg),INTENT(IN)                    :: out_int(*)
character(*),INTENT(IN),OPTIONAL                :: frm
integer(kind=irg),INTENT(IN),OPTIONAL           :: num
character(*),INTENT(IN),OPTIONAL                :: advance

if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") (out_int(i),i=1,num)
  else
    write (stdout, fmt=frm) (out_int(i),i=1,num)
  end if
 else
  write (stdout,*) (out_int(i),i=1,num)
 end if
else
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") out_int(1)
  else
    write (stdout, fmt=frm) out_int(1)
  end if
 else
  write (stdout,*) out_int(1)
 end if
end if

end subroutine WriteValueIntLong

! ###################################################################
!
!  subroutine WriteValueIntLongLong
!
!> @author Saransh, Carnegie Mellon University
!
!> @brief write one or more 8-byte integers
!
!> @param Qstring question string
!> @param out_int output string
!> @param num optional number of integers
!> @param frm optional formatting argument
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueIntLongLong(Qstring, out_int, num, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueIntLongLong

character(*), INTENT(IN)                        :: Qstring
integer(kind=ill),INTENT(IN)                    :: out_int(*)
character(*),INTENT(IN),OPTIONAL                :: frm
integer(kind=irg),INTENT(IN),OPTIONAL           :: num
character(*),INTENT(IN),OPTIONAL                :: advance

if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") (out_int(i),i=1,num)
  else
    write (stdout, fmt=frm) (out_int(i),i=1,num)
  end if
 else
  write (stdout,*) (out_int(i),i=1,num)
 end if
else
 if (PRESENT(frm)) then
    if (present(advance)) then
    write (stdout, fmt=frm, advance="no") out_int(1)
  else
    write (stdout, fmt=frm) out_int(1)
  end if
 else
  write (stdout,*) out_int(1)
 end if
end if

end subroutine WriteValueIntLongLong

! ###################################################################
!
!  subroutine WriteValueRealSingle
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief write one or more single precision reals
!
!> @param Qstring question string
!> @param out_real output string
!> @param num optional number of integers
!> @param frm optional formatting argument
!> @param stdout optional output unit identifier
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueRealSingle(Qstring, out_real, num, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueRealSingle

character(*), INTENT(IN)                        :: Qstring
real(kind=sgl),INTENT(IN)                       :: out_real(*)
character(*),INTENT(IN),OPTIONAL                :: frm
integer(kind=irg),INTENT(IN),OPTIONAL           :: num
character(*),INTENT(IN),OPTIONAL                :: advance

if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") (out_real(i),i=1,num)
  else
    write (stdout, fmt=frm) (out_real(i),i=1,num)
  end if
 else
  write (stdout,*) (out_real(i),i=1,num)
 end if
else
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") out_real(1)
  else
    write (stdout, fmt=frm) out_real(1)
  end if
 else
  write (stdout,*) out_real(1)
 end if
end if

end subroutine WriteValueRealSingle

! ###################################################################
!
!  subroutine WriteValueRealDouble
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief write one or more double precision reals
!
!> @param Qstring question string
!> @param out_real output string
!> @param num optional number of integers
!> @param frm optional formatting argument
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueRealDouble(Qstring, out_real, num, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueRealDouble

character(*), INTENT(IN)                        :: Qstring
real(kind=dbl),INTENT(IN)                       :: out_real(*)
character(*),INTENT(IN),OPTIONAL                :: frm
integer(kind=irg),INTENT(IN),OPTIONAL           :: num
character(*),INTENT(IN),OPTIONAL                :: advance

if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
 if (PRESENT(frm)) then
    if (present(advance)) then
    write (stdout, fmt=frm, advance="no") (out_real(i),i=1,num)
  else
    write (stdout, fmt=frm) (out_real(i),i=1,num)
  end if
 else
  write (stdout,*) (out_real(i),i=1,num)
 end if
else
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") out_real(1)
  else
    write (stdout, fmt=frm) out_real(1)
  end if
 else
  write (stdout,*) out_real(1)
 end if
end if

end subroutine WriteValueRealDouble

! ###################################################################
!
!  subroutine WriteValueRealComplex
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief write one or more single precision complex numbers
!
!> @param Qstring question string
!> @param out_real output string
!> @param num optional number of integers
!> @param frm optional formatting argument
!
!> @date 03/19/13 MDG 1.0 new routine
!> @date 06/05/14 MDG 2.0 changed io handling
!> @date 03/29/18 MDG 4.1 removed stdout argument
! ###################################################################
recursive subroutine WriteValueRealComplex(Qstring, out_cmplx, num, frm, advance)
!DEC$ ATTRIBUTES DLLEXPORT :: WriteValueRealComplex

character(*), INTENT(IN)                        :: Qstring
complex(kind=sgl),INTENT(IN)                    :: out_cmplx(*)
character(*),INTENT(IN),OPTIONAL                :: frm
integer(kind=irg),INTENT(IN),OPTIONAL           :: num
character(*),INTENT(IN),OPTIONAL                :: advance

if (len(Qstring).ne.0) call Message(Qstring, frm = "(A)",advance="no")

! one or more than one values expected ?
if (PRESENT(num)) then
 if (PRESENT(frm)) then
    if (present(advance)) then
    write (stdout, fmt=frm, advance="no") (out_cmplx(i),i=1,num)
  else
    write (stdout, fmt=frm) (out_cmplx(i),i=1,num)
  end if
 else
  write (stdout,*) (out_cmplx(i),i=1,num)
 end if
else
 if (PRESENT(frm)) then
  if (present(advance)) then
    write (stdout, fmt=frm, advance="no") out_cmplx(1)
  else
    write (stdout, fmt=frm) out_cmplx(1)
  end if
 else
  write (stdout,*) out_cmplx(1)
 end if
end if

end subroutine WriteValueRealComplex

recursive subroutine PrintMatrixd(s,a)
!DEC$ ATTRIBUTES DLLEXPORT :: PrintMatrixd

use local

IMPLICIT NONE

real(kind=dbl)   :: a(3,3)
integer(kind=irg):: i,j
character(4)     :: s

write (stdout,"(A/)") s
do i=1,3
  write (stdout,"(3(F12.5,2x))") (a(i,j),j=1,3)
end do
write (stdout,"(/)")

end subroutine

recursive subroutine PrintMatrixcd(s,a)
!DEC$ ATTRIBUTES DLLEXPORT :: PrintMatrixcd

use local

IMPLICIT NONE

complex(kind=dbl)   :: a(3,3)
integer(kind=irg):: i,j
character(4)     :: s

write (stdout,"(A/)") s
do i=1,3
  write (stdout,*) (a(i,j),j=1,3)
end do
write (stdout,"(/)")

end subroutine
!

end module io

max_stars_repo_name: josephtessmer/EMsoft

max_stars_repo_licenses: Intel Unlicense MIT-CMU BSD-3-Clause

max_issues_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_forks_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense


text:
C  *********************************************************
C  *                                                       *
C  *    SUBROUTINE 06.02.02/tstvip                         *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      SUBROUTINE TSTVIP  (SPCPOS, REFPOS, RELPRI, NDVI)

C  TSTVIP tests the PSVTIP function by computing the correct
C  effect of the function call on the current list of view
C  indices, and comparing it with the actual result.
C  Pass or fail is issued as a result.

C  ------- INPUT PARAMETERS
C  spcpos     : position in current list of specified index
C  refpos     : position in current list of reference index
C  relpri     : relative priority of specified index
C  ndvi       : number of defined view indices

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C relative input priority
      INTEGER    PHIGHR,     PLOWER
      PARAMETER (PHIGHR = 0, PLOWER = 1)

      INTEGER   SPCPOS, REFPOS, RELPRI, NDVI, SPCVW, REFVW, IVW
      INTEGER   ARSIZ
      PARAMETER (ARSIZ = 500)
      INTEGER   CURORD(ARSIZ), EXPORD(ARSIZ), ACTORD(ARSIZ)
      INTEGER   ISEND, IRECV, NSPOS, IDUM1

      LOGICAL   IAREQ

C check array size is sufficient
      IF (ARSIZ .LT. NDVI) THEN
         CALL UNMSG ('Allocated array size in TSTVIP is too small.')
      ENDIF

C  Use <inquire list of view indices> to determine
C        curord = current order
      DO 100 IVW = 1,NDVI
         CALL PQEVWI (WKID, IVW, ERRIND, IDUM1, CURORD(IVW))
         CALL CHKINQ ('pqevwi', ERRIND)
100   CONTINUE

C from parameters, calculate expord = expected order
C     nspos is new position of specified view index
      NSPOS = REFPOS
      IF     (SPCPOS .LT. REFPOS .AND. RELPRI .EQ. PHIGHR) THEN
         NSPOS = NSPOS-1
      ELSEIF (SPCPOS .GT. REFPOS .AND. RELPRI .EQ. PLOWER) THEN
         NSPOS = NSPOS+1
      ENDIF

C irecv and isend are receiving and sending positions in curord
C and expord, respectively.
      IRECV=0
      ISEND=0

150   CONTINUE
      ISEND = ISEND+1
      IRECV = IRECV+1
C skip sending of old position of specified view
      IF (ISEND .EQ. SPCPOS) ISEND = ISEND+1
C skip reception of new position of specified view
      IF (IRECV .EQ. NSPOS)  IRECV = IRECV+1
C check if done with loop
      IF (ISEND .GT. NDVI .OR. IRECV .GT. NDVI) GOTO 200
      EXPORD(IRECV) = CURORD(ISEND)
      GOTO 150

200   CONTINUE
C fill in last slot
      EXPORD(NSPOS) = CURORD(SPCPOS)

C  establish identity (not just position) of specified and reference view
      SPCVW = CURORD (SPCPOS)
      REFVW = CURORD (REFPOS)

C  invoke <set view transformation input priority> using
C        spcvw, refvw, relpri
      CALL PSVTIP (WKID, SPCVW, REFVW, RELPRI)

C  Use <inquire list of view indices> to determine
C        actord = actual order
      DO 300 IVW = 1,NDVI
         CALL PQEVWI (WKID, IVW, ERRIND, IDUM1, ACTORD(IVW))
         CALL CHKINQ ('pqevwi', ERRIND)
300   CONTINUE

C  pass/fail depending on (actord = expord)
      CALL IFPF (IAREQ (NDVI, ACTORD, EXPORD) )

      END

ext: f


text:
! Licensed under a 3-clause BSD style license - see LICENSE

module counterpart_pairing_fortran
! This module provides the Fortran code for the final act of cross-matching, in which
! sources are paired within their respective islands, and counterparts are established.

use, intrinsic :: ieee_arithmetic, only: IEEE_VALUE, IEEE_QUIET_NAN

implicit none

integer, parameter :: dp = kind(0.0d0)  ! double precision
real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp)

contains

subroutine find_island_probabilities(a_astro, a_photo, b_astro, b_photo, alist, alist_, blist, blist_, agrplen, bgrplen, c_array, &
    fa_array, fb_array, c_priors, fa_priors, fb_priors, amagref, bmagref, amodrefind, bmodrefind, abinsarray, abinlengths, &
    bbinsarray, bbinlengths, afrac_grids, aflux_grids, bfrac_grids, bflux_grids, afourier_grids, bfourier_grids, a_sky_inds, &
    b_sky_inds, rho, drho, n_fracs, large_len, mini_crpts_len, acountinds, bcountinds, afieldinds, bfieldinds, acontamprob, &
    bcontamprob, etaarray, xiarray, acontamflux, bcontamflux, probcarray, crptseps, probfaarray, afieldfluxs, afieldseps, &
    afieldetas, afieldxis, probfbarray, bfieldfluxs, bfieldseps, bfieldetas, bfieldxis)
    ! Calculate the final cross-matches. For each island, unique permutations of match and non-match are
    ! compared, with the highest probability combination stored. The matches and non-matches then have
    ! various secondary parameters, such as log-likelihood-ratio of match or their contamination probability
    ! derived and returned.
    integer, parameter :: dp = kind(0.0d0)  ! double precision
    ! Small and full-index versions of the island group numbers, and their lengths in the two catalogues.
    integer, intent(in) :: alist(:, :), alist_(:, :), blist(:, :), blist_(:, :), agrplen(:), bgrplen(:)
    ! Magnitude bin lengths for filters in the two catalogues.
    integer, intent(in) :: abinlengths(:, :), bbinlengths(:, :)
    ! Indices into filter, specific perturbation AUF model, and sky slice, for the respective catalogues.
    integer, intent(in) :: amagref(:), bmagref(:), amodrefind(:, :), bmodrefind(:, :), a_sky_inds(:), b_sky_inds(:)
    ! Number of contaminant fractions simulated in the perturbation AUF, and total lengths of small and
    ! large catalogues.
    integer, intent(in) :: n_fracs, large_len, mini_crpts_len
    ! TODO: *binsarray, *binslength, *_priors and *_array are all memmapped. Does this matter in practice?
    ! Should be able to take unique-sky-slice cutouts, but this'll need convenience functions for the slicing.
    ! Astrometry and photometry for the two respective catalogues.
    real(dp), intent(in) :: a_astro(:, :), a_photo(:, :), b_astro(:, :), b_photo(:, :)
    ! Photometric "counterpart" and "field" likelihood and prior arrays.
    real(dp), intent(in) :: c_array(:, :, :, :, :), fa_array(:, :, :, :), fb_array(:, :, :, :), c_priors(:, :, :), &
    fa_priors(:, :, :), fb_priors(:, :, :)
    ! Magnitude bins for both catalogues.
    real(dp), intent(in) :: abinsarray(:, :, :), bbinsarray(:, :, :)
    ! Perturbation AUF derivates: fraction of contaminated objects, average flux contamination, and AUF
    ! component fourier transforms for both catalogues.
    real(dp), intent(in) :: afrac_grids(:, :, :, :), aflux_grids(:, :, :), afourier_grids(:, :, :, :), bfrac_grids(:, :, :, :), &
    bflux_grids(:, :, :), bfourier_grids(:, :, :, :)
    ! Fourier-space grids for perturbation AUF.
    real(dp), intent(in) :: rho(:), drho(:)
    ! Indices of matches and non-matches in the two catalogues.
    integer, intent(out) :: acountinds(mini_crpts_len), bcountinds(mini_crpts_len), &
    afieldinds(size(a_astro, 1)), bfieldinds(size(b_astro, 1))
    ! Contamination probabilities of each source in the pairs of matches.
    real(dp), intent(out) :: acontamprob(n_fracs, mini_crpts_len), bcontamprob(n_fracs, mini_crpts_len)
    ! Eta and Xi -- log-likelihood-ratios in astrometry/photometry -- of each match, based on a one-to-one match.
    ! In the non-match case, nearest opposing island object's position and brightness are used.
    real(dp), intent(out) :: etaarray(mini_crpts_len), xiarray(mini_crpts_len), afieldetas(size(a_astro, 1)), &
    afieldxis(size(a_astro, 1)), bfieldetas(size(b_astro, 1)), bfieldxis(size(b_astro, 1))
    ! Average simulated contaminating flux for each match or non-match.
    real(dp), intent(out) :: acontamflux(mini_crpts_len), bcontamflux(mini_crpts_len), afieldfluxs(size(a_astro, 1)), &
    bfieldfluxs(size(b_astro, 1))
    ! Probability of (non-)match.
    real(dp), intent(out) :: probcarray(mini_crpts_len), probfaarray(size(a_astro, 1)), probfbarray(size(b_astro, 1))
    ! Separations between pair of matches, or of nearest opposing island object if a non-match.
    real(dp), intent(out) :: crptseps(mini_crpts_len), afieldseps(size(a_astro, 1)), bfieldseps(size(b_astro, 1))
    ! Indices into islands, or to map each island to a unique part of the larger array to avoid index clashes.
    integer :: i, j, match_cumulat_ind(size(alist, 2)), afield_cumulat_ind(size(alist, 2)), bfield_cumulat_ind(size(alist, 2))
    ! NaN, and temporary non-normalised posterior, and evidence, holders.
    real(dp) :: nan, prob, integral
    ! Smaller arrays holding one islands' worth of each output described above.
    real(dp), allocatable :: acrptscontp(:, :), bcrptscontp(:, :), etacrpts(:), xicrpts(:), acrptflux(:), &
    bcrptflux(:), crptsep(:), afieldflux(:), bfieldflux(:), afieldsep(:), afieldeta(:), afieldxi(:), &
    bfieldsep(:), bfieldeta(:), bfieldxi(:)
    integer, allocatable :: acrpts(:), bcrpts(:), afield(:), bfield(:)

    nan = IEEE_VALUE(nan, IEEE_QUIET_NAN)

    ! Start by calculating the cumulative maximum index of each island, assuming either 0%
    ! or 100% match rate.
    match_cumulat_ind(1) = 1
    afield_cumulat_ind(1) = 1
    bfield_cumulat_ind(1) = 1
    do i = 2, size(alist, 2)
        match_cumulat_ind(i) = match_cumulat_ind(i-1) + min(agrplen(i-1), bgrplen(i-1))
        afield_cumulat_ind(i) = afield_cumulat_ind(i-1) + agrplen(i-1)
        bfield_cumulat_ind(i) = bfield_cumulat_ind(i-1) + bgrplen(i-1)
    end do

    acountinds(:) = large_len + 1
    bcountinds(:) = large_len + 1
    etaarray(:) = -100.0_dp
    xiarray(:) = -100.0_dp
    acontamprob(:, :) = -100.0_dp
    bcontamprob(:, :) = -100.0_dp
    acontamflux(:) = -100.0_dp
    bcontamflux(:) = -100.0_dp
    probcarray(:) = -100.0_dp
    crptseps(:) = -100.0_dp

    afieldinds(:) = large_len + 1
    afieldfluxs(:) = -100.0_dp
    afieldseps(:) = -100.0_dp
    afieldetas(:) = -100.0_dp
    afieldxis(:) = -100.0_dp
    probfaarray(:) = -100.0_dp

    bfieldinds(:) = large_len + 1
    bfieldfluxs(:) = -100.0_dp
    bfieldseps(:) = -100.0_dp
    bfieldetas(:) = -100.0_dp
    bfieldxis(:) = -100.0_dp
    probfbarray(:) = -100.0_dp

!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i, j, acrpts, bcrpts, acrptscontp, bcrptscontp, etacrpts, xicrpts, acrptflux, bcrptflux, &
!$OMP& crptsep, afield, bfield, afieldflux, bfieldflux, afieldsep, afieldeta, afieldxi, bfieldsep, bfieldeta, bfieldxi, prob, &
!$OMP& integral) SHARED(a_astro, a_photo, b_astro, b_photo, c_array, fa_array, fb_array, c_priors, fa_priors, fb_priors, &
!$OMP& abinsarray, bbinsarray, abinlengths, bbinlengths, afrac_grids, aflux_grids, afourier_grids, bfrac_grids, bflux_grids, &
!$OMP& bfourier_grids, rho, drho, n_fracs, large_len, alist, blist, alist_, blist_, amagref, a_sky_inds, bmagref, b_sky_inds, &
!$OMP& amodrefind, bmodrefind, acountinds, bcountinds, afieldinds, bfieldinds, acontamprob, bcontamprob, etaarray, xiarray, &
!$OMP& acontamflux, bcontamflux, probcarray, crptseps, probfaarray, afieldfluxs, afieldseps, afieldetas, afieldxis, probfbarray, &
!$OMP& bfieldfluxs, bfieldseps, bfieldetas, bfieldxis, agrplen, bgrplen, match_cumulat_ind, afield_cumulat_ind, &
!$OMP& bfield_cumulat_ind, nan)
    do i = 1, size(alist, 2)
        if (agrplen(i) == 0) then
            bfieldinds(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = blist_(:bgrplen(i), i)
            probfbarray(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = 1
            do j = 1, bgrplen(i)
                bfieldfluxs(bfield_cumulat_ind(i)+j-1) = bflux_grids( &
                    bmodrefind(1, blist(j, i)+1)+1, bmodrefind(2, blist(j, i)+1)+1, bmodrefind(3, blist(j, i)+1)+1)
            end do
            bfieldseps(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = nan
            bfieldetas(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = nan
            bfieldxis(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = nan
        else if (bgrplen(i) == 0) then
            afieldinds(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = alist_(:agrplen(i), i)
            probfaarray(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = 1
            do j = 1, agrplen(i)
                afieldfluxs(afield_cumulat_ind(i)+j-1) = aflux_grids( &
                    amodrefind(1, alist(j, i)+1)+1, amodrefind(2, alist(j, i)+1)+1, amodrefind(3, alist(j, i)+1)+1)
            end do
            afieldseps(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = nan
            afieldetas(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = nan
            afieldxis(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = nan
        else
            ! Allocate all of the outputs from find_single_island_prob...
            allocate( &
                acrpts(min(agrplen(i), bgrplen(i))), bcrpts(min(agrplen(i), bgrplen(i))), &
                acrptscontp(n_fracs, min(agrplen(i), bgrplen(i))), bcrptscontp(n_fracs, min(agrplen(i), bgrplen(i))), &
                etacrpts(min(agrplen(i), bgrplen(i))), xicrpts(min(agrplen(i), bgrplen(i))), &
                acrptflux(min(agrplen(i), bgrplen(i))), bcrptflux(min(agrplen(i), bgrplen(i))), &
                crptsep(min(agrplen(i), bgrplen(i))), afield(agrplen(i)), bfield(bgrplen(i)), afieldflux(agrplen(i)), &
                bfieldflux(bgrplen(i)), afieldsep(agrplen(i)), afieldeta(agrplen(i)), afieldxi(agrplen(i)), &
                bfieldsep(bgrplen(i)), bfieldeta(bgrplen(i)), bfieldxi(bgrplen(i)))
            ! Lots of +1 indexing here as the pointers were all generated in Python. Note that *list_ does NOT
            ! get +1 indexed because it's never used and only ends up documenting the full-scale indexes of
            ! (non-)matches, and might as well stay in Python-indexing throughout.
            call find_single_island_prob(a_astro, a_photo, b_astro, b_photo, c_array, fa_array, fb_array, c_priors, fa_priors, &
                fb_priors, abinsarray, bbinsarray, abinlengths, bbinlengths, afrac_grids, aflux_grids, afourier_grids, &
                bfrac_grids, bflux_grids, bfourier_grids, rho, drho, n_fracs, large_len, alist(:agrplen(i), i)+1, &
                blist(:bgrplen(i), i)+1, alist_(:agrplen(i), i), blist_(:bgrplen(i), i), amagref(alist(:agrplen(i), i)+1)+1, &
                a_sky_inds(alist(:agrplen(i), i)+1)+1, bmagref(blist(:bgrplen(i), i)+1)+1, b_sky_inds(blist(:bgrplen(i), i)+1)+1, &
                amodrefind(:, alist(:agrplen(i), i)+1)+1, bmodrefind(:, blist(:bgrplen(i), i)+1)+1, acrpts, bcrpts, acrptscontp, &
                bcrptscontp, etacrpts, xicrpts, acrptflux, bcrptflux, crptsep, afield, bfield, &
                afieldflux, bfieldflux, afieldsep, afieldeta, afieldxi, bfieldsep, bfieldeta, bfieldxi, prob, integral)
            ! Similar to in the main counterpart_pairing loop, we can assume either 100% or 0%
            ! match rates, and handle the missing values later. This assignment is therefore just
            ! a version of "know how large each island is, and jump the N/M a/b sources for each
            ! island".
            afieldinds(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = afield
            probfaarray(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = prob / integral
            afieldfluxs(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = afieldflux
            afieldseps(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = afieldsep
            afieldetas(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = afieldeta
            afieldxis(afield_cumulat_ind(i):afield_cumulat_ind(i)+agrplen(i)-1) = afieldxi

            bfieldinds(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = bfield
            probfbarray(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = prob / integral
            bfieldfluxs(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = bfieldflux
            bfieldseps(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = bfieldsep
            bfieldetas(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = bfieldeta
            bfieldxis(bfield_cumulat_ind(i):bfield_cumulat_ind(i)+bgrplen(i)-1) = bfieldxi

            acountinds(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = acrpts
            bcountinds(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = bcrpts
            etaarray(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = etacrpts
            xiarray(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = xicrpts
            acontamprob(:, match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = acrptscontp
            bcontamprob(:, match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = bcrptscontp
            acontamflux(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = acrptflux
            bcontamflux(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = bcrptflux
            probcarray(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = prob / integral
            crptseps(match_cumulat_ind(i):match_cumulat_ind(i)+size(acrpts)-1) = crptsep

            deallocate(acrpts, bcrpts, acrptscontp, bcrptscontp, etacrpts, xicrpts, acrptflux, bcrptflux, crptsep, afield, bfield, &
            afieldflux, bfieldflux, afieldsep, afieldeta, afieldxi, bfieldsep, bfieldeta, bfieldxi)
        end if
    end do
!$OMP END PARALLEL DO

end subroutine find_island_probabilities

subroutine find_single_island_prob(a_astro, a_photo, b_astro, b_photo, c_array, fa_array, fb_array, c_priors, fa_priors, &
    fb_priors, abinsarray, bbinsarray, abinlengths, bbinlengths, afrac_grids, aflux_grids, afourier_grids, bfrac_grids, &
    bflux_grids, bfourier_grids, rho, drho, n_fracs, large_len, aperm, bperm, aperm_, bperm_, aused, qa, bused, qb, arefinds, &
    brefinds, acrpts, bcrpts, acrptscontp, bcrptscontp, etacrpts, xicrpts, acrptflux, bcrptflux, crptseps, afield, bfield, &
    afieldflux, bfieldflux, afieldseps, afieldeta, afieldxi, bfieldseps, bfieldeta, bfieldxi, prob, integral)
    ! Calculate a single island match and non-match permutation. For each combination of N catalogue "a"
    ! and M catalogue "b" sources, find the pairings and non-pairings that maximise the likelihood. Also
    ! calculates secondary parameters such as likelihood ratios for astrometry and photometry and
    ! contamination probabilities for matches.
    integer, parameter :: dp = kind(0.0d0)  ! double precision
    integer, parameter :: int64 = selected_int_kind(15)  ! 64-bit integer
    ! Small- and large-index versions of the indices of the objects in each catalogue in this island.
    integer, intent(in) :: aperm(:), bperm(:), aperm_(:), bperm_(:)
    ! Magnitude and sky filter indices for catalogues 'a' and 'b'.
    integer, intent(in) :: aused(:), qa(:), bused(:), qb(:)
    ! Perturbation AUF grid indices for each source in the two catalogues.
    integer, intent(in) :: arefinds(:, :), brefinds(:, :)
    ! Lengths of the magnitude bin arrays for each catalogue.
    integer, intent(in) :: abinlengths(:, :), bbinlengths(:, :)
    ! Number of contamination fractions recorded in each simulated perturbation AUF, and largest index value.
    integer, intent(in) :: n_fracs, large_len
    ! Astrometry and photometry for the two respective catalogues.
    real(dp), intent(in) :: a_astro(:, :), a_photo(:, :), b_astro(:, :), b_photo(:, :)
    ! Photometric "counterpart" and "field" likelihood and prior arrays.
    real(dp), intent(in) :: c_array(:, :, :, :, :), fa_array(:, :, :, :), fb_array(:, :, :, :), c_priors(:, :, :), &
    fa_priors(:, :, :), fb_priors(:, :, :)
    ! Magnitude bins for both catalogues.
    real(dp), intent(in) :: abinsarray(:, :, :), bbinsarray(:, :, :)
    ! Perturbation AUF derivates: fraction of contaminated objects, average flux contamination, and AUF
    ! component fourier transforms for both catalogues.
    real(dp), intent(in) :: afrac_grids(:, :, :, :), aflux_grids(:, :, :), afourier_grids(:, :, :, :), bfrac_grids(:, :, :, :), &
    bflux_grids(:, :, :), bfourier_grids(:, :, :, :)
    ! Fourier-space grids for perturbation AUF.
    real(dp), intent(in) :: rho(:), drho(:)
    ! Contamination probabilities of matches and non-matches from the island.
    real(dp), intent(out) :: acrptscontp(n_fracs, min(size(aperm), size(bperm))), &
    bcrptscontp(n_fracs, min(size(aperm), size(bperm)))
    ! Log-likelihood-ratios for astrometry and photometry for matches and non-matches. For non-matches,
    ! the nearest neighbour is chosen.
    real(dp), intent(out) :: etacrpts(min(size(aperm), size(bperm))), xicrpts(min(size(aperm), size(bperm))), &
    afieldeta(size(aperm)), afieldxi(size(aperm)), bfieldeta(size(bperm)), bfieldxi(size(bperm))
    ! Average contaminating simulated flux in the perturbation AUF generated for these (non-)matches.
    real(dp), intent(out) :: acrptflux(min(size(aperm), size(bperm))), bcrptflux(min(size(aperm), size(bperm))), &
    afieldflux(size(aperm)), bfieldflux(size(bperm))
    ! Separations between matches and non-matches. For non-matches, the nearest neighbour is picked.
    real(dp), intent(out) :: crptseps(min(size(aperm), size(bperm))), afieldseps(size(aperm)), bfieldseps(size(bperm))
    ! Non-normalised posterior and accompanying normalising evidence for the maximum-likelihood match.
    real(dp), intent(out) :: prob, integral
    ! Pairs of matched indices, and individual non-match indices, for sources in this island.
    integer, intent(out) :: acrpts(min(size(aperm), size(bperm))), bcrpts(min(size(aperm), size(bperm))), &
    afield(size(aperm)), bfield(size(bperm))
    ! Loop counter and indices.
    integer :: i, j, k, l, ff, N, t_ind
    ! Number of permutations or combinations for a given number of chosen items from an array.
    integer :: a_n_combs, b_n_combs, b_n_perms_per_comb
    ! Indices into magnitude bins for each source in the islands.
    integer :: bina(size(aperm)), binb(size(bperm))
    ! Counterpart and "field" probability for all combinations of match between opposing catalogue objects.
    real(dp) :: counterpartgrid(size(aperm), size(bperm)), afieldarray(size(aperm)), bfieldarray(size(bperm))
    ! Log-likelihood-ratio storage arrays for astrometry and photometry.
    real(dp) :: etagrid(size(aperm), size(bperm)), xigrid(size(aperm), size(bperm))
    ! Average simulated contaminating flux for each source in the island.
    real(dp) :: acontamfluxgrid(size(aperm)), bcontamfluxgrid(size(bperm))
    ! "Field" source prior and likelihood for each object in the island.
    real(dp) :: Nfa(size(aperm)), Nfb(size(bperm)), fa(size(aperm)), fb(size(bperm))
    ! Storage for average flux contamination and probability of flux contamination of potential matches.
    real(dp) :: aF(n_fracs), bF(n_fracs), acontamprobgrid(n_fracs, size(aperm), size(bperm)), &
    bcontamprobgrid(n_fracs, size(aperm), size(bperm))
    ! Fourier-space representations of centroiding component of AUF, and its corresponding centroid
    ! uncertainty, and temporary storage of a particular simulated perturbation AUF component.
    real(dp) :: afourgauss(size(drho)), bfourgausses(size(rho)-1, size(bperm)), aoffs(size(drho)), asig, boffs(size(drho)), bsig
    ! Real-space probabilities of match in various contamination combinations, and fourier-space
    ! representations of AUFs.
    real(dp) :: G, G0cc(size(drho)), G0cn(size(drho)), G0nc(size(drho)), G0nn(size(drho)), Gcc, Gcn, Gnc, Gnn
    ! Temporary storage: counterpart prior and likelihood, match evidence, posterior probability,
    ! and sky separation; NaN storage.
    real(dp) :: Nc, cdmdm, pr, tempprob, sep, tempsep, seps(size(aperm), size(bperm)), nan
    ! Arrays to keep track of which island indices are matches and which non-matches for a given permutation.
    integer, allocatable :: ya(:), yb(:), ta(:), tb(:)
    ! Flag for removing the ya or yb indices from aperm or bperm to get ta or tb.
    logical :: flag
    ! Combinations of catalogue a indices, and permutations of catalogue b indices, to provide the overall
    ! set of all a-b matches and non-matches.
    integer, allocatable :: a_combinations(:, :), b_permutations(:, :)
    ! Temporary storage of 64-bit integer factorial numbers.
    integer(int64) :: fac1, fac2, fac3

    nan = IEEE_VALUE(nan, IEEE_QUIET_NAN)

    do j = 1, size(aperm)
        bina(j) = maxloc(a_photo(aperm(j), aused(j)) - abinsarray(:abinlengths(aused(j), qa(j)), aused(j), qa(j)), &
                         mask=(a_photo(aperm(j), aused(j)) - abinsarray(:abinlengths(aused(j), qa(j)), aused(j), qa(j)) >= 0), &
                         dim=1)
        ! For the field sources we don't know which other filter to use, so we
        ! just default to using the first filter in the opposing catalogue,
        ! but it shouldn't matter since it ought to be independent.
        Nfa(j) = fa_priors(1, aused(j), qa(j))
        fa(j) = fa_array(bina(j), 1, aused(j), qa(j))
        acontamfluxgrid(j) = aflux_grids(arefinds(1, j), arefinds(2, j), arefinds(3, j))
    end do

    do j = 1, size(bperm)
        binb(j) = maxloc(b_photo(bperm(j), bused(j)) - bbinsarray(:bbinlengths(bused(j), qb(j)), bused(j), qb(j)), &
                         mask=(b_photo(bperm(j), bused(j)) - bbinsarray(:bbinlengths(bused(j), qb(j)), bused(j), qb(j)) >= 0), &
                         dim=1)
        Nfb(j) = fb_priors(bused(j), 1, qb(j))
        fb(j) = fb_array(binb(j), bused(j), 1, qb(j))
        bcontamfluxgrid(j) = bflux_grids(brefinds(1, j), brefinds(2, j), brefinds(3, j))
    end do

    do j = 1, size(bperm)
        bsig = b_astro(bperm(j), 3)
        bfourgausses(:, j) = exp(-2.0_dp * pi**2 * (rho(:size(rho)-1)+drho/2.0_dp)**2 * bsig**2)
    end do

    afieldarray = Nfa*fa
    bfieldarray = Nfb*fb

    do j = 1, size(aperm)
        aF = afrac_grids(:, arefinds(1, j), arefinds(2, j), arefinds(3, j))
        aoffs = afourier_grids(:, arefinds(1, j), arefinds(2, j), arefinds(3, j))
        asig = a_astro(aperm(j), 3)
        afourgauss = exp(-2.0_dp * pi**2 * (rho(:size(rho)-1)+drho/2.0_dp)**2 * asig**2)
        do k = 1, size(bperm)
            call haversine(a_astro(aperm(j), 1), b_astro(bperm(k), 1), a_astro(aperm(j), 2), b_astro(bperm(k), 2), sep)
            ! sep comes out of haversine in degrees, but contam_match_prob
            ! assumes everything is in arcseconds, so convert sep here.
            sep = sep * 3600.0_dp
            bF = bfrac_grids(:, brefinds(1, k), brefinds(2, k), brefinds(3, k))
            boffs = bfourier_grids(:, brefinds(1, k), brefinds(2, k), brefinds(3, k))

            ! Calculate the probability densities of all four combinations
            ! of perturbation and non-perturbation AUFs.
            G0nn = afourgauss*bfourgausses(:, k)
            G0cc = aoffs*boffs*G0nn
            G0cn = aoffs*G0nn
            G0nc = boffs*G0nn
            call contam_match_prob(G0cc, G0cn, G0nc, G0nn, rho(:size(rho)-1)+drho/2.0_dp, drho, sep, Gcc, Gcn, Gnc, Gnn)
            ! G would be in units of per square arcseconds, but we need it
            ! in units of per square degree to compare to Nf.
            G = Gcc * 3600.0_dp**2

            do ff = 1, n_fracs
                pr = aF(ff)*bF(ff)*Gcc + aF(ff)*(1-bF(ff))*Gcn + (1-aF(ff))*bF(ff)*Gnc + (1-aF(ff))*(1-bF(ff))*Gnn
                ! Marginalise over the opposite source contamination probability
                ! to calculate specific source's contamination chance.
                acontamprobgrid(ff, j, k) = min(1.0_dp, max(0.0_dp, (aF(ff)*bF(ff)*Gcc + aF(ff)*(1-bF(ff))*Gcn)/pr))
                bcontamprobgrid(ff, j, k) = min(1.0_dp, max(0.0_dp, (aF(ff)*bF(ff)*Gcc + (1-aF(ff))*bF(ff)*Gnc)/pr))
            end do
            Nc = c_priors(bused(k), aused(j), qb(k))
            cdmdm = c_array(binb(k), bina(j), bused(k), aused(j), qb(k))
            counterpartgrid(j, k) = Nc * G * cdmdm

            if (abs(fa(j) * fb(k)) < 1e-7_dp) then
                etagrid(j, k) = 10.0_dp
            else if (abs(cdmdm) < 1e-7_dp) then
                etagrid(j, k) = -10.0_dp
            else
                etagrid(j, k) = log10(cdmdm / (fa(j) * fb(k)))
            end if

            if (abs(Nfa(j) * Nfb(k)) < 1e-7_dp) then
                xigrid(j, k) = 10.0_dp
            else if (abs(Nc * G) < 1e-7_dp) then
                xigrid(j, k) = -10.0_dp
            else
                xigrid(j, k) = log10(Nc * G / (Nfa(j) * Nfb(k)))
            end if

            seps(j, k) = sep
        end do
    end do

    ! Start with the case of no matches between any island objects. Assume
    ! by default that no sources match in the island. Add a small normalising constant just in case.
    tempprob = product(afieldarray) * product(bfieldarray) + 1e-10_dp
    integral = tempprob
    prob = tempprob

    ! Default values for arrays need assigning to counterpart grids, and the field arrays fully
    ! filling in to initialise everything.
    acrpts(:) = large_len + 1
    bcrpts(:) = large_len + 1
    acrptscontp(:, :) = -100.0_dp
    bcrptscontp(:, :) = -100.0_dp
    etacrpts(:) = -100.0_dp
    xicrpts(:) = -100.0_dp
    acrptflux(:) = -100.0_dp
    bcrptflux(:) = -100.0_dp
    crptseps(:) = -100.0_dp

    afield(:) = aperm_(:)
    bfield(:) = bperm_(:)
    afieldflux(:) = acontamfluxgrid(:)
    bfieldflux(:) = bcontamfluxgrid(:)

    do j = 1, size(aperm)
        tempsep = 1e10_dp
        do k = 1, size(bperm)
            if (seps(j, k) < tempsep) then
                afieldseps(j) = seps(j, k)
                afieldeta(j) = etagrid(j, k)
                afieldxi(j) = xigrid(j, k)
                tempsep = seps(j, k)
            end if
        end do
    end do

    do k = 1, size(bperm)
        tempsep = 1e10_dp
        do j = 1, size(aperm)
            if (seps(j, k) < tempsep) then
                bfieldseps(k) = seps(j, k)
                bfieldeta(k) = etagrid(j, k)
                bfieldxi(k) = xigrid(j, k)
                tempsep = seps(j, k)
            end if
        end do
    end do

    do N = 1, min(size(aperm), size(bperm))
        ! Combinations of m-pick-n are m! / n! / (m - n)!
        call factorial(int(size(aperm), int64), fac1)
        call factorial(int(N, int64), fac2)
        call factorial(int(size(aperm)-N, int64), fac3)
        a_n_combs = int(fac1 / fac2 / fac3)
        call factorial(int(size(bperm), int64), fac1)
        call factorial(int(N, int64), fac2)
        call factorial(int(size(bperm)-N, int64), fac3)
        b_n_combs = int(fac1 / fac2 / fac3)
        ! Permutations of an n-sized array (n-pick-n) are n!
        call factorial(int(N, int64), fac1)
        b_n_perms_per_comb = int(fac1)
        allocate(a_combinations(N, a_n_combs), b_permutations(N, b_n_combs*b_n_perms_per_comb))
        call calc_combs(size(aperm), a_n_combs, N, a_combinations)
        call calc_permcombs(size(bperm), N, b_n_perms_per_comb, b_n_combs, b_permutations)
        allocate(ya(size(a_combinations, 1)), yb(size(b_permutations, 1)))
        allocate(ta(size(aperm, 1) - size(a_combinations, 1)), tb(size(bperm, 1) - size(b_permutations, 1)))
        do i = 1, size(a_combinations, 2)
            do j = 1, size(b_permutations, 2)
                ! For paired sources, order matters, so we have to find the
                ! index of the array holding the source's overall catalogue
                ! index that matches that index in the permutation list.
                ya = a_combinations(:, i)
                yb = b_permutations(:, j)
                ! This should be the same as np.delete(np.arange(0, len(aperm)), ya), in which we look
                ! at each aperm index (np.arange(0, len(aperm))), ask if it's in ya at all, and if it
                ! isn't then we keep it in ta.
                t_ind = 1
                do k = 1, size(aperm)
                    flag = .True.
                    do l = 1, size(ya)
                        if (k == ya(l)) then
                            flag = .False.
                        end if
                    end do
                    if (flag) then
                        ta(t_ind) = k
                        t_ind = t_ind + 1
                    end if
                end do
                t_ind = 1
                do k = 1, size(bperm)
                    flag = .True.
                    do l = 1, size(yb)
                        if (k == yb(l)) then
                            flag = .False.
                        end if
                    end do
                    if (flag) then
                        tb(t_ind) = k
                        t_ind = t_ind + 1
                    end if
                end do
                tempprob = 1
                do k = 1, size(ya)
                    tempprob = tempprob * counterpartgrid(ya(k), yb(k))
                end do
                do k = 1, size(ta)
                    tempprob = tempprob * afieldarray(ta(k))
                end do
                do k = 1, size(tb)
                    tempprob = tempprob * bfieldarray(tb(k))
                end do
                integral = integral + tempprob
                if (tempprob > prob) then
                    prob = tempprob
                    ! Reset the grids to include "default" values:
                    acrpts(:) = large_len + 1
                    bcrpts(:) = large_len + 1
                    acrptscontp(:, :) = -100.0_dp
                    bcrptscontp(:, :) = -100.0_dp
                    etacrpts(:) = -100.0_dp
                    xicrpts(:) = -100.0_dp
                    acrptflux(:) = -100.0_dp
                    bcrptflux(:) = -100.0_dp
                    crptseps(:) = -100.0_dp

                    afield(:) = large_len + 1
                    bfield(:) = large_len + 1
                    afieldflux(:) = -100.0_dp
                    bfieldflux(:) = -100.0_dp
                    do k = 1, size(ya)
                        acrpts(k) = aperm_(ya(k))
                        acrptflux(k) = acontamfluxgrid(ya(k))
                        bcrpts(k) = bperm_(yb(k))
                        bcrptflux(k) = bcontamfluxgrid(yb(k))
                        acrptscontp(:, k) = acontamprobgrid(:, ya(k), yb(k))
                        bcrptscontp(:, k) = bcontamprobgrid(:, ya(k), yb(k))
                        etacrpts(k) = etagrid(ya(k), yb(k))
                        xicrpts(k) = xigrid(ya(k), yb(k))
                        crptseps(k) = seps(ya(k), yb(k))
                    end do
                    do k = 1, size(ta)
                        afield(k) = aperm_(ta(k))
                        afieldflux(k) = acontamfluxgrid(ta(k))
                    end do
                    do k = 1, size(tb)
                        bfield(k) = bperm_(tb(k))
                        bfieldflux(k) = bcontamfluxgrid(tb(k))
                    end do
                    afieldseps(:) = nan
                    afieldeta(:) = nan
                    afieldxi(:) = nan
                    ! TODO: Do we want these to be nearest of all sources, or nearest non-paired object?
                    do k = 1, size(ta)
                        tempsep = 1e10_dp
                        do l = 1, size(tb)
                            if (seps(ta(k), tb(l)) < tempsep) then
                                afieldseps(k) = seps(ta(k), tb(l))
                                afieldeta(k) = etagrid(ta(k), tb(l))
                                afieldxi(k) = xigrid(ta(k), tb(l))
                                tempsep = seps(ta(k), tb(l))
                            end if
                        end do
                    end do
                    bfieldseps(:) = nan
                    bfieldeta(:) = nan
                    bfieldxi(:) = nan
                    do l = 1, size(tb)
                        tempsep = 1e10_dp
                        do k = 1, size(ta)
                            if (seps(ta(k), tb(l)) < tempsep) then
                                bfieldseps(l) = seps(ta(k), tb(l))
                                bfieldeta(l) = etagrid(ta(k), tb(l))
                                bfieldxi(l) = xigrid(ta(k), tb(l))
                                tempsep = seps(ta(k), tb(l))
                            end if
                        end do
                    end do
                end if
            end do
        end do
        deallocate(a_combinations, b_permutations)
        deallocate(ya, yb, ta, tb)
    end do

end subroutine find_single_island_prob

subroutine factorial(N, g)
    ! Calculates N! = N * (N-1) * (N-2) * ... * 1, the factorial of N.
    integer, parameter :: int64 = selected_int_kind(15)  ! 64-bit integer
    ! The number to calculate the factorial of.
    integer(int64), intent(in) :: N
    ! The returned value g = N!.
    integer(int64), intent(out) :: g
    ! Loop counter.
    integer(int64) :: k

    g = 1
    do k = 1, N
        g = g * k
    end do

end subroutine factorial

subroutine contam_match_prob(Fcc, Fcn, Fnc, Fnn, rho, drho, sep, Gcc, Gcn, Gnc, Gnn)
    ! Calculate the probability density of two sources being two detections of a single sky object
    ! given the separation between them, including the combinations of hypotheses that the sources
    ! independently suffer, or do not suffer, from a statistical description of blended source
    ! perturbation.
    integer, parameter :: dp = kind(0.0d0)  ! double precision
    ! Combinations of fourier-space representations of convolutions of AUFs. "c" and "n" represent
    ! "contaminated" and "non-contaminated", for catalogue "a" then "b" for F[a][b]. Thus the four
    ! arrays are the four combinations of AUFs that suffer perturbation or not for the two sources.
    real(dp), intent(in) :: Fcc(:), Fcn(:), Fnc(:), Fnn(:)
    ! Fourier-space representation of the F** arrays.
    real(dp), intent(in) :: rho(:), drho(:)
    ! Sky separation, in arcseconds, between the two objects in question.
    real(dp), intent(in) :: sep
    ! Output probability densities -- again, per combination of (non-)contamination -- that the
    ! two objects are a single source detected twice given their on-sky separation.
    real(dp), intent(out) :: Gcc, Gcn, Gnc, Gnn
    ! Loop counter.
    integer :: j
    ! Hankel transform parameters.
    real(dp) :: j0, z
    Gcc = 0.0_dp
    Gcn = 0.0_dp
    Gnc = 0.0_dp
    Gnn = 0.0_dp

    do j = 1, size(rho)
        z = rho(j)*sep*2.0_dp*pi
        call jy01a_j0(z, j0)
        Gcc = Gcc + rho(j) * Fcc(j) * j0 * drho(j)
        Gcn = Gcn + rho(j) * Fcn(j) * j0 * drho(j)
        Gnc = Gnc + rho(j) * Fnc(j) * j0 * drho(j)
        Gnn = Gnn + rho(j) * Fnn(j) * j0 * drho(j)
    end do

    Gcc = Gcc * 2.0_dp * pi
    Gcn = Gcn * 2.0_dp * pi
    Gnc = Gnc * 2.0_dp * pi
    Gnn = Gnn * 2.0_dp * pi

end subroutine contam_match_prob

recursive subroutine perm(i, a, perm_grid, n, ind)
    ! Iterate over a to find all permutations of the array. Adapted from
    ! http://rosettacode.org/wiki/Permutations#Fortran.
    ! Iteration counter, and size of array a.
    integer, intent(in) :: i, n
    ! Temporary storage array, array into which to save all permutations, and recursion counter.
    integer, intent(inout) :: a(:), perm_grid(:, :), ind
    ! Loop counter and temporary storage variable.
    integer :: j, t
    if (i == n) then
        perm_grid(:, ind) = a
        ind = ind + 1
    else
        do j = i, n
            t = a(i)
            a(i) = a(j)
            a(j) = t
            call perm(i + 1, a, perm_grid, n, ind)
            t = a(i)
            a(i) = a(j)
            a(j) = t
        end do
    end if

end subroutine

subroutine calc_combs(n, n_combs, m_pick, comb_grid)
    ! Calculate all combinations of length m_pick from an array of length n.
    ! Size of array, total number of combinations, and size of sub-set array to select combinations of.
    integer, intent(in) :: n, n_combs, m_pick
    ! Output grid containing all n-pick-m combinations.
    integer, intent(out) :: comb_grid(m_pick, n_combs)
    ! Temporary array to store m_pick-combinations to shuffle, and a recursive index counter.
    integer :: a(m_pick), ind
    ind = 1
    call comb(1, a, comb_grid, n, ind, m_pick)

end subroutine calc_combs

recursive subroutine comb(m, a, comb_grid, n_max, ind, m_pick)
    ! Iterate over all m_pick-combinations of a (of length n_max), storing them in comb_grid.
    ! Adapted from https://rosettacode.org/wiki/Combinations#Fortran.
    ! Recursive counter, size of combinations to select, and size of array to pick items from.
    integer, intent (in) :: m, m_pick, n_max
    ! Temporary storage for swapping items, array to save combinations to, and recursion counter.
    integer, intent(inout) :: a(:), comb_grid(:, :), ind
    ! Loop counter.
    integer :: n

    if (m > m_pick) then
      comb_grid(:, ind) = a
      ind = ind + 1
    else
      do n = 1, n_max
        if (m == 1) then
          a(m) = n
          call comb(m + 1, a, comb_grid, n_max, ind, m_pick)
        else if (n > a(m - 1)) then
          a(m) = n
          call comb(m + 1, a, comb_grid, n_max, ind, m_pick)
        end if
      end do
    end if

end subroutine comb

subroutine calc_permcombs(n, m_pick, n_perms, n_combs, combperm_grid)
    ! Calculate all permutations of selecting M items from an array of length N.
    ! Length of input array, length of output permutations array, and the number of combinations from
    ! n-pick-m and the number of permutations of each unique combination for the array.
    integer, intent(in) :: n, m_pick, n_perms, n_combs
    ! Grid iterating over every permutation of m_pick items from an n-length array.
    integer, intent(out) :: combperm_grid(m_pick, n_perms*n_combs)
    ! Loop counters.
    integer :: i, ind
    ! Temporary storage to shuffle combinations, and storage for all unique combinations of n-pick-m.
    integer :: a(m_pick), comb_grid(m_pick, n_combs)
    ind = 1
    call comb(1, a, comb_grid, n, ind, m_pick)
    ind = 1
    do i = 1, n_combs
        call perm(1, comb_grid(:, i), combperm_grid, m_pick, ind)
    end do

end subroutine calc_permcombs

end module counterpart_pairing_fortran

max_stars_repo_name: lsst-uk/macauff

max_issues_repo_name: lsst-uk/macauff

max_forks_repo_name: lsst-uk/macauff


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
cdis
cdis
cdis

        subroutine read_profiles(i4time_sys,heights_3d,             ! I
     1                   lat_pr,lon_pr,obstype,c5_name_a,           ! O
     1                   lat,lon,topo,i4time_ob_pr,                 ! I
     1                   MAX_PR,MAX_PR_LEVELS,                      ! I
     1                   l_use_raob,l_use_all_nontower_lvls,        ! I
     1                   ob_pr_u , ob_pr_v ,                        ! O
     1                   max_obs,obs_point,nobs_point,weight_prof,  ! I/O
     1                   nlevels_obs_pr, n_profiles,                ! O
     1                   rlat_radar,rlon_radar,rheight_radar,       ! I
     1                   n_vel_grids,                               ! I
     1                   u_mdl_bkg_4d,v_mdl_bkg_4d,NTMIN,NTMAX,     ! I
     1                   ilaps_cycle_time,r_missing_data,           ! I
     1                   imax,jmax,kmax,                            ! I
     1                   istatus                )                   ! O

!       1992 Steve Albers
!       Note that the profiler data in the .PRO files are in knots...
!       1994 Keith Brewster   Added reading of sounding data
c       1995 Keith Brewster   Re-Added reading of sounding data, improved
c                             error handling
c       1996 Steve Albers     Added read of ob times from PRO files
c       1996 Steve Albers     Read nearest PRO file, even if its time does
c                             not exactly match the LAPS analysis time.
c                             Accept only those profiler obs whose
c                             mid-window times are within one LAPS cycle
c                             time of the current LAPS analysis time.
c	2006 Yuanfu Xie	      Use of the fraction grid values of obs_point.

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

        use mem_namelist, ONLY: iwrite_output

        include 'barnesob.inc'
        type (barnesob) :: obs_point(max_obs)

!       LAPS Grid Dimensions

        real lat(imax,jmax)
        real lon(imax,jmax)
        real topo(imax,jmax)

!       Profile Stuff
        real lat_pr(MAX_PR)
        real lon_pr(MAX_PR)
        real elev_pr(MAX_PR)
        real rcycles_pr(MAX_PR)
        integer i4time_ob_pr(MAX_PR)

!       Profiler Observations

        integer nlevels_obs_pr(MAX_PR)
        character*8 obstype(MAX_PR)
        real ob_pr_ht_obs(MAX_PR,MAX_PR_LEVELS)                             ! L
        real ob_pr_pr_obs(MAX_PR,MAX_PR_LEVELS)                             ! L
!       real ob_pr_di_obs(MAX_PR_LEVELS)                                    ! L
!       real ob_pr_sp_obs(MAX_PR_LEVELS)                                    ! L
        real ob_pr_u_obs(MAX_PR,MAX_PR_LEVELS)                              ! L
        real ob_pr_v_obs(MAX_PR,MAX_PR_LEVELS)                              ! L
        real ob_pr_rms_obs(MAX_PR,MAX_PR_LEVELS)                            ! L
        real ob_pr_t_obs(MAX_PR,MAX_PR_LEVELS)                              ! L
        real ob_pr_td_obs(MAX_PR,MAX_PR_LEVELS)                             ! L
        real ob_pr_ht(MAX_PR,kmax)                                          ! L
        real ob_pr_di(MAX_PR,kmax)                                          ! L
        real ob_pr_sp(MAX_PR,kmax)                                          ! L
        real ob_pr_u (MAX_PR,kmax) ! Vertically interpolated Profiler wind  ! O
        real ob_pr_v (MAX_PR,kmax) ! Vertically interpolated Profiler wind  ! O

!       Profiler surface data (not currently used)
        real sfc_t(MAX_PR), sfc_p(MAX_PR), sfc_rh(MAX_PR)                   ! L
        real sfc_u(MAX_PR), sfc_v(MAX_PR)                                   ! L

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

        real heights_3d(imax,jmax,kmax)

        dimension u_mdl_bkg_4d(imax,jmax,kmax,NTMIN:NTMAX)
        dimension v_mdl_bkg_4d(imax,jmax,kmax,NTMIN:NTMAX)

        character*255 c_filespec
        character ext*31
        character*5 c5_name, c5_name_a(MAX_PR)
        character*9 a9time_ob
        character*12 c_obstype

        logical l_use_raob, l_use_all_nontower_lvls

        r_mspkt = .518

        write(6,*)' Subroutine read_profiles: i4time = ',i4time_sys

!       Initialize

        do i_pr = 1,MAX_PR
            nlevels_obs_pr(i_pr) = 0
        enddo

        do i_pr = 1,MAX_PR
            do level = 1,kmax

                ob_pr_ht(i_pr,level) = r_missing_data
                ob_pr_di(i_pr,level) = r_missing_data
                ob_pr_sp(i_pr,level) = r_missing_data
                ob_pr_u(i_pr,level)  = r_missing_data
                ob_pr_v(i_pr,level)  = r_missing_data

            enddo
        enddo

        i4time_prg = i4time_sys

        ext = 'prg'
        if(iwrite_output .ge. 1)then
            call open_lapsprd_file(32,i4time_prg,ext,istatus)
            if(istatus .ne. 1)go to 880
        endif

! ***   Read in profiler data    ***************************************

!       Open nearest PRO file to the LAPS analysis time
        ext = 'pro'
        call get_filespec(ext,2,c_filespec,istatus)
        call get_file_time(c_filespec,i4time_sys,i4time_prof)

        lun = 12
        call read_pro_data(lun,i4time_prof,ext                        ! I
     1                         ,MAX_PR,MAX_PR_LEVELS                  ! I
     1                         ,n_profiles                            ! O
     1                         ,nlevels_obs_pr,lat_pr,lon_pr,elev_pr  ! O
     1                         ,c5_name_a,i4time_ob_pr,obstype        ! O
     1                         ,ob_pr_ht_obs                          ! O
     1                         ,ob_pr_u_obs,ob_pr_v_obs               ! O
     1                         ,ob_pr_rms_obs                         ! O
     1                         ,sfc_t,sfc_p,sfc_rh,sfc_u,sfc_v        ! O
     1                         ,istatus)                              ! O

        n_profilers = n_profiles

c ***   Read in sonde data    ***************************************
c

      write(6,*)

      if(.not. l_use_raob)then
          write(6,*)' Not using raobs, l_use_raob = ',l_use_raob
      endif

      i4time_raob_file_window = 0

      ext = 'snd'
      call get_filespec(ext,2,c_filespec,istatus)
      call get_file_time(c_filespec,i4time_sys,i4time_nearest)

      i4time_diff = abs(i4time_sys - i4time_nearest)
      if(i4time_diff .le. i4time_raob_file_window)then
          write(6,*)' Nearest SND file is within time window'
     1                ,i4time_diff,i4time_raob_file_window
      else
          write(6,*)' Warning: Nearest SND file is outside time window'
     1                ,i4time_diff,i4time_raob_file_window
          go to 600
      endif

      i4time_snd = i4time_nearest

      lun = 12
      mode = 1 ! key levels off of wind data
      call read_snd_data(lun,i4time_snd,ext                             ! I
     1                         ,MAX_PR,MAX_PR_LEVELS                    ! I
     1                         ,lat,lon,topo,imax,jmax,kmax             ! I
     1                         ,heights_3d,.true.                       ! I
     1                         ,mode                                    ! I
     1                         ,n_profiles                              ! I/O
     1                         ,nlevels_obs_pr,lat_pr,lon_pr,elev_pr    ! O
     1                         ,c5_name_a,i4time_ob_pr,obstype          ! O
     1                         ,ob_pr_ht_obs,ob_pr_pr_obs               ! O
     1                         ,ob_pr_u_obs,ob_pr_v_obs                 ! O
     1                         ,ob_pr_t_obs,ob_pr_td_obs                ! O
     1                         ,istatus)                                ! O

 600  continue

      n_snd=n_profiles-n_profilers

      write(6,*)
      write(6,*) ' Read ',n_profilers,' wind profiler(s).'
      write(6,*) ' Read ',n_snd,' sounding(s).'
      write(6,*)
c
c     Process all wind profiles.  Interpolate heights to LAPS levels.
c
      DO i_pr=1,n_profiles

            if(i_pr .le. 10)then
                idebug = 1
            else
                idebug = 0
            endif

            rcycles_pr(i_pr) = float(i4time_sys - i4time_ob_pr(i_pr))
     1                                      / float(ilaps_cycle_time)

            if(i_pr .le. 200 .or. i_pr .eq. (i_pr/10)*10)then
                iwrite = 1
                istatus = 1
            else
                iwrite = 0
                istatus = 100 ! Suppress write statement in 'latlon_to_rlapsgrid'
            endif

!           Determine if profile is in the LAPS domain

            call latlon_to_rlapsgrid(lat_pr(i_pr),lon_pr(i_pr),lat,lon
     1                              ,imax,jmax,ri,rj,istatus)

            i_ob = nint(ri)
            j_ob = nint(rj)
            if(i_ob .ge. 1 .and. i_ob .le. imax .and.
     1         j_ob .ge. 1 .and. j_ob .le. jmax      )then
                if(iwrite .eq. 1)
     1             write(6,*)'Profile  # ',i_pr,' In Bounds - Doing '
     1                   ,'Vertical Interpolation'
            else
                if(iwrite .eq. 1)
     1             write(6,*)'Profile  # ',i_pr,' Out of Domain Bounds'
                nlevels_obs_pr(i_pr)=0 ! This effectively throws out the profile
            endif

            call get_windob_time_window(obstype(i_pr),i4_window_ob
     1                                               ,istatus)

            rcyc_thresh = float(i4_window_ob)
     1                   /float(ilaps_cycle_time)

            rcyc_thresh = min(1.0,rcyc_thresh)

!           Determine if profile was obtained close enough in time....
            if(abs(rcycles_pr(i_pr)) .gt. rcyc_thresh)then
                if(iwrite .eq. 1)
     1             write(6,*)'Profile  # ',i_pr,' Out of time bounds:'
     1                                    ,rcycles_pr(i_pr)
                nlevels_obs_pr(i_pr)=0 ! This effectively throws out the profile
            endif

!  ***  Interpolate the profiles to the LAPS grid levels  *******

            if(nlevels_obs_pr(i_pr) .gt. 0)then

              if(l_use_all_nontower_lvls .OR.
     1           obstype(i_pr)(1:5) .eq. 'TOWER' .OR.
     1           obstype(i_pr)(1:5) .eq. 'SODAR'
     1                                           )then ! remap all levels

                write(6,311)i_pr,i_ob,j_ob,nlevels_obs_pr(i_pr)
     1                     ,obstype(i_pr),c5_name_a(i_pr)
 311            format(1x,' Remapping profile ',4i6,1x,a8,1x,a5
     1                ,' (all levels)')

                rklaps_min = r_missing_data
                rklaps_max = 0.
                n_good_lvl = 0

                do lvl = 1,nlevels_obs_pr(i_pr)
                    ob_height = ob_pr_ht_obs(i_pr,lvl)
                    ob_u      = ob_pr_u_obs(i_pr,lvl)
                    ob_v      = ob_pr_v_obs(i_pr,lvl)
                    rklaps = height_to_zcoord2(ob_height
     1                                        ,heights_3d,imax,jmax,kmax
     1                                        ,i_ob,j_ob,istatus)
                    klaps = nint(rklaps)

                    if(istatus .eq. 1)then
                        rklaps_min = min(rklaps,rklaps_min)
                        rklaps_max = max(rklaps,rklaps_max)
                        n_good_lvl = n_good_lvl + 1

!                       Obtain time terms
                        call get_time_term(u_mdl_bkg_4d,imax,jmax,kmax
     1                                    ,NTMIN,NTMAX
     1                                    ,i_ob,j_ob,klaps
     1                                    ,i4time_sys,i4time_ob_pr(i_pr)
     1                                    ,u_time_interp,u_diff_term
     1                                    ,istatus)

!                       u_diff_term = du/dt * [t(ob) - t(anal)]
!                       u_diff      = du/dt * [t(anal) - t(ob)]
                        u_diff = -u_diff_term

                        call get_time_term(v_mdl_bkg_4d,imax,jmax,kmax
     1                                    ,NTMIN,NTMAX
     1                                    ,i_ob,j_ob,klaps
     1                                    ,i4time_sys,i4time_ob_pr(i_pr)
     1                                    ,v_time_interp,v_diff_term
     1                                    ,istatus)
!                       v_diff_term = dv/dt * [t(ob) - t(anal)]
!                       v_diff      = dv/dt * [t(anal) - t(ob)]
                        v_diff = -v_diff_term

!                       Add to data structure (full sampling)
                        nobs_point = nobs_point + 1
                        obs_point(nobs_point)%i = i_ob
                        obs_point(nobs_point)%j = j_ob
                        obs_point(nobs_point)%k = klaps
                        obs_point(nobs_point)%ri = ri	! Yuanfu
                        obs_point(nobs_point)%rj = rj	! Yuanfu
                        obs_point(nobs_point)%rk = rklaps
                        obs_point(nobs_point)%valuef(1) = ob_u + u_diff
                        obs_point(nobs_point)%valuef(2) = ob_v + v_diff
                        obs_point(nobs_point)%weight = weight_prof
                        obs_point(nobs_point)%vert_rad_rat = 1.0
                        if(obstype(i_pr)(1:5) .eq. 'TOWER')then
                            call downcase(obstype(i_pr),c_obstype)
                            obs_point(nobs_point)%type = c_obstype
                        else
                            obs_point(nobs_point)%type = 'prof'
                        endif
                        obs_point(nobs_point)%file = 'pro'
                        obs_point(nobs_point)%i4time =
     1                                        i4time_ob_pr(i_pr)
                        if(obstype(i_pr)(1:4) .eq. 'RAOB')then
                          if(.not. l_use_raob)then
                            obs_point(nobs_point)%l_withhold = .true.
                          endif
                        endif
                    endif ! istatus

                    call uv_to_disp(ob_u,ob_v,ob_di,ob_sp)

                    if(iwrite_output .ge. 1)then
312                     write(32,313,err=314)ri,rj,rklaps,ob_di,ob_sp
     1                                      ,obstype(i_pr)
313                     format(1x,3f10.5,2f10.3,1x,a8)
314                     continue
                    endif

                enddo ! lvl

!               Calculate mean # of laps levels between sounding levels
                if(n_good_lvl .gt. 1)then
                    rklaps_range = rklaps_max - rklaps_min
                    rklaps_mean = rklaps_range / float(n_good_lvl-1)
                    n_cross = int(rklaps_range)
                    wt_factor = (min(rklaps_mean/0.5,1.0))
                    write(6,*)
     1                ' n_good_lvl,n_cross,rk-range/mean,wt_factor'
     1                 ,n_good_lvl,n_cross,rklaps_range,rklaps_mean
     1                 ,wt_factor

!                   Adjust the weights for this profile?
                endif

              else ! interpolate from levels to LAPS grid
                do level = 1,kmax

                    ht = heights_3d(i_ob,j_ob,level)

                    ob_pr_ht(i_pr,level) = ht

                    call get_time_term(u_mdl_bkg_4d,imax,jmax,kmax
     1                                ,NTMIN,NTMAX
     1                                ,i_ob,j_ob,level
     1                                ,i4time_sys,i4time_ob_pr(i_pr)
     1                                ,u_time_interp,u_diff_term
     1                                ,istatus)

!                   u_diff_term = du/dt * [t(ob) - t(anal)]
!                   u_diff      = du/dt * [t(anal) - t(ob)]
                    u_diff = -u_diff_term

                    call get_time_term(v_mdl_bkg_4d,imax,jmax,kmax
     1                                ,NTMIN,NTMAX
     1                                ,i_ob,j_ob,level
     1                                ,i4time_sys,i4time_ob_pr(i_pr)
     1                                ,v_time_interp,v_diff_term
     1                                ,istatus)
!                   v_diff_term = dv/dt * [t(ob) - t(anal)]
!                   v_diff      = dv/dt * [t(anal) - t(ob)]
                    v_diff = -v_diff_term

                    call interp_prof(ob_pr_ht_obs,ob_pr_u_obs,     ! I
     1                               ob_pr_v_obs,                  ! I
     1                               u_diff,                       ! I
     1                               v_diff,                       ! I
     1                               ob_pr_u(i_pr,level),          ! O
     1                               ob_pr_v(i_pr,level),          ! O
     1                               ob_pr_di(i_pr,level),         ! O
     1                               ob_pr_sp(i_pr,level),         ! O
     1                               i_pr,ht,level,nlevels_obs_pr, ! I
     1                               lat_pr,lon_pr,i_ob,j_ob,      ! I
     1                               r_missing_data,               ! I
     1                               heights_3d,imax,jmax,kmax,    ! I
     1                               MAX_PR,MAX_PR_LEVELS,         ! I
     1                               n_vel_grids,istatus)          ! I/O

                    if(idebug .eq. 1)write(6,411,err=412)i_pr,level
     1                              ,ob_pr_ht(i_pr,level)
     1                              ,ob_pr_di(i_pr,level)
     1                              ,ob_pr_sp(i_pr,level)
     1                              ,ob_pr_u(i_pr,level)
     1                              ,ob_pr_v(i_pr,level)
     1                              ,u_diff
     1                              ,v_diff
411                 format(1x,2i4,f8.1,6f7.1)

412                 if(iwrite_output .ge. 1)then
                        write(32,313,err=414)ri,rj,float(level)
     1                        ,ob_pr_di(i_pr,level),ob_pr_sp(i_pr,level)
     1                        ,obstype(i_pr)
414                     continue
                    endif

                enddo ! level

              endif ! use all levels

            else
              if(iwrite .eq. 1)then
                  write(6,*)' This profile is set to 0 levels',i_pr
     1                      ,obstype(i_pr)
              endif

            endif ! # levels > 0

        enddo  ! i_pr

        close(32)
        istatus=1
        return

  880   CONTINUE
        write(6,*) ' Error opening PRG file'
        istatus=0
        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
SUBROUTINE start
!
!   Start or restart the run
!
  IMPLICIT NONE
!
!   Local vars and arrays
!________________________________________________________________________________
  WRITE(*,'(a/)') '=== Start or restart the run ==='
!________________________________________________________________________________
!
END SUBROUTINE start

max_stars_repo_name: alexbovet/fastpart

max_issues_repo_name: alexbovet/fastpart

max_forks_repo_name: alexbovet/fastpart