text:
module restraints_constraints_xml_module
use system_module, only : dp, system_abort, operator(//)
use constraints_module
use dynamicalsystem_module
use quip_common_module
implicit none
private

type(DynamicalSystem), pointer, private :: parse_ds
logical, private :: parse_in_restraints = .false., parse_in_constraints = .false.

public :: init_restraints_constraints

contains

   subroutine restraints_constraints_startElement_handler(URI, localname, name, attributes)
      character(len=*), intent(in)   :: URI
      character(len=*), intent(in)   :: localname
      character(len=*), intent(in)   :: name
      type(dictionary_t), intent(in) :: attributes

      character(len=1024) :: value
      integer :: status, di_status, tau_status, t0_status
      integer :: n
      integer :: atom_1, atom_2, atom_3, Z
      real(dp) :: k
      character(len=1024) :: bound_str
      integer :: bound
      real(dp) :: c
      real(dp) :: d, plane_n(3), di, p, q(3), SF
      real(dp) :: t0, tau
      real(dp) :: egap
      real(dp) :: tol
      character(len=20) :: type_str
      logical :: print_summary

      if (parse_in_restraints) type_str = "restraint"
      if (parse_in_constraints) type_str = "constraint"

      if (name == 'restraints' .or. name == 'constraints') then ! new restraints stanze
	 if (name == 'restraints') then
	    parse_in_restraints = .true.
	    parse_in_constraints = .false.
	    call QUIP_FoX_get_value(attributes, "N", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read N in restraints stanza")
	    read (value, *) n
	    call print("allocating array for "//n//" restraints")
	    if (allocated(parse_ds%restraint)) deallocate(parse_ds%restraint)
	    allocate(parse_ds%restraint(n))
	 else ! constraints
	    parse_in_constraints = .true.
	    parse_in_restraints = .false.
	    call QUIP_FoX_get_value(attributes, "N", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read N in constraints stanza")
	    read (value, *) n
	    call print("allocating array for "//n//" constraints")
	    if (allocated(parse_ds%constraint)) deallocate(parse_ds%constraint)
	    allocate(parse_ds%constraint(n))
	 endif

      else if (parse_in_restraints .or. parse_in_constraints) then
	 call QUIP_FoX_get_value(attributes, "print_summary", value, status)
	 if (status == 0) then
	    read (value, *) print_summary
	 else
	    print_summary=.true.
	 endif

	 if (parse_in_restraints) then
	    call QUIP_FoX_get_value(attributes, "k", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read k in bond_length "//trim(type_str))
	    read (value, *) k
	    call QUIP_FoX_get_value(attributes, "bound", value, status)
	    if (status == 0) then
	       read (value, *) bound_str
	       call interpret_bound_string(trim(bound_str),bound)
	    else
	       bound = BOTH_UPPER_AND_LOWER_BOUNDS
	    endif
	 else
	    call QUIP_FoX_get_value(attributes, "tol", value, status) !constraint tolerance
	    if (status == 0) then
	       read (value, *) tol
	    else
	       tol=-1._dp !default
	    endif
	 endif

	 if (name == 'bond_length') then
	    call QUIP_FoX_get_value(attributes, "atom_1", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_1 in bond_length "// &
	      trim(type_str))
	    read (value, *) atom_1
	    call QUIP_FoX_get_value(attributes, "atom_2", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_2 in bond_length "// &
	      trim(type_str))
	    read (value, *) atom_2
	    call QUIP_FoX_get_value(attributes, "d", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read d in bond_length "//trim(type_str))
	    read (value, *) d

	    tau = -1.0_dp
	    di = -1.0_dp
	    call QUIP_FoX_get_value(attributes, "tau", value, status)
	    if (status == 0) then
	       read (value, *) tau
	    endif
	    call QUIP_FoX_get_value(attributes, "di", value, status)
	    if (status == 0) then
	       read (value, *) di
	    endif
	    call QUIP_FoX_get_value(attributes, "t0", value, status)
	    if (status == 0) then
	       read (value, *) t0
	    endif
	    if ((status == 0 .and. (tau <= 0.0_dp .or. di < 0.0_dp)) .or. &
	        (status /= 0 .and. (tau > 0.0_dp .or. di >= 0.0_dp))) then
	      call system_abort("got t0, but tau or di is missing or invalid")
	    endif

	    if (tau > 0.0_dp) then
	       if (parse_in_restraints) then
		  call constrain_bondlength(parse_ds, atom_1, atom_2, d, di=di, t0=t0, tau=tau, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength(parse_ds, atom_1, atom_2, d, di=di, t0=t0, tau=tau, tol=tol, print_summary=print_summary)
	       endif
	    else
	       if (parse_in_restraints) then
		  call constrain_bondlength(parse_ds, atom_1, atom_2, d, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength(parse_ds, atom_1, atom_2, d, tol=tol, print_summary=print_summary)
	       endif
	    endif

	 else if (name == 'bond_length_dev_pow') then
	    call QUIP_FoX_get_value(attributes, "atom_1", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_1 in bond_length_dev_pow "// &
	      trim(type_str))
	    read (value, *) atom_1
	    call QUIP_FoX_get_value(attributes, "atom_2", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_2 in bond_length_dev_pow "// &
	      trim(type_str))
	    read (value, *) atom_2
	    call QUIP_FoX_get_value(attributes, "p", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read p in bond_length_dev_pow "//trim(type_str))
	    read (value, *) p
	    call QUIP_FoX_get_value(attributes, "d", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read d in bond_length_dev_pow "//trim(type_str))
	    read (value, *) d

	    tau = -1.0_dp
	    di = -1.0_dp
	    call QUIP_FoX_get_value(attributes, "tau", value, status)
	    if (status == 0) then
	       read (value, *) tau
	    endif
	    call QUIP_FoX_get_value(attributes, "di", value, status)
	    if (status == 0) then
	       read (value, *) di
	    endif
	    call QUIP_FoX_get_value(attributes, "t0", value, status)
	    if (status == 0) then
	       read (value, *) t0
	    endif
	    if ((status == 0 .and. (tau <= 0.0_dp .or. di < 0.0_dp)) .or. &
	        (status /= 0 .and. (tau > 0.0_dp .or. di >= 0.0_dp))) then
	      call system_abort("got t0, but tau or di is missing or invalid")
	    endif

	    if (tau > 0.0_dp) then
	       if (parse_in_restraints) then
		  call constrain_bondlength_dev_pow(parse_ds, atom_1, atom_2, p, d, di=di, t0=t0, tau=tau, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength_dev_pow(parse_ds, atom_1, atom_2, p, d, di=di, t0=t0, tau=tau, tol=tol, print_summary=print_summary)
	       endif
	    else
	       if (parse_in_restraints) then
		  call constrain_bondlength_dev_pow(parse_ds, atom_1, atom_2, p, d, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength_dev_pow(parse_ds, atom_1, atom_2, p, d, tol=tol, print_summary=print_summary)
	       endif
	    endif

	 else if (name == 'bond_length_sq') then

	    call QUIP_FoX_get_value(attributes, "atom_1", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_1 in bond_length_sq "//trim(type_str))
	    read (value, *) atom_1
	    call QUIP_FoX_get_value(attributes, "atom_2", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_2 in bond_length_sq "//trim(type_str))
	    read (value, *) atom_2
	    call QUIP_FoX_get_value(attributes, "d", value, status)
	    if (status == 0) then
	       read (value, *) d
	       if (parse_in_restraints) then
		  call constrain_bondlength_sq(parse_ds, atom_1, atom_2, d, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength_sq(parse_ds, atom_1, atom_2, d, tol=tol, print_summary=print_summary)
	       endif
	    else
	       if (parse_in_restraints) then
		  call constrain_bondlength_sq(parse_ds, atom_1, atom_2, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength_sq(parse_ds, atom_1, atom_2, tol=tol, print_summary=print_summary)
	       endif
	    endif

	 else if (name == 'bond_angle_cos') then

	    call QUIP_FoX_get_value(attributes, "atom_1", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_1 in bond_angle_cos "//trim(type_str))
	    read (value, *) atom_1
	    call QUIP_FoX_get_value(attributes, "atom_2", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_2 in bond_angle_cos "//trim(type_str))
	    read (value, *) atom_2
	    call QUIP_FoX_get_value(attributes, "atom_3", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_3 in bond_angle_cos "//trim(type_str))
	    read (value, *) atom_3
	    call QUIP_FoX_get_value(attributes, "c", value, status)
	    if (status == 0) then
	       read (value, *) c
	       if (parse_in_restraints) then
		  call constrain_bondanglecos(parse_ds, atom_1, atom_2, atom_3, c, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondanglecos(parse_ds, atom_1, atom_2, atom_3, c, tol=tol, print_summary=print_summary)
	       endif
	    else
	       if (parse_in_restraints) then
		  call constrain_bondanglecos(parse_ds, atom_1, atom_2, atom_3, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondanglecos(parse_ds, atom_1, atom_2, atom_3, tol=tol, print_summary=print_summary)
	       endif
	    endif

	 else if (name == 'bond_length_diff') then

	    call QUIP_FoX_get_value(attributes, "atom_1", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_1 in bond_length_diff "//trim(type_str))
	    read (value, *) atom_1
	    call QUIP_FoX_get_value(attributes, "atom_2", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_2 in bond_length_diff "//trim(type_str))
	    read (value, *) atom_2
	    call QUIP_FoX_get_value(attributes, "atom_3", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom_3 in bond_length_diff "//trim(type_str))
	    read (value, *) atom_3
	    call QUIP_FoX_get_value(attributes, "d", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read d in bond_length_diff "//trim(type_str))
	    read (value, *) d

	    call QUIP_FoX_get_value(attributes, "di", value, di_status)
	    if (di_status == 0) read (value, *) di
	    call QUIP_FoX_get_value(attributes, "tau", value, tau_status)
	    if (tau_status == 0) read (value, *) tau
	    call QUIP_FoX_get_value(attributes, "t0", value, t0_status)
	    if (t0_status == 0) read (value, *) t0
	    if (count( (/ di_status, tau_status, t0_status /) == 0) /= 0 .and. &
	        count( (/ di_status, tau_status, t0_status /) == 0) /= 3) then
	       call system_abort("restraint_startElement_handler needs either all or none of di, tau, t0")
	    endif

	    if (di_status == 0) then ! relax target value with time
	       if (parse_in_restraints) then
		  call constrain_bondlength_diff(parse_ds, atom_1, atom_2, atom_3, d, di=di, tau=tau, t0=t0, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength_diff(parse_ds, atom_1, atom_2, atom_3, d, di=di, tau=tau, t0=t0, tol=tol, print_summary=print_summary)
	       endif
	    else ! fixed target value
	       if (parse_in_restraints) then
		  call constrain_bondlength_diff(parse_ds, atom_1, atom_2, atom_3, d, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_bondlength_diff(parse_ds, atom_1, atom_2, atom_3, d, tol=tol, print_summary=print_summary)
	       endif
	    endif

	 else if (name == 'gap_energy') then

	    call QUIP_FoX_get_value(attributes, "egap", value, status)
	    if (status == 0) then
	       read (value, *) egap
	       if (parse_in_restraints) then
		  call constrain_gap_energy(parse_ds, egap, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_gap_energy(parse_ds, egap, tol=tol, print_summary=print_summary)
	       endif
	    else
	       call system_abort("restraint_startElement_handler failed to read egap in gap_energy "//trim(type_str))
	       !if (parse_in_restraints) then
	       !  call constrain_gap_energy(parse_ds, restraint_k=k, print_summary=print_summary)
	       !else
	       !  call constrain_gap_energy(parse_ds, tol=tol, print_summary=print_summary)
	       !endif
	    endif

	 else if (name == 'atom_plane') then

	    call QUIP_FoX_get_value(attributes, "atom", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read atom in atom_plane "//trim(type_str))
	    read (value, *) atom_1
	    call QUIP_FoX_get_value(attributes, "normal", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read normal in atom_plane "//trim(type_str))
	    read (value, *) plane_n
	    call QUIP_FoX_get_value(attributes, "d", value, status)
	    if (status == 0) then
	       read (value, *) d
	       if (parse_in_restraints) then
		  call constrain_atom_plane(parse_ds, atom_1, plane_n, d, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_atom_plane(parse_ds, atom_1, plane_n, d, tol=tol, print_summary=print_summary)
	       endif
	    else
	       if (parse_in_restraints) then
		  call constrain_atom_plane(parse_ds, atom_1, plane_n, restraint_k=k, bound=bound, print_summary=print_summary)
	       else
		  call constrain_atom_plane(parse_ds, atom_1, plane_n, tol=tol, print_summary=print_summary)
	       endif
	    endif

	 else if (name == 'struct_factor_like_mag') then

	    call QUIP_FoX_get_value(attributes, "Z", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read Z in struct_factor_like_mag "//trim(type_str))
	    read (value, *) Z
	    call QUIP_FoX_get_value(attributes, "q", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read q in struct_factor_like_mag "//trim(type_str))
	    read (value, *) q
	    call QUIP_FoX_get_value(attributes, "SF", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read SF in struct_factor_like_mag "//trim(type_str))
	    read (value, *) SF
	    if (parse_in_restraints) then
	       call constrain_struct_factor_like_mag(parse_ds, Z, q, SF, restraint_k=k, bound=bound, print_summary=print_summary)
	    else
	       call constrain_struct_factor_like_mag(parse_ds, Z, q, SF, tol=tol, print_summary=print_summary)
	    endif

	 else if (name == 'struct_factor_like_r') then

	    call QUIP_FoX_get_value(attributes, "Z", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read Z in struct_factor_like_r "//trim(type_str))
	    read (value, *) Z
	    call QUIP_FoX_get_value(attributes, "q", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read q in struct_factor_like_r "//trim(type_str))
	    read (value, *) q
	    call QUIP_FoX_get_value(attributes, "SF", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read SF in struct_factor_like_r "//trim(type_str))
	    read (value, *) SF
	    if (parse_in_restraints) then
	       call constrain_struct_factor_like_r(parse_ds, Z, q, SF, restraint_k=k, bound=bound, print_summary=print_summary)
	    else
	       call constrain_struct_factor_like_r(parse_ds, Z, q, SF, tol=tol, print_summary=print_summary)
	    endif

	 else if (name == 'struct_factor_like_i') then

	    call QUIP_FoX_get_value(attributes, "Z", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read Z in struct_factor_like_i "//trim(type_str))
	    read (value, *) Z
	    call QUIP_FoX_get_value(attributes, "q", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read q in struct_factor_like_i "//trim(type_str))
	    read (value, *) q
	    call QUIP_FoX_get_value(attributes, "SF", value, status)
	    if (status /= 0) call system_abort("restraint_startElement_handler failed to read SF in struct_factor_like_i "//trim(type_str))
	    read (value, *) SF
	    if (parse_in_restraints) then
	       call constrain_struct_factor_like_i(parse_ds, Z, q, SF, restraint_k=k, bound=bound, print_summary=print_summary)
	    else
	       call constrain_struct_factor_like_i(parse_ds, Z, q, SF, tol=tol, print_summary=print_summary)
	    endif

	 else
	    call system_abort("unknown " // trim(type_str) //" type '"//name//"' in " // trim(type_str) //"s xml stanza")
	 endif
      endif

   end subroutine

   subroutine restraints_constraints_endElement_handler(URI, localname, name)
      character(len=*), intent(in)   :: URI
      character(len=*), intent(in)   :: localname
      character(len=*), intent(in)   :: name

      if (name == 'restraints') then ! end of restraints stanza
	 parse_in_restraints = .false.
      else if (name == 'constraints') then ! end of constraints stanza
	 parse_in_constraints = .false.
      endif

   end subroutine

   !% Initialise restraints and constraints for a dynamical system from an XML containing string.
   !%
   !% XML format for constraints
   !%
   !% \begin{verbatim}
   !% <constraints N="n">
   !%   <bond_length atom_1="i" atom_2="j" d="r" />
   !%   <bond_length_sq atom_1="i" atom_2="j" d="r" />
   !%   <bond_length_diff atom_1="i" atom_2="j" atom_3="k" d="r" />
   !%   <bond_angle_cos atom_1="i" atom_2="j" atom_3="k" c="r" />
   !% </constraints N="n">
   !% \end{verbatim}
   !%
   !% In all types, constraint value (c or d) is optional.  If not
   !% specified, defaults to value at initialization time, just like
   !% underlying constraint adding routines.
   !%
   !% Format for restraints is similar, except surrounding stanza name
   !% is `restraints', and the additional {\tt k} attribute
   !% specifying spring constant for each restraint
   !
   subroutine init_restraints_constraints(this, param_str)
      type(DynamicalSystem), intent(inout), target :: this
      character(len=*), intent(in) :: param_str

      type(xml_t) :: fxml

      if (len_trim(param_str) == 0) return

      parse_ds => this

      call open_xml_string(fxml, param_str)

      call parse(fxml, startElement_handler = restraints_constraints_startElement_handler, &
	 endElement_handler = restraints_constraints_endElement_handler)

      call close_xml_t(fxml)

   end subroutine init_restraints_constraints

   subroutine interpret_bound_string(bound_str,bound)
     character(len=*), intent(in) :: bound_str
     integer, intent(out) :: bound

     if ( trim(bound_str) == trim(BOUND_STRING(UPPER_BOUND)) ) then
        bound = UPPER_BOUND
     elseif ( trim(bound_str) == trim(BOUND_STRING(LOWER_BOUND)) ) then
        bound = LOWER_BOUND
     elseif ( trim(bound_str) == trim(BOUND_STRING(BOTH_UPPER_AND_LOWER_BOUNDS)) ) then
        bound = BOTH_UPPER_AND_LOWER_BOUNDS
     else
        call system_abort("bound must take a value of "//trim(BOUND_STRING(UPPER_BOUND))//", "// &
           trim(BOUND_STRING(LOWER_BOUND))//" and "//trim(BOUND_STRING(BOTH_UPPER_AND_LOWER_BOUNDS)))
     endif

   end subroutine interpret_bound_string

end module restraints_constraints_xml_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
*     ------------------------------------------------------------------
*       O U T L S J
*     ------------------------------------------------------------------
*
      SUBROUTINE OUTLSJ(NC,C,IPACK,JA,JB,IPTR,IPT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,NBUG6,NBUG7,IFULL
      COMMON /FOUT/NF,NG,NR,NL,NZ,NN,NV,NS,IFLAG,NIJ,IDIM
      COMMON/INFORM/IREAD,IWRITE,IOUT,ISC(4),JSC(3),IALL,ISCW
      POINTER(QNOC,NOCCSH(1)),(QNELCSH,NELCSH(8,1)),
     :       (QNOCORB,NOCORB(8,1)),(QJ1,J1QNRD(15,1))
      POINTER(QIAJCMP,IAJCMP(1)),(QLJCOMP,LJCOMP(1)),
     :       (QNJCOMP,NJCOMP(1)),(QIAJCLD,IAJCLD(1)),
     :       (QLJCLSD,LJCLSD(1))
      COMMON /NDIMS/ QNOC,QNELCSH,QNOCORB,QJ1,NCFG
      COMMON /NON30/ QIAJCMP,QNJCOMP,QLJCOMP,QIAJCLD,QLJCLSD,MAXORB
      PARAMETER (LSTACK=128)
      DIMENSION C(NC),IPACK(NC),JA(NC),JB(NC),
     :          IPTR(NC),IPT(NC)
      INTEGER NCOUNT(8),ISTACK(LSTACK),II(4),IEL(4)
      EQUIVALENCE (NCOUNT(1),NF),(II(1),I1),(II(2),I2),(II(3),I3),
     :            (II(4),I4)
      CHARACTER*1 INT(8)
      DATA INT/'F','G','R','L','Z','N','V','S'/
*
*     Formats for Integrals
*
   10 FORMAT(1X,A1,I2,'(',A3,',',A3,')',I6)
   11 FORMAT(1X,A1,I2,'(',2A3,',',2A3,')',I6)
   12 FORMAT(1X,A1,2X,'(',A3,',',A3,')',I6)
*
*     Format for Coefficients
*
   20 FORMAT(F14.8,A1,2I4,I3)
*
*     Format for Coefficient terminator
*
   30 FORMAT(14X,'*',16X,' ')
   31 FORMAT(1X,'*')
*
* --- TERMINATE  INTEGERAL LISTS and Rewind
*
      DO 1 J = 1,8
         ENDFILE(UNIT=ISC(J))
         REWIND(UNIT=ISC(J))
    1 CONTINUE
*
*===  Begin processing data
*
      NINT = 0
      DO 100 ICASE = 1,8
         N = NCOUNT(ICASE)
         IF (ICASE .NE. 3 .AND. ICASE .NE. 4) THEN
            DO 102 J = 1,N
              READ(ISC(ICASE)) C(J),IPACK(J),JA(J),JB(J)
  102       CONTINUE
         ELSE
            DO 104 J = 1,N
               READ(ISC(ICASE)) C(J),IPACK(J),JA(J),JB(J),IPTR(J)
  104       CONTINUE
         END IF
         CALL QSORT(N,IPACK,IPT,ISTACK,LSTACK,IERR)
         IF (IERR .EQ. 1) THEN
            WRITE(ISCW,*) ' Stack dimension not large enough for sort'
            CALL EXIT(1)
         END IF
*
*        Output the list of integrals with pointers to the data
*
         LAST = 0
  110    J = LAST +1
         LAST = J
         IF (J .LE. N) THEN
*
*          Unpack electron data
*
           K = IPACK(IPT(J))
           I4 = MOD(K,64)
           K = K/64
           IF (ICASE.LE.2 .OR. ICASE.EQ.4 .OR. ICASE.EQ.5) THEN
             I2 = MOD(K,64)
             K = K/64
           ELSE
              I3 = MOD(K,64)
              K = K/64
              I2 = MOD(K,64)
              K = K/64
              I1 = MOD(K,64)
              K = K/64
              IF (ICASE .GT. 5) K = K - 1
           END IF
*
*           Find  last item in the list with this integral
*
  120      LAST = LAST + 1
           IF (LAST .LE. N) THEN
             IF (IPACK(IPT(J)) .EQ. IPACK(IPT(LAST))) GO TO 120
           END IF
           LAST = LAST -1
           NINT = NINT + 1
           IF (ICASE .LE. 2) THEN
             WRITE(IOUT,10) INT(ICASE),K,IAJCMP(I2),IAJCMP(I4),LAST
           ELSE IF (ICASE .EQ. 4 .OR. ICASE .EQ. 5) THEN
             WRITE(IOUT,12) INT(ICASE),IAJCMP(I2),IAJCMP(I4),LAST
           ELSE
             DO 140 J = 1,4
               IF (II(J) .GE. 44 .AND. IALL .EQ. 0) THEN
                 IEL(J) = IAJCLD(64-II(J))
               ELSE
                 IEL(J) = IAJCMP(II(J))
               END IF
  140        CONTINUE
             WRITE(IOUT,11) INT(ICASE),K,(IEL(J),J=1,4),LAST
           END IF
           GO TO 110
        END IF
        WRITE(IOUT,31)
*
*             Write out the data for the integrals
*
        DO 150 J = 1,N
          K = IPT(J)
          WRITE(IOUT,20) C(K),INT(ICASE),JA(K),JB(K)
  150   CONTINUE
        WRITE(IOUT,30)
        IF (ICASE .EQ. 2) THEN
          WRITE(IOUT,31)
          WRITE(IOUT,31)
        END IF
  100 CONTINUE
      WRITE(ISCW,*) 'The total number of integrals =',NINT
      IDIM = NINT
      RETURN
      END

ext: f


text:
SUBROUTINE interp_bilinear(di,dj,var,velp)
    !== give 4 corner points of a square, interpolate point values inside
    !== di is the distance of the particle to the left face
    !== dj is the distance of the particle to the southern face

    IMPLICIT NONE
    REAL*8, INTENT(in) :: di,dj
    REAL*4, INTENT(in), DIMENSION( 2, 2 ) :: var
    REAL*8, INTENT(out) :: velp
    REAL*8 :: i1,i4

    ! calcuate the bilinear interpolation
    i1 = (var(2,1) - var(1,1))*di + var(1,1)
    i4 = (var(2,2) - var(1,2))*di + var(1,2)
    velp = (i4 - i1)*dj + i1

END SUBROUTINE interp_bilinear

max_stars_repo_name: jinbow/Octupus

max_issues_repo_name: jinbow/Octupus

max_forks_repo_name: jinbow/Octupus


text:
! Copyright (c) 2003, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!   Internal subprogram.
!   The main program contains only the call to the internal subprogram,
!   and there is no program name (i.e., 'program p').

  call f()

contains
  subroutine f()
    parameter(n=1)
    integer result(n), expect(n)
    data expect/2/
    b = 2
    result(1) = b
    call check( result, expect, n )
  end subroutine
end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
! { dg-do compile }
! PR fortran/29403
program p
   character(len=10) a, b, c
   integer i
   i = 1
   print ('(I0)'), i
   a = '(I0,'
   b = 'I2,'
   c = 'I4)'
   call prn(a, b, c, i)
   print (1,*), i       ! { dg-error "in PRINT statement" }
end program p

subroutine prn(a, b, c, i)
  integer i
  character(len=*) a, b, c
  print (a//(b//c)), i, i, i
  print trim(a//trim(b//c)), i, i, i
end subroutine prn

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
c
c ----------------------------------------------------------------------
c
      subroutine fsigld(fi,sc,gc,fmass,fwidth,smass,mNLSP,idecay , fsi)
c
c This subroutine computes an off-shell fermion wavefunction from a
c flowing-in external fermion and a scalar boson, for the NLSP-boson-
c Goldstino vertex. The h.c. of the NLSP decay is handled via the
c input parameter idecay.
c
c input:
c       complex fi(6)          : flow-in  fermion                   |fi>
c       complex sc(3)          : input    scalar                      s
c       complex gc(2)          : coupling constants                  gsf
c       real    fmass          : mass  of output fermion f'
c       real    fwidth         : width of output fermion f'
c       integer idecay         :  1 for NLSP decay to Goldstino
c                              : -1 for Goldstino to NLSP (h.c. of above)
c
c output:
c       complex fsi(6)         : off-shell fermion             |f',s,fi>
c
      implicit none
      double complex  fi(6), sc(3), gc(2), fsi(6), s1, s2, s3, s4, ds
      double complex  p14p, p14m, p23p, p23m
      double precision  pf(0:3), fmass, fwidth, mNLSP, smass, pf2
      integer idecay

      double complex ci, cZero
      parameter( ci = ( 0.0d0, 1.0d0 ), cZero = ( 0.0d0, 0.0d0 ) )
c
      fsi(5) = fi(5) - sc(2)
      fsi(6) = fi(6) - sc(3)

      pf(0) = dble( fsi(5))
      pf(1) = dble( fsi(6))
      pf(2) = dimag(fsi(6))
      pf(3) = dimag(fsi(5))
      pf2 = pf(0)**2 - pf(1)**2 - pf(2)**2 - pf(3)**2

      if ( idecay.ne.1 .or. idecay.ne.-1 ) then
         write(6,*) 'error in idecay of FSIGLD'
         stop
      end if

      p14p = dble(sc(2)) + dimag(sc(2))
      p14m = dble(sc(2)) - dimag(sc(2))
      p23p = dble(sc(3)) + dimag(sc(3))*ci
      p23m = dble(sc(3)) - dimag(sc(3))*ci

      ds = -sc(1)/dcmplx( pf2-fmass**2, fmass*fwidth )

      s1 = -idecay*gc(1)*fi(1)*smass**2
      s2 = -idecay*gc(1)*fi(2)*smass**2
      s3 = gc(1)*mNLSP*( fi(1)*p14p + fi(2)*p23m )
      s4 = gc(1)*mNLSP*( fi(1)*p23p + fi(2)*p14m )

      if ( gc(2).ne.cZero ) then
         s1 = s1 + gc(2)*mNLSP*( fi(3)*p14m - fi(4)*p23m )
         s2 = s2 + gc(2)*mNLSP*(-fi(3)*p23p + fi(4)*p14p )
         s3 = s3 - gc(2)*idecay*fi(3)*smass**2
         s4 = s4 - gc(2)*idecay*fi(4)*smass**2
      end if

      fsi(1) = ( (pf(0)-pf(3))*s3 - dconjg(fsi(6))*s4 + fmass*s1 )*ds
      fsi(2) = (       -fsi(6)*s3 +  (pf(0)+pf(3))*s4 + fmass*s2 )*ds
      fsi(3) = ( (pf(0)+pf(3))*s1 + dconjg(fsi(6))*s2 + fmass*s3 )*ds
      fsi(4) = (        fsi(6)*s1 +  (pf(0)-pf(3))*s2 + fmass*s4 )*ds
c
      return
      end


text:
module mGf90
implicit none

private
public srand_mtGaus
contains

real function srand_mtGaus(N, mean, cvm) result(y)
implicit none

!The dimension of the multivariate vector
integer, intent(in)   :: N

real,    intent(in)   :: mean(N)  !Mean values of the multivariate distribution

real,    intent(in)   :: cvm(N,N) !Covariance matrix of the multivariate distribution
dimension             :: y(N)
integer    :: nullty  ! i/o for cholesky subroutine
integer    :: error
integer    :: i, j
real       :: cvm_(N*(N+1)/2)
real       :: chol(N*(N+1)/2)

!Convert cvm to a vector with length N*(N+1)/2
do i = 1, N
  do j = 1,i
    cvm_(i*(i-1)/2+j) = cvm(i,j)
  enddo
enddo

call cholesky(cvm_,N,N*(N+1)/2,chol,nullty,error)

y = multiGauss(chol,mean,N)
end function srand_mtGaus

real function gasdev() result(gval)
!This function generates an array of independent Gaussian variates
implicit none
real :: FAC, R, V1, V2, X
real, save :: GSET
integer, save :: ISET = 0

IF (ISET.EQ.0) THEN
   R = 99
   do while( R .ge. 1.0d0 )
      call random_number(V1)
      V1= 2d0*V1 - 1d0

      call random_number(V2)
      V2= 2d0*V2 - 1d0

      R = V1**2 + V2**2
   end do
   FAC  = SQRT( -2.0d0*LOG(R)/R )
   GSET = V1*FAC
   gval = V2*FAC
   ISET = 1
ELSE
   gval=GSET
   ISET=0
ENDIF

RETURN
END function gasdev

subroutine cholesky ( a, n, nn, u, nullty, ifault )
!*****************************************************************************80
!
!! CHOLESKY computes the Cholesky factorization of a PDS matrix.
!
!  Discussion:
!
!    For a positive definite symmetric matrix A, the Cholesky factor U
!    is an upper triangular matrix such that A = U' * U.
!
!    This routine was originally named "CHOL", but that conflicted with
!    a built in MATLAB routine name.
!
!  Modified:
!
!    01 February 2008
!
!  Author:
!
!    Michael Healy
!    Modifications by AJ Miller.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Michael Healy,
!    Algorithm AS 6:
!    Triangular decomposition of a symmetric matrix,
!    Applied Statistics,
!    Volume 17, Number 2, 1968, pages 195-197.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A((N*(N+1))/2), a positive definite matrix
!    stored by rows in lower triangular form as a one dimensional array,
!    in the sequence
!    A(1,1),
!    A(2,1), A(2,2),
!    A(3,1), A(3,2), A(3,3), and so on.
!
!    Input, integer ( kind = 4 ) N, the order of A.
!
!    Input, integer NN, the dimension of the array used to store A,
!    which should be at least (N*(N+1))/2.
!
!    Output, real ( kind = 8 ) U((N*(N+1))/2), an upper triangular matrix,
!    stored by columns, which is the Cholesky factor of A.  The program is
!    written in such a way that A and U can share storage.
!    It is also the same as lower triangular matrix stored by rows (by B. Chen on 12 Aug 2019)
!
!    Output, integer ( kind = 4 ) NULLTY, the rank deficiency of A.  If NULLTY is zero,
!    the matrix is judged to have full rank.
!
!    Output, integer ( kind = 4 ) IFAULT, an error indicator.
!    0, no error was detected;
!    1, if N < 1;
!    2, if A is not positive semi-definite.
!    3, NN < (N*(N+1))/2.
!
!  Local Parameters:
!
!    Local, real ( kind = 8 ) ETA, should be set equal to the smallest positive
!    value such that 1.0 + ETA is calculated as being greater than 1.0 in the
!    accuracy being used.
!
  implicit none

  integer ( kind = 4 ), intent(in) :: n
  integer ( kind = 4 ), intent(in) :: nn

  real    ( kind = 8 ), intent(in) :: a(nn)
  real    ( kind = 8 ), parameter :: eta = 1.0D-19
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icol
  integer ( kind = 4 ), intent(out) :: ifault
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) irow
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kk
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ), intent(out) :: nullty
  real    ( kind = 8 ), intent(out) :: u(nn)
  real    ( kind = 8 ) w
  real    ( kind = 8 ) x

  ifault = 0
  nullty = 0

  if ( n <= 0 ) then
    ifault = 1
    return
  end if

  if ( nn < ( n * ( n + 1 ) ) / 2 ) then
    ifault = 3
    return
  end if

  j = 1
  k = 0
  ii = 0
!
!  Factorize column by column, ICOL = column number.
!
  do icol = 1, n

    ii = ii + icol
    x = eta * eta * a(ii)
    l = 0
    kk = 0
!
!  IROW = row number within column ICOL.
!
    do irow = 1, icol

      kk = kk + irow
      k = k + 1
      w = a(k)
      m = j

      do i = 1, irow - 1
        l = l + 1
        w = w - u(l) * u(m)
        m = m + 1
      end do

      l = l + 1

      if ( irow == icol ) then
        exit
      end if

      if ( u(l) /= 0.0D+00 ) then

        u(k) = w / u(l)

      else

        u(k) = 0.0D+00

        if ( abs ( x * a(k) ) < w * w ) then
          ifault = 2
          return
        end if

      end if

    end do
!
!  End of row, estimate relative accuracy of diagonal element.
!
    if ( abs ( w ) <= abs ( eta * a(k) ) ) then

      u(k) = 0.0D+00
      nullty = nullty + 1

    else

      if ( w < 0.0D+00 ) then
        ifault = 2
        return
      end if

      u(k) = sqrt ( w )

    end if

    j = j + icol

  end do

  return
end subroutine cholesky

!-----------------------------------------------------------------------------
subroutine subchl ( a, b, n, u, nullty, ifault, ndim, det )

!*****************************************************************************80
!
!! SUBCHL computes the Cholesky factorization of a (subset of a) PDS matrix.
!
!  Modified:
!
!    11 February 2008
!
!  Author:
!
!    FORTRAN77 version by Michael Healy, PR Freeman
!    FORTRAN90 version by  John Burkardt
!
!  Reference:
!
!    PR Freeman,
!    Remark AS R44:
!    A Remark on AS 6 and AS7: Triangular decomposition of a symmetric matrix
!    and Inversion of a positive semi-definite symmetric matrix,
!    Applied Statistics,
!    Volume 31, Number 3, 1982, pages 336-339.
!
!    Michael Healy,
!    Algorithm AS 6:
!    Triangular decomposition of a symmetric matrix,
!    Applied Statistics,
!    Volume 17, Number 2, 1968, pages 195-197.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A((M*(M+1))/2), a positive definite matrix
!    stored by rows in lower triangular form as a one dimensional array,
!    in the sequence
!    A(1,1),
!    A(2,1), A(2,2),
!    A(3,1), A(3,2), A(3,3), and so on.
!    In the simplest case, M, the order of A, is equal to N.
!
!    Input, integer ( kind = 4 ) B(N), indicates the order in which the
!    rows and columns of A are to be used.  In the simplest case,
!    B = (1,2,3...,N).
!
!    Input, integer ( kind = 4 ) N, the order of the matrix, that is,
!    the matrix formed by using B to select N rows and columns of A.
!
!    Output, real ( kind = 8 ) U((N*(N+1))/2), an upper triangular matrix,
!    stored by columns, which is the Cholesky factor of A.  The program is
!    written in such a way that A and U can share storage.
!
!    Output, integer ( kind = 4 ) NULLTY, the rank deficiency of A.
!    If NULLTY is zero, the matrix is judged to have full rank.
!
!    Output, integer ( kind = 4 ) IFAULT, an error indicator.
!    0, no error was detected;
!    1, if N < 1;
!    2, if A is not positive semi-definite.
!
!    Input, integer ( kind = 4 ) NDIM, the dimension of A and U, which might
!    be presumed to be (N*(N+1))/2.
!
!    Output, real ( kind = 8 ) DET, the determinant of the matrix.
!
  implicit none

  integer ( kind = 4 ) n
  integer ( kind = 4 ) ndim

  real    ( kind = 8 ) a(ndim)
  integer ( kind = 4 ) b(n)
  real    ( kind = 8 ) det
  real    ( kind = 8 ), parameter :: eta = 1.0D-09
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icol
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) ij
  integer ( kind = 4 ) irow
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kk
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) nullty
  real    ( kind = 8 ) u(ndim)
  real    ( kind = 8 ) w
  real    ( kind = 8 ) x

  ifault = 0
  nullty = 0
  det = 1.0D+00

  if ( n <= 0 ) then
    ifault = 1
    return
  end if

  ifault = 2
  j = 1
  k = 0

  do icol = 1, n

    ij = ( b(icol) * ( b(icol) - 1 ) ) / 2
    ii = ij + b(icol)
    x = eta * eta * a(ii)
    l = 0

    do irow = 1, icol

      kk = ( b(irow) * ( b(irow) + 1 ) ) / 2
      k = k + 1
      jj = ij + b(irow)
      w = a(jj)
      m = j

      do i = 1, irow - 1
        l = l + 1
        w = w - u(l) * u(m)
        m = m + 1
      end do

      l = l + 1

      if ( irow == icol ) then
        exit
      end if

      if ( u(l) /= 0.0D+00 ) then

        u(k) = w / u(l)

      else

        if ( abs ( x * a(kk) ) < w * w ) then
          ifault = 2
          return
        end if

        u(k) = 0.0D+00

      end if

    end do

    if ( abs ( eta * a(kk) ) <= abs ( w ) ) then

      if ( w < 0.0D+00 ) then
        ifault = 2
        return
      end if

      u(k) = sqrt ( w )

    else

      u(k) = 0.0D+00
      nullty = nullty + 1

    end if

    j = j + icol
    det = det * u(k) * u(k)

  end do

  return
end subroutine subchl

real function multigauss(R,mu,n) result(G)
!!$ Generates a one dimensional array of length n containing multivariate Gaussian pseudo-random numbers
!!$ where R is the lower-triangular Cholesky factor of the covariance matrix of the desired distribution
!!$ and mu is a one dimensional array of length n, containing the mean value for each random variable.
!!$ R must be a one dimensional array, containing the lower-triangular matrix stored by row, starting from the
!!$ uppermost, leftmost entry (first row, first column).
   implicit none
   integer, intent(in) :: n
   real, intent(in) :: R(n*(n+1)/2)
   real, intent(in) :: mu(n)
   real :: Nu(n)
   dimension :: G(n)
   integer :: i, j

   if( n*(n+1)/2 .ne. size(R) ) then
      write(6,*) ' n*(n+1)/2 != size(R) in multiGauss '
      write(6,*) ' Cholesky factor matrix has size ', size(R)
      write(6,*) ' but you are requesting a vector of size ',n
      stop
   else

!!$ generate an array of independent Gaussian variates
      do j = 1, n
         Nu(j) = gasdev()
      end do

!!$ start with the desired mean, and add the product of
!!$ [the lower-triangular Cholesky factor of the covariance matrix] x [ the vector of iid Gaussian variates]
      do i = 1, n
         G(i) = mu(i)
         do j =  1,i
            G(i) = G(i) + R( i*(i-1)/2 + j ) * Nu(j)
         end do
      end do

   end if
   return
 end function multigauss
end module mGf90

max_stars_repo_name: BingzhangChen/IBM

max_issues_repo_name: BingzhangChen/IBM

max_forks_repo_name: BingzhangChen/IBM


text:
!************************************************************
!
!  This example shows how to iterate over group members using
!  H5Literate.
!
!  This file is intended for use with HDF5 Library version 1.8
!  with --enable-fortran2003
!
!
!************************************************************
MODULE g_iterate

  USE HDF5
  USE ISO_C_BINDING
  IMPLICIT NONE

CONTAINS

!************************************************************
!
!  Operator function.  Prints the name and type of the object
!  being examined.
!
! ************************************************************

  INTEGER FUNCTION op_func(loc_id, name, info, operator_data) bind(C)

    USE HDF5
    USE ISO_C_BINDING
    IMPLICIT NONE

    INTEGER(HID_T), VALUE :: loc_id
    CHARACTER(LEN=1), DIMENSION(1:10) :: name ! must have LEN=1 for bind(C) strings
    TYPE(C_PTR) :: info
    TYPE(C_PTR) :: operator_data

    INTEGER   :: status, i, len

    TYPE(H5O_info_t), TARGET :: infobuf
    TYPE(C_PTR) :: ptr
    CHARACTER(LEN=10) :: name_string

    !
    ! Get type of the object and display its name and type.
    ! The name of the object is passed to this FUNCTION by
    ! the Library.
    !

    DO i = 1, 10
       name_string(i:i) = name(i)(1:1)
    ENDDO

    CALL H5Oget_info_by_name_f(loc_id, name_string, infobuf, status)

    ! Include the string up to the C NULL CHARACTER
    len = 0
    DO
       IF(name_string(len+1:len+1).EQ.C_NULL_CHAR.OR.len.GE.10) EXIT
       len = len + 1
    ENDDO

    IF(infobuf%type.EQ.H5O_TYPE_GROUP_F)THEN
       WRITE(*,*) "Group: ", name_string(1:len)
    ELSE IF(infobuf%type.EQ.H5O_TYPE_DATASET_F)THEN
       WRITE(*,*) "Dataset: ", name_string(1:len)
    ELSE IF(infobuf%type.EQ.H5O_TYPE_NAMED_DATATYPE_F)THEN
       WRITE(*,*) "Datatype: ", name_string(1:len)
    ELSE
       WRITE(*,*) "Unknown: ", name_string(1:len)
    ENDIF

    op_func = 0 ! return successful

  END FUNCTION op_func

END MODULE g_iterate

PROGRAM main

  USE HDF5
  USE ISO_C_BINDING
  USE g_iterate

  IMPLICIT NONE

  CHARACTER(LEN=17), PARAMETER :: filename  = "h5ex_g_iterate.h5"
  INTEGER(HID_T) :: file ! Handle
  INTEGER :: status
  TYPE(C_FUNPTR) :: funptr
  TYPE(C_PTR) :: ptr
  INTEGER(hsize_t) :: idx
  INTEGER :: ret_value

  !
  ! Open file.
  !
  CALL H5Fopen_f(filename, H5F_ACC_RDONLY_F, file, status)

  !
  ! Begin iteration.
  !
  WRITE(*,'(A)') "Objects in root group:"

  idx = 0
  funptr = C_FUNLOC(op_func) ! call back function
  ptr    = C_NULL_PTR

  CALL H5Literate_f(file, H5_INDEX_NAME_F, H5_ITER_NATIVE_F, idx, funptr, ptr, ret_value, status)

  !
  ! Close and release resources.
  !
  CALL H5Fclose_f(file, status)

END PROGRAM main

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
*
*
*     --------------------------------------------------------------
*     R L S P 3
*     --------------------------------------------------------------
*                                                                  *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vanderbilt University,  Nashville             October 1996   *
*
      SUBROUTINE RLSP3(K,JA1,JA2,JA3,K1,K2,K3,K4,KA,REC)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/MEDEFN/IHSH,NJ(16),LJ(16),NOSH1(16),NOSH2(16),J1QN1(31,3),
     :     J1QN2(31,3),IJFUL(16)
      REC=ZERO
*
      IF((JA3.GT.JA1).AND.(JA3.GT.JA2)) THEN
        IF(JA1-JA2.LT.0) THEN
* 1 2 3
          CALL RLSP31(K,JA1,JA2,JA3,K1,K2,K3,K4,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          CALL RLSP32(K,JA1,JA2,JA3,K1,K2,K3,K4,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          CALL RLSP31(K,JA1,JA2,JA3,K1,K2,K3,K4,KA,1,IAT,REC)
          CALL RLSP32(K,JA1,JA2,JA3,K1,K2,K3,K4,KA,1,IAT,R1)
	  REC=REC*R1
        ELSEIF(JA1-JA2.GT.0) THEN
* 2 1 3
          CALL RLSP31(K,JA2,JA1,JA3,K2,K1,K3,K4,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          CALL RLSP32(K,JA2,JA1,JA3,K2,K1,K3,K4,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          CALL RLSP31(K,JA2,JA1,JA3,K2,K1,K3,K4,KA,1,IAT,REC)
          CALL RLSP32(K,JA2,JA1,JA3,K2,K1,K3,K4,KA,1,IAT,R1)
	  REC=REC*R1
          IF(MOD(K1+K2-K3,4).NE.0)REC=-REC
	ELSE
	  STOP
	ENDIF
      ELSEIF((JA3.LT.JA1).AND.(JA3.LT.JA2)) THEN
        IF(JA1-JA2.LT.0) THEN
* 3 1 2
          IPL=ITREXG(K1,K4,K2,KA,IKKL)+1
          IF(IKKL.EQ.0) RETURN
          IGL=IPL+IKKL-1
          CALL RLSP31(K,JA3,JA1,JA2,K4,K1,J12,K2,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          DO 3 I=IPL,IGL,2
            J12=I-1
            CALL DLSA6(K,K4,K3,KA,K1,K2,J12,0,IAT,R)
            IF(IAT.NE.0) THEN
              CALL RLSP32(K,JA3,JA1,JA2,K4,K1,J12,K2,KA,0,IAT,R)
              IF(IAT.NE.0) THEN
                CALL RLSP32(K,JA3,JA1,JA2,K4,K1,J12,K2,KA,1,IAT,R1)
                CALL DLSA6(K,K4,K3,KA,K1,K2,J12,1,IAT,R2)
                IF(MOD(2*K1+K2+K4-J12-K3,4).NE.0)R1=-R1
                REC=REC+R1*R2
              ENDIF
            ENDIF
    3     CONTINUE
	  IF(DABS(REC).LT.EPS)RETURN
          CALL RLSP31(K,JA3,JA1,JA2,K4,K1,J12,K2,KA,1,IAT,R1)
          REC=REC*R1
        ELSEIF(JA1-JA2.GT.0) THEN
* 3 2 1
          IPL=ITREXG(K2,K4,K1,KA,IKKL)+1
          IF(IKKL.EQ.0) RETURN
          IGL=IPL+IKKL-1
          CALL RLSP31(K,JA3,JA2,JA1,K4,K2,J12,K1,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          DO 4 I=IPL,IGL,2
            J12=I-1
            CALL DLSA6(K,K4,K3,KA,K2,K1,J12,0,IAT,R)
            IF(IAT.NE.0) THEN
              CALL RLSP32(K,JA3,JA2,JA1,K4,K2,J12,K1,KA,0,IAT,R)
              IF(IAT.NE.0) THEN
                CALL RLSP32(K,JA3,JA2,JA1,K4,K2,J12,K1,KA,1,IAT,R1)
                CALL DLSA6(K,K4,K3,KA,K2,K1,J12,1,IAT,R2)
                IF(MOD(3*K2+2*K1+K4-J12-2*K3,4).NE.0)R1=-R1
                REC=REC+R1*R2
              ENDIF
            ENDIF
    4     CONTINUE
	  IF(DABS(REC).LT.EPS)RETURN
          CALL RLSP31(K,JA3,JA2,JA1,K4,K2,J12,K1,KA,1,IAT,R1)
          REC=REC*R1
	ELSE
	  STOP
	ENDIF
      ELSE
        IF(JA1-JA2.LT.0) THEN
* 1 3 2
          IPL=ITREXG(K1,K4,K2,KA,IKKL)+1
          IF(IKKL.EQ.0) RETURN
          IGL=IPL+IKKL-1
          CALL RLSP31(K,JA1,JA3,JA2,K1,K4,J12,K2,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          DO 5 I=IPL,IGL,2
            J12=I-1
            CALL DLSA6(K,K4,K3,KA,K1,K2,J12,0,IAT,R)
            IF(IAT.NE.0) THEN
              CALL RLSP32(K,JA1,JA3,JA2,K1,K4,J12,K2,KA,0,IAT,R)
              IF(IAT.NE.0) THEN
                CALL RLSP32(K,JA1,JA3,JA2,K1,K4,J12,K2,KA,1,IAT,R1)
                CALL DLSA6(K,K4,K3,KA,K1,K2,J12,1,IAT,R2)
                REC=REC+R1*R2
              ENDIF
            ENDIF
    5     CONTINUE
	  IF(DABS(REC).LT.EPS)RETURN
          CALL RLSP31(K,JA1,JA3,JA2,K1,K4,J12,K2,KA,1,IAT,R1)
          REC=REC*R1
          IF(MOD(K1+K2-K3,4).NE.0)REC=-REC
        ELSEIF(JA1-JA2.GT.0) THEN
* 2 3 1
          IPL=ITREXG(K2,K4,K1,KA,IKKL)+1
          IF(IKKL.EQ.0) RETURN
          IGL=IPL+IKKL-1
          CALL RLSP31(K,JA2,JA3,JA1,K2,K4,J12,K1,KA,0,IAT,R)
          IF(IAT.EQ.0) RETURN
          DO 6 I=IPL,IGL,2
            J12=I-1
            CALL DLSA6(K,K4,K3,KA,K2,K1,J12,0,IAT,R)
            IF(IAT.NE.0) THEN
              CALL RLSP32(K,JA2,JA3,JA1,K2,K4,J12,K1,KA,0,IAT,R)
              IF(IAT.NE.0) THEN
                CALL RLSP32(K,JA2,JA3,JA1,K2,K4,J12,K1,KA,1,IAT,R1)
                CALL DLSA6(K,K4,K3,KA,K2,K1,J12,1,IAT,R2)
                REC=REC+R1*R2
              ENDIF
            ENDIF
    6     CONTINUE
	  IF(DABS(REC).LT.EPS)RETURN
          CALL RLSP31(K,JA2,JA3,JA1,K2,K4,J12,K1,KA,1,IAT,R1)
          REC=REC*R1
	ELSE
	  STOP
	ENDIF
      ENDIF
      REC=REC/DSQRT(DBLE(J1QN1(JA1,K)*J1QN1(JA2,K)*J1QN1(JA3,K)))
      RETURN
      END


text:
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      SUBROUTINE CFFTMB (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV,
     1                  WORK, LENWRK, IER)
      INTEGER  LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER
      COMPLEX       C(LENC)
      REAL     WSAVE(LENSAV)     ,WORK(LENWRK)
      LOGICAL XERCON
C
      IER = 0
C
      IF (LENC .LT. (LOT-1)*JUMP + INC*(N-1) + 1) THEN
        IER = 1
        CALL XERFFT ('CFFTMB ', 6)
      ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
        IER = 2
        CALL XERFFT ('CFFTMB ', 8)
      ELSEIF (LENWRK .LT. 2*LOT*N) THEN
        IER = 3
        CALL XERFFT ('CFFTMB ', 10)
      ELSEIF (.NOT. XERCON(INC,JUMP,N,LOT)) THEN
        IER = 4
        CALL XERFFT ('CFFTMB ', -1)
      ENDIF
C
      IF (N .EQ. 1) RETURN
C
      IW1 = N+N+1
      CALL CMFM1B (LOT,JUMP,N,INC,C,WORK,WSAVE,WSAVE(IW1),
     1                           WSAVE(IW1+1))
      RETURN
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: thkrz/fstpack

max_forks_repo_name: thkrz/fstpack


text:
! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
!
! --- End License

!
!////////////////////////////////////////////////////////////////////////
!
!      FTLinkedListUnitTests.f90
!      Created: January 15, 2013 5:25 PM
!      By: David Kopriva
!
!      Create a linked list of FTValue objects. The class can store
!      any objects inheriting from FTObject. For instance, one could
!      have a linked list whose first record is a linked list and second
!      is a value.
!
!      The FTLinkedListIterator class, used for stepping through a linked
!      list, is also and tested demonstrated here.
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE FTLinkedListClassTests
         IMPLICIT NONE
!
!        ------------
!        Basic tests
!        ------------
!
         CALL basicTests
!
!        --------------------------------------
!        Do more with deleting objects in lists
!        --------------------------------------
!
         CALL TestDeletingObjects
!
!        ---------------------------------
!        Now test appending lists to lists
!        ---------------------------------
!
         CALL testAppendingLists

      END SUBROUTINE FTLinkedListClassTests
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE basicTests
         USE FTAssertions
         USE FTValueClass
         USE FTLinkedListClass
         USE FTLinkedListIteratorClass
         IMPLICIT NONE
!
!        ------------------------------------------------------
!        r1,r2,r3 are three object pointers to be stored in the
!        linked list
!        ------------------------------------------------------
!
         TYPE(FTValue), POINTER :: r1, r2, r3
!
!        ------------------------------------------------------------
!        We need a pointer of the base class type to pass to the list
!        ------------------------------------------------------------
!
         CLASS(FTObject), POINTER :: objectPtr
!
!        -------------------------------------------------
!        Here we will define the list as a pointer, but it
!        can be a non-pointer, too, like the iterator.
!        -------------------------------------------------
!
         CLASS (FTLinkedList)        , POINTER :: list
         TYPE (FTLinkedListIterator), POINTER :: iterator

         INTEGER                      :: i
         REAL                         :: singleTol = 2*EPSILON(1.0e0)
!
!        -------------------------------------------------------
!        Allocate and initialize the linked list. Remember that
!        init initializes the object with reference count of one
!        and ownership to this calling subroutine.
!        -------------------------------------------------------
!
         ALLOCATE(list)
         CALL list % init()
!
!        --------------------------------------------------------------------------
!        Test Reference counting.
!           Retaining an object implies that this subroutine wants to
!           share ownership of an object. (In fact, it already does, so the retain
!           is redundant.)
!           Releasing an object implies that the caller releases its share of
!           the object. Note that, because of the init call, we still own the
!           object after the retain+release. Finally, remember to balance init/retains
!           and releases: init + # retains = # releases
!        --------------------------------------------------------------------------
!
         CALL FTAssertEqual(1,list % refCount(),"Reference counting: Initial object refCount")
         CALL list % retain()
         CALL FTAssertEqual(2,list % refCount(), "Reference counting: Test retain")
         CALL releaseFTLinkedList(list)
         CALL FTAssertEqual(1,list % refCount(),"Reference counting: test release")
!
!        ---------------------------------------------------------------------------------
!        A linked list that is just initialized has no items in it, so its COUNT is zero.
!        ---------------------------------------------------------------------------------
!
         CALL FTAssertEqual(0,list % COUNT(),"Initial list size")
!
!        ------------------------------------------------------------------
!        Add some items to the linked list.
!        As we add items the list's count increases.
!        Also, the list takes a share of ownership of the object added
!        so the reference count of the object is increased. Once we
!        add an object to the list, we may not need it anymore, in which case
!        we reliquish our share of the ownership by releasing it. We should
!        check to see if the object needs to be deallocated for safety.
!
!        The first item that we add is an FTValueObject that represents an
!        integer.
!        ------------------------------------------------------------------
!
         ALLOCATE(r1)
         CALL r1 % initWithValue(1)
         objectPtr => r1
         CALL list % add(objectPtr)
         CALL FTAssertEqual(1,list % COUNT(),"List size after adding one object")
!
         CALL FTAssertEqual(2,r1 % refCount(),&
         "Reference counting: Stored object should have reference count increased")
         CALL releaseFTValue(r1)

         CALL FTAssertEqual(1,objectPtr % refCount(),&
         "Reference counting: Stored object should have reference count decreased")
!
!        ---------------------------------------------------------------------
!        The second item in the list is an FTValueObject that stores a string.
!        We will use r2 later in the subroutine, so we don't give up ownership
!        at this time.
!        ---------------------------------------------------------------------
!
         ALLOCATE(r2)
         CALL r2 % initWithValue("r2 is a string")
         objectPtr => r2
         CALL list % add(objectPtr)
         CALL FTAssertEqual(2,list % COUNT(),"List size after adding second object")
!
!        -------------------------------------
!        The final object we add stores a real
!        -------------------------------------
!
         ALLOCATE(r3)
         CALL r3 % initWithValue(3.14)
         objectPtr => r3
         CALL list % add(objectPtr)
         CALL FTAssertEqual(3,list % COUNT(),"List size after adding third object")
         CALL releaseFTValue(r3)
!
!        ---------------------------------------------------------------------------------
!        Check integrity of stored objects. We iterate
!        through the linked list with an iterator. The
!        iterator is initialized with the list, which
!        is why we have defined the list as a pointer.
!        Note that the init function gives us ownership
!        of the iterator, and the iterator takes an ownership
!        stake in the list. When we get to the end of this subroutine,
!        since we are the sole owner of the iterator (we could check with  % shouldDestruct),
!        we must destruct the iterator. Since the iterator is the sole owner of
!        the list after we release it, it will destruct the list.
!        ---------------------------------------------------------------------------------
!
         ALLOCATE(iterator)
         CALL iterator % initWithFTLinkedList(list)
         CALL FTAssertEqual(2,list % refCount(),"Ref count increase on addition of list to iterator")
!
!        -------------------------------------------
!        Iterate through the list from the beginning
!        -------------------------------------------
!
         CALL iterator % setToStart() !(This is automatically done in the init)
         i = 1

         DO WHILE (.NOT.iterator % isAtEnd() )
!
!           --------------------------------------------------------------------------
!           Get a pointer to the iterator's current object. Just pointing does not
!           imply ownership, so if we wanted to save the reference to the object after
!           destructing the iterator and list, we'd have to retain it.
!           --------------------------------------------------------------------------
!
            objectPtr => iterator % object()

            SELECT TYPE(v=>objectPtr)
               TYPE IS (FTValue)
                 SELECT CASE(i)
                     CASE(1)
                        CALL FTAssertEqual(1,v % integerValue(),"First item is integer value")
                     CASE(2)
                        CALL FTAssertEqual("r2 is a string",v % stringValue(14),"Second item is string value")
                     CASE(3)
                        CALL FTAssertEqual(3.14,v % realValue(),singleTol,"Third item in list is real value")
                  END SELECT
               CLASS DEFAULT
                  CALL FTAssert(.false.,"Unknown type stored in linked list")
            END SELECT
            CALL iterator % moveToNext()
            i = i + 1
         END DO
!
!        ---------------------------------------------------
!        Delete the object r2 from the list
!        The number of items should be decreased by one, and
!        what were r1 and r3 should still be there.
!        ---------------------------------------------------
!
         objectPtr => r2
         CALL list % remove(objectPtr)
         CALL FTAssertEqual(2,list % COUNT(),"List has two objects after removing one")
         CALL FTAssertEqual(1,r2 % refCount(),msg = "Refcount after removing object")

         CALL iterator % setToStart()
         i = 1

         DO WHILE (.NOT.iterator % isAtEnd())
            objectPtr => iterator % object()
            SELECT TYPE(v=>objectPtr)
               TYPE IS (FTValue)
                 SELECT CASE(i)
                     CASE(1)
                        CALL FTAssertEqual(1,v % integerValue(),"First item in list doesn't have proper value")
                     CASE(3)
                        CALL FTAssertEqual(3.14,v % realValue(),singleTol,"third item in list doesn't have proper value")
                  END SELECT
               CLASS DEFAULT
                  CALL FTAssert(.false.,"Known type stored in linked list")
            END SELECT
            CALL iterator % moveToNext()
            i = i + 1
         END DO
!
!        ------------------------------------------------------
!        By releasing an object we promise not to reference it.
!        Otherwise, it is possible to get an undefined pointer.
!        ------------------------------------------------------
!
         CALL releaseFTLinkedList(list)
         CALL FTAssertEqual(1,list % refCount(),"Ref count decrease on release")
!
!        -------------------------------------------------------------------
!        Normally we would now check if the list should be deallocated. But
!        since we know that it's refCount = 1, we won't and simply nullify
!        the pointer. Note that the list itself is retained by the iterator.
!        -------------------------------------------------------------------
!
         list => NULL()
!
!        ------------------------------------------------------------------------------
!        Clean up iterator. Its refCount, since we know we are the only owners, should
!        cause it to be destructed. On release, the iterator will deallocate the list
!        since it is the last owner.
!        ------------------------------------------------------------------------------
!
         CALL releaseFTLinkedListIterator(iterator)

      END SUBROUTINE basicTests
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE testAppendingLists
         USE FTAssertions
         USE FTValueClass
         USE FTLinkedListClass
         USE FTLinkedListIteratorClass
         IMPLICIT NONE
!
         TYPE (FTValue)             , POINTER :: v
         CLASS(FTObject)            , POINTER :: objectPtr
         CLASS(FTLinkedList)        , POINTER :: list1, list2
         TYPE (FTMutableObjectArray), POINTER :: array

         TYPE (FTLinkedListIterator), POINTER :: iterator
         INTEGER                              :: j, N

         ALLOCATE(iterator)
!
!        ----------------------------------------------
!        Create the two lists that will be concatenated
!        ----------------------------------------------
!
         ALLOCATE(list1,list2)
         CALL list1 % init()
         CALL list2 % init()
!
!        -----------------------------
!        Add some objects to the lists
!        -----------------------------
!
         DO j = 1, 5
            ALLOCATE(v)
            CALL v % initWithValue(j)
            objectPtr => v
            CALL list1 % add(objectPtr)
            CALL releaseFTValue(v)
         END DO

         DO j = 6, 10
            ALLOCATE(v)
            CALL v % initWithValue(j)
            objectPtr => v
            CALL list2 % add(objectPtr)
            CALL releaseFTValue(v)
         END DO
!
!        -----------------------------------
!        Add the elements of list2 to list 1
!        -----------------------------------
!
         CALL list1 % addObjectsFromList(list2)
         CALL FTAssertEqual(10, list1 % COUNT(),"Append list increases list size")
!
!        -------------------------------------------
!        See that the new list contains the old one.
!        Note that objects are owned by both lists.
!        -------------------------------------------
!
         CALL iterator % initWithFTLinkedList(list1)
         j = 1
         DO WHILE (.NOT.iterator % isAtEnd())
            v => valueFromObject(iterator % object())
            CALL FTAssertEqual(j, v % integerValue(),"Item value stored properly")
            IF ( j >= 6 )     THEN
               objectPtr => iterator % object()
               CALL FTAssertEqual(2, objectPtr % refCount(), "Records owned by two lists")
            END IF
            CALL iterator % moveToNext()
            j = j + 1
         END DO
!
!        --------------------------------------------------
!        Now delete the list2. Its contents should still be
!        in list1
!        --------------------------------------------------
!
         CALL releaseFTLinkedList(list2)
         CALL FTAssertEqual(.TRUE., .NOT. ASSOCIATED(list2),"List has only one owner and should deallocate on release")
!
!        --------------------------------------------------
!        List1 should have its contents plus the other list
!        with refCount 1
!        --------------------------------------------------
!
         CALL iterator % setToStart()
         j = 1
         DO WHILE (.NOT.iterator % isAtEnd())

            v         => valueFromObject(iterator % object())
            CALL FTAssertEqual(j, v % integerValue(),"Item value stored properly afer release of added list")
            CALL FTAssertEqual(1, v % refCount(),"Item value pointed to by list refCount")

            objectPtr => iterator % object()
            CALL FTAssertEqual(1, objectPtr % refCount(), "Objects owned by one list")
            CALL iterator % moveToNext()
            j = j + 1
         END DO
!
!        -------------------------------------
!        Create an object array from the list.
!        -------------------------------------
!
         array => list1 % allObjects()
         DO j = 1, array % COUNT()
            v => valueFromObject(array % objectAtIndex(j))
            CALL FTAssertEqual(j, v % integerValue(),"Item value stored in array created from list")
            objectPtr => array % objectAtIndex(j)
            CALL FTAssertEqual(2, objectPtr % refCount(), "Objects owned by one list and one array")
         END DO
         CALL releaseFTMutableObjectArray(array)
         CALL FTAssert(test = .NOT. ASSOCIATED(array),msg = "Array unreferenced")
!
!        -------------------------------------------
!        Now reverse the list and iterate through it
!        -------------------------------------------
!
         CALL list1 % reverse()
         N = list1 % COUNT()
         j = N
         CALL iterator % setToStart
         DO WHILE (.NOT.iterator % isAtEnd())

            v => valueFromObject(iterator % object())
            CALL FTAssertEqual(j, v % integerValue(),"Item value stored properly afer release of added list")
            CALL FTAssertEqual(1, v % refCount(),"Item value pointed to by list refCount")

            objectPtr => iterator % object()
            CALL FTAssertEqual(1, objectPtr % refCount(), "Objects owned by one list")
            CALL iterator % moveToNext()
            j = j - 1
         END DO
!
!        --------
!        Clean up
!        --------
!
         CALL releaseFTLinkedList(list1)
         CALL releaseFTLinkedListIterator(iterator)

      END SUBROUTINE testAppendingLists
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE TestDeletingObjects
         USE FTLinkedListClass
         USE FTLinkedListIteratorClass
         USE FTValueClass
         USE FTAssertions
         IMPLICIT NONE
!
!        -----------------------------------------------
!        Add a number of objects to a list and then test
!        out removing them
!        -----------------------------------------------
!
         TYPE (FTValue)           , POINTER   :: v
         CLASS(FTObject)          , POINTER   :: obj
         CLASS(FTLinkedList)      , POINTER   :: list
         CLASS(FTLinkedListRecord), POINTER   :: recordPtr

         TYPE(FTLinkedListIterator), POINTER :: iterator
         INTEGER                             :: j
!
!        ----------------------------------------------
!        Create the two lists that will be concatenated
!        ----------------------------------------------
!
         ALLOCATE(list)
         CALL list % init()
!
!        -----------------------------
!        Add some objects to the lists
!        -----------------------------
!
         DO j = 1, 6
            ALLOCATE(v)
            CALL v % initWithValue(j)
            obj => v
            CALL list % add(obj)
            CALL releaseFTValue(v)
         END DO
!
!        ------------------------------------------------------------
!        Create an iterator on the list and remove an "interior" item
!        ------------------------------------------------------------
!
         ALLOCATE(iterator)
         CALL iterator % initwithFTLinkedList(list)
!
!        ---------------
!        Delete the tail
!        ---------------
!
         CALL iterator % setToStart()
         j = 1
         DO WHILE( .NOT.iterator % isAtEnd() )
            IF ( j == 6 )     THEN
               v => valueFromObject(iterator % object())
               CALL FTAssertEqual(6,v % integerValue(),"Value of object to be deleted")
               CALL iterator % removeCurrentRecord()
            END IF
            CALL iterator % moveToNext()
            j = j + 1
         END DO

         CALL FTAssertEqual(5,list % COUNT(),"Count after deletion of tail object")
!!
!!        -----------------------
!!        Delete the third record
!!        -----------------------
!!
         CALL iterator % setToStart()
         j = 1
         DO WHILE( .NOT.iterator % isAtEnd() )
            IF ( j == 3 )     THEN
               v   => valueFromObject(iterator % object())
               CALL FTAssertEqual(3,v % integerValue(),"Value of object to be deleted")
               CALL iterator % removeCurrentRecord()
               EXIT
            END IF
            CALL iterator % moveToNext()
            j = j + 1
         END DO

         CALL FTAssertEqual(4,list % COUNT(),"Count after deletion of middle object")
         v   => valueFromObject(iterator % object())
         CALL FTAssertEqual(4,v % integerValue(),"Value of current object after deleting")
!
!        ---------------------------------
!        Make sure connections are correct
!        ---------------------------------
!
         recordPtr => iterator % currentRecord()
         v         => valueFromObject(recordPtr % previous % recordObject)
         CALL FTAssertEqual(2,v % integerValue(), "Value of previous object after deleting")
         v         => valueFromObject(recordPtr % next % recordObject)
         CALL FTAssertEqual(5,v % integerValue(), "Value of next object after deleting")
!
!        -----------------------
!        Delete the first record
!        -----------------------
!
         CALL iterator % setToStart()
         CALL iterator % removeCurrentRecord()
         CALL FTAssertEqual(3, list % COUNT(), "count after deleting head")
         v   => valueFromObject(iterator % object())
         CALL FTAssertEqual(2,v % integerValue(),"Value of current head after deleting")
!
!        --------
!        Clean up
!        --------
!
         CALL releaseFTLinkedList(list)
         CALL releaseFTLinkedListIterator(iterator)

      END SUBROUTINE TestDeletingObjects

max_stars_repo_name: HigherOrderMethods/FTObjectLibrary

max_issues_repo_name: HigherOrderMethods/FTObjectLibrary

max_forks_repo_name: HigherOrderMethods/FTObjectLibrary


text:
! /******************* C O N F O R T *******************
!     (c) 2017 Dominik Gronkiewicz <gronki@gmail.com>
!     Distributed under MIT License.
!
!     Example file:
!     --------------------------
!     key1 value1
!
!     # comment
!     key2 value2 value3   #  another comment
!
!     key3 "very long # text
!     with newline"
!
!     key4   5.0
!     key5   6.0  7.5  # hoorayy it's the end
! ******************************************************/

module confort

  use iso_c_binding
  implicit none

  type, bind(C) :: config
    type(c_ptr) :: buffer
    integer(c_size_t) ::  buffer_sz
    integer(c_size_t) :: n_records
    integer(c_size_t) :: records_sz
    type(c_ptr) :: records
  end type

  integer, parameter :: MINCF_OK =                0
  integer, parameter :: MINCF_ERROR =             1
  integer, parameter :: MINCF_ARGUMENT_ERROR =    ishft(1,1)
  integer, parameter :: MINCF_MEMORY_ERROR =      ishft(1,2)
  integer, parameter :: MINCF_FILE_NOT_FOUND =    ishft(1,3)
  integer, parameter :: MINCF_SYNTAX_ERROR =      ishft(1,4)
  integer, parameter :: MINCF_NOT_FOUND =         ishft(1,5)

  ! read the config

  interface mincf_read
    module procedure :: mincf_read_stdin
    module procedure :: mincf_read_file
  end interface

  ! get entry

  interface mincf_get
    module procedure :: mincf_get_or_error
    module procedure :: mincf_get_default
  end interface

  ! Other procedures

  interface
    subroutine mincf_free(cfg) &
          & bind(C,name='mincf_free')
      use iso_c_binding
      import :: config
      type(config), intent(inout) :: cfg
    end subroutine
  end interface

contains

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

  subroutine mincf_read_stdin(cfg,errno)

    type(config), intent(inout) :: cfg
    integer, intent(inout), optional :: errno
    integer :: errno_local

    interface
      subroutine c_mincf_read_stdin(cfg,errno) &
            & bind(C, name='fort_mincf_read_stdin')
        use iso_c_binding
        import :: config
        type(config), intent(inout) :: cfg
        integer(c_int), intent(out) :: errno
      end subroutine
    end interface

    if ( present(errno) ) then
      call c_mincf_read_stdin(cfg,errno)
    else
      call c_mincf_read_stdin(cfg,errno_local)
      if ( errno_local .ne. 0 ) then
        call mincf_free(cfg)
        error stop "confort: fatal error while reading the configuration"
      end if
    end if
  end subroutine

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

  subroutine mincf_read_file(cfg,fn,errno)

    type(config), intent(inout) :: cfg
    character(len=*), intent(in) :: fn
    integer, intent(inout), optional :: errno
    integer :: errno_local

    interface
      subroutine c_mincf_read_file(cfg,fn,sz,errno) &
            & bind(C, name='fort_mincf_read_file')
        use iso_c_binding
        import :: config
        type(config), intent(inout) :: cfg
        integer(c_size_t), intent(in), value :: sz
        character(c_char), intent(in) :: fn(sz)
        integer(c_int), intent(out) :: errno
      end subroutine
    end interface

    if ( present(errno) ) then
      call c_mincf_read_file(cfg, fn, len(fn,c_size_t), errno)
    else
      call c_mincf_read_file(cfg, fn, len(fn,c_size_t), errno_local)
      if ( errno_local .ne. 0 ) then
        call mincf_free(cfg)
        error stop "confort: fatal error while reading the configuration"
      end if
    end if

  end subroutine

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

  subroutine mincf_get_or_error(cfg,key,buf,errno)

    interface
      subroutine c_mincf_get(cfg,key,key_sz,buf,sz,errno) &
            & bind(C, name='fort_mincf_get')
        use iso_c_binding
        import config
        type(config), intent(out) :: cfg
        integer(c_size_t), intent(in), value :: key_sz
        character(c_char), intent(in) :: key(key_sz)
        integer(c_size_t), intent(in), value :: sz
        character(c_char), intent(inout) :: buf(sz)
        integer(c_int), intent(out) :: errno
      end subroutine
    end interface

    type(config), intent(out) :: cfg
    character(len=*), intent(in) :: key
    character(len=*), intent(out) :: buf
    integer, intent(inout), optional :: errno
    integer :: errno_local

    if ( present(errno) ) then
      call c_mincf_get(cfg, &
          & key, len(key,c_size_t), &
          & buf, len(buf,c_size_t), &
          & errno)
    else
      call c_mincf_get(cfg, &
          & key, len(key,c_size_t), &
          & buf, len(buf,c_size_t), &
          & errno_local)
      if ( errno_local .ne. 0 ) then
        call mincf_free(cfg)
        error stop "confort: fatal error while reading the configuration"
      end if
    end if

  end subroutine

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

  logical function mincf_exists(cfg,key)
    type(config), intent(out) :: cfg
    character(len=*), intent(in) :: key
    integer :: errno

    interface
      subroutine c_mincf_get_exists(cfg,key,key_sz,errno) &
            & bind(C, name='fort_mincf_get_exists')
        use iso_c_binding
        import config
        type(config), intent(out) :: cfg
        integer(c_size_t), intent(in), value :: key_sz
        character(c_char), intent(in) :: key(key_sz)
        integer(c_int), intent(out) :: errno
      end subroutine
    end interface

    call c_mincf_get_exists(cfg, key, len(key,c_size_t), errno)

    mincf_exists = ( errno .eq. MINCF_OK )
  end function

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

  subroutine mincf_get_default(cfg,key,buf,defvalue,errno)
    type(config), intent(out) :: cfg
    character(len=*), intent(in) :: key, defvalue
    character(len=*), intent(out) :: buf
    integer, intent(inout), optional :: errno
    integer :: errno_local

    interface
      subroutine c_mincf_get_default(cfg,key,key_sz,buf,sz,defvalue,defvalue_sz,errno) &
            & bind(C, name='fort_mincf_get_default')
        use iso_c_binding
        import config
        type(config), intent(out) :: cfg
        integer(c_size_t), intent(in), value :: key_sz, defvalue_sz
        character(c_char), intent(in) :: key(key_sz), defvalue(defvalue_sz)
        integer(c_size_t), intent(in), value :: sz
        character(c_char), intent(inout) :: buf(sz)
        integer(c_int), intent(out) :: errno
      end subroutine
    end interface

    call c_mincf_get_default(cfg, &
        & key, len(key,c_size_t), &
        & buf, len(buf,c_size_t), &
        & defvalue, len(defvalue,c_size_t), &
        & errno_local)

    if ( present(errno) ) then
      errno = errno_local
    end if
  end subroutine

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

  subroutine mincf_print_error(errno,file,line)
    integer, intent(in) :: errno
    character(len=*), intent(in), optional :: file
    integer, intent(in), optional :: line

    character(len=128) :: prefix = "confort"

    if ( present(file) .and. present(line) ) then
      write (prefix,"(A,' line ',I0)") file,line
    end if

    if ( iand(errno,MINCF_ERROR) .ne. 0 ) &
    & write (0,"(A,': ',A)") trim(prefix),"There was an error."
    if ( iand(errno,MINCF_NOT_FOUND) .ne. 0 ) &
    & write (0,"(A,': ',A)") trim(prefix),"The entry was not found."
    if ( iand(errno,MINCF_FILE_NOT_FOUND) .ne. 0 ) &
    & write (0,"(A,': ',A)") trim(prefix),"File not found."
    if ( iand(errno,MINCF_SYNTAX_ERROR) .ne. 0 ) &
    & write (0,"(A,': ',A)") trim(prefix),"Syntax error in configuration file."
    if ( iand(errno,MINCF_ARGUMENT_ERROR) .ne. 0 ) &
    & write (0,"(A,': ',A)") trim(prefix),"Incorrect argument."
    if ( iand(errno,MINCF_MEMORY_ERROR) .ne. 0 ) &
    & write (0,"(A,': ',A)") trim(prefix),"Memory error"
  end subroutine

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

end module

max_stars_repo_name: gronki/libconfort

max_issues_repo_name: gronki/libconfort

max_forks_repo_name: gronki/libconfort


text:
C$Procedure  RDENCI  ( Read encoded integers from text file )

      SUBROUTINE RDENCI ( UNIT, N, DATA )

C$ Abstract
C
C     Read N encoded integers from a text file, decoding them into
C     their equivalent integers.
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     None.
C
C$ Keywords
C
C     CONVERSION
C     NUMBERS
C     UTILITY
C
C$ Declarations

      INTEGER               UNIT
      INTEGER               N
      INTEGER               DATA(*)

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C      UNIT      I    Fortran unit number of input text file.
C      N         I    Number of integers to read and decode.
C      DATA      I    List of decoded integers.
C
C$ Detailed_Input
C
C     UNIT     The Fortran unit number for a previously opened text
C              file. All reading will begin at the CURRENT POSITION
C              in the text file.
C
C     N        The number of encoded integers to be read from the
C              text file attached to UNIT.
C
C$ Detailed_Output
C
C     DATA     List of decoded integers read from the text file
C              attached to UNIT.
C
C$ Parameters
C
C      None.
C
C$ Exceptions
C
C     1)   If N, the number of data items, is not positive, the error
C          SPICE(INVALIDARGUMENT) will be signalled.
C
C     2)   If an error occurs while reading from the text file
C          attached to  UNIT, the error SPICE(FILEREADFAILED) will
C          be signalled.
C
C     3)   If an error occurs while decoding a number, the error
C          SPICE(DECODINGERROR) will be signalled.
C
C$ Files
C
C     See the description of UNIT in Detailed_Input.
C
C$ Particulars
C
C     This routine will read N encoded integers beginning at the
C     current position in a previously opened text file. The current
C     position in a file is defined to be the text line immediately
C     following the last text line that was written or read. The
C     integers will be decoded and placed into a list of integers
C     which will be passed back to the caller. The encoded integers
C     are represented as quoted character strings so that a Fortran
C     list directed read may be used to read the encoded values,
C     rather than a formatted read with the format specifier
C     FMT = '(A)'.
C
C     This routine is one of a pair of routines which are used to
C     encode and decode integers:
C
C           WRENCI -- Encode and write integers to a file.
C           RDENCI -- Read and decode integers from a file.
C
C     The encoding/decoding of integers is performed to provide a
C     portable means for transferring data values.
C
C     Currently the encoded integers are represented as signed
C     hexadecimal numbers See INT2HX.FOR and HX2INT.FOR for details.
C
C$ Examples
C
C     Suppose we have the following input file which contains the values
C     1 - 100 encoded, and that the input file has already been opened
C     for reading. The arrow, '-->', indicates the current position in
C     the file.
C
C        -->'1' '2' '3' '4' '5' '6' '7' '8' '9' 'A' 'B' 'C' 'D' 'E'
C           'F' '10' '11' '12' '13' '14' '15' '16' '17' '18' '19'
C           '1A' '1B' '1C' '1D' '1E' '1F' '20' '21' '22' '23' '24'
C           '25' '26' '27' '28' '29' '2A' '2B' '2C' '2D' '2E' '2F'
C           '30' '31' '32' '33' '34' '35' '36' '37' '38' '39' '3A'
C           '3B' '3C' '3D' '3E' '3F' '40'
C           '41' '42' '43' '44' '45' '46' '47' '48' '49' '4A' '4B'
C           '4C' '4D' '4E' '4F' '50' '51' '52' '53' '54' '55' '56'
C           '57' '58' '59' '5A' '5B' '5C' '5D' '5E' '5F' '60' '61'
C           '62' '63' '64'
C
C     Then, the following code fragment would read and decode these
C     100 values.
C
C           N = 100
C           CALL RDENCI( UNIT, N, DATA )
C
C      Upon returning, the array data would contain the values 1 - 100.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     K.R. Gehringer (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 19-OCT-1992 (KRG)
C
C-&

C$ Index_Entries
C
C      read and decode encoded integers from a text file
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN
C
C     Local parameters
C
      INTEGER               ERRLEN
      PARAMETER           ( ERRLEN = 80 )

      INTEGER               MAXENC
      PARAMETER           ( MAXENC = 64 )

      INTEGER               WRKSIZ
      PARAMETER           ( WRKSIZ = 64 )
C
C     Local variables
C
      CHARACTER*(ERRLEN)    ERRMSG
      CHARACTER*(MAXENC)    WORK(WRKSIZ)

      INTEGER               I
      INTEGER               IOSTAT
      INTEGER               ITMBEG
      INTEGER               NITMS

      LOGICAL               ERROR

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'RDENCI' )
      END IF
C
C     Check to see if the number of data items is less than or equal
C     to zero. If it is, signal an error.
C
      IF ( N .LT. 1 ) THEN

         CALL SETMSG ( 'The number of data items to be read was' //
     .                 ' not positive: #.'                        )
         CALL ERRINT ( '#', N                                     )
         CALL SIGERR ( 'SPICE(INVALIDARGUMENT)'                   )
         CALL CHKOUT ( 'RDENCI'                                   )
         RETURN

      END IF
C
C     Initialize the beginning location to place the decoded data
C     items.
C
      ITMBEG = 1
C
C     We read in the encoded numbers in blocks of size WRKSIZ, and if
C     there was not a read error we will attempt to decode the numbers.
C     We signal an error if either:
C
C                (1) there is a read error
C                (2) there is an error decoding a number.
C
C     NOTE: EOF is interpreted as a read error because we know a priori
C           exactly how many data items we need to read: N.
C
C     Begin decoding the encoded data items read from the input file
C     in blocks of size NITMS. Each time the number of data items
C     NITMS is reached, decode the encoded numbers into the data array.
C
      DO WHILE ( ITMBEG .LE. N )
C
C        The number of items is either the size of the workspace, or
C        the number of data items which remain to be processed, which
C        should always be less than or equal to the size of the
C        workspace.
C
         NITMS = MIN ( WRKSIZ, N - ITMBEG + 1 )
C
C        Read in a block of data items to be decoded.
C
         READ (UNIT,*,IOSTAT=IOSTAT) ( WORK(I), I = 1, NITMS )
C
C        Check to see if we got a read error: IOSTAT .NE. 0. If we did,
C        then signal an error. EOF is considered to be a read error,
C        since we know exactly how many data items we expect to read.
C
         IF ( IOSTAT .NE. 0 ) THEN

            CALL SETMSG ( 'Error reading from logical unit #,' //
     .                    ' IOSTAT = #.'                        )
            CALL ERRINT ( '#', UNIT                             )
            CALL ERRINT ( '#', IOSTAT                           )
            CALL SIGERR ( 'SPICE(FILEREADFAILED)'               )
            CALL CHKOUT ( 'RDENCI'                              )
            RETURN

         END IF
C
C        Begin to decode the data items into the data array. Signal an
C        error if we cannot decode a data item.
C
         DO I = 1, NITMS

            CALL HX2INT ( WORK(I), DATA(ITMBEG + I - 1), ERROR, ERRMSG )
            IF ( ERROR ) THEN

               CALL SETMSG ( 'Decoding error occurred while'       //
     .                       ' attempting to decode item #: #. #'   )
               CALL ERRINT ( '#', I                                 )
               CALL ERRCH  ( '#', WORK(I)                           )
               CALL ERRCH  ( '#', ERRMSG                            )
               CALL SIGERR ( 'SPICE(DECODINGERROR)'                 )
               CALL CHKOUT ( 'RDENCI'                               )
               RETURN

            END IF

         END DO
C
C        Position the data item pointer at the next location to begin
C        placing the decoded items in the array DATA, and continue
C        processing the until done.
C
         ITMBEG = ITMBEG + NITMS

      END DO

      CALL CHKOUT ( 'RDENCI' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
 SUBROUTINE POLATES6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, &
                     NO,RLAT,RLON,IBO,LO,GO,IRET)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! $Revision: 71314 $
!
! SUBPROGRAM:  POLATES6   INTERPOLATE SCALAR FIELDS (BUDGET)
!   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM PERFORMS BUDGET INTERPOLATION
!           FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
!           IT REQUIRES A GRID FOR THE OUTPUT FIELDS (KGDSO(1)>=0).
!           THE ALGORITHM SIMPLY COMPUTES (WEIGHTED) AVERAGES
!           OF NEIGHBOR POINTS ARRANGED IN A SQUARE BOX
!           CENTERED AROUND EACH OUTPUT GRID POINT AND STRETCHING
!           NEARLY HALFWAY TO EACH OF THE NEIGHBORING GRID POINTS.
!           OPTIONS ALLOW CHOICES OF NUMBER OF POINTS IN EACH RADIUS
!           FROM THE CENTER POINT (IPOPT(1)) WHICH DEFAULTS TO 2
!           (IF IPOPT(1)=-1) MEANING THAT 25 POINTS WILL BE AVERAGED;
!           FURTHER OPTIONS ARE THE RESPECTIVE WEIGHTS FOR THE RADIUS
!           POINTS STARTING AT THE CENTER POINT (IPOPT(2:2+IPOPT(1))
!           WHICH DEFAULTS TO ALL 1 (IF IPOPT(1)=-1 OR IPOPT(2)=-1).
!           ANOTHER OPTION IS THE MINIMUM PERCENTAGE FOR MASK,
!           I.E. PERCENT VALID INPUT DATA REQUIRED TO MAKE OUTPUT DATA,
!           (IPOPT(3+IPOPT(1)) WHICH DEFAULTS TO 50 (IF -1).
!           ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
!           THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
!           (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
!           THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
!             (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
!             (KGDS(1)=001) MERCATOR CYLINDRICAL
!             (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
!             (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
!             (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
!             (KGDS(1)=203) ROTATED EQUIDISTANT CYLINDRICAL (E-STAGGER)
!             (KGDS(1)=205) ROTATED EQUIDISTANT CYLINDRICAL (B-STAGGER)
!           WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
!           AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
!           AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
!           INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
!           OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
!           EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
!           THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
!
! PROGRAM HISTORY LOG:
!   96-04-10  IREDELL
!   96-10-04  IREDELL  NEIGHBOR POINTS NOT BILINEAR INTERPOLATION
! 1999-04-08  IREDELL  SPLIT IJKGDS INTO TWO PIECES
! 2001-06-18  IREDELL  INCLUDE MINIMUM MASK PERCENTAGE OPTION
! 2015-01-27  GAYNO    REPLACE CALLS TO GDSWIZ WITH NEW MERGED
!                      VERSION OF GDSWZD.
!
! USAGE:    CALL POLATES6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
!    &                    NO,RLAT,RLON,IBO,LO,GO,IRET)
!
!   INPUT ARGUMENT LIST:
!     IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
!                IPOPT(1) IS NUMBER OF RADIUS POINTS
!                (DEFAULTS TO 2 IF IPOPT(1)=-1);
!                IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
!                (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
!                IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
!                (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
!     KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
!     KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
!     MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
!                OR DIMENSION OF INPUT GRID FIELDS IF KM=1
!     MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
!                OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
!     KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
!     IBI      - INTEGER (KM) INPUT BITMAP FLAGS
!     LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
!     GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
!
!   OUTPUT ARGUMENT LIST:
!     NO       - INTEGER NUMBER OF OUTPUT POINTS
!     RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES
!     RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES
!     IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
!     LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
!     GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
!     IRET     - INTEGER RETURN CODE
!                0    SUCCESSFUL INTERPOLATION
!                2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
!                3    UNRECOGNIZED OUTPUT GRID
!                31   INVALID UNDEFINED OUTPUT GRID
!                32   INVALID BUDGET METHOD PARAMETERS
!
! SUBPROGRAMS CALLED:
!   GDSWZD       GRID DESCRIPTION SECTION WIZARD
!   IJKGDS0      SET UP PARAMETERS FOR IJKGDS1
!   (IJKGDS1)    RETURN FIELD POSITION FOR A GIVEN GRID POINT
!   POLFIXS      MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 USE GDSWZD_MOD
!
 IMPLICIT NONE
!
 INTEGER,    INTENT(IN   ) :: IBI(KM), IPOPT(20), KM, MI, MO
 INTEGER,    INTENT(IN   ) :: KGDSI(200), KGDSO(200)
 INTEGER,    INTENT(  OUT) :: IBO(KM), IRET, NO
!
 LOGICAL*1,  INTENT(IN   ) :: LI(MI,KM)
 LOGICAL*1,  INTENT(  OUT) :: LO(MO,KM)
!
 REAL,       INTENT(IN   ) :: GI(MI,KM)
 REAL,       INTENT(  OUT) :: GO(MO,KM), RLAT(MO), RLON(MO)
!
 REAL,       PARAMETER     :: FILL=-9999.
!
 INTEGER                   :: IB, I1, IJKGDS1, IJKGDSA(20)
 INTEGER                   :: JB, J1, K, LB, LSW, MP, N
 INTEGER                   :: N11(MO), NB, NB1, NB2, NB3, NB4, NV
!
 REAL                      :: PMP,RLOB(MO),RLAB(MO)
 REAL                      :: WB, WO(MO,KM), XI, YI
 REAL                      :: XPTB(MO),YPTB(MO),XPTS(MO),YPTS(MO)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
 IRET=0
 IF(KGDSO(1).GE.0) THEN
   CALL GDSWZD(KGDSO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
   IF(NO.EQ.0) IRET=3
 ELSE
   IRET=31
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  SET PARAMETERS
 NB1=IPOPT(1)
 IF(NB1.EQ.-1) NB1=2
 IF(IRET.EQ.0.AND.NB1.LT.0) IRET=32
 LSW=1
 IF(IPOPT(1).EQ.-1.OR.IPOPT(2).EQ.-1) LSW=0
 IF(IRET.EQ.0.AND.LSW.EQ.1.AND.NB1.GT.15) IRET=32
 MP=IPOPT(3+IPOPT(1))
 IF(MP.EQ.-1.OR.MP.EQ.0) MP=50
 IF(MP.LT.0.OR.MP.GT.100) IRET=32
 PMP=MP*0.01
 IF(IRET.EQ.0) THEN
   NB2=2*NB1+1
   NB3=NB2*NB2
   NB4=NB3
   IF(LSW.EQ.1) THEN
     NB4=IPOPT(2)
     DO IB=1,NB1
       NB4=NB4+8*IB*IPOPT(2+IB)
     ENDDO
   ENDIF
 ELSE
   NB2=0
   NB3=0
   NB4=0
 ENDIF
 DO K=1,KM
   DO N=1,NO
     GO(N,K)=0.
     WO(N,K)=0.
   ENDDO
 ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
 CALL IJKGDS0(KGDSI,IJKGDSA)
 DO NB=1,NB3
!  LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS
   JB=(NB-1)/NB2-NB1
   IB=NB-(JB+NB1)*NB2-NB1-1
   LB=MAX(ABS(IB),ABS(JB))
   WB=1
   IF(LSW.EQ.1) WB=IPOPT(2+LB)
   IF(WB.NE.0) THEN
     DO N=1,NO
       XPTB(N)=XPTS(N)+IB/REAL(NB2)
       YPTB(N)=YPTS(N)+JB/REAL(NB2)
     ENDDO
     CALL GDSWZD(KGDSO, 1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV)
     CALL GDSWZD(KGDSI,-1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV)
     IF(IRET.EQ.0.AND.NV.EQ.0.AND.LB.EQ.0) IRET=2
     DO N=1,NO
       XI=XPTB(N)
       YI=YPTB(N)
       IF(XI.NE.FILL.AND.YI.NE.FILL) THEN
         I1=NINT(XI)
         J1=NINT(YI)
         N11(N)=IJKGDS1(I1,J1,IJKGDSA)
       ELSE
         N11(N)=0
       ENDIF
     ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  INTERPOLATE WITH OR WITHOUT BITMAPS
     DO K=1,KM
       DO N=1,NO
         IF(N11(N).GT.0) THEN
           IF(IBI(K).EQ.0.OR.LI(N11(N),K)) THEN
             GO(N,K)=GO(N,K)+WB*GI(N11(N),K)
             WO(N,K)=WO(N,K)+WB
           ENDIF
         ENDIF
       ENDDO
     ENDDO
   ENDIF
 ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  COMPUTE OUTPUT BITMAPS AND FIELDS
 DO K=1,KM
   IBO(K)=IBI(K)
   DO N=1,NO
     LO(N,K)=WO(N,K).GE.PMP*NB4
     IF(LO(N,K)) THEN
       GO(N,K)=GO(N,K)/WO(N,K)
     ELSE
       IBO(K)=1
       GO(N,K)=0.
     ENDIF
   ENDDO
 ENDDO
 IF(KGDSO(1).EQ.0) CALL POLFIXS(NO,MO,KM,RLAT,RLON,IBO,LO,GO)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE POLATES6

max_stars_repo_name: AerodyneLabs/Stratocast

max_issues_repo_name: AerodyneLabs/Stratocast

max_forks_repo_name: AerodyneLabs/Stratocast


program test
  use iso_c_binding
  interface
     subroutine foo(val) bind(C)
       use iso_c_binding
       integer, VALUE :: val
     end subroutine foo
  end interface
  integer :: i
  i = 5
  call foo (i)
end program test


text: #redirect Yoshi Akaha

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search