text:
      double complex function amp_2gam2g_mppppm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      -       +         +         +
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6

c--- Equation (2.9) with r=2, m=2 and labels permuted appropriately
c--- negative helicity gluon i=1
      amp_2gam2g_mppppm=
     & za(j5,j1)*za(j6,j1)**3/(za(j5,j6)*za(j6,j1)*za(j1,j2)*za(j2,j5))
     & *za(j6,j5)/(za(j6,j3)*za(j3,j5))
     & *za(j6,j5)/(za(j6,j4)*za(j4,j5))

      return
      end

      double complex function amp_2gam2g_pmpppm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      +       -         +         +
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6

c--- Equation (2.9) with r=2, m=2 and labels permuted appropriately
c--- negative helicity gluon i=2
      amp_2gam2g_pmpppm=
     & za(j5,j2)*za(j6,j2)**3/(za(j5,j6)*za(j6,j1)*za(j1,j2)*za(j2,j5))
     & *za(j6,j5)/(za(j6,j3)*za(j3,j5))
     & *za(j6,j5)/(za(j6,j4)*za(j4,j5))

      return
      end

      double complex function amp_2gam2g_ppmppm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      +       +         -         +
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6

c--- Equation (2.9) with r=2, m=2 and labels permuted appropriately
c--- negative helicity photon i=3
      amp_2gam2g_ppmppm=
     & za(j5,j3)*za(j6,j3)**3/(za(j5,j6)*za(j6,j1)*za(j1,j2)*za(j2,j5))
     & *za(j6,j5)/(za(j6,j3)*za(j3,j5))
     & *za(j6,j5)/(za(j6,j4)*za(j4,j5))

      return
      end

      double complex function amp_2gam2g_pppmpm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      +       +         +         -
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6

c--- Equation (2.9) with r=2, m=2 and labels permuted appropriately
c--- negative helicity photon i=4
      amp_2gam2g_pppmpm=
     & za(j5,j4)*za(j6,j4)**3/(za(j5,j6)*za(j6,j1)*za(j1,j2)*za(j2,j5))
     & *za(j6,j5)/(za(j6,j3)*za(j3,j5))
     & *za(j6,j5)/(za(j6,j4)*za(j4,j5))

      return
      end

      double complex function amp_2gam2g_mmpppm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      -       -         +         +
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6
      double precision t

c--- Equation (4.1)
c--- note: I have interpreted, e.g. <5|1+2|6> as [5|1+2|6> in order
c--- to get the correct <| and |] required by the helicities
      amp_2gam2g_mmpppm=
     &+(zb(j5,j1)*za(j1,j6)+zb(j5,j2)*za(j2,j6))*za(j5,j6)*t(j1,j2,j6)
     & /(zb(j5,j2)*zb(j1,j6)*zb(j2,j1)*za(j6,j3)
     &  *za(j6,j4)*za(j5,j3)*za(j5,j4))
     &-za(j1,j6)*zb(j4,j5)*(zb(j3,j1)*za(j1,j2)+zb(j3,j6)*za(j6,j2))
     & /(zb(j1,j6)*za(j4,j5)*zb(j5,j2)*za(j6,j3)*t(j1,j3,j6))
     &-za(j1,j6)*zb(j3,j5)*(zb(j4,j1)*za(j1,j2)+zb(j4,j6)*za(j6,j2))
     & /(zb(j1,j6)*za(j3,j5)*zb(j5,j2)*za(j6,j4)*t(j1,j4,j6))
      return
      end

      double complex function amp_2gam2g_pmmppm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      +       -         -         +
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6
      double precision t

c--- Equation (4.2)
c--- note: I have interpreted, e.g. <4|3+6|2> as [4|3+6|2> in order
c--- to get the correct <| and |] required by the helicities
      amp_2gam2g_pmmppm=
     &+zb(j4,j5)*za(j2,j6)**2*(zb(j1,j2)*za(j2,j3)+zb(j1,j6)*za(j6,j3))
     & /(za(j4,j5)*zb(j5,j3)*za(j1,j6)*s(j1,j2)*t(j1,j2,j6))
     &+za(j3,j6)*zb(j5,j4)*(zb(j1,j3)*za(j3,j2)+zb(j1,j6)*za(j6,j2))
     & /(zb(j3,j6)*za(j5,j4)*zb(j5,j2)*za(j1,j6)*t(j1,j3,j6))
     &+za(j3,j6)*zb(j1,j5)**2*(zb(j4,j3)*za(j3,j2)+zb(j4,j6)*za(j6,j2))
     & /(zb(j3,j6)*zb(j2,j5)*za(j6,j4)*s(j1,j2)*t(j1,j2,j5))
     &+zb(j5,j1)*za(j6,j2)*(zb(j5,j1)*za(j1,j6)+zb(j5,j4)*za(j4,j6))
     & *t(j1,j4,j5)/(zb(j3,j6)*zb(j5,j2)*zb(j5,j3)
     &              *za(j6,j1)*za(j6,j4)*za(j4,j5)*s(j1,j2))
     &+(zb(j1,j3)*za(j3,j2)+zb(j1,j6)*za(j6,j2))
     & *(zb(j5,j1)*za(j1,j6)*(zb(j5,j2)*za(j2,j6)+zb(j5,j3)*za(j3,j6))
     &  +zb(j5,j4)*za(j4,j6)*zb(j5,j2)*za(j2,j6))
     & /(zb(j3,j6)*zb(j5,j2)*zb(j5,j3)*za(j6,j1)
     &  *za(j6,j4)*za(j4,j5)*s(j1,j2))

      return
      end

      double complex function amp_2gam2g_mpmppm(j1,j2,j3,j4,j5,j6,za,zb)
c--- -i*Amplitude for the process
c---       0  -->  qb(j5) + q(j6) + g(j1) + g(j2) + gam(j3) + gam(j4)
c--- helicities:     +       -      -       +         -         +
c---
c--- Taken from "Multi-Photon Amplitudes for Next-to-Leading Order QCD"
c---  V. Del Duca, W. Kilgore and F. Maltoni, hep-ph/9910253
c---
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6
      double precision t

c--- Equation (4.3)
c--- note: I have interpreted, e.g. <4|2+5|3> as [4|2+5|3> in order
c--- to get the correct <| and |] required by the helicities
      amp_2gam2g_mpmppm=
     &+za(j1,j6)*zb(j2,j5)*(zb(j4,j2)*za(j2,j3)+zb(j4,j5)*za(j5,j3))
     & /(zb(j1,j6)*za(j2,j5)*zb(j5,j3)*za(j6,j4)*t(j1,j4,j6))
     &+za(j1,j6)*zb(j4,j5)*zb(j6,j2)
     & *(zb(j2,j4)*za(j4,j3)+zb(j2,j5)*za(j5,j3))
     & /(zb(j1,j6)*za(j4,j5)*zb(j5,j3)*s(j1,j2)*t(j1,j2,j6))
     &+za(j3,j6)*zb(j2,j5)*za(j5,j1)
     & *(zb(j4,j2)*za(j2,j1)+zb(j4,j5)*za(j5,j1))
     & /(zb(j3,j6)*za(j2,j5)*za(j6,j4)*s(j1,j2)*t(j3,j4,j6))
     &+za(j1,j6)*zb(j2,j5)*(zb(j2,j4)*za(j4,j3)+zb(j2,j5)*za(j5,j3))
     & /(zb(j1,j6)*zb(j3,j5)*za(j5,j4)*za(j4,j6)*s(j1,j2))
     &+t(j1,j3,j6)*(za(j5,j1)*zb(j2,j6)*zb(j5,j3)*za(j3,j6)
     &             +zb(j1,j2)*zb(j6,j5)*za(j1,j6)*za(j5,j1))
     & /(zb(j1,j6)*zb(j5,j3)*zb(j3,j6)*za(j2,j5)
     &  *za(j5,j4)*za(j4,j6)*s(j1,j2))

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
Describe Users/WarrenJones here.

20080202 12:39:11 nbsp Welcome to the Wiki. Thanks for updating Warren Jones Integrative Massage and Bodywork. As a business owner you may find Welcome to the Wiki/Business Owner useful. Please let Wiki Community us know if you have any questions. Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    %W% %G%
      real function cpu_secs( cput0 )
      real cput0

c     This routine returns a delta value, to get cpu time in seconds
c     from the beginning of the process, pass in a value of 0.0

c     To time a section of code do the following:
c        cpu_init = cpu_secs( 0.0 )
c           ...   code to be timed
c        cpu_used = cpu_secs( cpu_init )

c     This is a base level routine that calls a "C" routine for
c     portability. If you do not have a "C" compiler and your system
c     supports a FORTRAN call to get the cpu seconds, replace the
c     following code with your system's equivalent call.

      real cput1

c*** this is the "C" routine
      call cputmr( cput1 )
c***********************************

      cpu_secs = cput1 - cput0

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

      SUBROUTINE prp_xcfft_box(gmaxxc,bmat,kxc1d,kxc2d,kxc3d)
c*********************************************************************
c     Determine the dimensions kxc1d, kxc2d, kxc3d
c     of the dimension of the xc-energy and potential fft-box
c     needed for the fast calculation in interstitial region
c
c     dimensions kxc(i)d for xc-potential/energy FFT in interstitial
c     added.
c        s. bluegel, IFF, Oct. 97
c
c     subroutine boxdim added to subroutine prp_xcfft
c        s.bluegel, IFF, 18.Nov.97
c*********************************************************************

      USE m_constants
      USE m_ifft, ONLY : ifft235
      USE m_boxdim
      IMPLICIT NONE

      INTEGER kxc1d,kxc2d,kxc3d
      REAL bmat(3,3),gmaxxc

      REAL arltv1,arltv2,arltv3
      INTEGER mk1,mk2,mk3,iofile,ksfft
C ..

C ... Intrinsic Functions
      INTRINSIC int
C
C------->          ABBREVIATIONS
C
C   ksfft=(0,1) : KEY OF SELECTING FFT-PRDOGRAM AND RADIX-TYPE
C                      0  PROGRAM, RADIX-2 ONLY
C                      1  PROGRAM, RADIX-2, RADIX-3,RADIX-5
c   iofile      : device number for in and output
c   gmax        : cut-off wavevector for charge density
c   gmaxxc      : cut-off wavevector for xc-potential and energy
c                 in interstitial
c   arltv(i)    : length of reciprical lattice vector along
c                 direction (i)
c
C---> Determine rkmax box of size mk1, mk2, mk3,
c     for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < gmaxxc
C
      CALL boxdim(bmat,arltv1,arltv2,arltv3)

c     (add 1+1 due to integer rounding, strange k_vector in BZ)
      mk1 = int(gmaxxc/arltv1) + 2
      mk2 = int(gmaxxc/arltv2) + 2
      mk3 = int(gmaxxc/arltv3) + 2
c
      kxc1d = 2*mk1
      kxc2d = 2*mk2
      kxc3d = 2*mk3
c
c---> fft's are usually fastest for low primes
c     (restrict kxcid to: kxcid=  (2**P) * (3**Q) * (5**R)
c
      iofile = oUnit
      ksfft = 1
      !WRITE (oUnit,*) 'minimum: kxc1d,kxc2d,kxc3d',kxc1d,kxc2d,kxc3d
      kxc1d = ifft235(ksfft,kxc1d,2.0)
      kxc2d = ifft235(ksfft,kxc2d,2.0)
      kxc3d = ifft235(ksfft,kxc3d,2.0)
      !WRITE (oUnit,*) 'ifft235: kxc1d,kxc2d,kxc3d',kxc1d,kxc2d,kxc3d

      RETURN
      END

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


module mod_fenwick_tree
  use, intrinsic :: iso_fortran_env
  implicit none
  integer, private, parameter :: intkind = int32
  type fenwick_tree
    private
    integer(int32) :: n, p
    integer(intkind), allocatable :: body(:)
  contains
    procedure :: init => init
    procedure :: add => add
    procedure :: sum0 => sum0
    procedure :: sum1 => sum1
    generic :: sum => sum0, sum1
    procedure :: get => get
    procedure :: lower_bound => lower_bound
    procedure :: update => update
    procedure :: max => max0
  end type
  interface fenwick_tree
    module procedure :: newft0
  end interface
contains
  type(fenwick_tree) function newft0(n) result(res)
    integer(int32), intent(in) :: n
    integer(int32) :: p
    p = 1
    do while (lshift(p, 1) <= n)
      p = lshift(p, 1)
    end do
    res%n = n
    res%p = p
    allocate(res%body(n))
    res%body = 0
  end
  subroutine init(this, n)
    class(fenwick_tree), intent(out) :: this
    integer(int32), intent(in) :: n
    integer(int32) :: p
    p = 1
    do while (lshift(p, 1) <= n)
      p = lshift(p, 1)
    end do
    this%n = n
    this%p = p
    allocate(this%body(n))
    this%body = 0
  end
  subroutine add(this, i, v)
    class(fenwick_tree), intent(inout) :: this
    integer(int32), intent(in) :: i
    integer(intkind), intent(in) :: v
    integer(int32) :: x
    x = i
    do while (x <= this%n)
      this%body(x) = this%body(x) + v
      x = x + and(x, -x)
    end do
  end
  integer(intkind) function sum0(this, i) result(res)
    class(fenwick_tree), intent(in) :: this
    integer(int32), intent(in) :: i
    integer(int32) :: x
    res = 0
    x = i
    do while (x > 0)
      res = res + this%body(x)
      x = x - and(x, -x)
    end do
  end
  integer(intkind) function sum1(this, l, r) result(res)
    class(fenwick_tree), intent(in) :: this
    integer(int32), intent(in) :: l, r
    integer(int32) :: i, j
    res = 0
    i = l - 1
    j = r
    do while (j > i)
      res = res + this%body(j)
      j = j - and(j, -j)
    end do
    do while (i > j)
      res = res - this%body(i)
      i = i - and(i, -i)
    end do
  end
  integer(intkind) function get(this, i) result(res)
    class(fenwick_tree), intent(in) :: this
    integer(int32), intent(in) :: i
    integer(int32) :: x
    res = sum1(this, i, i)
  end
  integer(int32) function lower_bound(this, v) result(res)
    class(fenwick_tree), intent(in) :: this
    integer(intkind), intent(in) :: v
    integer(int32) :: k
    integer(intkind) :: w
    res = 0
    if (v <= 0) return
    k = this%p
    w = v
    do while (k > 0)
      if (res + k <= this%n .and. this%body(res + k) < w) then
        w = w - this%body(res + k)
        res = res + k
      end if
      k = rshift(k, 1)
    end do
    res = res + 1
  end
  subroutine update(this, i, v)
    class(fenwick_tree), intent(inout) :: this
    integer(int32), intent(in) :: i
    integer(intkind), intent(in) :: v
    integer(int32) :: x
    x = i
    do while (x <= this%n)
      if (this%body(x) < v) this%body(x) = v
      x = x + and(x, -x)
    end do
  end
  integer(intkind) function max0(this, i) result(res)
    class(fenwick_tree), intent(in) :: this
    integer(int32), intent(in) :: i
    integer(int32) :: x
    x = i
    res = 0
    do while (x > 0)
      res = max(res, this%body(x))
      x = x - and(x, -x)
    end do
  end
  integer(intkind) function lis_length(a) result(res)
    integer, intent(in) :: a(:)
    type(fenwick_tree) :: bit
    integer(intkind) :: tmp
    integer(int32) :: n, i
    n = size(a)
    call init(bit, n)
    res = 0
    do i = 1, n
      tmp = max0(bit, a(i) - 1) + 1
      res = max(res, tmp)
      call update(bit, a(i), tmp)
    end do
  end
  integer(intkind) function inversion_num(a) result(res)
    integer(int32), intent(in) :: a(:)
    type(fenwick_tree) :: bit
    integer(int32) :: n, i
    n = size(a)
    call init(bit, n)
    res = 0
    do i = 1, n
      call add(bit, a(i), 1_intkind)
      res = res + i - sum0(bit, a(i))
    end do
  end
end module mod_fenwick_tree


program main
    use SoybeanClass
    implicit none

    type(Soybean_) :: soy

    ! let us deform soybean!!
    call soy%init(config="soy.json")
    !call soy%init(config="Tutorial/obj/realSoybeanConfig.json")

    ! (1) 設定事項の確認
    call soy%checkProperties(Simulation=PF_DEFORMATION_ANALYSIS)

    ! この状態では，ただ空のオブジェクトがあるだけ．
    ! (2)材料条件をどう設定する?

    ! 範囲，属性，ID指定できるようにする．
    ! 材料情報を設定する．
    call soy%setProperties(density=.true.)
    call soy%setProperties(YoungModulus=.true.,default_value=1.0d0)
    call soy%setProperties(PoissonRatio=.true.,default_value=0.2d0)
    ! 初期条件を設定する．
    call soy%setProperties(InitialDisplacement=.true.,default_value=0.0d0)
    call soy%setProperties(InitialStress=.true.,default_value=0.0d0)
    ! 境界条件を設定する．
    call soy%setProperties(BoundaryTractionForce=.true.,default_value=0.0d0)
    call soy%setProperties(BoundaryDisplacement=.true.,default_value=0.0d0)
    ! 重力加速度を設定する．
    call soy%setProperties(Gravity=.true.,default_value=0.0d0)

    ! 必要な情報が全て設定されているか確認する．
    call soy%checkProperties(Simulation=PF_DEFORMATION_ANALYSIS)

    ! 必要な情報が全て設定されていれば，戻り値は.true.
    print *, soy%readyFor(Simulation=PF_DEFORMATION_ANALYSIS)

    ! (6) シミュレーションの実行
    !call soy%DeformationAnalysis()

    ! (7) 結果の確認
    !call soy%vtk("soy_deformed/soy")

end program main


text:
      PROGRAM MPEX08
C
C Produce a Mercator projection of the whole globe, using a
C version of MAPUSR which dots the grid lines and dashes the
C continental outlines.
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
        PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C Define the label for the top of the map.
C
      CHARACTER*30 PLBL
C
      DATA PLBL / 'ILLUSTRATING THE USE OF MAPUSR' /
C
C Open GKS.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C Weird up the projection a little.
C
      CALL SUPMAP (9,0.,0.,90.,0.,0.,0.,0.,1,15,2,0,IERR)
C
C Put the label at the top of the plot.
C
      CALL SET (0.,1.,0.,1.,0.,1.,0.,1.,1)
      CALL PWRIT (.5,.975,PLBL,30,2,0,0)
C
C Draw a boundary around the edge of the plotter frame.
C
      CALL BNDARY
C
C Advance the frame.
C
      CALL FRAME
C
C Close GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
C
C Done.
C
      STOP
C
      END
      SUBROUTINE MAPUSR (IPRT)
C
C This version of MAPUSR forces the grid lines to be dotted and
C the outlines to be dashed.
C
C Certain local parameters must be saved from call to call.
C
      SAVE IDTF,IDBD
C
C If IPRT is positive, a part is about to be drawn.  Save the
C dotted/solid flag and/or the distance between dots and then
C reset them and/or the dash pattern.
C
      IF (IPRT.GT.0) THEN
        IF (IPRT.EQ.2) THEN
          CALL MAPGTI ('DL',IDTF)
          CALL MAPGTI ('DD',IDBD)
          CALL MAPSTI ('DL',1)
          CALL MAPSTI ('DD',24)
        ELSE IF (IPRT.EQ.5) THEN
          CALL MAPGTI ('DL',IDTF)
          CALL MAPSTI ('DL',0)
          CALL DASHDB (21845,0,0,0)
        END IF
C
      ELSE
C
C Otherwise, a part has just been drawn.  Restore saved settings
C and/or select a solid dash pattern.
C
        IF (IPRT.EQ.-2) THEN
          CALL MAPSTI ('DL',IDTF)
          CALL MAPSTI ('DD',IDBD)
        ELSE IF (IPRT.EQ.-5) THEN
          CALL MAPSTI ('DL',IDTF)
          CALL DASHDB (65535,0,0,0)
        END IF
C
      END IF
C
C Done.
C
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
       program chkfil
       character*80 card
       read(5,901) card
901    format(a80)
       inptp = index( card, 'NPTP' )
       iplt2 = index( card, 'PLT'  )
       iexit = 0
       if ( inptp .ne. 0 ) iexit = 1
       if ( iplt2 .ne. 0 ) iexit = iexit + 10
       if ( iexit .eq. 0 ) go to 700
       if ( iexit .eq. 1 ) open ( 77, file='nogood1', status='unknown' )
       if ( iexit .eq. 10) open ( 77, file='nogood2', status='unknown' )
       if ( iexit .eq. 11) open ( 77, file='nogood3', status='unknown' )
       write ( 77, * ) ' iexit=',iexit
       close ( 77 )
700    call exit( iexit )
       end


!##############################################################################
! PROGRAM household2
!
! ## The life cycle model with wage risk
!
! This code is published under the GNU General Public License v3
!                         (https://www.gnu.org/licenses/gpl-3.0.en.html)
!
! Authors: Hans Fehr and Fabian Kindermann
!          contact@ce-fortran.com
!
! #VC# VERSION: 1.0  (23 January 2018)
!
!##############################################################################
include "prog05_02m.f90"

program household2

    use globals
    use toolbox

    implicit none
    real*8 :: x(1+n_w)
    real*8 :: fret, low(1+n_w), up(1+n_w)
    integer :: j, iw

    ! discretize w
    call log_normal_discrete(w, weight_w, mu_w, sig_w)

    ! lower and upper border and initial guess
    low = 0d0
    up(1) = mu_w
    up(2:1+n_w) = R*mu_w+w
    x = up/2d0

    ! minimization routine
    call fminsearch(x, fret, low, up, utility)

    ! set up data for output
    do iw = 1, n_w
        wag(1,iw) = mu_w
        inc(1,iw) = mu_w
        sav(1,iw) = a(2,1)

        wag(2,iw) = w(iw)
        inc(2,iw) = w(iw)+R*a(2,1)
        sav(2,iw) = a(3,iw)

        wag(3,iw) = 0d0
        inc(3,iw) = R*a(3,iw)
        sav(3,iw) = 0d0
    enddo

    ! output
    write(*,'(/a/)')' AGE   CONS   WAGE    INC    SAV'
    do j = 1, 3
        write(*,'(i4,4f7.2,a)')j,E(c(j,:)),E(wag(j,:)),&
            E(inc(j,:)),E(sav(j,:)),' (MEAN)'
        write(*,'(4x,4f7.2,a/)')Std(c(j,:)),Std(wag(j,:)),&
            Std(inc(j,:)),Std(sav(j,:)),' (STD)'
    enddo

    write(*,'(/2(a,f6.2))')' E(w) = ',sum(weight_w*w),'   Var(w) = ', &
        sum(weight_w*w**2)-sum(weight_w*w)**2

end program


text: #redirect Health Education and Promotion


text:
program compare
integer a, b
c        fortran 77 I/O statements, for simplicity
read(*,*) a, b

if (a .lt. b) write(*, *) a, ' is less than ', b
if (a .eq. b) write(*, *) a, ' is equal to ', b
if (a .gt. b) write(*, *) a, ' is greater than ', b
end

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
program P18221
    implicit none
    integer :: size, i, j, sungyu_x, sungyu_y, gyosu_x, gyosu_y, top, left, bottom, right, barrier
    integer, dimension(1000, 1000) :: room
    read *, size
    do i = 1, size
        read *, (room(i, j), j = 1, size)
    end do
    do i = 1, size
        do j = 1, size
            if (room(i, j) == 2) then
                sungyu_x = i
                sungyu_y = j
            else if (room(i, j) == 5) then
                gyosu_x = i
                gyosu_y = j
            end if
        end do
    end do
    left = min(sungyu_x, gyosu_x)
    right = max(sungyu_x, gyosu_x)
    top = min(sungyu_y, gyosu_y)
    bottom = max(sungyu_y, gyosu_y)
    barrier = 0
    do i = left, right
        do j = top, bottom
            if (room(i, j) == 1) then
                barrier = barrier + 1
            end if
        end do
    end do
    if (((sungyu_x - gyosu_x) ** 2 + (sungyu_y - gyosu_y) ** 2) >= 25 .and. barrier >= 3) then
        write (*, '(I0)') 1
    else
        write (*, '(I0)') 0
    end if
end program P18221

max_stars_repo_name: daily-boj/kiwiyou

max_issues_repo_name: daily-boj/kiwiyou

max_forks_repo_name: daily-boj/kiwiyou


text:
	INCLUDE 'VICMAIN_FOR'
	SUBROUTINE MAIN44
C
C---- VICAR/IBIS PROGRAM ROWOP.
C     PURPOSE:  TO DELETE OR SELECT SELECTED ROWS OF THE INPUT IBIS
C	INTERFACE FILE WHILE COPYING THE FILE
C
C	REVISION 2	FRANK EVANS	MARCH 1986
C			ADDED SELECT OPTION, PUT IN NCOPY PARAMETER,
C			MADE THE PDF AGREE WITH THE CODE
C
C	REVISION 3	FRANK EVANS	JUNE 1987
C			ADDED PICK OPTION, IMPROVED HELP FILE
C
C	REVISION 4	NILES RITTER	MAY 1994
C			PORTED TO UNIX; USES IBIS2 CALLS
c
c       REVISION 5      BARBARA MCGUFFIE  OCTOBER 26, 1998
c			MODIFIED FOR DOUBLE PRECISION DATA
C
	IMPLICIT NONE

	REAL	BUF(100), PREC, RANGE(200), LOWER(100), UPPER(100)
        real*8  dbuf(100)
	INTEGER	KEYCOL(100), NCOPY(2)
	INTEGER STATUS,IBISIN,IBISOUT,USIZE,OUTSIZE
        INTEGER INRECORD,INRECORDD
	BYTE	COL(4000000)
	INTEGER*4 FULLCOL(1000000)
	INTEGER NKEY,DEF,NRANGE,COUNT,NUMCOPIES,INDEXCOL
	INTEGER RUNIT,WUNIT,CLEN,NCOL,I,OUTLEN,IREC
	INTEGER ICOL,PTR,N,INPTR,J,OFFSET,GR1DIM
	LOGICAL*1 KEEP(1000000)
	LOGICAL  DELETE, SELECT, PICK, XVPTST
	CHARACTER*12 ORG,CFORM
	CHARACTER*5 FORMATS(1024)
	EQUIVALENCE(COL,FULLCOL)  !Only used to save space.

        integer doubflag

C	Get the input parameters

	DELETE = XVPTST('DELETE')
	SELECT = XVPTST('SELECT')
	PICK = XVPTST('PICK')
	CALL XVPARM ('KEYCOL', KEYCOL, NKEY, DEF, 100)
	CALL XVPARM ('RANGE', RANGE, NRANGE, DEF, 200)
	CALL XVP ('PREC', PREC, COUNT)
	CALL XVP ('NCOPY', NCOPY, COUNT)
	NUMCOPIES = NCOPY(1)
	INDEXCOL = NCOPY(2)
	IF (COUNT .EQ. 1)  INDEXCOL = 0

c		Open the input interface file
	CALL XVUNIT(RUNIT,'INP',1,STATUS,' ')
	CALL XVP('GR1DIM',GR1DIM,DEF)

	CALL IBIS_FILE_OPEN( RUNIT, IBISIN, 'READ',
     +     GR1DIM, 0, ' ', ' ',STATUS )
     	IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(RUNIT,STATUS,1)

C    Get the File description
	CALL IBIS_FILE_GET( IBISIN, 'NR', CLEN, 1, 1)
	CALL IBIS_FILE_GET( IBISIN, 'NC', NCOL, 1, 1)
	CALL IBIS_FILE_GET( IBISIN, 'ORG', ORG, 1, 1)
	CALL IBIS_FILE_GET( IBISIN, 'FORMATS', FORMATS, 1, 1024)
C
C    Check for double precision columns
C
         doubflag =  0
         do i = 1, ncol
             if ( formats(i) .eq. 'DOUB' ) then
                 doubflag = 1
                 go to 10
             end if
         end do

C    Create the Record
 10     continue
        if ( doubflag .eq. 0 ) then
  	    CALL IBIS_RECORD_OPEN( IBISIN, INRECORD, ' ',
     +            KEYCOL, NKEY, 'REAL', STATUS)
        else
  	    CALL IBIS_RECORD_OPEN( IBISIN, INRECORDD, ' ',
     +            KEYCOL, NKEY, 'DOUB', STATUS)
        end if
    	IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)

	IF (2*NKEY .NE. NRANGE) THEN
	    CALL XVMESSAGE ('The number of ranges must equal ' //
     * 'the number of key columns.', ' ')
	    CALL ABEND
	ENDIF

	DO I = 1, NKEY
	    LOWER(I) = RANGE(2*I-1) - PREC
	    UPPER(I) = RANGE(2*I) + PREC
	ENDDO

C		Go through each row and flag acceptable rows
	IF (DELETE) THEN
	    OUTLEN = 0
	    DO IREC = 1, CLEN
    		IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)
                 if ( doubflag .eq. 0 ) then
	    	   CALL IBIS_RECORD_READ(INRECORD,BUF,IREC,STATUS )
                else
	    	   CALL IBIS_RECORD_READ(INRECORDD,DBUF,IREC,STATUS )
                end if
		KEEP(IREC) = .TRUE.
		DO I = 1, NKEY
                  if ( doubflag .eq. 0 ) then
		    IF ( (BUF(I) .GE. LOWER(I)) .AND.
     +		         (BUF(I) .LE. UPPER(I)) )  KEEP(IREC) = .FALSE.
                  else
		    IF ( (DBUF(I) .GE. LOWER(I)) .AND.
     +		         (DBUF(I) .LE. UPPER(I)) ) KEEP(IREC) = .FALSE.
                  end if
		ENDDO
		IF (KEEP(IREC)) OUTLEN = OUTLEN+1
	    ENDDO
	ELSE IF (SELECT) THEN
	    OUTLEN = 0
	    DO IREC = 1, CLEN
                 if ( doubflag .eq. 0 ) then
	    	   CALL IBIS_RECORD_READ(INRECORD,BUF,IREC,STATUS )
                else
	    	   CALL IBIS_RECORD_READ(INRECORDD,DBUF,IREC,STATUS )
                end if
    		IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)
		KEEP(IREC) = .TRUE.
		DO I = 1, NKEY
                  if ( doubflag .eq. 0 ) then
		    IF (.NOT. ((BUF(I) .GE. LOWER(I)) .AND.
     +		         (BUF(I) .LE. UPPER(I))) )  KEEP(IREC) = .FALSE.
                  else
		    IF ( .NOT. ((DBUF(I) .GE. LOWER(I)) .AND.
     +		         (DBUF(I) .LE. UPPER(I))) ) KEEP(IREC) = .FALSE.
                  end if
		ENDDO
		IF (KEEP(IREC)) OUTLEN = OUTLEN+1
	    ENDDO
	ELSE IF (PICK) THEN
	    OUTLEN = 0
	    DO IREC = 1, CLEN
                 if ( doubflag .eq. 0 ) then
	    	   CALL IBIS_RECORD_READ(INRECORD,BUF,IREC,STATUS )
                else
	    	   CALL IBIS_RECORD_READ(INRECORDD,DBUF,IREC,STATUS )
                end if
    		IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)
		KEEP(IREC) = .FALSE.
		DO I = 1, NKEY
                  if ( doubflag .eq. 0 ) then
		    IF (((BUF(I) .GE. LOWER(I)) .AND.
     +		         (BUF(I) .LE. UPPER(I))) )  KEEP(IREC) = .TRUE.
                  else
		    IF (((DBUF(I) .GE. LOWER(I)) .AND.
     +		         (DBUF(I) .LE. UPPER(I))) ) KEEP(IREC) = .TRUE.
                  end if
		ENDDO
		IF (KEEP(IREC)) OUTLEN = OUTLEN+1
	    ENDDO
	ENDIF

        if ( doubflag .eq. 0 ) then
  	     CALL IBIS_RECORD_CLOSE( INRECORD, STATUS )
        else
  	     CALL IBIS_RECORD_CLOSE( INRECORDD, STATUS )
        end if
    	IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)

c***************************************************************
c

C		Copy from input to output one column at a time
c		    keeping only the desired rows.
c		If necessary make duplicate copies of the rows
c		    and put out the index numbers.

	CALL XVUNIT(WUNIT,'OUT',1,STATUS,' ')

	CALL IBIS_FILE_OPEN(WUNIT,IBISOUT,'WRITE',
     +       NCOL,NUMCOPIES*OUTLEN, FORMATS, ORG, STATUS)
    	IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(WUNIT,STATUS,1)

	DO ICOL = 1, NCOL
	  IF (ICOL .NE. INDEXCOL) THEN
	     CALL IBIS_COLUMN_GET(IBISOUT,'FORMAT',CFORM ,ICOL,STATUS)
	     IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)

	     IF (CFORM(1:1).EQ.'A') THEN !dont format ASCII
	       CALL IBIS_COLUMN_SET(IBISIN,'U_FORMAT','NONE',ICOL,
     +                              STATUS)
    	       IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)
	       CALL IBIS_COLUMN_SET(IBISOUT,'U_FORMAT','NONE',ICOL,
     +                              STATUS)
	       IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)
	     ENDIF

	     CALL IBIS_COLUMN_GET(IBISOUT,'U_SIZE', USIZE ,ICOL,STATUS)
	     IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)
	     CALL IBIS_COLUMN_READ(IBISIN,COL,ICOL,1,CLEN,STATUS)
	     IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)

	     OUTSIZE = OUTLEN * USIZE

	     PTR = 1
	     INPTR = 1

	     DO I = 1, CLEN
               IF (KEEP(I)) THEN
		   DO J = 1, USIZE
			COL(PTR) = COL(INPTR)
			PTR=PTR+1
			INPTR=INPTR+1
		   ENDDO
	       ELSE
	    	   INPTR = INPTR+USIZE
	       ENDIF
  	     ENDDO

	     OFFSET=OUTSIZE
	     DO N = 1, NUMCOPIES-1
		DO I = 1, OUTSIZE
		   COL(OFFSET+I) = COL(I)
		 ENDDO
		 OFFSET=OFFSET+OUTSIZE
	     ENDDO
	     CALL IBIS_COLUMN_WRITE(IBISOUT,COL,
     +          ICOL,1,NUMCOPIES*OUTLEN,STATUS)
	     IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)
	  ELSE
	     CALL IBIS_COLUMN_SET(IBISOUT,'U_FORMAT','FULL',ICOL,STATUS)
	     IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)
	     OFFSET=0
	     DO N = 1, NUMCOPIES
	        DO I = 1, OUTLEN
		   FULLCOL(OFFSET+I) = N
	        ENDDO
		OFFSET=OFFSET+OUTLEN
	     ENDDO
	     CALL IBIS_COLUMN_WRITE(IBISOUT,FULLCOL,
     +               ICOL,1,NUMCOPIES*OUTLEN,STATUS)
	     IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)
	  ENDIF
	ENDDO

	CALL IBIS_FILE_CLOSE(IBISIN,' ',STATUS)
	IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISIN,STATUS,1)
	CALL IBIS_FILE_CLOSE(IBISOUT,' ',STATUS)
	IF (STATUS.NE.1) CALL IBIS_SIGNAL(IBISOUT,STATUS,1)

	RETURN
	END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
!  cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  L. Zhang
!  Northwestern University
!  This subroutine solves for the residual for all degrees of freedom
!____________________________________________________________________
!  L. Zhang, 06/24/2004
!  Tulane University
!  Revised the subroutine to array
!  cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine block(xloc, dloc, doloc, p, q_p, hk, ien, f_fluids,rngface, f_stress,&
                 ne_local,ien_local,node_local,nn_local,fden,fvis,I_fluid)
    use global_constants
    use run_variables
    use fluid_variables
    use solid_variables, only: nn_solid
    use r_common, only: density_solid, vis_solid
    use mpi_variables
    use pml_variables
    use lumpedmass
    implicit none

    integer ien(nen,ne)
    real* 8 xloc(nsd,nn)
    real* 8 dloc(ndf,nn),doloc(ndf,nn)
    real* 8 p(ndf,nn+nn_PML),q(ndf,nn),hk(ne)
    integer flagpmlelm
    real(8) betaPML

    real* 8 x(nsd,nen)
    real* 8 d(ndf,nen),d_old(ndf,nen)
    real* 8 qvel(ndf,nen),qvoldel(ndf,nen)

    real* 8 eft0,det,effd,effm,effc
    real* 8 sh(0:nsd,nen),ph(0:nsd,nen)
    real* 8 xr(nsd,nsd),cf(nsd,nsd),sx(nsd,nsd)
    real* 8 f_stress(nsd,nsd,nn)

    real* 8 qrs(ndf), qrt(ndf)  ! interpolation and derivatives of qv
    real* 8 drt(ndf),drs(ndf)
    real* 8 dr(nsd,ndf)
    real* 8 u,v,w,pp,ug
    real* 8 tau(nsd,nsd)
    real* 8 hg,taum,tauc,vel,ree, taul
    real* 8 res_c,res_a(nsd),res_t(nsd)
    real* 8 prs_c,prs_t(nsd),p_vec(3),prs_cc(nsd)
    real* 8 mu,nu,ro,g(nsd),sigmaPe(nsd)
    real* 8 tempc(ndf),temp
    real* 8 dtinv,oma,ama
    integer inl, ie, isd, iq, node,jsd
    integer ieface,irng, rngface(neface,ne)
    real* 8 res_pml_c, res_pml_a(nsd), res_pml_t(nsd)

    real* 8 f_fluids(nsd,nn)
    real* 8 I_fluid(nn)
    real* 8 fnode(nsd,nen),fq(nsd)
    integer nn_local
    integer node_local(nn_local)
    !======================================
    ! Defined by Chu
    ! updated by Jack, 05/29/2014, for PML
    real* 8 q_d(ndf,nen)
    real* 8 q_res_c(nen), q_pml_c(nen)
    real* 8 q_p(ndf,nn+nn_PML)
    real* 8 q_res_a(nsd,nen), q_pml_a(nsd,nen)
    real* 8 diag(12)
    !======================================
    ! Define by Xingshi
    ! MPI varialbes & implicit FSI force
    real(8) fden(nn)
    real(8) local_den(nen)
    real(8) fvis(nn)
    real(8) local_vis(nen)
    !======================================
    ! varibles for mpi implementation
    integer ne_local ! # of element on each processor
    integer ien_local(ne_local) ! subregion-->whole region element index
    integer ie_local ! loop parameter
    integer icount
    !--------------------------------------------------
    !  real*8 TC,ZC,RC,P0,dens0          ! defined in global_variables now
    real*8 kappa !compressibility for artificial fluids
    !==================================================
    q_res_a(1:nsd,1:nen) = 0
    q_p(1:ndf,1:nn+nn_PML) = 0
    q_d(1:ndf,1:nen) = 1 !set each ndf for each node as 1
    q_res_c(1:nen) = 0
    !---------------------------------------------------
    !  TC=273.15
    !  ZC=1.4
    !  RC=2.87058e6
    !  P0=1.01325e6
    !  dens0=1.2922e-3
    !................................
    ! Jack 05/20/2014
    ! Defined in global_variables now
    !---------------------------------------------------
    kappa=1.0e4
    !---------------------------------------------------
    dtinv = 1.0/dt
    if(steady) dtinv = 0.0
    oma   = 1.0 - alpha
    ama   = 1.0 - oma
    if (nn_PML_local > 0) betaPML=( uvwpAvgPML(udf)/sqrt(ZC*P0/dens0) )/( 1 - uvwpAvgPML(udf)**2/(ZC*P0/dens0) )
    !betaPML=0.0d0
    !=================================================
    do icount=1, nn_local
        node=node_local(icount)
        p(1:nsd,node)=p(1:nsd,node)+f_fluids(1:nsd,node)
    enddo

    include "reconstruct_tauij.f90"
    !===================================================
    ! f_fluids is actually the FSI force at fluid nodes,
    ! then we will just subscrib it from p(!:nsd) which is the
    ! residuals for the momentum equations.
    ! 1 let fnod=0 then fndoe and fq 's effects will be canceled out
    ! 2 do the subscribition after the elements loop
    ! Xingshi 09/15/2008
    !===================================================
    do ie_local=1,ne_local              ! loop over subregion elements
        ie=ien_local(ie_local)
        flagpmlelm=0
        do inl=1,nen
            x(1:nsd,inl) = xloc(1:nsd,ien(inl,ie))
            fnode(:,inl) = 0.0
            d(1:ndf,inl) =  dloc(1:ndf,ien(inl,ie))
            d_old(1:ndf,inl) = doloc(1:ndf,ien(inl,ie))
            f_stress(1:nsd,1:nsd,ien(inl,ie)) = 0.0
            !----------------------------------------------
            node=ien(inl,ie)
            local_den(inl)=(dloc(ndf,node)*dens0/(ZC*P0)+dens0)*(1.0-I_fluid(node))+&
                           (density_solid+den_liq)*I_fluid(node)
            local_vis(inl)=fvis(ien(inl,ie))
            !------------------- check PML nodes ------------------------------
            if (seqcPML(node) .ne. 0) flagpmlelm=flagpmlelm+1
            !------------------------------------------------------------------
        enddo
        hg = hk(ie)
        !======================= LOOP THRU QUADRATURE POINTS ==========================
        do iq=1,nquad  ! loop over the quadrature points in each element
        !...  calculate the shape function and the weight at quad point
            if (nsd==2) then
                if (nen.eq.3) then !calculate shape function at quad point
                    include "sh2d3n.h"
                elseif (nen.eq.4) then
                    include "sh2d4n.h"
                endif
            elseif (nsd==3) then
                if (nen.eq.4) then !calculate shape function at quad point
                    include "sh3d4n.h"
                elseif (nen.eq.8) then
                    include "sh3d8n.h"
                endif
            endif
            eft0 = abs(det) * wq(iq) ! calculate the weight at each quad pt
            ph(0:nsd,1:nen) = sh(0:nsd,1:nen)*eft0
            !...  initialize d, dd/dx, dd/dy, dd/dz, dd/dt
            drs(:) = 0.0
            drt(:) = 0.0
            dr(1:nsd,1:ndf)=0.0
            !-----PML-------------
            qrs(1:ndf) = 0.0
            qrt(1:ndf) = 0.0
            !---------------------
            fq(:)=0.0
            !... calculate vi, dvi/dxj
            do inl=1,nen
                node=ien(inl,ie)
                tempc(1:nsd) = ama*d(1:nsd,inl)+oma*d_old(1:nsd,inl)
                drs(1:nsd) = drs(1:nsd)+sh(0,inl)*tempc(1:nsd)
                do isd=1,nsd
                    dr(isd,1:nsd) = dr(isd,1:nsd)+sh(isd,inl)*tempc(1:nsd)
                enddo
                fq(:) = fq(:) + sh(0,inl)*fnode(:,inl)
            enddo

            ro=0.0
            mu=0.0
            sigmaPe(1:nsd)=0.0
            !... calculate dvi/dt, p, dp/dxi
            do inl=1,nen
                node=ien(inl,ie)
                drt(1:ndf)=drt(1:ndf)+sh(0,inl)*(d(1:ndf,inl)-d_old(1:ndf,inl))*dtinv
                !------PML---------
                if (seqcPML(node) > 0) then
                    qrs(1:ndf) = qrs(1:ndf)+sh(0,inl)*qv(1:ndf,node)
                    qrt(1:ndf) = qrt(1:ndf)+sh(0,inl)*(qv(1:ndf,node)-qvold(1:ndf,node))*dtinv
                endif
                sigmaPe(1:nsd)=sigmaPe(1:nsd)+sh(0,inl)*sigmaPML(1:nsd,node)
                !------------------
                drs(pdf)=drs(pdf)+sh(0,inl)*d(pdf,inl)
                dr(1:nsd,pdf)=dr(1:nsd,pdf)+sh(1:nsd,inl)*d(pdf,inl)
                !----------------------------------------
                ro=ro+sh(0,inl)*local_den(inl)  ! local fluid density
                mu=mu+sh(0,inl)*local_vis(inl)  ! local fluid viscosity
            enddo
            !... define u=v1, v=v2, w=v3, pp=p
            if (nsd==2) then
                u = drs(udf)
                v = drs(vdf)
                w = 0.0
            elseif (nsd==3) then
                u = drs(udf)
                v = drs(vdf)
                w = drs(wdf)
            endif
            pp= drs(pdf)

            if(stokes) then ! if stokes flow
                u = 0.0
                v = 0.0
                w = 0.0
            endif
            !....  calculate liquid constant and gravity
            g(1:nsd)  = gravity(1:nsd)  ! gravatitional force

            ! believe nu is calculated only for turbulent model
            if (nsd==2) then
                nu = delta(4)*turb_kappa**2*hg**2 * sqrt(2*dr(1,1)**2+(dr(2,1)+dr(1,2))**2 &
                                                        +2*dr(2,2)**2)
            elseif (nsd==3) then
                nu = delta(4)*turb_kappa**2*hg**2 * sqrt(2*dr(1,1)**2+(dr(2,1)+dr(1,2))**2 &
                                                        +2*dr(2,2)**2+(dr(3,1)+dr(1,3))**2 &
                                                        +2*dr(3,3)**2+(dr(3,2)+dr(2,3))**2)
            endif
            mu = mu + nu*ro

            !....  calculate each term in the residual equation
            res_c = 0.0                ! Continuity residual
            res_pml_c = 0.0
            res_pml_a(1:nsd) = 0.0
            q_pml_c(1:nen) = 0.0
            q_pml_a(1:nsd,1:nen) = 0.0
            if (nsd==2) then
                do inl=1,nen
                    node=ien(inl,ie)
                    res_c = res_c + (ZC+pp/P0)*(sh(xsd,inl)*d(udf,inl)+ sh(ysd,inl)*d(vdf,inl))+ &
                            (u*sh(xsd,inl)*d(pdf,inl)+v*sh(ysd,inl)*d(pdf,inl))/P0*(1.0-I_fluid(node))
                    q_res_c(inl) = (sh(0,inl)*dtinv+(ZC+pp/P0)*(u*sh(xsd,inl)+v*sh(ysd,inl)))/P0*(1.0-I_fluid(node)) &
                                 +sh(0,inl)*dtinv/kappa*I_fluid(node)
                enddo
            elseif (nsd==3) then
                do inl=1,nen
                    node = ien(inl,ie)
                    res_c = res_c + ZC*(sh(xsd,inl)*d(udf,inl) &
                                  + sh(ysd,inl)*d(vdf,inl) &
                                  + sh(zsd,inl)*d(wdf,inl)) &
                                  + (u*sh(xsd,inl)*d(pdf,inl)+v*sh(ysd,inl)*d(pdf,inl)+w*sh(zsd,inl)*d(pdf,inl))&
                                  / P0*(1.0-I_fluid(node))
                    q_res_c(inl) = (sh(0,inl)*dtinv+ZC*(u*sh(xsd,inl)+v*sh(ysd,inl)+w*sh(zsd,inl)) )/P0*(1.0-I_fluid(node)) &
                                   + sh(0,inl)*dtinv/kappa*I_fluid(node)
                enddo
            endif

            do inl=1,nen
                node=ien(inl,ie)
                res_c=res_c+sh(0,inl)*( d(ndf,inl)-d_old(ndf,inl) )*dtinv* &
                      (1.0/P0*(1.0-I_fluid(node))+1.0/kappa*I_fluid(node))
            enddo  ! add dp/dt term

            do isd = 1, nsd
                if (nsd==2) then
                    res_a(isd)=ro*( drt(isd) + u*dr(1,isd) + v*dr(2,isd) &
                                  + drs(isd)*(dr(1,1) + dr(2,2)) - g(isd) ) &
                               + ( drs(isd)*(u*dr(1,pdf) + v*dr(2,pdf)) + drt(pdf)*drs(isd) )*dens0/(ZC*P0) &
                               - fq(isd)
                    do inl = 1, nen
                        q_res_a(isd,inl)=ro*(sh(0,inl)*q_d(isd,inl)*dtinv+u*sh(1,inl)*q_d(isd,inl)+v*sh(2,inl)*q_d(isd,inl) &
                                             + q_d(isd,inl)*(dr(1,1) + dr(2,2))) &
                                         + q_d(isd,inl)*(u*dr(1,pdf) + v*dr(2,pdf) + drt(pdf))*dens0/(ZC*P0)
                    enddo ! get res_a for u and v for momentum equation
                elseif (nsd==3) then
                    res_a(isd)=ro*(drt(isd)+u*dr(1,isd)+v*dr(2,isd)+w*dr(3,isd)-g(isd))-fq(isd)
                    do inl = 1, nen
                        q_res_a(isd,inl)=ro*(sh(0,inl)*dtinv+u*sh(1,inl)+v*sh(2,inl)+w*sh(3,inl))
                    enddo ! get res_a for u v and w for momentum equation
                endif
            enddo

            !c....   Calculate stress tau(ij) and pressure
            do isd = 1,nsd
                do jsd = 1,nsd
                    tau(isd,jsd) = mu*(dr(isd,jsd) + dr(jsd,isd))
                    do inl=1,nen
                        ! Compute fluid stress
                        f_stress(isd,jsd,ien(inl,ie))= f_stress(isd,jsd,ien(inl,ie)) + tau(isd,jsd)
                    enddo
                enddo
                if (nsd==2) then
                    tau(isd,isd)=tau(isd,isd)-2.0*mu*(dr(1,1)+dr(2,2))/3.0
                elseif (nsd==3) then
                    tau(isd,isd)=tau(isd,isd)-2.0*mu*(dr(1,1)+dr(2,2)+dr(3,3))/3.0
                endif
            enddo

!--------------- PML terms in Continuity/Momentum eqns --------------------------
            if (flagpmlelm > 0) then
                !------------ residual correction for NS eqns ----------------------
                if (nsd==2) then
                    do inl=1,nen
                        node=ien(inl,ie)
                        !q_pml_c(inl)=(sigmaPe(xsd)+sigmaPe(ysd))/P0
                        q_pml_a(xsd,inl)=(sigmaPe(xsd)+sigmaPe(ysd))*ro*sh(0,inl) + betaPML*sigmaPe(xsd)*ro*u*sh(0,inl)
                        q_pml_a(ysd,inl)=(sigmaPe(xsd)+sigmaPe(ysd))*ro*sh(0,inl) + betaPML*sigmaPe(xsd)*ro*u*sh(0,inl)
                    enddo
                    res_pml_c=( sigmaPe(xsd)*qrs(pdf) + sigmaPe(ysd)*(pp-qrs(pdf)) )/P0 + &
                                betaPML*sigmaPe(xsd)*((ZC+pp/P0)*u-ZC*uvwpAvgPML(udf))
                    res_pml_a(xsd)=ro*( sigmaPe(xsd)*qrs(udf)+sigmaPe(ysd)*(u-qrs(udf)) ) + betaPML*sigmaPe(xsd)*ro*u*u
                    res_pml_a(ysd)=ro*( sigmaPe(xsd)*qrs(vdf)+sigmaPe(ysd)*(v-qrs(vdf)) ) + betaPML*sigmaPe(xsd)*ro*u*v
                endif
                res_c=res_c+res_pml_c
                res_a(1:nsd)=res_a(1:nsd)+res_pml_a(1:nsd)
                !q_res_c(1:nen)=q_res_c(1:nen)+q_pml_c(1:nen)
                q_res_a(1:nsd,1:nen)=q_res_a(1:nsd,1:nen)+q_pml_a(1:nsd,1:nen)
                !------------ residual for auxiliary variables ----------------------
                ! res_pml_a & res_pml_c have been recycled to calculate dq/dt+\sigma*q+dF1/dx
                res_pml_c=0.0
                res_pml_a(:)=0.0
                q_pml_a(:,:)=0.0
                q_pml_c(:)=0.0
                if (nsd==2) then
                    ! PML Continuity
                    res_pml_c = ( qrt(pdf) + sigmaPe(xsd)*qrs(pdf) + (ZC*P0+pp)*dr(xsd,udf) + u*dr(xsd,pdf) )/P0 + &
                                betaPML*sigmaPe(xsd)*( (ZC+pp/P0)*u - ZC*uvwpAvgPML(udf) )
                    do inl=1,nen
                        node=ien(inl,ie)
                        q_pml_c(inl) = ( sh(0,inl)*dtinv + (ZC+pp/P0)*(u*sh(xsd,inl)) )/P0 + &
                                       betaPML*sigmaPe(xsd)*ZC*sh(0,inl)
                    enddo
                    ! PML Momentum
                    do isd = 1, nsd
                        res_pml_a(isd) = ro*( qrt(isd) + sigmaPe(xsd)*qrs(isd) + u*dr(1,isd) + drs(isd)*dr(1,1) ) + &
                                         ( drs(isd)*u*dr(1,pdf) + qrt(pdf)*drs(isd) )*dens0/(ZC*P0) + &
                                         betaPML*sigmaPe(xsd)*ro*u*drs(isd)
                        do inl = 1, nen
                            q_pml_a(isd,inl)=ro*( sh(0,inl)*q_d(isd,inl)*dtinv+u*sh(1,inl)*q_d(isd,inl) + q_d(isd,inl)*dr(1,1) ) &
                                             + q_d(isd,inl)*(u*dr(1,pdf) + qrt(pdf))*dens0/(ZC*P0) &
                                             + betaPML*sigmaPe(xsd)*ro*u*sh(0,inl)
                        enddo
                    enddo
                    res_pml_a(xsd) = res_pml_a(xsd) + betaPML*sigmaPe(xsd)*pp
                endif

                do inl=1,nen
                    node = ien(inl,ie)
                    if (seqcPML(node) > 0) then
                        p(pdf,nn+seqcPML(node)) = p(pdf,nn+seqcPML(node)) - ph(0,inl)*res_pml_c
                        q_p(pdf,nn+seqcPML(node)) = q_p(pdf,nn+seqcPML(node)) - ph(0,inl)*q_pml_c(inl)
                        p(1:nsd,nn+seqcPML(node)) = p(1:nsd,nn+seqcPML(node)) - ph(0,inl)*res_pml_a(1:nsd)
                        q_p(1:nsd,nn+seqcPML(node)) = q_p(1:nsd,nn+seqcPML(node)) - ph(0,inl)*q_pml_a(1:nsd,inl)
                        p(udf,nn+seqcPML(node)) = p(udf,nn+seqcPML(node)) + ph(xsd,inl)*(pp-tau(1,1))
                        p(vdf,nn+seqcPML(node)) = p(vdf,nn+seqcPML(node)) - ph(xsd,inl)*tau(2,1)
                        q_p(udf,nn+seqcPML(node))=q_p(1,nn+seqcPML(node)) + ph(xsd,inl)*mu*(sh(1,inl)*q_d(1,inl)*(2.0-2.0/3.0))
                        q_p(vdf,nn+seqcPML(node))=q_p(2,nn+seqcPML(node)) + ph(xsd,inl)*mu*(sh(1,inl)+sh(2,inl))*q_d(2,inl)
                    endif
                enddo
                !--------------------------------------------------------------------
            endif
!---------------------------------------------------------------------------------------
            res_t(1:nsd) = dr(1:nsd,pdf) + res_a(1:nsd)

            include "add_tauijcj_resmomentum.f90"

            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            ! Mickael 02/01/2005
            ! TAUm, TAUc and TAUl (l for  lsic), See Tezduyar and Sathe, 2003
            vel  = sqrt(u*u+v*v+w*w)     ! magnitude of the velocity
            ree  = vel*hg*ro/(2.0*mu)    ! Reynolds number
            if(steady.or.(.not.taudt)) then !stablization, taum
                taum = 1.0/sqrt((2.0*vel/hg)**2+(4.0*mu/hg/hg)**2)
            else
                taum = 1.0/sqrt((2.0/dt)**2+(2.0*vel/hg)**2+(4.0*mu/hg/hg)**2)
            endif
            tauc = taum/ro

            if (ree.le.3.0) then
                taul = hg*vel*ree/6.0
            else
                taul = hg*vel/2.0
            endif

            prs_t(1:nsd) = res_t(1:nsd)*taum
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            ! Mickael 02/01/2005
            ! TAUl (l for  lsic), See Tezduyar and Sathe, 2003
            prs_c = ro*res_c*taul
            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            prs_cc(1:nsd) = res_t(1:nsd)*tauc
            !.... calculate the residual at each degree of freedom
            do inl=1,nen ! loop over number of nodes in an element
                node = ien(inl,ie)
                if (nsd==2) then
                    temp = ro*(u*ph(xsd,inl)+v*ph(ysd,inl))
                elseif (nsd==3) then
                    temp = ro*(u*ph(xsd,inl)+v*ph(ysd,inl)+w*ph(zsd,inl))
                endif
!=============================================================
                ! Continuty Equation
                p(pdf,node) = p(pdf,node)-ph(0,inl)*res_c
                q_p(pdf,node) = q_p(pdf,node)+ph(0,inl)*q_res_c(inl) ! get diag for continuity
                ! Momentum Equation (Euler Residual)
                ! why originally it is minus?
                p(1:nsd,node) = p(1:nsd,node)-ph(0,inl)*res_a(1:nsd)
                q_p(1:nsd,node) = q_p(1:nsd,node)+ph(0,inl)*q_res_a(1:nsd,inl)
!=============================================================
                ! Momentum Equation (C: Diffusion terms)
                if (nsd==2) then
                    do isd=1,nsd
                        p(isd,node)=p(isd,node) + ph(isd,inl)*pp -   &
                                                  ph(1,inl)*tau(1,isd) -  &
                                                  ph(2,inl)*tau(2,isd)
                        !p(isd,node)=p(isd,node) + mu*ph(isd,inl)*(dr(1,1)+dr(2,2))*2.0/3.0
                    enddo
                    q_p(1,node)=q_p(1,node)+ph(1,inl)*mu*(sh(1,inl)*q_d(1,inl)*(2.0-2.0/3.0)) + &
                                            ph(2,inl)*mu*sh(2,inl)*q_d(1,inl)
                    q_p(2,node)=q_p(2,node)+ph(1,inl)*mu*sh(1,inl)*q_d(2,inl) + &
                                            ph(2,inl)*(mu*sh(2,inl)*q_d(2,inl)*(2.0-2.0/3.0))
                elseif (nsd==3) then
                    do isd=1,nsd
                        p(isd,node)=p(isd,node) + ph(isd,inl)*pp -   &
                                                  ph(1,inl)*tau(1,isd) -  &
                                                  ph(2,inl)*tau(2,isd) -  &
                                                  ph(3,inl)*tau(3,isd)
                        !p(isd,node)=p(isd,node)+mu*ph(isd,inl)*(dr(1,1)+dr(2,2)+dr(3,3))*2.0/3.0
                    enddo
                    q_p(1,node)=q_p(1,node)+ph(1,inl)*mu*(sh(1,inl)*(2.0 - 2.0/3.0))+&
                                            ph(2,inl)*mu*sh(2,inl)+ph(3,inl)*mu*sh(3,inl)
                    q_p(2,node)=q_p(2,node)+ph(1,inl)*mu*sh(1,inl)+ph(2,inl)*mu*(sh(2,inl)*(2.0 -2.0/3.0))+&
                                            ph(3,inl)*mu*sh(3,inl)
                    q_p(3,node)=q_p(3,node)+ph(1,inl)*mu*sh(1,inl)+ph(2,inl)*mu*sh(2,inl)+&
                                            ph(3,inl)*mu*(sh(3,inl)*(2.0-2.0/3.0))
                endif
                ! Stablization with Tau_moment
                if (nsd==2) then
                    p(pdf,node) = p(pdf,node) - ph(xsd,inl)*prs_cc(udf)  &
                                              - ph(ysd,inl)*prs_cc(vdf)
                    q_p(pdf,node) = q_p(pdf,node)+tauc*(sh(1,inl)*sh(1,inl)+sh(2,inl)*sh(2,inl))*eft0
                elseif (nsd==3) then
                    p(pdf,node) = p(pdf,node) - ph(xsd,inl)*prs_cc(udf)  &
                                              - ph(ysd,inl)*prs_cc(vdf)  &
                                              - ph(zsd,inl)*prs_cc(wdf)
                    q_p(pdf,node) = q_p(pdf,node)+tauc*(sh(1,inl)**2+sh(2,inl)**2+sh(3,inl)**2)*eft0
                endif             ! Stablization with Tau_cont
                p(1:nsd,node) = p(1:nsd,node) - prs_t(1:nsd)*temp - ph(1:nsd,inl)*prs_c
                q_p(1:nsd,node) = q_p(1:nsd,node)+taum*temp*q_res_a(1:nsd,inl)+ &
                                         ph(1:nsd,inl)*taul*ro*sh(1:nsd,inl)
            enddo
        enddo ! end of qudrature pts loop
    enddo ! end of element loop
continue
!=====================================
! Apply boundary condition du/dx=0 on outedge
!call out2d4n(rngface,dloc,ien,xloc,p)
!======================================
end subroutine block

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


module sfldy
  use ode

  type SFluidQty
    real(dp), allocatable :: w(:) ! Cell width and height and depth.
    real(dp), allocatable :: o(:) ! Offset from top-left of the grid, or top-left-front in 3D.
    real(dp) :: cs                ! Cell size.
  end type SFluidQty

  interface interpolate
    module procedure :: lint
  end interface interpolate
contains

  subroutine sfldyEuler(u, x, h)
    ! Forward Euler for the integration of the velocities.
    implicit none
    class(SFluidQty), intent(in) :: u(:)
    real(dp), intent(in)         :: h ! Fluid vectors (u, v, w), time step
    real(dp), intent(inout)      :: x(:)
    real(dp)                     :: lint
    real(dp)                     :: du(size(u))
    integer                      :: i

    ! Obtain terivatives.
    do i = 1, size(u)
      du(i) = interp(u(i), x) / u(i) % cs
    end do

    ! Integrate x.
    x = x + du * h

  end subroutine sfldyEuler

  function interp(u, x) result(interpolation)
    implicit none
    class(SFluidQty), intent(in)     :: u
    real(dp), intent(in)             :: x(:)
    real(dp)                         :: interpolation
    real(dp), dimension(size(x))     :: xt
    integer,  dimension(size(x))     :: ix
    real(dp), dimension(2**size(x))  :: grid
    real(dp), parameter              :: offset = 1+1d-30

    ! Clamping vertex coordinates.
    xt = min(max(x - u % o, 0._dp), u % w - offset)
    ! Making vertices integers.
    ix = nint(x)
    ! Setting up origin to the top-left or top-left-front corner.
    xt = xt - ix

    ! Find grid positions.
    grid = FindCoord(u, x)

    ! This is for 2D only.
    interpolation = lint( [ lint(grid(1:2), x(1)), lint(grid(3:4), x(1)) ], x(2) )
  end function interp

  function lint(c, x)
    implicit none
    real(dp), intent(in) :: c(:), x ! Coeffcients and x
    real(dp)             :: lint

    lint = ( c(1) - c(2) ) * x - c(1)
  end function lint

  function FindCoord(u, x) result(grid)
    ! Scales the offset to the cell width.
    implicit none
    real(dp), intent(in)         :: x(:)
    class(SFluidQty), intent(in) :: u
    real(dp)                     :: grid(2**size(x)), tmp
    integer                      :: i, j, k, sx, sg
    integer, parameter :: disp(24) = [0,1,0,1,0,1,0,1,0,0,1,1,&
                                      0,0,1,1,0,0,0,0,1,1,1,1]

    sx  = size(x)
    sg  = size(grid)

    ! Loop over Grid Displacements.
    gd: do i = 1, sg
      ! Sum Only Relevant Quantities.
      sorq1: do j = 1, sx
        ! Save the value of x(i) in tmp.
        tmp = x(j) + disp(i + 8*(j-1))
        sorq2: do k = 1, sx
          ! If j = i we need to skip to the next value of j.
          if (j == k) cycle sorq2
          ! Check whether this is correct, do not be an idiot!
          tmp = tmp + ( x(k) + disp(i + 8*(k-1)) ) * u % w(j)
          ! Check whether this is correct, do not be an idiot!
        end do sorq2
        ! Assign tmp to the coordinate i.
        grid(i + 8*(j-1)) = tmp
      end do sorq1
    end do gd
  end function FindCoord

end module sfldy


text:
      BLOCK DATA IFX1BD
CIFX1BD
C     DEFINITION OF VARIABLES IN /IFPX1/ AND /IFPX0/
C*****
C
C     COMMON /IFPX1/
C     --------------
C
C     N          = TOTAL NUMBER OF PAIRED ENTRIES IN THE IBD AND
C                  IPR ARRAYS
C                = (TOTAL DIMENSION OF IBD + ACTIVE IPR ARRAYS)/2
C
C     IBD ARRAYS = ARRAYS CONTAINING PAIRED ENTRIES OF BULK DATA
C                  CARD NAMES
C
C     IPR ARRAYS = ARRAYS CONTAINING PAIRED ENTRIES OF BULK DATA
C                  PARAMETER NAMES
C
C     CAUTION 1 -- THE TOTAL DIMENSION OF THE IBD AND IPR ARRAYS
C                  MUST BE A MULTIPLE OF 62 (OR, IN OTHER WORDS,
C                  AN EVEN MULTIPLE OF 31)
C
C                  SEE NOTES 1 AND 2 BELOW
C
C     ICC ARRAYS = ARRAYS CONTAINING PAIRED ENTRIES OF CASE CONTROL
C                  FLAG NAMES FOR USE IN RESTART RUNS
C
C     CAUTION 2 -- THE TOTAL DIMENSION OF THE ICC ARRAYS MUST BE A
C                  MULTIPLE OF 62 (OR, IN OTHER WORDS, AN EVEN
C                  MULTIPLE OF 31)
C
C                  SEE NOTE 3 BELOW
C
C     NOTES
C     -----
C
C              1.  IF NEW BULK DATA CARD NAMES ARE TO BE ADDED,
C                  USE THE EXISTING PADDING WORDS (OF THE 4H****
C                  TYPE) IN THE IBD ARRAYS.  IF NECESSARY, EXPAND
C                  THE IBD ARRAYS KEEPING CAUTION 1 IN MIND.
C
C              2.  IF NEW BULK DATA PARAMETER NAMES ARE TO BE ADDED,
C                  USE THE EXISTING PADDING WORDS (OF THE 4H****
C                  TYPE) IN THE IPR ARRAYS.  IF NECESSARY, EXPAND
C                  THE IPR ARRAYS KEEPING CAUTION 1 IN MIND.
C
C              3.  IF NEW CASE CONTROL FLAG NAMES ARE TO BE ADDED,
C                  USE THE EXISTING PADDING WORDS (OF THE 4H****
C                  TYPE) IN THE ICC ARRAYS.  IF NECESSARY, EXPAND
C                  THE ICC ARRAYS KEEPING CAUTION 2 IN MIND.
C
C              4.  THE IBD ARRAYS ARE IN SYSCHRONIZTION WITH THE I ARRAY
C                  IN IFX2BD, IFX3BD, IFX4BD, IFX5BD, AND IFX6BD
C                  (E.G. CONM1 POSITIONS IN IBD2, CONTINUATION 3, THE DA
C                  FOR CONM1 IN IFX2BD IS IN I2, CONTINUATION 3 CARD)
C*****
C
C     COMMON /IFPX0/
C     --------------
C
C     LBDPR      = (TOTAL DIMENSION OF THE IBD AND IPR ARRAYS)/62
C
C     LCC        = (TOTAL DIMENSION OF THE ICC ARRAYS)/62
C
C     IWRDS      = ARRAY WHOSE DIMENSION IS EQUAL TO (LBDPR + LCC).
C                  ALL (LBDPR + LCC) WORDS IN THE ARRAY INITIALLY
C                  SET TO ZERO.
C
C     IPARPT     = POINTER THAT POINTS TO THE PAIRED ENTRY IN THE
C                  IBD AND IPR ARRAYS THAT CONTAINS THE FIRST
C                  BULK DATA PARAMETER NAME.  AS PER THE DIFINITIONS
C                  OF THE VARIABLES IN COMMON /IFPX1/, THIS POINTS
C                  TO THE FIRST WORD OF THE IPR1 ARRAY.  HENCE, WE
C                  HAVE --
C
C                  IPARPT = 1 + (TOTAL DIMENSION OF IBD ARRAYS)/2
C
      COMMON /IFPX0 /    LBDPR    , LCC      , IWRDS(18), IPARPT
      COMMON /IFPX1 / N, IBD1(100), IBD2(100), IBD3(100), IBD4(100),
     1                   IBD5(100), IBD6(100), IBD7(100), IBD8(100),
     2                   IPR1(100), IPR2( 92),
     3                   ICC1(100), ICC2( 24)
C
C*****
C     INITIALIZATION OF VARIABLES IN COMMON /IFPX1/
C*****
C
      DATA N / 496 /
C*****
C     THE IBD ARRAYS CONTAIN PAIRED ENTRIES OF BULK DATA CARD NAMES
C*****
      DATA IBD1 /
     1     4HGRID,4H    , 4HGRDS,4HET  , 4HADUM,4H1   , 4HSEQG,4HP   ,
     *     4HCORD,4H1R  , 4HCORD,4H1C  , 4HCORD,4H1S  , 4HCORD,4H2R  ,
     *     4HCORD,4H2C  , 4HCORD,4H2S  , 4HPLOT,4HEL  , 4HSPC1,4H    ,
     3     4HSPCA,4HDD  , 4HSUPO,4HRT  , 4HOMIT,4H    , 4HSPC ,4H    ,
     *     4HMPC ,4H    , 4HFORC,4HE   , 4HMOME,4HNT  , 4HFORC,4HE1  ,
     *     4HMOME,4HNT1 , 4HFORC,4HE2  , 4HMOME,4HNT2 , 4HPLOA,4HD   ,
     5     4HSLOA,4HD   , 4HGRAV,4H    , 4HTEMP,4H    , 4HGENE,4HL   ,
     *     4HPROD,4H    , 4HPTUB,4HE   , 4HPVIS,4HC   , 4HADUM,4H2   ,
     *     4HPTRI,4HA1  , 4HPTRI,4HA2  , 4HPTRB,4HSC  , 4HPTRP,4HLT  ,
     7     4HPTRM,4HEM  , 4HPQUA,4HD1  , 4HPQUA,4HD2  , 4HPQDP,4HLT  ,
     *     4HPQDM,4HEM  , 4HPSHE,4HAR  , 4HPTWI,4HST  , 4HPMAS,4HS   ,
     *     4HPDAM,4HP   , 4HPELA,4HS   , 4HCONR,4HOD  , 4HCROD,4H    ,
     9     4HCTUB,4HE   , 4HCVIS,4HC   /
      DATA IBD2 /
     1     4HADUM,4H3   , 4HCTRI,4HA1  , 4HCTRI,4HA2  , 4HCTRB,4HSC  ,
     *     4HCTRP,4HLT  , 4HCTRM,4HEM  , 4HCQUA,4HD1  , 4HCQUA,4HD2  ,
     *     4HCQDP,4HLT  , 4HCQDM,4HEM  , 4HCSHE,4HAR  , 4HCTWI,4HST  ,
     3     4HCONM,4H1   , 4HCONM,4H2   , 4HCMAS,4HS1  , 4HCMAS,4HS2  ,
     *     4HCMAS,4HS3  , 4HCMAS,4HS4  , 4HCDAM,4HP1  , 4HCDAM,4HP2  ,
     *     4HCDAM,4HP3  , 4HCDAM,4HP4  , 4HCELA,4HS1  , 4HCELA,4HS2  ,
     5     4HCELA,4HS3  , 4HCELA,4HS4  , 4HMAT1,4H    , 4HMAT2,4H    ,
     *     4HCTRI,4HARG , 4HCTRA,4HPRG , 4HDEFO,4HRM  , 4HPARA,4HM   ,
     *     4HMPCA,4HDD  , 4HLOAD,4H    , 4HEIGR,4H    , 4HEIGB,4H    ,
     7     4HEIGC,4H    , 4HADUM,4H4   , 4H    ,4H    , 4HMATS,4H1   ,
     *     4HMATT,4H1   , 4HOMIT,4H1   , 4HTABL,4HEM1 , 4HTABL,4HEM2 ,
     *     4HTABL,4HEM3 , 4HTABL,4HEM4 , 4HTABL,4HES1 , 4HTEMP,4HD   ,
     9     4HADUM,4H5   , 4HADUM,4H6   /
      DATA IBD3 /
     1     4HADUM,4H7   , 4HMATT,4H2   , 4HADUM,4H8   , 4HCTOR,4HDRG ,
     *     4HSPOI,4HNT  , 4HADUM,4H9   , 4HCDUM,4H1   , 4HCDUM,4H2   ,
     *     4HCDUM,4H3   , 4HCDUM,4H4   , 4HCDUM,4H5   , 4HCDUM,4H6   ,
     3     4HCDUM,4H7   , 4HCDUM,4H8   , 4HCDUM,4H9   , 4HPDUM,4H1   ,
     *     4HPDUM,4H2   , 4HPDUM,4H3   , 4HDMI ,4H    , 4HDMIG,4H    ,
     *     4HPTOR,4HDRG , 4HMAT3,4H    , 4HDLOA,4HD   , 4HEPOI,4HNT  ,
     5     4HFREQ,4H1   , 4HFREQ,4H    , 4HNOLI,4HN1  , 4HNOLI,4HN2  ,
     *     4HNOLI,4HN3  , 4HNOLI,4HN4  , 4HRLOA,4HD1  , 4HRLOA,4HD2  ,
     *     4HTABL,4HED1 , 4HTABL,4HED2 , 4HSEQE,4HP   , 4HTF  ,4H    ,
     7     4HTIC ,4H    , 4HTLOA,4HD1  , 4HTLOA,4HD2  , 4HTABL,4HED3 ,
     *     4HTABL,4HED4 , 4HTSTE,4HP   , 4HDSFA,4HCT  , 4HAXIC,4H    ,
     *     4HRING,4HAX  , 4HCCON,4HEAX , 4HPCON,4HEAX , 4HSPCA,4HX   ,
     9     4HMPCA,4HX   , 4HOMIT,4HAX  /
      DATA IBD4 /
     1     4HSUPA,4HX   , 4HPOIN,4HTAX , 4HSECT,4HAX  , 4HPRES,4HAX  ,
     *     4HTEMP,4HAX  , 4HFORC,4HEAX , 4HMOMA,4HX   , 4HEIGP,4H    ,
     *     4HPDUM,4H4   , 4HPDUM,4H5   , 4HPDUM,4H6   , 4HTABD,4HMP1 ,
     3     4HPDUM,4H7   , 4HPDUM,4H8   , 4HPDUM,4H9   , 4HFREQ,4H2   ,
     *     4HCONC,4HT1  , 4HCONC,4HT   , 4HTRAN,4HS   , 4HRELE,4HS   ,
     *     4HLOAD,4HC   , 4HSPCS,4HD   , 4HSPCS,4H1   , 4HSPCS,4H    ,
     5     4HBDYC,4H    , 4HMPCS,4H    , 4HBDYS,4H    , 4HBDYS,4H1   ,
     *     4HBARO,4HR   , 4HCBAR,4H    , 4HPBAR,4H    , 4HDARE,4HA   ,
     *     4HDELA,4HY   , 4HDPHA,4HSE  , 4HPLFA,4HCT  , 4HGNEW,4H    ,
     7     4HGTRA,4HN   , 4HTABR,4HNDG , 4HMATT,4H3   , 4HRFOR,4HCE  ,
     *     4HTABR,4HND1 , 4HPLOA,4HD4  , 4HUSET,4H    , 4HUSET,4H1   ,
     *     4HRAND,4HPS  , 4HRAND,4HT1  , 4HRAND,4HT2* , 4HPLOA,4HD1  ,
     9     4HPLOA,4HD2  , 4HDTI ,4H    /
      DATA IBD5 /
     1     4HTEMP,4HP1  , 4HTEMP,4HP2  , 4HTEMP,4HP3  , 4HTEMP,4HRB  ,
     *     4HGRID,4HB   , 4HFSLI,4HST  , 4HRING,4HFL  , 4HPRES,4HPT  ,
     *     4HCFLU,4HID2 , 4HCFLU,4HID3 , 4HCFLU,4HID4 , 4HAXIF,4H    ,
     3     4HBDYL,4HIST , 4HFREE,4HPT  , 4HASET,4H    , 4HASET,4H1   ,
     *     4HCTET,4HRA  , 4HCWED,4HGE  , 4HCHEX,4HA1  , 4HCHEX,4HA2  ,
     *     4HDMIA,4HX   , 4HFLSY,4HM   , 4HAXSL,4HOT  , 4HCAXI,4HF2  ,
     5     4HCAXI,4HF3  , 4HCAXI,4HF4  , 4HCSLO,4HT3  , 4HCSLO,4HT4  ,
     *     4HGRID,4HF   , 4HGRID,4HS   , 4HSLBD,4HY   , 4HCHBD,4HY   ,
     *     4HQHBD,4HY   , 4HMAT4,4H    , 4HMAT5,4H    , 4HPHBD,4HY   ,
     7     4HMATT,4H4   , 4HMATT,4H5   , 4HQBDY,4H1   , 4HQBDY,4H2   ,
     *     4HQVEC,4HT   , 4HQVOL,4H    , 4HRADL,4HST  , 4HRADM,4HTX  ,
     *     4HSAME,4H    , 4HNOSA,4HME  , 4HINPU,4HT   , 4HOUTP,4HUT  ,
     9     4HCQDM,4HEM1 , 4HPQDM,4HEM1 /
      DATA IBD6 /
     1     4HCIHE,4HX1  , 4HCIHE,4HX2  , 4HCIHE,4HX3  , 4HPIHE,4HX   ,
     *     4HPLOA,4HD3  , 4HSPCD,4H    , 4HCYJO,4HIN  , 4HCNGR,4HNT  ,
     *     4HCQDM,4HEM2 , 4HPQDM,4HEM2 , 4HCQUA,4HD4  , 4HMAT8,4H    ,
     3     4HCAER,4HO1  , 4HPAER,4HO1  , 4HAERO,4H    , 4HSPLI,4HNE1 ,
     *     4HSPLI,4HNE2 , 4HSET1,4H    , 4HSET2,4H    , 4HMKAE,4HRO2 ,
     *     4HMKAE,4HRO1 , 4HFLUT,4HTER , 4HAEFA,4HCT  , 4HFLFA,4HCT  ,
     5     4HCBAR,4HAO  , 4HPLIM,4HIT  , 4HPOPT,4H    , 4HPLOA,4HDX  ,
     *     4HCRIG,4HD1  , 4HPCOM,4HP   , 4HPCOM,4HP1  , 4HPCOM,4HP2  ,
     *     4HPSHE,4HLL  , 4HCRIG,4HD2  , 4HCTRI,4HAAX , 4HPTRI,4HAAX ,
     7     4HCTRA,4HPAX , 4HPTRA,4HPAX , 4HVIEW,4H    , 4HVARI,4HAN  ,
     *     4HCTRI,4HM6  , 4HPTRI,4HM6  , 4HCTRP,4HLT1 , 4HPTRP,4HLT1 ,
     *     4HTEMP,4HG   , 4HTEMP,4HP4  , 4HCRIG,4HDR  , 4HCRIG,4HD3  ,
     9     4HCTRS,4HHL  , 4HPTRS,4HHL  /
      DATA IBD7 /
     1     4HCAER,4HO2  , 4HCAER,4HO3  , 4HCAER,4HO4  , 4HPAER,4HO2  ,
     *     4HPAER,4HO3  , 4HPAER,4HO4  , 4HSPLI,4HNE3 , 4HGUST,4H    ,
     *     4HCAER,4HO5  , 4HPAER,4HO5  , 4HDARE,4HAS  , 4HDELA,4HYS  ,
     3     4HDPHA,4HSES , 4HTICS,4H    , 4HMATP,4HZ1  , 4HMATP,4HZ2  ,
     *     4HMTTP,4HZ1  , 4HMTTP,4HZ2  , 4HMAT6,4H    , 4HMATT,4H6   ,
     *     4HCEML,4HOOP , 4HSPCF,4HLD  , 4HCIS2,4HD8  , 4HPIS2,4HD8  ,
     5     4HGEML,4HOOP , 4HREMF,4HLUX , 4HBFIE,4HLD  , 4HMDIP,4HOLE ,
     *     4HPROL,4HATE , 4HPERM,4HBDY , 4HCFFR,4HEE  , 4HCFLS,4HTR  ,
     *     4HCFHE,4HX1  , 4HCFHE,4HX2  , 4HCFTE,4HTRA , 4HCFWE,4HDGE ,
     7     4HMATF,4H    , 4HCELB,4HOW  , 4HPELB,4HOW  , 4HNOLI,4HN5  ,
     *     4HNOLI,4HN6  , 4HCFTU,4HBE  , 4HPFTU,4HBE  , 4HNFTU,4HBE  ,
     *     4HSTRE,4HAML1, 4HSTRE,4HAML2, 4HCRRO,4HD   , 4HCRBA,4HR   ,
     9     4HCRTR,4HPLT , 4HCRBE,4H1   /
      DATA IBD8 /
     1     4HCRBE,4H2   , 4HCRBE,4H3   , 4HCRSP,4HLINE, 4HCTRI,4HA3  ,
     *     4HTABL,4HEM5 , 4HCPSE,4H2   , 4HCPSE,4H3   , 4HCPSE,4H4   ,
     *     4HPPSE,4H    , 82*4H****    /
C*****
C     THE IPR ARRAYS CONTAIN PAIRED ENTRIES OF BULK DATA PARAMETER NAMES
C*****
      DATA IPR1 /
     1     4HGRDP,4HNT  , 4HWTMA,4HSS  , 4HIRES,4H    , 4HLFRE,4HQ   ,
     5     4HHFRE,4HQ   , 4HLMOD,4HES  , 4HG   ,4H    , 4HW3  ,4H    ,
     9     4HW4  ,4H    , 4HMODA,4HCC  , 4HCOUP,4HMASS, 4HCPBA,4HR   ,
     3     4HCPRO,4HD   , 4HCPQU,4HAD1 , 4HCPQU,4HAD2 , 4HCPTR,4HIA1 ,
     7     4HCPTR,4HIA2 , 4HCPTU,4HBE  , 4HCPQD,4HPLT , 4HCPTR,4HPLT ,
     1     4HCPTR,4HBSC , 4HMAXI,4HT   , 4HEPSH,4HT   , 4HTABS,4H    ,
     5     4HSIGM,4HA   , 4HBETA,4H    , 4HRADL,4HIN  , 4HBETA,4HD   ,
     9     4HNT  ,4H    , 4HEPSI,4HO   , 4HCTYP,4HE   , 4HNSEQ,4HS   ,
     3     4HNLOA,4HD   , 4HCYCI,4HO   , 4HCYCS,4HEQ  , 4HKMAX,4H    ,
     7     4HKIND,4HEX  , 4HNODJ,4HE   , 4HP1  ,4H    , 4HP2  ,4H    ,
     1     4HP3  ,4H    , 4HVREF,4H    , 4HPRIN,4HT   , 4HISTA,4HRT  ,
     5     4HKDAM,4HP   , 4HGUST,4HAERO, 4HIFTM,4H    , 4HMACH,4H    ,
     9     4HQ   ,4H    , 4HHOPT,4H    /
      DATA IPR2 /
     1     4HGRDE,4HQ   , 4HSTRE,4HSS  , 4HSTRA,4HIN  , 4HNINT,4HPTS ,
     5     4HASET,4HOUT , 4HAUTO,4HSPC , 4HVOLU,4HME  , 4HSURF,4HACE ,
     9     4HKTOU,4HT   , 4HAPRE,4HSS  , 4HATEM,4HP   , 4HSTRE,4HAML ,
     3     4HPGEO,4HM   , 4HSIGN,4H    , 4HZORI,4HGN  , 4HFXCO,4HOR  ,
     7     4HFYCO,4HOR  , 4HFZCO,4HOR  , 4HKGGI,4HN   , 4HIREF,4H    ,
     1     4HMINM,4HACH , 4HMAXM,4HACH , 4HMTYP,4HE   , 4H****,4H****,
     5     44*4H****    /
C*****
C     THE ICC ARRAYS CONTAIN PAIRED ENTRIES OF CASE CONTROL FLAG NAMES
C     FOR USE IN RESTART RUNS
C*****
      DATA ICC1 /
     1     4HMPC$,4H    , 4HSPC$,4H    , 4HLOAD,4H$   , 4HMETH,4HOD$ ,
     5     4HDEFO,4HRM$ , 4HTEMP,4HLD$ , 4HTEMP,4HMT$ , 4HIC$ ,4H    ,
     9     4HAOUT,4H$   , 4HLOOP,4H$   , 4HLOOP,4H1$  , 4HDLOA,4HD$  ,
     3     4HFREQ,4H$   , 4HTF$ ,4H    , 4HPLOT,4H$   , 4HTSTE,4HP$  ,
     7     4HPOUT,4H$   , 4HTEMP,4HMX$ , 4HDSCO,4H$   , 4HK2PP,4H$   ,
     1     4HM2PP,4H$   , 4HB2PP,4H$   , 4HCMET,4HHOD$, 4HSDAM,4HP$  ,
     5     4HPLCO,4H$   , 4HNLFO,4HRCE$, 4HXYOU,4HT$  , 4HFMET,4HHOD$,
     9     4HRAND,4HOM$ , 4HAXYO,4HUT$ , 4HNOLO,4HOP$ , 4HGUST,4H$   ,
     3     4HQOUT,4H$   , 4HBOUT,4H$   ,
     7     32*4H****    /
      DATA ICC2 /
     1     24*4H****    /
C
C*****
C     INITIALIZATION OF VARIABLES IN COMMON /IFPX0/
C*****
C
C     THE VALUES ASSIGNED BELOW TO THE VARIABLES IN COMMON /IFPX0/
C     ARE AS PER THEIR DEFINITIONS GIVEN EARLIER IN THE COMMENTS
C     AND ARE DERIVED FROM THE COMMON /IFPX1/ INFORMATION
C*****
      DATA LBDPR, LCC, IWRDS, IPARPT /16, 2, 18*0, 401/
C
      END


text:
      subroutine gf_unpack2(cgrib,lcgrib,iofst,lencsec2,csec2,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    gf_unpack2
!   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2002-04-09
!
! ABSTRACT: This subroutine unpacks Section 2 (Local Use Section)
!           as defined in GRIB Edition 2.
!
! PROGRAM HISTORY LOG:
! 2002-04-09  Gilbert
!
! USAGE:    CALL gf_unpack2(cgrib,lcgrib,iofst,lencsec2,csec2,ierr)
!   INPUT ARGUMENT LIST:
!     cgrib    - Character array containing Section 2 of the GRIB2 message
!     lcgrib   - Length (in bytes) of GRIB message array cgrib.
!     iofst    - Bit offset of the beginning of Section 2.
!
!   OUTPUT ARGUMENT LIST:
!     iofst    - Bit offset at the end of Section 2, returned.
!     lencsec2 - Length (in octets) of Local Use data
!     csec2()  - Pointer to a character*1 array containing local use data
!     ierr     - Error return code.
!                0 = no error
!                2 = Array passed is not section 2
!                6 = memory allocation error
!
! REMARKS: None
!
! ATTRIBUTES:
!   LANGUAGE: Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      character(len=1),intent(in) :: cgrib(lcgrib)
      integer,intent(in) :: lcgrib
      integer,intent(inout) :: iofst
      integer,intent(out) :: lencsec2
      integer,intent(out) :: ierr
      character(len=1),pointer,dimension(:) :: csec2

      ierr=0
      lencsec2=0
      nullify(csec2)

      call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
      iofst=iofst+32
      lencsec2=lensec-5
      call gbyte(cgrib,isecnum,iofst,8)         ! Get Section Number
      iofst=iofst+8
      ipos=(iofst/8)+1

      if ( isecnum.ne.2 ) then
         ierr=6
         print *,'gf_unpack2: Not Section 2 data. '
         return
      endif

      allocate(csec2(lencsec2),stat=istat)
      if (istat.ne.0) then
         ierr=6
         nullify(csec2)
         return
      endif

      csec2(1:lencsec2)=cgrib(ipos:ipos+lencsec2-1)
      iofst=iofst+(lencsec2*8)

      return    ! End of Section 2 processing
      end

max_stars_repo_name: martinremy/wps

max_issues_repo_name: martinremy/wps

max_forks_repo_name: WRF-CMake/WPS


text:
C$Procedure      INCNSG ( Intersection of cone and line segment )

      SUBROUTINE INCNSG ( APEX,   AXIS,  ANGLE, ENDPT1,
     .                    ENDPT2, NXPTS, XPT1,  XPT2   )

C$ Abstract
C
C     Compute the points of intersection of a specified nappe of a cone
C     and a line segment.
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     CONE
C     GEOMETRY
C     INTERSECTION
C     LINE
C     MATH
C     SEGMENT
C
C$ Declarations

      IMPLICIT NONE

      DOUBLE PRECISION      APEX   ( 3 )
      DOUBLE PRECISION      AXIS   ( 3 )
      DOUBLE PRECISION      ANGLE
      DOUBLE PRECISION      ENDPT1 ( 3 )
      DOUBLE PRECISION      ENDPT2 ( 3 )
      INTEGER               NXPTS
      DOUBLE PRECISION      XPT1   ( 3 )
      DOUBLE PRECISION      XPT2   ( 3 )

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     APEX       I   Apex of cone.
C     AXIS       I   Axis of cone.
C     ANGLE      I   Angle of cone.
C     ENDPT1,
C     ENDPT2     I   Endpoints of line segment.
C     NXPTS      O   Number of intersection points.
C     XPT1       O   First intersection point, if it exists.
C     XPT2       O   Second intersection point, if it exists.
C
C$ Detailed_Input
C
C     APEX       is the apex (tip) of the cone. In this routine's
C                documentation, we'll consider the cone to be a
C                semi-infinite pyramid with circular cross-section. In
C                some contexts, this object is called one "nappe" of
C                the complete cone.
C
C     AXIS       is an axis vector of the cone.
C
C     ANGLE      is the angular separation from AXIS of the rays
C                comprising the cone. Let the notation
C
C                   < A, B >
C
C                denote the dot product of vectors A and B, and let
C
C                   ||A||
C
C                denote the norm of vector A. Then the cone is the set
C                of points
C
C                             X-APEX       AXIS
C                   { X:  < ----------,  -------- >  =  cos(ANGLE) }
C                           ||X-APEX||   ||AXIS||
C
C
C     ENDPT1,
C     ENDPT2     are endpoints of a line segment. These points
C                must be distinct.
C
C$ Detailed_Output
C
C     NXPTS      is the number of points of intersection of the input
C                line segment and cone.
C
C     XPT1       is the point of intersection of the segment and cone
C                that is closest to ENDPT1, if an intersection exists.
C                If there are no intersections, XPT1 is undefined.
C
C     XPT2       is the point of intersection of the segment and cone
C                that is farthest from ENDPT1, if two points of
C                intersection exist. If there are not two
C                intersections, XPT2 is undefined.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If AXIS is the zero vector, the error SPICE(ZEROVECTOR)
C         will be signaled.
C
C     2)  If ANGLE is less than zero, the error SPICE(INVALIDANGLE)
C         will be signaled.
C
C     3)  If ENDPT1 and ENDPT2 coincide, the error
C         SPICE(ENDPOINTSMATCH) will be signaled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine is used by the SPICELIB DSK subsystem. In
C     particular, it is used to determine whether a ray contacts a
C     latitude boundary of a volume element in either planetocentric
C     latitudinal or planetodetic coordinates.
C
C$ Examples
C
C     The numerical results shown for these examples may differ across
C     platforms. The results depend on the SPICE kernels used as input
C     (if any), the compiler and supporting libraries, and the machine
C     specific arithmetic implementation.
C
C     1) Compute the intersection of a line segment and cone in
C        a simple case for which the results can easily be checked.
C
C        Let the apex of the cone be at the origin. Let the axis
C        of the cone lie on the +X axis. Let the angle of the cone
C        be 45 degrees. Let the line segment have endpoints
C
C           ENDPT1 = ( 1,   -2, sqrt(3)/2 )
C           ENDPT2 = ( 1,    2, sqrt(3)/2 )
C
C        We expect there to be two points of intersection:
C
C           XPT1   = ( 1, -1/2, sqrt(3)/2 )
C           XPT2   = ( 1,  1/2, sqrt(3)/2 )
C
C
C        Example code begins here.
C
C
C              PROGRAM EX1
C              IMPLICIT NONE
C        C
C        C     SPICELIB functions
C        C
C              DOUBLE PRECISION      RPD
C        C
C        C     Local parameters
C        C
C              CHARACTER*(*)         FMT1
C              PARAMETER           ( FMT1 = '(A,3F13.8)' )
C
C              CHARACTER*(*)         FMT2
C              PARAMETER           ( FMT2 = '(A,I2)' )
C        C
C        C     Local variables
C        C
C              DOUBLE PRECISION      ANGLE
C              DOUBLE PRECISION      APEX   ( 3 )
C              DOUBLE PRECISION      AXIS   ( 3 )
C              DOUBLE PRECISION      ENDPT1 ( 3 )
C              DOUBLE PRECISION      ENDPT2 ( 3 )
C              DOUBLE PRECISION      SQ3
C              DOUBLE PRECISION      XPT1   ( 3 )
C              DOUBLE PRECISION      XPT2   ( 3 )
C
C              INTEGER               NXPTS
C
C        C
C        C     Set up the cone's geometric attributes.
C        C
C              CALL VPACK ( 0.D0, 0.D0, 0.D0, APEX )
C              CALL VPACK ( 1.D0, 0.D0, 0.D0, AXIS )
C
C              ANGLE = 45.D0 * RPD()
C        C
C        C     Initialize the line segment's endpoints.
C        C
C              SQ3 = SQRT( 3.D0  )
C
C              CALL VPACK ( 1.D0, -2.D0, SQ3/2, ENDPT1 )
C              CALL VPACK ( 1.D0,  2.D0, SQ3/2, ENDPT2 )
C        C
C        C     Find the points of intersection.
C        C
C              CALL INCNSG ( APEX,   AXIS,  ANGLE, ENDPT1,
C             .              ENDPT2, NXPTS, XPT1,  XPT2   )
C
C              WRITE (*,*) ' '
C              WRITE (*,FMT1) 'Apex:        ', APEX
C              WRITE (*,FMT1) 'Axis:        ', AXIS
C              WRITE (*,FMT1) 'Angle (deg): ', ANGLE/RPD()
C              WRITE (*,FMT1) 'Endpoint 1:  ', ENDPT1
C              WRITE (*,FMT1) 'Endpoint 2:  ', ENDPT2
C              WRITE (*,*) ' '
C              WRITE (*,FMT2) 'Number of intersection points: ',
C             .            NXPTS
C              WRITE (*,*) ' '
C              WRITE (*,FMT1) 'Point 1:    ', XPT1
C              WRITE (*,FMT1) 'Point 2:    ', XPT2
C              WRITE (*,*) ' '
C
C              END
C
C
C     When this program was executed on a PC/Linux/gfortran/64-bit
C     platform, the output was:
C
C
C        Apex:           0.00000000   0.00000000   0.00000000
C        Axis:           1.00000000   0.00000000   0.00000000
C        Angle (deg):   45.00000000
C        Endpoint 1:     1.00000000  -2.00000000   0.86602540
C        Endpoint 2:     1.00000000   2.00000000   0.86602540
C
C        Number of intersection points:  2
C
C        Point 1:       1.00000000  -0.50000000   0.86602540
C        Point 2:       1.00000000   0.50000000   0.86602540
C
C
C$ Restrictions
C
C     1)  This routine is designed to avoid arithmetic overflow in
C         normal cases, such as those in which the line segment is
C         nearly parallel to the cone. However, it is possible to cause
C         arithmetic overflow by using input vectors with extremely
C         large magnitudes.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0 26-OCT-2016 (NJB)
C
C-&

C$ Index_Entries
C
C     intersection of line segment and cone
C
C-&

C
C     SPICELIB functions
C
      DOUBLE PRECISION      HALFPI
      DOUBLE PRECISION      PI
      DOUBLE PRECISION      VDIST
      DOUBLE PRECISION      VDOT

      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
      DOUBLE PRECISION      COSTOL
      PARAMETER           ( COSTOL = 1.D-10 )

      DOUBLE PRECISION      HSPTOL
      PARAMETER           ( HSPTOL = 1.D-14 )

      INTEGER               UBPL
      PARAMETER           ( UBPL   = 4 )

C
C     Local variables
C
      DOUBLE PRECISION      A
      DOUBLE PRECISION      AXMAG
      DOUBLE PRECISION      B
      DOUBLE PRECISION      C
      DOUBLE PRECISION      CA2
      DOUBLE PRECISION      COLAT
      DOUBLE PRECISION      COSANG
      DOUBLE PRECISION      COSERR
      DOUBLE PRECISION      DIR    ( 3 )
      DOUBLE PRECISION      DMAG
      DOUBLE PRECISION      DP1
      DOUBLE PRECISION      DP2
      DOUBLE PRECISION      LOCANG
      DOUBLE PRECISION      MAXLAT
      DOUBLE PRECISION      MAXP   ( 3 )
      DOUBLE PRECISION      MINLAT
      DOUBLE PRECISION      MINP   ( 3 )
      DOUBLE PRECISION      NRMPLN ( UBPL )
      DOUBLE PRECISION      OFF1   ( 3 )
      DOUBLE PRECISION      OFF2   ( 3 )
      DOUBLE PRECISION      ORIGIN ( 3 )
      DOUBLE PRECISION      PLNX   ( 3 )
      DOUBLE PRECISION      S1
      DOUBLE PRECISION      S2
      DOUBLE PRECISION      UAXIS  ( 3 )
      DOUBLE PRECISION      UDIR   ( 3 )
      DOUBLE PRECISION      UOFF1  ( 3 )
      DOUBLE PRECISION      UOFF2  ( 3 )
      DOUBLE PRECISION      UUAX
      DOUBLE PRECISION      UV1    ( 3 )
      DOUBLE PRECISION      UV2    ( 3 )
      DOUBLE PRECISION      V1     ( 3 )
      DOUBLE PRECISION      V1MAG
      DOUBLE PRECISION      V2     ( 3 )
      DOUBLE PRECISION      V2MAG
      DOUBLE PRECISION      VTEMP  ( 3 )
      DOUBLE PRECISION      VTEMP2 ( 3 )
      DOUBLE PRECISION      W2
      DOUBLE PRECISION      WU
      DOUBLE PRECISION      WUAX
      DOUBLE PRECISION      X      ( 3 )
      DOUBLE PRECISION      XOFF1  ( 3 )
      DOUBLE PRECISION      XOFF2  ( 3 )
      DOUBLE PRECISION      XFORM  ( 3, 3 )
      DOUBLE PRECISION      Y      ( 3 )
      DOUBLE PRECISION      Z      ( 3 )

      INTEGER               I
      INTEGER               N
      INTEGER               NPLNX

      LOGICAL               IN1
      LOGICAL               IN2
      LOGICAL               ISBRCK
      LOGICAL               NEG1
      LOGICAL               NEG2

C
C     Saved values
C
      SAVE                  ORIGIN
      SAVE                  Y
      SAVE                  Z

C
C     Initial values
C
      DATA                  ORIGIN / 3 * 0.D0         /
      DATA                  Y      / 0.D0, 1.D0, 0.D0 /
      DATA                  Z      / 0.D0, 0.D0, 1.D0 /

C
C     Use quasi-discovery check-in. We'll check in before
C     code sections that can generate SPICE errors, and check
C     out afterward. When those code sections are skipped,
C     we avoid traceback participation.
C
      IF ( RETURN() ) THEN
         RETURN
      END IF

C
C     No intersection was found so far.
C
      NXPTS = 0

C
C     The cone's axis vector must be non-zero.
C
      CALL UNORM ( AXIS, UAXIS, AXMAG )

      IF ( AXMAG .EQ. 0.D0 ) THEN

         CALL CHKIN  ( 'INCNSG'                                       )
         CALL SETMSG ( 'The cone''s axis vector must be non-zero '
     .   //            'but sadly, it failed to meet this criterion.' )
         CALL SIGERR ( 'SPICE(ZEROVECTOR)'                            )
         CALL CHKOUT ( 'INCNSG'                                       )
         RETURN

      END IF

C
C     The cone's angular radius must be non-negative.
C
      IF ( ANGLE .LT. 0.D0 ) THEN

         CALL CHKIN  ( 'INCNSG'                             )
         CALL SETMSG ( 'The cone''s angular radius must be '
     .   //            ' non-negative but was # (radians).' )
         CALL ERRDP  ( '#', ANGLE                           )
         CALL SIGERR ( 'SPICE(INVALIDANGLE)'                )
         CALL CHKOUT ( 'INCNSG'                             )
         RETURN

      END IF

C
C     The endpoints of the segment must be distinct. Check this after
C     computing a unit direction vector for the line segment.
C
      CALL VSUB  ( ENDPT2, ENDPT1, DIR )

      CALL UNORM ( DIR, UDIR, DMAG )

      IF ( DMAG .EQ. 0.D0 ) THEN

         CALL CHKIN  ( 'INCNSG'                              )
         CALL SETMSG ( 'The distance between the segment''s '
     .   //            'endpoints was zero. First endpoint: '
     .   //            '(# # #).'                            )
         CALL ERRDP  ( '#', ENDPT1(1)                        )
         CALL ERRDP  ( '#', ENDPT1(2)                        )
         CALL ERRDP  ( '#', ENDPT1(3)                        )
         CALL SIGERR ( 'SPICE(ENDPOINTSMATCH)'               )
         CALL CHKOUT ( 'INCNSG'                              )
         RETURN

      END IF

C
C     Store the cosine of the cone's angular radius. We'll treat all
C     cases with COSANG equal to 0 as though the cone is actually a
C     plane normal to the axis and containing the apex.
C
      COSANG = COS( ANGLE )
      LOCANG = ANGLE
C
C     We'll work with a local axis that has angular separation of
C     no more than pi/2 from the nappe.
C
      IF ( COSANG .LT. 0.D0 ) THEN

         COSANG   = - COSANG
         LOCANG   =   PI() - ANGLE

         UAXIS(1) = - UAXIS(1)
         UAXIS(2) = - UAXIS(2)
         UAXIS(3) = - UAXIS(3)

      END IF

C
C     Compute the offsets of the endpoints of the segment from
C     the cone's apex.
C
      CALL VSUB ( ENDPT1, APEX, OFF1 )
      CALL VSUB ( ENDPT2, APEX, OFF2 )

C
C     Deal with some of the simple cases first.
C
      CALL VHAT ( OFF1, UOFF1 )
      CALL VHAT ( OFF2, UOFF2 )

      DP1  = VDOT( UOFF1, UAXIS )
      DP2  = VDOT( UOFF2, UAXIS )

C
C     The given axis is inside the nappe defined by the angular radius.
C
C     There's no intersection if both endpoints are in the interior of
C     the nappe of the cone (since the nappe is convex).
C
      IN1  = DP1 .GE. COSANG
      IN2  = DP2 .GE. COSANG
C
C     If the line segment lies on the far side of the plane that
C     contains the apex and is orthogonal to the axis, there's no
C     intersection.
C
      NEG1 = DP1 .LT. 0.D0
      NEG2 = DP2 .LT. 0.D0

      IF (      ( IN1  .AND. IN2  )
     .     .OR. ( NEG1 .AND. NEG2 ) ) THEN
C
C        The segment is in the interior of the cone or
C        on the far side of the plane.
C
         NXPTS = 0

         RETURN

      END IF

C
C     Here's where we handle the half-space case.
C
      IF ( ABS(COSANG) .LT. HSPTOL ) THEN
C
C        See whether the ray emanating from the first endpoint and
C        having direction UDIR hits the plane normal to the axis and
C        containing the apex. We'll call this plane NRMPLN.
C
C        NVP2PL can signal an error only if the input axis is the
C        zero vector. We've ensured that it isn't.
C
         CALL NVP2PL ( UAXIS,  APEX, NRMPLN )
         CALL INRYPL ( ENDPT1, UDIR, NRMPLN, NPLNX, PLNX )

C
C        If the ray doesn't hit the plane, we're done. Otherwise,
C        check the intercept.
C
         IF ( NPLNX .EQ. 1 ) THEN
C
C           The ray does hit the plane. If the intersection is on the
C           line segment, we have a solution.
C
            IF ( VDIST(PLNX, ENDPT1) .LE. DMAG ) THEN
C
C              The intercept is not further along the ray than the
C              second endpoint. It's a valid solution.
C
               NXPTS = 1
               CALL VEQU ( PLNX, XPT1 )

            END IF

         END IF
C
C        This is the end of the half-space case.
C
         RETURN

      END IF

C
C     At this point we've disposed of the trivial cases. We'll
C     set up a quadratic equation for the intersection of the
C     line segment with the surface of the cone's nappe.
C
C     Due to round-off errors, the solution of the quadratic may
C     either be inaccurate or may not be found at all. We'll
C     examine the solutions we find and solve the problem by
C     an alternate method if necessary. However, the quadratic
C     method is fast, so we give it priority.
C
C     The equation of a ray starting at ENDPT1 and having unit
C     direction vector UDIR is
C
C        RAY  = { ENDPT1 + s*UDIR, s >= 0 }                          (1)
C
C     The equation of the nappe of the cone is
C
C        CONE = { X: < X - APEX, UAXIS > = ||X-APEX|| * cos(ANGLE) } (2)
C
C     where ANGLE is the angular radius of the cone and UAXIS is the
C     unit axis vector. Substituting the right hand side expression of
C     (1) for X in equation (2) and squaring both sides yields a
C     quadratic equation for S. We'll derive the coefficients of the
C     equation below.
C
C     Let
C
C        Q  = X - APEX
C        W  = ENDPT1 - APEX
C        U  = UDIR
C        CA = cos(ANGLE)
C
C     We can translate the cone and ray by -APEX, and (1) and (2)
C     can be re-written as
C
C        RAY  = { W + s*U, s >= 0 }                                  (3)
C
C        CONE = { Q: < Q, UAXIS > = ||Q|| * cos(ANGLE) }             (4)
C
C
C        Substituting the ray expression for Q, we obtain
C
C           < W + s*U, UAXIS > = ||W+s*U|| * CA                      (5)
C
C        and squaring both sides yields
C
C                      2                                      2   2
C             <W,UAXIS>  + 2*<W,UAXIS>*<U,UAXIS>*s + <U,UAXIS> * s
C
C                    2                2       2
C           = ( ||W||  + 2*<W,U>*s + s  ) * CA                       (6)
C
C
C       Collecting coefficients of powers of s, we have
C
C                         2     2     2
C              ( <U,UAXIS>  - CA ) * s
C
C                                            2
C            + 2 * ( <W,UAXIS>*<U,UAXIS> - CA * <W,U> ) * s
C
C                       2        2    2
C            + <W,UAXIS>  - ||W|| * CA
C
C
C         =  0                                                       (7)
C
C
C      Before continuing, we observe that the only non-unit vector
C      in (7) is W. So the coefficients in (7) have no possibility
C      of overflowing unless the vertex of the ray is very far from
C      the apex of the cone.
C
C      W has been computed above as OFF1.
C
C
C         [ Consider adding check on OFF1 here. ]
C
C
C     Intermediate values:
C
      UUAX = VDOT( UDIR, UAXIS )
      WUAX = VDOT( OFF1, UAXIS )
      WU   = VDOT( OFF1, UDIR  )
      W2   = VDOT( OFF1, OFF1  )
      CA2  = COSANG * COSANG

C
C     Quadratic coefficients:
C
      A    = ( UUAX * UUAX ) - CA2

      B    = 2 * (  ( WUAX * UUAX ) - ( CA2 * WU )  )

      C    = ( WUAX * WUAX ) - ( W2 * CA2 )

C
C     We're not interested in solutions that lie outside
C     of the line segment. The length of the segment is
C     DMAG.
C
C     Solve the equation, using DMAG as an upper bound
C     on the magnitude of the roots.
C
      CALL ZZCNQUAD ( A, B, C, DMAG, N, S1, S2 )

C
C     Compute the possible intersection points and test them
C     to make sure they really are solutions.
C
      IF ( N .GT. 0 ) THEN
C
C        Start with the solution closest to the ray's vertex.
C        Compute XPT1 and make sure it's on the correct nappe
C        of the cone.
C
         IF ( S1 .GE. 0.D0 ) THEN

            XPT1(1) = ENDPT1(1) + S1*UDIR(1)
            XPT1(2) = ENDPT1(2) + S1*UDIR(2)
            XPT1(3) = ENDPT1(3) + S1*UDIR(3)

            CALL VSUB  ( XPT1, APEX, V1 )

C
C           See whether V1 is on the cone.
C
            CALL UNORM ( V1, UV1, V1MAG )

            IF ( V1MAG .GT. 0.D0 ) THEN

               COSERR = ABS( VDOT(UV1,UAXIS) - COSANG )
            ELSE
               COSERR = 0.D0
            END IF

            IF (      ( V1MAG  .EQ. 0.D0   )
     .           .OR. ( COSERR .LT. COSTOL ) ) THEN
C
C              The root is on the cone (on the apex if V1MAG is zero).
C
C              We accept this root. Update NXPTS. Note that this is
C              not necessarily the final value of NXPTS; that
C              depends on the validity of the second root.

               NXPTS = 1

            END IF

         END IF

         IF ( N .EQ. 2 ) THEN
C
C           Check the second root.
C
            IF ( S2 .GE. 0.D0 ) THEN

               XPT2(1) = ENDPT1(1) + S2*UDIR(1)
               XPT2(2) = ENDPT1(2) + S2*UDIR(2)
               XPT2(3) = ENDPT1(3) + S2*UDIR(3)

               CALL VSUB  ( XPT2, APEX, V2 )

C
C              See whether V2 is on the cone.
C
               CALL UNORM ( V2, UV2, V2MAG )

               IF ( V2MAG .GT. 0.D0 ) THEN

                  COSERR = ABS( VDOT(UV2,UAXIS) - COSANG )
               ELSE
                  COSERR = 0.D0
               END IF

               IF (      ( V2MAG  .EQ. 0.D0   )
     .              .OR. ( COSERR .LT. COSTOL ) ) THEN
C
C                 The root is on the cone (on the apex if V2MAG is
C                 zero).
C
C                 We accept this root.
C
                  NXPTS = NXPTS + 1

                  IF ( NXPTS .EQ. 1 ) THEN
C
C                    This is the only valid root; overwrite XPT1.
C
                     CALL VEQU ( XPT2, XPT1 )

                  END IF

               END IF

            END IF

         END IF

      END IF

C
C     We're not done yet. If we have fewer roots than we should, we'll
C     need to solve the problem by an alternate method.
C
C     If we have two roots, we're in good shape. Otherwise we must
C     determine how many roots should be found.
C
      IF ( NXPTS .LT. 2 ) THEN
C
C        We must determine the expected number of roots, and if
C        we didn't come up with them, we must find the roots
C        by an alternate method.
C
C        We'll examine the containment of the endpoints within the
C        cone.
C
C        The case where both endpoints are inside the cone was handled
C        earlier.
C
C        If one endpoint is inside the cone and one is outside,
C        we expect to have one root.
C
         IF (     ( IN1 .AND. ( .NOT. IN2 ) )
     .       .OR. ( IN2 .AND. ( .NOT. IN1 ) )  ) THEN
C
C           There's supposed to be one root. If we found none, find one
C           now.
C
            IF ( NXPTS .EQ. 0 ) THEN
C
C              ZZCXBRUT signals an error if the axis is the zero
C              vector, but not otherwise. We've already ruled out this
C              situation. Therefore, we don't check in before the
C              following call.
C
               CALL ZZCXBRUT ( APEX,   UAXIS,  LOCANG,
     .                         ENDPT1, ENDPT2, XPT1,  ISBRCK )

               IF ( ISBRCK ) THEN
C
C                 As long as the root was bracketed, XPT1 is a
C                 solution.
C
                  NXPTS = 1

               END IF

            END IF

         ELSE

            CALL CHKIN ( 'INCNSG' )
C
C           Both endpoints are outside the cone. We could have zero to
C           two roots. If the minimum angular separation of the segment
C           from the axis is less than ANGLE, we expect to find two
C           roots; if it's equal to ANGLE, we expect to find one, and
C           if it's greater than ANGLE, none.
C
C           We'll transform OFF1 and OFF2 into a reference frame in
C           which angular separation from the axis is equivalent to
C           colatitude. Then we'll find the maximum latitude attained
C           on the segment.
C
C           We'll count the roots we find, so we'll start at zero.
C
            NXPTS = 0

            CALL FRAME ( UAXIS, X, Y )

            DO I = 1, 3

               XFORM(1,I) = X(I)
               XFORM(2,I) = Y(I)
               XFORM(3,I) = UAXIS(I)

            END DO

            CALL MXV ( XFORM, OFF1, XOFF1 )
            CALL MXV ( XFORM, OFF2, XOFF2 )

            CALL ZZSGLATX ( XOFF1, XOFF2, MINLAT, MINP, MAXLAT, MAXP )

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

C
C           COLAT is the colatitude of the point of maximum latitude.
C
            COLAT = HALFPI() - MAXLAT

            IF ( COLAT .LT. LOCANG ) THEN
C
C              MAXP is inside the cone. There should be an intersection
C              on the segment between XOFF1 and MAXP and another
C              between MAXP and XOFF2.
C
               CALL ZZCXBRUT ( ORIGIN, Z,    LOCANG,
     .                         XOFF1,  MAXP, VTEMP,  ISBRCK )

               IF ( ISBRCK ) THEN
C
C                 Convert VTEMP to the original frame, then translate
C                 it so that it's represented as an offset from the
C                 origin.
C
                  CALL MTXV ( XFORM,  VTEMP, VTEMP2 )
                  CALL VADD ( VTEMP2, APEX,  XPT1   )

                  NXPTS = 1

               END IF

               CALL ZZCXBRUT ( ORIGIN, Z,     LOCANG,
     .                         MAXP,   XOFF2, VTEMP,  ISBRCK )

               IF ( ISBRCK ) THEN
C
C                 Convert VTEMP to the original frame, then translate
C                 it so that it's represented as an offset from the
C                 origin.
C
                  CALL MTXV ( XFORM,  VTEMP, VTEMP2 )
                  CALL VADD ( VTEMP2, APEX,  XPT2   )

                  IF ( NXPTS .EQ. 1 ) THEN
C
C                    Both roots are valid.
C
                     NXPTS = 2

                  ELSE
C
C                    The second root is the only valid root. Move it
C                    into XPT1.
C
                     CALL VEQU ( XPT2, XPT1 )

                     NXPTS = 1

                  END IF

               END IF

            ELSE IF ( COLAT .EQ. LOCANG ) THEN
C
C              The root corresponds to a point of tangency of
C              the segment and cone. This occurs at the point
C              having maximum latitude: MAXP.
C
               CALL VEQU ( MAXP, XPT1 )

               NXPTS = 1

C
C           Note that if COLAT > LOCANG, there are no roots.
C
            END IF

            CALL CHKOUT ( 'INCNSG' )

         END IF
C
C        This is the end of portion of the "brute force" branch in
C        which both endpoints are outside the cone.
C
      END IF

C
C     NXPTS  has been set.
C
C     If NXPTS is 1, then XPT1 is set.
C
C     If NXPTS is 2, then both XPT1 and XPT2 are set.
C

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework