text:
      double precision function Ampvirtsq_AQaq_nonid(p1,p2,p3,p4)
c--- this routine is a wrapper to the 1-loop/Born interference
c--- of D.&S. 0906.0008; it performs the interference and sum over
c--- helicities and returns the ME in the same fashion as Hqarbvsq
      implicit none
      include 'constants.f'
      include 'zprods_com.f'
      integer p1,p2,p3,p4,j1,j2,j3,j4,h1,h2
      double precision sum,ren,H4prenorm
      double complex amplo(2,2),ampv(2,2),
     . A0Hqarbmppm,A0Hqarbmpmp,A41Hqarbmppm,A41Hqarbmpmp

c--- the routine expects momenta to be labelled:
c---                                 a(p1)+q(p2) --> Q(p3)+A(p4)
c---   and in this routine we use:   0 --> A(j1)+Q(j2)+a(j3)+q(j4)
c---
c--- so we permute appropriately (NOTE: different from LO routine)
      j1=p4
      j2=p3
      j3=p2
      j4=p1

c--- Born terms
c--- coded amplitudes
      amplo(1,2)=A0Hqarbmppm(j1,j2,j3,j4,za,zb)
      amplo(1,1)=A0Hqarbmpmp(j1,j2,j3,j4,za,zb)
c--- charge conjugation
c      amplo(2,1)=A0Hqarbmppm(j2,j1,j4,j3,za,zb)
c      amplo(2,2)=A0Hqarbmpmp(j2,j1,j4,j3,za,zb)
c--- equivalent to complex conjugation for Born amplitudes
      amplo(2,1)=dconjg(amplo(1,2))
      amplo(2,2)=dconjg(amplo(1,1))

c--- 1-loop terms
c--- coded amplitudes
      ampv(1,2) =A41Hqarbmppm(j1,j2,j3,j4,za,zb)
      ampv(1,1) =A41Hqarbmpmp(j1,j2,j3,j4,za,zb)
c--- charge conjugation
      ampv(2,1) =A41Hqarbmppm(j2,j1,j4,j3,za,zb)
      ampv(2,2) =A41Hqarbmpmp(j2,j1,j4,j3,za,zb)

c--- get renormalization factor
      ren=H4prenorm()

c--- square, with appropriate overall factor (in D.&S. notation,
c---   overall 4*C^2*gs^4*ason2pi is applied in gg_hgg_v.f)
      sum=zip
      do h1=1,2
      do h2=1,2
      ampv(h1,h2)=ampv(h1,h2)+ren*amplo(h1,h2)

      sum=sum+dble(ampv(h1,h2)*dconjg(amplo(h1,h2)))
      enddo
      enddo

      Ampvirtsq_AQaq_nonid=sum*V*xn/4d0

      return
      end

      double precision function Ampvirtsq_AQaq_ident(p1,p2,p3,p4)
c--- this routine is a wrapper to the 1-loop/Born interference
c--- of D.&S. 0906.0008; it performs the interference and sum over
c--- helicities and returns the ME in the same fashion as Hqarbvsq
      implicit none
      include 'constants.f'
      include 'zprods_com.f'
      integer p1,p2,p3,p4,j1,j2,j3,j4,h1,h2
      double precision sum,ren,H4prenorm
      double complex amplo_a(2,2),amplo_b(2,2),
     . amp42_a(2,2),amp42_b(2,2),
     . A0Hqarbmpmp,A42Hqarbmpmp

c--- the routine expects momenta to be labelled:
c---                                 a(p1)+q(p2) --> Q(p3)+A(p4)
c---   and in this routine we use:   0 --> A(j1)+Q(j2)+a(j3)+q(j4)
c---
c--- so we permute appropriately (NOTE: different from LO routine)
      j1=p4
      j2=p3
      j3=p2
      j4=p1

c--- Note that the A42 amplitude only contributes for identical
c--- quarks and, therefore, we only need the case of equal helicities

c--- First calculate with the normal ordering, (1,2,3,4)
c--- Born terms
c--- coded amplitudes
c      amplo_a(1,2)=A0Hqarbmppm(j1,j2,j3,j4,za,zb)
      amplo_a(1,1)=A0Hqarbmpmp(j1,j2,j3,j4,za,zb)
c--- charge conjugation
c      amplo_a(2,1)=A0Hqarbmppm(j2,j1,j4,j3,za,zb)
c      amplo_a(2,2)=A0Hqarbmpmp(j2,j1,j4,j3,za,zb)
c--- equivalent to complex conjugation for Born amplitudes
      amplo_a(2,2)=dconjg(amplo_a(1,1))

c--- 1-loop terms
c--- coded amplitudes
c      amp41_a(1,2) =A41Hqarbmppm(j1,j2,j3,j4,za,zb)
c      amp41_a(1,1) =A41Hqarbmpmp(j1,j2,j3,j4,za,zb)
c--- charge conjugation
c      amp41_a(2,1) =A41Hqarbmppm(j2,j1,j4,j3,za,zb)
c      amp41_a(2,2) =A41Hqarbmpmp(j2,j1,j4,j3,za,zb)

c--- coded amplitudes
c      amp42_a(1,2) =A42Hqarbmppm(j1,j2,j3,j4,za,zb)
      amp42_a(1,1) =A42Hqarbmpmp(j1,j2,j3,j4,za,zb)
c--- charge conjugation
c      amp42_a(2,1) =A42Hqarbmppm(j2,j1,j4,j3,za,zb)
      amp42_a(2,2) =A42Hqarbmpmp(j2,j1,j4,j3,za,zb)

c--- Now calculate with the exhanged order,    (1,4,3,2)
c--- Born terms
c--- coded amplitudes
c      amplo_b(1,2)=A0Hqarbmppm(j1,j4,j3,j2,za,zb)
      amplo_b(1,1)=A0Hqarbmpmp(j1,j4,j3,j2,za,zb)
c--- charge conjugation
c      amplo_b(2,1)=A0Hqarbmppm(j4,j1,j2,j3,za,zb)
c      amplo_b(2,2)=A0Hqarbmpmp(j4,j1,j2,j3,za,zb)
c--- equivalent to complex conjugation for Born amplitudes
      amplo_b(2,2)=dconjg(amplo_b(1,1))

c--- 1-loop terms
c--- coded amplitudes
c      amp41_b(1,2) =A41Hqarbmppm(j1,j4,j3,j2,za,zb)
c      amp41_b(1,1) =A41Hqarbmpmp(j1,j4,j3,j2,za,zb)
c--- charge conjugation
c      amp41_b(2,1) =A41Hqarbmppm(j4,j1,j2,j3,za,zb)
c      amp41_b(2,2) =A41Hqarbmpmp(j4,j1,j2,j3,za,zb)

c--- coded amplitudes
c      amp42_b(1,2) =A42Hqarbmppm(j1,j4,j3,j2,za,zb)
      amp42_b(1,1) =A42Hqarbmpmp(j1,j4,j3,j2,za,zb)
c--- charge conjugation
c      amp42_b(2,1) =A42Hqarbmppm(j4,j1,j2,j3,za,zb)
      amp42_b(2,2) =A42Hqarbmpmp(j4,j1,j2,j3,za,zb)

c--- get renormalization factor
      ren=H4prenorm()

c--- square, with appropriate overall factor (in D.&S. notation,
c---   overall 4*C^2*gs^4*ason2pi is applied in gg_hgg_v.f)
      sum=zip
      do h1=1,2
c      do h2=1,2
      h2=h1

c--- note: A42 is renormalized with the opposite sign
c      amp41_a(h1,h2)=amp41_a(h1,h2)+ren*amplo_a(h1,h2)
c      amp41_b(h1,h2)=amp41_b(h1,h2)+ren*amplo_b(h1,h2)
      amp42_a(h1,h2)=amp42_a(h1,h2)-ren*amplo_a(h1,h2)
      amp42_b(h1,h2)=amp42_b(h1,h2)-ren*amplo_b(h1,h2)

c--- We commment out these two contributions in order that we can
c--- take over the existing structure in gg_hgg_v.f. This routine
c--- then returns the same quantity as Hqaqavsq.f
c      sum=sum+dble(amp41_a(h1,h2)*dconjg(amplo_a(h1,h2)))
c     .       +dble(amp41_b(h1,h2)*dconjg(amplo_b(h1,h2)))
c      if (h1 .eq. h2) then
      sum=sum-1d0/xn*dble(amp42_a(h1,h2)*dconjg(amplo_b(h1,h2)))
     .       -1d0/xn*dble(amp42_b(h1,h2)*dconjg(amplo_a(h1,h2)))
c      endif
c      enddo
      enddo

      Ampvirtsq_AQaq_ident=sum*V*xn/4d0

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
subroutine banfac ( w, nroww, nrow, nbandl, nbandu, iflag )

!*****************************************************************************80
!
!! BANFAC factors a banded matrix without pivoting.
!
!  Discussion:
!
!    BANFAC returns in W the LU-factorization, without pivoting, of
!    the banded matrix A of order NROW with (NBANDL+1+NBANDU) bands
!    or diagonals in the work array W.
!
!    Gauss elimination without pivoting is used.  The routine is
!    intended for use with matrices A which do not require row
!    interchanges during factorization, especially for the totally
!    positive matrices which occur in spline calculations.
!
!    The matrix storage mode used is the same one used by LINPACK
!    and LAPACK, and results in efficient innermost loops.
!
!    Explicitly, A has
!
!      NBANDL bands below the diagonal
!      1     main diagonal
!      NBANDU bands above the diagonal
!
!    and thus, with MIDDLE=NBANDU+1,
!    A(I+J,J) is in W(I+MIDDLE,J) for I=-NBANDU,...,NBANDL, J=1,...,NROW.
!
!    For example, the interesting entries of a banded matrix
!    matrix of order 9, with NBANDL=1, NBANDU=2:
!
!      11 12 13  0  0  0  0  0  0
!      21 22 23 24  0  0  0  0  0
!       0 32 33 34 35  0  0  0  0
!       0  0 43 44 45 46  0  0  0
!       0  0  0 54 55 56 57  0  0
!       0  0  0  0 65 66 67 68  0
!       0  0  0  0  0 76 77 78 79
!       0  0  0  0  0  0 87 88 89
!       0  0  0  0  0  0  0 98 99
!
!    would appear in the first 1+1+2=4 rows of W as follows:
!
!       0  0 13 24 35 46 57 68 79
!       0 12 23 34 45 56 67 78 89
!      11 22 33 44 55 66 77 88 99
!      21 32 43 54 65 76 87 98  0
!
!    All other entries of W not identified in this way with an
!    entry of A are never referenced.
!
!    This routine makes it possible to solve any particular linear system
!    A*X=B for X by the call
!
!      call banslv ( w, nroww, nrow, nbandl, nbandu, b )
!
!    with the solution X contained in B on return.
!
!    If IFLAG=2, then one of NROW-1, NBANDL, NBANDU failed to be nonnegative,
!    or else one of the potential pivots was found to be zero
!    indicating that A does not have an LU-factorization.  This
!    implies that A is singular in case it is totally positive.
!
!  Modified:
!
!    14 February 2007
!
!  Author:
!
!    Carl DeBoor
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer, 2001,
!    ISBN: 0387953663,
!    LC: QA1.A647.v27.
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) W(NROWW,NROW).
!    On input, W contains the "interesting" part of a banded
!    matrix A, with the diagonals or bands of A stored in the
!    rows of W, while columns of A correspond to columns of W.
!    On output, W contains the LU-factorization of A into a unit
!    lower triangular matrix L and an upper triangular matrix U
!    (both banded) and stored in customary fashion over the
!    corresponding entries of A.
!
!    Input, integer ( kind = 4 ) NROWW, the row dimension of the work array W.
!    NROWW must be at least NBANDL+1 + NBANDU.
!
!    Input, integer ( kind = 4 ) NROW, the number of rows in A.
!
!    Input, integer ( kind = 4 ) NBANDL, the number of bands of A below
!    the main diagonal.
!
!    Input, integer ( kind = 4 ) NBANDU, the number of bands of A above
!    the main diagonal.
!
!    Output, integer ( kind = 4 ) IFLAG, error flag.
!    1, success.
!    2, failure, the matrix was not factored.
!
  implicit none

  integer ( kind = 4 ) nrow
  integer ( kind = 4 ) nroww

  real ( kind = 8 ) factor
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) middle
  integer ( kind = 4 ) nbandl
  integer ( kind = 4 ) nbandu
  real ( kind = 8 ) pivot
  real ( kind = 8 ) w(nroww,nrow)

  iflag = 1

  if ( nrow < 1 ) then
    iflag = 2
    return
  end if
!
!  W(MIDDLE,*) contains the main diagonal of A.
!
  middle = nbandu + 1

  if ( nrow == 1 ) then
    if ( w(middle,nrow) == 0.0D+00 ) then
      iflag = 2
    end if
    return
  end if
!
!  A is upper triangular.  Check that the diagonal is nonzero.
!
  if ( nbandl <= 0 ) then

    do i = 1, nrow-1
      if ( w(middle,i) == 0.0D+00 ) then
        iflag = 2
        return
      end if
    end do

    if ( w(middle,nrow) == 0.0D+00 ) then
      iflag = 2
    end if

    return
!
!  A is lower triangular.  Check that the diagonal is nonzero and
!  divide each column by its diagonal.
!
  else if ( nbandu <= 0 ) then

    do i = 1, nrow-1

      pivot = w(middle,i)

      if ( pivot == 0.0D+00 ) then
        iflag = 2
        return
      end if

      do j = 1, min ( nbandl, nrow-i )
        w(middle+j,i) = w(middle+j,i) / pivot
      end do

    end do

    return

  end if
!
!  A is not just a triangular matrix.
!  Construct the LU factorization.
!
  do i = 1, nrow - 1
!
!  W(MIDDLE,I) is the pivot for the I-th step.
!
    if ( w(middle,i) == 0.0D+00 ) then
      iflag = 2
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'BANFAC - Fatal error!'
      write ( *, '(a,i8)' ) '  Zero pivot encountered in column ', i
      stop
    end if
!
!  Divide each entry in column I below the diagonal by PIVOT.
!
    do j = 1, min ( nbandl, nrow-i )
      w(middle+j,i) = w(middle+j,i) / w(middle,i)
    end do
!
!  Subtract A(I,I+K)*(I-th column) from (I+K)-th column (below row I).
!
    do k = 1, min ( nbandu, nrow-i )
      factor = w(middle-k,i+k)
      do j = 1, min ( nbandl, nrow-i )
        w(middle-k+j,i+k) = w(middle-k+j,i+k) - w(middle+j,i) * factor
      end do
    end do

  end do
!
!  Check the last diagonal entry.
!
  if ( w(middle,nrow) == 0.0D+00 ) then
    iflag = 2
  end if

  return
end

max_stars_repo_name: horvatm/b-spline-vector-function

max_issues_repo_name: horvatm/b-spline-vector-function

max_forks_repo_name: horvatm/b-spline-vector-function


text:
C$Procedure     EDTCMD ( Edit a file using a specified text editor )

      SUBROUTINE EDTCMD ( CMD, FILE )

C$ Abstract
C
C     Edit a file using a specified editor.
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     SYSTEM
C     UTILITY
C
C$ Declarations

      CHARACTER*(*)         CMD
      CHARACTER*(*)         FILE

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     CMD        I   Command string used to invoke editor.
C     FILE       I   Name of file to edit.
C
C$ Detailed_Input
C
C     CMD            is a character string containing the command
C                    used to invoke a text editor available on the
C                    system under which the calling program is running.
C                    This routine will invoke the specified editor
C                    using FILE as the target file to edit.  The name
C                    of the file to be edited is not included in the
C                    command; this name is input as a separate argument.
C
C                    Case sensitivity of CMD varies with the system on
C                    which the calling program is run.
C
C                    Trailing white space in CMD is not significant.
C
C
C     FILE           is the name of a file that is to be edited.  FILE
C                    need not exist at the time this routine is called.
C
C                    Case sensitivity of FILE varies with the system on
C                    which the calling program is run.
C
C                    Trailing white space in FILE is not significant.
C
C$ Detailed_Output
C
C     None.  See $Particulars for further information on the action of
C     this routine.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If the specified edit command fails, the error will be
C         diagnosed by routines called by this routine.
C
C     2)  If the editing session started by this routine is terminated
C         abnormally, the effect on the operation of the calling program
C         is unspecified.
C
C$ Files
C
C     See $Particulars.
C
C$ Particulars
C
C     This routine should be used with caution; calls to this routine
C     may have unintended side effects on the operation of the calling
C     program.  A solid understanding of the global operation of the
C     calling program is a prerequisite for programmers wishing to
C     use this routine.
C
C     The input argument FILE should unambiguously designate a file
C     that can be edited by the specified editor on the system under
C     which the calling program is being run.  The calling program
C     should have read or write privileges consistent with the intended
C     mode of access to FILE.
C
C     This routine may fail to recover in a predictable fashion from
C     abnormal termination of an editing session.
C
C$ Examples
C
C     1)   On a VAX/VMS system, the EDT editor could be invoked by
C          the calls
C
C             CALL EDTCMD ( 'EDIT/EDT',  FILE  )
C
C          or
C
C             CALL EDTCMD ( 'EDIT/EDT/COMMAND = <command file>',  FILE )
C
C
C     2)   On a Unix system, the emacs editor could be invoked
C          (normally) by the calls
C
C              CALL EDTCMD ( 'emacs', FILE )
C
C          or
C
C              CALL EDTCMD ( '/usr/bin/emacs', FILE )
C
C
C$ Restrictions
C
C     1)   The means by which this routine invokes an editor are system-
C          dependent; invoking the editor may have side effects that
C          affect the operation of the calling program.  For example,
C          on Unix systems, this routine may start a new shell in which
C          to run the editor; starting a new shell may interfere with
C          any sequential file I/O in progress at the time the shell is
C          started.
C
C          See the code for implementation details.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    Beta Version 2.27.0, 10-MAR-2014 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-INTEL.
C
C-    Beta Version 2.26.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-LINUX-64BIT-IFORT.
C
C-    Beta Version 2.25.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-GFORTRAN.
C
C-    Beta Version 2.24.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-64BIT-GFORTRAN.
C
C-    Beta Version 2.23.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-64BIT-GCC_C.
C
C-    Beta Version 2.22.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL.
C
C-    Beta Version 2.21.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL-CC_C.
C
C-    Beta Version 2.20.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL-64BIT-CC_C.
C
C-    Beta Version 2.19.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-NATIVE_C.
C
C-    Beta Version 2.18.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-WINDOWS-64BIT-IFORT.
C
C-    Beta Version 2.17.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-LINUX-64BIT-GFORTRAN.
C
C-    Beta Version 2.16.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-64BIT-MS_C.
C
C-    Beta Version 2.15.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-INTEL_C.
C
C-    Beta Version 2.14.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-IFORT.
C
C-    Beta Version 2.13.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-GFORTRAN.
C
C-    Beta Version 2.12.0, 18-MAR-2009 (BVS)
C
C        Updated for PC-LINUX-GFORTRAN.
C
C-    Beta Version 2.11.0, 18-MAR-2009 (BVS)
C
C        Updated for MAC-OSX-GFORTRAN.
C
C-    Beta Version 2.10.0, 19-FEB-2008 (BVS)
C
C        Updated for PC-LINUX-IFORT.
C
C-    Beta Version 2.9.0, 14-NOV-2006 (BVS)
C
C        Updated for PC-LINUX-64BIT-GCC_C.
C
C-    Beta Version 2.8.0, 14-NOV-2006 (BVS)
C
C        Updated for MAC-OSX-INTEL_C.
C
C-    Beta Version 2.7.0, 14-NOV-2006 (BVS)
C
C        Updated for MAC-OSX-IFORT.
C
C-    Beta Version 2.6.0, 14-NOV-2006 (BVS)
C
C        Updated for PC-WINDOWS-IFORT.
C
C-    Beta Version 2.5.0, 26-OCT-2005 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-GCC_C.
C
C-    Beta Version 2.4.0, 03-JAN-2005 (BVS)
C
C        Updated for PC-CYGWIN_C.
C
C-    Beta Version 2.3.0, 03-JAN-2005 (BVS)
C
C        Updated for PC-CYGWIN.
C
C-    Beta Version 2.2.5, 17-JUL-2002 (BVS)
C
C        Added MAC-OSX environments.
C
C-    Beta Version 2.2.4, 08-OCT-1999 (WLT)
C
C        The environment lines were expanded so that the supported
C        environments are now explicitely given.  New
C        environments are WIN-NT
C
C-    Beta Version 2.2.3, 20-SEP-1999 (NJB)
C
C        CSPICE and PC-LINUX environment lines were added.  Some
C        typos were corrected.
C
C-    Beta Version 2.2.2, 28-JUL-1999 (WLT)
C
C        The environment lines were expanded so that the supported
C        environments are now explicitly given.  New
C        environments are PC-DIGITAL, SGI-O32 and SGI-N32.
C
C-    Beta Version 2.2.1, 18-MAR-1999 (WLT)
C
C        The environment lines were expanded so that the supported
C        environments are now explicitly given.  Previously,
C        environments such as SUN-SUNOS and SUN-SOLARIS were implied
C        by the environment label SUN.
C
C-    Beta Version 2.2.0, 12-AUG-1996 (WLT)
C
C        Added DEC-OSF1 to the list of supported environments
C
C-    Beta Version 2.1.0, 10-JAN-1996 (WLT)
C
C        Added PC-LAHEY to the list of supported environments.
C
C-    Beta Version 2.0.0, 16-JUN-1995 (WLT)(HAN)
C
C        Created master file from collection of machine dependent
C        routines.  Copyright notice added.
C
C-    Beta Version 1.0.0, 16-AUG-1994 (NJB)
C
C-&

C$ Index_Entries
C
C     invoke a text editor within a program
C
C-&

C
C     SPICELIB functions
C
      INTEGER               RTRIM
      LOGICAL               RETURN

C
C     Local parameters
C
      INTEGER               CMDLEN
      PARAMETER           ( CMDLEN = 255 )

C
C     Local variables
C
      CHARACTER*(CMDLEN)    LOCCMD

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

C
C
C     SUN:
C
C        Computer:         Sun SPARCstation 2
C        Operating System: Sun OS 4.1.2
C        Fortran:          Sun FORTRAN 1.3.1
C
C     HP:
C
C        Computer:         HP 715/50
C        Operating System: HP-UX 9.01
C        Fortran:          HP-UX.09.00.24
C                          HP-UX FORTRAN/9000
C                             Series 700 B2408A.09.00
C                             Series 800 B2409B.09.00
C
C     NEXT:
C
C        Computer:         NeXT
C        Operating System: NeXtStep 3.0, 3.2
C        Fortran:          Absoft Fortran V3.2
C        NEXT (NeXT 3.0, Absoft Fortran 3.2):
C
C        Computer:         Alpha/OSF1
C        Operating System: OSF1 V3.2
C        Fortran:        : DEC Fortran Compiler Driver V3.5-053
C
C
C     Build the edit command to be passed to the system.
C
      LOCCMD = CMD
      CALL SUFFIX ( FILE, 1, LOCCMD )

C
C     For safety, append a null to the command.  If the user has
C     linked to a C version of the "system" function, this may
C     save us much grief.
C
      CALL SUFFIX ( CHAR(0), 0, LOCCMD )

C
C     Invoke the editor.
C
      CALL EXESYS (   LOCCMD( :RTRIM(LOCCMD) )   )

      CALL CHKOUT ( 'EDTCMD' )
      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
C PR rtl-optimization/58968.f
C { dg-do compile { target powerpc*-*-*} }
C { dg-options "-mcpu=power7 -O3 -w -ffast-math  -funroll-loops" }
      SUBROUTINE MAKTABS(IW,SOME,LBOX1,LBOX2,LBOX3,NSPACE,NA,NB,
     *            LBST,X,
     *            NX,IAMA,IAMI,IBMA,IBMI,MNUM,IDIM,MSTA,IBO,
     *            IDSYM,ISYM1,NSYM,
     *            NACT,LWRK,KTAB,LGMUL,
     *            LCON,LCOA,LCOB,
     *            LANDET,LBNDET,NAST,NBST,LSYMA,LSYMB,LGCOM,
     *            MINI,MAXI,LSPA,LSPB,LDISB,
     *            LSAS,LSBS,LSAC,LSBC,
     *            ITGA,ITGB,IAST,IBST,NCI,NA1EX,NB1EX,FDIRCT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL SOME
      DIMENSION LBOX1(NSPACE),LBOX2(NSPACE),LBOX3(NSPACE),LBST(NSPACE)
      DIMENSION X(NX)
      DIMENSION IAMA(NSPACE),IAMI(NSPACE),IBMA(NSPACE),IBMI(NSPACE)
      DIMENSION MNUM(NSPACE),IDIM(NSPACE),MSTA(NSPACE+1),IBO(NACT)
      DIMENSION LWRK(43),KTAB(NSYM),LGMUL(NSYM,NSYM)
      DIMENSION LCON(NA)
      DIMENSION LCOA(NSYM,ITGA),LCOB(NSYM,ITGB)
      DIMENSION LANDET(NSPACE,ITGA),LBNDET(NSPACE,ITGB)
      DIMENSION NAST(ITGA+1),NBST(ITGB+1)
      DIMENSION LSYMA(IAST),LSYMB(IBST)
      DIMENSION LGCOM(ITGB,ITGA)
      DIMENSION MINI(NSPACE),MAXI(NSPACE)
      DIMENSION LSPA(IAST),LSPB(IBST)
      DIMENSION LDISB(NSYM,ITGB,ITGA)
      DIMENSION LSAS(NSYM+1,ITGA),LSBS(NSYM+1,ITGB)
      DIMENSION LSAC(IAST),LSBC(IBST)
      LOGICAL FDIRCT
      LCOA = 0
      LCOB = 0
      ISTA1 = LBST(1)
      CALL RESETCO(LBOX1,NSPACE,NB,IBMA,IBMI,LBOX2)
      NAST(1) = 0
      NBST(1) = 0
      DO II=1,ITGA
         ITOT = 1
         DO JJ=1,NSPACE
            ITOT = ITOT * LANDET(JJ,II)
         ENDDO
         NAST(II+1) = NAST(II) + ITOT
      ENDDO
      DO II=1,ITGB
         ITOT = 1
         DO JJ=1,NSPACE
            ITOT = ITOT * LBNDET(JJ,II)
         ENDDO
         NBST(II+1) = NBST(II) + ITOT
      ENDDO
      ICOMP = 0
      CALL RESETCO(LBOX1,NSPACE,NA,IAMA,IAMI,LBOX3)
      NA1EX = 0
      NB1EX = 0
      CALL RESETCO(LBOX1,NSPACE,NB,IBMA,IBMI,LBOX3)
      DO IIB = 1,ITGB
         CALL RESETDE(LBOX1,NSPACE,NB,MSTA,LCON)
         DO KKB=NBST(IIB)+1,NBST(IIB+1)
            DO II=1,NSPACE
               LBOX2(II) = LBOX1(II)
            ENDDO
            IEBS = NB+1
            DO ISPB1=NSPACE,1,-1
               IOC1 = LBOX1(ISPB1)
               IEBE = IEBS - 1
               IEBS = IEBS - IOC1
               LBOX2(ISPB1) = LBOX2(ISPB1)-1
               DO IB1=IEBE,IEBS,-1
                  IO1 = LCON(IB1)
                  IGBE = IEBE - LBOX1(ISPB1)
                  DO ISPB2=ISPB1,NSPACE
                     IGBS = IGBE + 1
                     IGBE = IGBE + LBOX1(ISPB2)
                     LBOX2(ISPB2) = LBOX2(ISPB2) + 1
                     IGBA = MAX(IB1+1,IGBS)
                     DO IGAP=IGBA,IGBE+1
                        DO JJ=ISTA,IEND
                           NB1EX = NB1EX + 1
                        ENDDO
                        ISTA = LCON(IGAP)+1
                        IEND = LCON(IGAP+1)-1
                        IF (IGAP.EQ.IGBE) IEND=MSTA(ISPB2+1)-1
                     ENDDO
                     LBOX2(ISPB2) = LBOX2(ISPB2) - 1
                  ENDDO
               ENDDO
               LBOX2(ISPB1) = LBOX2(ISPB1) + 1
            ENDDO
            CALL MOVEUP2(LBOX1,NSPACE,NB,MSTA,LCON)
         ENDDO
         CALL PUSHCO(LBOX1,NSPACE,NB,IBMA,IBMI,LBOX3,IEND)
      ENDDO
      RETURN
      END


! { dg-do run }
!
! Test for the fix for PR34640. In this case, final testing of the
! patch revealed that in some cases the actual descriptor was not
! being passed to procedure dummy pointers.
!
! Contributed by Thomas Koenig  <tkoenig@netcologne.de>
!
module x
  use iso_c_binding
  implicit none
  type foo
     complex :: c
     integer :: i
  end type foo
contains
  subroutine printit(c, a)
    complex, pointer, dimension(:) :: c
    integer :: i
    integer(kind=c_intptr_t) :: a
    a = transfer(c_loc(c(2)),a)
  end subroutine printit
end module x

program main
  use x
  use iso_c_binding
  implicit none
  type(foo), dimension(5), target :: a
  integer :: i
  complex, dimension(:), pointer :: pc
  integer(kind=c_intptr_t) :: s1, s2, s3
  a%i = 0
  do i=1,5
     a(i)%c = cmplx(i**2,i)
  end do
  pc => a%c
  call printit(pc, s3)

  s1 = transfer(c_loc(a(2)%c),s1)
  if (s1 /= s3) STOP 1

  s2 = transfer(c_loc(pc(2)),s2)
  if (s2 /= s3) STOP 2

end program main


! ###############################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

   MODULE OFP3_ELFE_2D_Interface

   INTERFACE

      SUBROUTINE OFP3_ELFE_2D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW )

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  WRT_BUG, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, ELOUT_ELFE_BIT, FATAL_ERR, IBIT, INT_SC_NUM, MBUG, MOGEL,                   &
                                         WARN_ERR, NELE, NCQUAD4, NCQUAD4K, NCSHEAR, NCTRIA3, NCTRIA3K, SOL_NAME
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  OFP3_ELFE_2D_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO
      USE FEMAP_ARRAYS, ONLY          :  FEMAP_EL_NUMS, FEMAP_EL_VECS
      USE PARAMS, ONLY                :  OTMSKIP, POST
      use model_stuf, only            :  pcomp_props
      USE MODEL_STUF, ONLY            :  ANY_ELFE_OUTPUT, EDAT, EPNT, ETYPE, FCONV, EID, ELMTYP, ELOUT, METYPE, NUM_EMG_FATAL_ERRS,&
                                         PLY_NUM, TYPE, STRESS
      USE LINK9_STUFF, ONLY           :  EID_OUT_ARRAY, MAXREQ, OGEL
      USE OUTPUT4_MATRICES, ONLY      :  OTM_ELFE, TXT_ELFE

      IMPLICIT NONE

      CHARACTER( 1*BYTE), PARAMETER   :: IHDR      = 'Y'   ! An input to subr WRITE_GRID_OUTPUTS, called herein

      INTEGER(LONG), INTENT(IN)       :: FEMAP_SET_ID      ! Set ID for FEMAP output
      INTEGER(LONG), INTENT(IN)       :: ITE               ! Unit number for text files for OTM row descriptors
      INTEGER(LONG), INTENT(IN)       :: JVEC              ! Solution vector number
      INTEGER(LONG), INTENT(INOUT)    :: OT4_EROW          ! Row number in OT4 file for elem related OTM descriptors
      integer(long)                   :: num_pcomp_elems   ! number of elements that are composites (used to prevent output of engr
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = OFP3_ELFE_2D_BEGEND

      END SUBROUTINE OFP3_ELFE_2D

   END INTERFACE

   END MODULE OFP3_ELFE_2D_Interface


text:
module m_fftw_test

contains
   subroutine test_sub()
      use iso_c_binding
      include 'fftw3.f03'
      type(C_PTR)                                     :: plan
      complex(C_DOUBLE_COMPLEX), dimension(1024,1000) :: in, out

      plan = fftw_plan_dft_2d(1000,1024, in,out, FFTW_FORWARD,FFTW_ESTIMATE)

      call fftw_execute_dft(plan, in, out)
      call fftw_destroy_plan(plan)
   endsubroutine test_sub
end module

program test_prog
   use m_fftw_test

   call test_sub()
 end program test_prog

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
LOGICAL FUNCTION OPNLIST3( FID, PGNAME ) RESULT( OFLAG )

    !!***********************************************************************
    !! Version "$Id: opnlist3.f90 118 2019-06-15 20:50:32Z coats $"
    !! EDSS/Models-3 I/O API.
    !! Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr.,
    !! (C) 2003-2013 Baron Advanced Meteorological Systems,
    !! (C) 2007-2013 Carlie J. Coats, Jr., and
    !! (C) 2014 UNC Institute for the Environment.
    !! Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1
    !! See file "LGPL.txt" for conditions of use.
    !!.........................................................................
    !!  function body starts at line  83
    !!
    !!  FUNCTION:
    !!       Open a FILE-LIST multi-file data set for status FSREAD3
    !!       with logical name FLIST3( FID ).
    !!
    !!  RETURN VALUE:
    !!       TRUE iff it succeeds in opening the file, reading its
    !!        attributes, and storing the relevant ones in STATE3.EXT
    !!
    !!  PRECONDITIONS REQUIRED:
    !!       All files in the logical-name list already exist.
    !!
    !!  REVISION  HISTORY:
    !!       Prototype  03/2002 by CJC for I/O API V2.2
    !!       Modified   02/2015 by CJC for I/O API 3.2: USE M3UTILIO
    !!
    !!      Modified 10/2015 by CJC for I/O API 3.2: use NF_*() instead of NC*()
    !!      for netCDF-Fortran 4.x compatibility
    !!***********************************************************************

    USE M3UTILIO
    USE MODNCFIO

    IMPLICIT NONE

    !!...........   INCLUDES:

    INCLUDE 'STATE3.EXT'

    !!...........   ARGUMENTS and their descriptions:

    INTEGER      , INTENT(IN   ) :: FID     !  index into STATE3 tables
    CHARACTER*(*), INTENT(IN   ) :: PGNAME  !  name of calling program

    !!...........   EXTERNAL FUNCTIONS and their descriptions:

    LOGICAL, EXTERNAL :: OPNFIL3  !  does work of opening "old" files

    !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

    CHARACTER*512   EQNAME  !  environment value of FLIST3( FID )
    CHARACTER*16    FIL16   !  scratch file-name buffer

    CHARACTER*256   MESG    !  for m3msg2, m3warn
    INTEGER         IDUM, IERR
    INTEGER         F, I, J, L, V

    !!  Support for FILE-LIST multi-file input data sets:
    !!  LISTFLAG = true iff FLIST3( FID ) is a FILE-LIST
    !!  LISTCNT  = # of elements in the list
    !!  LISTNAME( 1:LISTCNT ) = list of the logical names

    LOGICAL         AFLAG, EFLAG
    INTEGER         LISTCNT
    CHARACTER*16    LISTNAME( MXFILE3 )

    !!...........   STATEMENT FUNCTION:  REAL, REAL*8 "definitely unequal"

    LOGICAL         DBLERR
    REAL*8          P, Q

    DBLERR( P, Q ) = ( (P - Q)**2  .GT.  1.0E-10*( P*P + Q*Q + 1.0E-5 ) )

    !!.............................................................................
    !!   begin body of subroutine  OPNFIL3

    IF ( .NOT. STRLIST( FLIST3( FID ), 'Logical name list',  MXFILE3, LISTCNT, LISTNAME ) ) THEN

         MESG = FIL16 // ':'//TRIM( EQNAME )
        CALL M3MSG2( MESG )
        CALL M3WARN( 'OPEN3', 0, 0, 'Bad FILE_LIST.' )
        FLIST3( FID ) = CMISS3
        OFLAG         = .FALSE.
        RETURN

    END IF

    EFLAG = .FALSE.

    L     = ILCNT3
    IFRST3( FID ) = L + 1
    NLIST3( FID ) = LISTCNT
    CDFID3( FID ) = LSTFIL3
    SDATE3( FID ) = 99999999
    STIME3( FID ) = 0

    DO  F = 1, LISTCNT

        I = INDEX1( CMISS3, MXFILE3, FLIST3 )
        FIL16 = LISTNAME( F )
        CALL NAMEVAL( FIL16, EQNAME )
        INQUIRE ( FILE = EQNAME, EXIST = AFLAG )

        IF ( I .LE. 0 ) THEN

            MESG = 'Could not open ' // FIL16 // ' Max # of files already opened.'
            CALL M3MSG2( MESG )
            EFLAG = .TRUE.

        ELSE IF ( AFLAG ) THEN

            IDUM = MAX ( INDEX( EQNAME, '-v' ) , INDEX( EQNAME, '-V' ) )
            IF ( IDUM .GT. 0 ) THEN
                EQNAME( IDUM:IDUM+1 ) = '  ' !  fix the '-v' (etc.)
                VOLAT3( I ) = .TRUE.         !  volatile file
            ELSE
                VOLAT3( I ) = .FALSE.
            END IF

            FLIST3( I ) = FIL16

            IF ( OPNFIL3( EQNAME, I, FSREAD3, PGNAME ) ) THEN

                L = L + 1
                ILIST3( L ) = I

                IF ( F .EQ. 1 ) THEN

                    SDATE3( FID ) = SDATE3( I )
                    STIME3( FID ) = STIME3( I )
                    TSTEP3( FID ) = TSTEP3( I )
                    FTYPE3( FID ) = FTYPE3( I )
                    GDTYP3( FID ) = GDTYP3( I )
                    NCOLS3( FID ) = NCOLS3( I )
                    NROWS3( FID ) = NROWS3( I )
                    NLAYS3( FID ) = NLAYS3( I )
                    NVARS3( FID ) = NVARS3( I )
                    NTHIK3( FID ) = NTHIK3( I )
                    P_ALP3( FID ) = P_ALP3( I )
                    P_BET3( FID ) = P_BET3( I )
                    P_GAM3( FID ) = P_GAM3( I )
                    XCENT3( FID ) = XCENT3( I )
                    YCENT3( FID ) = YCENT3( I )
                    XORIG3( FID ) = XORIG3( I )
                    YORIG3( FID ) = YORIG3( I )
                    XCELL3( FID ) = XCELL3( I )
                    YCELL3( FID ) = YCELL3( I )

                    DO  V = 1, NVARS3( FID )
                        VLIST3( V,FID ) = VLIST3( V,I )
                        VTYPE3( V,FID ) = VTYPE3( V,I )
                    END DO

                ELSE IF ( .NOT. EFLAG ) THEN

                    IF ( TSTEP3( I ) .NE. TSTEP3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad TSTEP for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( FTYPE3( I ) .NE. FTYPE3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad FTYPE for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( GDTYP3( I ) .NE. GDTYP3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad GDTYP for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( NCOLS3( I ) .NE. NCOLS3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad NCOLS for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( NROWS3( I ) .NE. NROWS3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad NROWS for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( NLAYS3( I ) .NE. NLAYS3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad NLAYS for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( NTHIK3( I ) .NE. NTHIK3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad NTHIK for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( NVARS3( I ) .NE. NVARS3( FID ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad NVARS for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( P_ALP3( I ), P_ALP3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad P_ALP for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( P_BET3( I ), P_BET3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad P_BET for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( P_GAM3( I ), P_GAM3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad P_GAM for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( XORIG3( I ), XORIG3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad XORIG for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( YORIG3( I ), YORIG3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad YORIG for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( XCENT3( I ), XCENT3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad XCENT for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( YCENT3( I ), YCENT3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad YCENT for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( XCELL3( I ), XCELL3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad XCELL for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    IF ( DBLERR( YCELL3( I ), YCELL3(FID) ) ) THEN
                        EFLAG = .TRUE.
                        MESG = 'Bad YCELL for ' //FIL16// ' in file list ' // FLIST3( FID )
                        CALL M3WARN( 'OPEN3', 0, 0, MESG )
                    END IF

                    DO  V = 1, NVARS3( FID )

                        IF ( VLIST3(V,FID) .NE. VLIST3(V,I) ) THEN
                            EFLAG = .TRUE.
                            MESG = 'Bad VNAME for ' //FIL16// ' in file list ' // FLIST3( FID )
                            CALL M3WARN( 'OPEN3', 0, 0, MESG )
                        END IF

                        IF ( VTYPE3(V,FID) .NE. VTYPE3(V,I) ) THEN
                            EFLAG = .TRUE.
                            MESG = 'Bad VNAME for ' //FIL16// ' in file list ' // FLIST3( FID )
                            CALL M3WARN( 'OPEN3', 0, 0, MESG )
                        END IF

                    END DO

                    IF ( EFLAG ) THEN
                        FLIST3( I ) = CMISS3
                    ELSE IF ( SECSDIFF( SDATE3( FID ),STIME3( FID ),     &
                                        SDATE3(   I ),STIME3(   I ) ) .LT. 0 ) THEN
                        SDATE3( FID ) = SDATE3( I )
                        STIME3( FID ) = STIME3( I )
                    END IF

                END IF              !  if F = 1 or not

            ELSE

                EFLAG = .TRUE.
                FLIST3( I ) = CMISS3

            END IF          !  if opnfil3() or not

        ELSE                !  i > 0 but not aflag:

            EFLAG = .TRUE.
            MESG  = FIL16//':'//TRIM(EQNAME)
            CALL M3MSG2( MESG )
            CALL M3WARN( 'OPNLIST3', 0, 0, 'File not available.' )
            FLIST3( I ) = CMISS3

        END IF      ! if i le 0; else if AFLAG; or not

    END DO                  !  end loop on names in this list

    IF ( EFLAG ) THEN       !  process error...

        L = ILCNT3
        DO  F = ILCNT3 + 1, L+1
            I = CDFID3( F )
            IF ( I .GT. 0 ) THEN
                IERR = NF_CLOSE( I )
                IF ( IERR .NE. 0 ) THEN
                    CALL M3ABORT( FLIST3( F ), I, IERR, 'Error NF_CLOSE-ing file' )
                END IF
            END IF
            FLIST3( I ) = CMISS3
            CDFID3( I ) = IMISS3
        END DO

        FLIST3( FID ) = CMISS3
        CDFID3( FID ) = IMISS3
        OFLAG         = .FALSE.
        RETURN

    ELSE                    !  process success

        DO  I = IFRST3(FID), IFRST3(FID) + NLIST3(FID) - 1
            F = ILIST3( I )
            J = JSTEP3( SDATE3(   F ), STIME3(   F ),       &
                        SDATE3( FID ), STIME3( FID ), TSTEP3( FID ) )
            BEGRC3( F ) = J
            ENDRC3( F ) = J + MXREC3( F ) - 1
            COUNT3      = MAX( COUNT3, ILIST3( F ) )
            MXREC3( FID ) = MAX( MXREC3( FID ), ENDRC3( F ) )
        END DO
        NLIST3( FID ) = LISTCNT
        IFRST3( FID ) = ILCNT3 + 1
        ILCNT3        = ILCNT3 + LISTCNT
        OFLAG         = .TRUE.

    END IF                  !  if error, or not

    RETURN

END FUNCTION OPNLIST3

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
	SUBROUTINE TA_DCOD ( cldt, bufrta_old, bufrta_new,
     +						gemfil, nhours, iret)
C************************************************************************
C* TA_DCOD								*
C*									*
C* This routine decodes bulletins containing TAMDAR BUFR messages      	*
C* received directly from AirDat into GEMPAK sounding format and 	*
C* writes the data into the output sounding file.			*
C*									*
C* TA_DCOD ( CLDT, BUFRTA_OLD, BUFRTA_NEW, GEMFIL, NHOURS, IRET )	*
C*									*
C* Input parameters:							*
C*	CLDT		CHAR*		Date-time from command line	*
C*	BUFRTA_OLD	CHAR*		Old TAMDAR BUFR table file   	*
C*	BUFRTA_NEW	CHAR*		New TAMDAR BUFR table file	*
C*	GEMFIL		CHAR*		Output file name template	*
C*	NHOURS		INTEGER		Max # of hours before run time	*
C*					for creating BUFR output	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code:			*
C*					  0 = normal return		*
C*									*
C**									*
C* Log:									*
C* C. Caruso Magee/NCEP	07/07 						*
C* C. Caruso Magee/NCEP	10/07  	Add pressure (007004) to input descrpor	*
C*                             	list.					*
C* m.gamazaychikov/SAIC	02/08	Change output stream to sounding format	*
C* m.gamazaychikov/SAIC	07/08	Cleaned up				*
C* T. Piper/SAIC	08/08	Changed '.dummy/' to '.dummy.'		*
C* T. Piper/SAIC	09/08	Completely restructured logic		*
C* T. Piper/SAIC	10/08	Create unique "Observation ID"		*
C* T. Piper/SAIC	12/08	Added capability to process new format	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'BRIDGE.PRM'
	INCLUDE		'tacmn.cmn'
C*
C*      Number of non-NCEP BUFR tables files.
C*
        PARAMETER       ( NBUFRT = 2 )
C
	PARAMETER	( NACRN=1000, NLV=50, NPARMS=6 )
C
	CHARACTER*(*)	cldt, bufrta_old, bufrta_new, gemfil
C
	CHARACTER	acrn(NACRN)*9, afdate(NACRN, NLV)*13, bbb*8,
     +			bfstyp*8, bufrtf(NBUFRT)*(DCMXLN), buhd*8,
     +			bull*(DCMXBF), bulldt*8, cbf,
     +			cborg*8, cbull*(DCMXBF), cdum*13,
     +                  cimnem(NCIMN)*8, coun(NACRN)*2, cracrn*9,
     +                  dattim*13, filnam*132, newfil*132,
     +			parms(NPARMS)*4, prmfil*72,
     +			rimnem(NRIMN)*8, rundt*12, seqnum*8,
     +			stat(NACRN)*2, sysdt*12
C
        REAL            corn
        REAL            alat(NACRN, NLV), alon(NACRN, NLV),
     +                  datai(NPARMS*NLV), datao(NPARMS*NLV),
     +			flat(NACRN), flon(NACRN),
     +                  poaf(NACRN, NLV),
     +                  wspd(NACRN, NLV), wdir(NACRN, NLV),
     +                  hmdt(NACRN, NLV), tmpk(NACRN, NLV),
     +			hght(NACRN, NLV), pres(NACRN, NLV),
     +                  selv(NACRN)
C
	INTEGER		iacrn(NACRN), iarcft(5), ibull(DCMXBF/4),
     +			irptdt(5), irundt(5), iubfmf(NBUFRT),
     +			iubftf(NBUFRT), nmlev(NACRN), nxdsc(NBUFRT)
C
	LOGICAL         addstn, amatch, bullok, cirflg, datflg,
     +			exist, mrgdat, pkflg, stmflg
C
	EQUIVALENCE	( cbull(1:4), ibull(1) )
C*
C*	Number of expected descriptors within Section 3 of each
C*	type of non-NCEP BUFR message.
C*
	PARAMETER	( NXDSC1 = 23 )
	PARAMETER	( NXDSC2 = 25 )
C*
C*	The following array will hold the list of expected descriptors
C*	within Section 3 of each type of non-NCEP BUFR message.
C*	The first dimension of this array must be at least as large
C*	as the largest of the above NXDSC values.
C*
	CHARACTER	cxdsc( NXDSC2, NBUFRT )*6
C*
C*	Expected descriptors within Section 3 of a TAMDAR BUFR message
C*	prior to October 2008.
C*
	DATA		( cxdsc ( ii, 1 ), ii = 1, NXDSC1 )
     +          / '301011', '301013', '301021', '007002',
     +            '007004', '007010', '010082', '008004', '001008',
     +		  '204002', '031021', '001013',
     +            '012001', '011001', '011002', '011037',
     +            '011038', '020041', '204000', '204007',
     +            '031021', '013009', '204000'/
C*
C*	Expected descriptors within Section 3 of a TAMDAR BUFR message
C*	after October 2008.
C*
	DATA		( cxdsc ( ii, 2 ), ii = 1, NXDSC2 )
     +          / '301011', '301013', '301021', '007002',
     +            '007004', '007010', '010082', '008004', '001008',
     +		  '001095', '001009',
     +		  '204002', '031021', '001013',
     +            '012001', '011001', '011002', '011037',
     +            '011038', '020041', '204000', '204007',
     +            '031021', '013009', '204000'/
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
C
        bufrtf(1) = bufrta_old
        nxdsc(1) = NXDSC1
C
        bufrtf(2) = bufrta_new
        nxdsc(2) = NXDSC2

C
C*  Set the pointers for the interface arrays.
C
	CALL TA_IFSP  ( rimnem, cimnem, ierfsp )
	IF  ( ierfsp .ne. 0 )  THEN
	    RETURN
	END IF

C
C++++	BUFR input setup
C
C*  Loop on each type of non-NCEP BUFR tables file.
C
	DO ii = 1, NBUFRT
C
C*  Open the tables file.
C
	    CALL FL_TBOP  ( bufrtf(ii), 'bufr', iubftf(ii), ierspn )
	    IF  ( ierspn .ne. 0 )  THEN
	        CALL DC_WLOG  ( 0, 'FL', ierspn, bufrtf(ii), ierwlg )
	        RETURN
	    END IF
C
C*  Open the messages file.
C
	    CALL FL_GLUN  ( iubfmf(ii), iergln )
	    IF  ( iergln .ne. 0 )  THEN
	        CALL DC_WLOG  ( 0, 'FL', iergln, ' ', ierwlg )
	        RETURN
	    END IF
	    WRITE ( UNIT = cbf, FMT = '(I1.1)')  ii
	    OPEN  ( UNIT = iubfmf(ii),
     +		    FILE = '.dummy.dctama_' // cbf,
     +		    FORM = 'UNFORMATTED' )
C
C*  Connect the tables and messages files.
C
	    CALL OPENBF  ( iubfmf(ii), 'IN', iubftf(ii) )
C
C*  Close the tables file.
C
            CALL FL_CLOS  ( iubftf(ii), iercls )
	    IF  ( iercls .ne. 0 )  THEN
		CALL DC_WLOG  ( 0, 'FL', iercls, ' ', ierwlg )
	    END IF
C
	END DO

C
C*  Constants for SN_CRFP
C
	iflsrc = 4
	maxstn = NACRN
	maxtim = 24
	stmflg = .true.
	addstn = .true.
	cirflg = .true.
	prmfil = 'sntdar.pack'

C
C======================================================================
C  Process bulltins...
C
	DO WHILE ( .true. )
C
C*  Get a new bulletin from the input pipe.
C
	  CALL DC_GBUL ( bull, lenb, ifdtyp, iergbl )
	  IF ( iergbl .ne. 0 ) THEN
C
C*  A time-out occurred while waiting for a new bulletin
C*  on the input pipe.  Shut down the decoder and exit.
C
	      CALL DC_WLOG  ( 0, 'DC', iergbl, ' ', ierwlg )
		DO ii = 1, NBUFRT
		    CALL CLOSBF  ( iubfmf (ii) )
		END DO
		CALL FL_CLAL ( iercal )
	      RETURN
	  END IF
C
C*  Decode the header information from this bulletin.
C
	  IF  ( ifdtyp .ne. 0 )  THEN
C
C*  Do not decode AFOS products.
C
	    bullok = .false.
	  ELSE
C
C*  Decode WMO products.  Initialize the arrays.
C
	    nrept = 0
	    naircft = 0
	    bullok = .true.
	    DO na = 1, NACRN
	      iacrn( na ) = IMISSD
	      selv( na ) = RMISSD
	      acrn( na ) = ' '
	      stat( na ) = '--'
	      coun( na ) = '--'
	      nmlev( na ) = 0
	      flat( na ) = RMISSD
	      flon( na ) = RMISSD
	      DO nl = 1, nlv
		poaf( na, nl ) = RMISSD
		alat( na, nl ) = RMISSD
		alon( na, nl ) = RMISSD
		wspd( na, nl ) = RMISSD
		wdir( na, nl ) = RMISSD
		hmdt( na, nl ) = RMISSD
		tmpk( na, nl ) = RMISSD
		hght( na, nl ) = RMISSD
		pres( na, nl ) = RMISSD
		afdate(na, nl) = ' '
	      END DO
	    END DO
C
C*  Decode WMO products.
C
	    CALL DC_GHDR  ( bull, lenb, seqnum, buhd, cborg,
     +                          bulldt, bbb, ibptr, ierghd )
	    IF  ( ierghd .ne. 0 )  THEN
	      CALL DC_WLOG  ( 2, 'DC', ierghd, ' ', ierwlg )
	      bullok = .false.
	    ELSE
C
C*  Start an entry for this bulletin in the decoder log.
C
	          logmsg = '####################' //
     +			   '####################'
	      CALL DC_WLOG  ( 2, 'DC', 2, logmsg, ierwlg )
	      logmsg = seqnum // buhd // cborg // bulldt // bbb
	      CALL DC_WLOG  ( 2, 'DC', 2, logmsg, ierwlg )
	    END IF
	  END IF
C
	  IF  ( bullok )  THEN
C
C*  Get the system time.
C
	    itype = 1
	    CALL CSS_GTIM  ( itype, sysdt, iergtm )
	    IF  ( iergtm .ne. 0 )  THEN
	      CALL DC_WLOG  ( 2, 'SS', iergtm, ' ', ierwlg )
	      bullok = .false.
	    END IF
	  END IF
C
	  IF  ( bullok )  THEN
C
C*  If a date-time was entered on the command line, then
C*  use it as the run date-time.  Otherwise, use the
C*  system time as the run date-time.
C
	    IF  ( cldt .eq. 'SYSTEM' )  THEN
	      rundt = sysdt
	    ELSE
	      CALL TI_STAN  ( cldt, sysdt, rundt, ierstn )
	      IF  ( ierstn .ne. 0 )  THEN
	        CALL DC_WLOG  ( 2, 'TI', ierstn, ' ', ierwlg )
	        bullok = .false.
	      END IF
	    END IF
	  END IF
C
	  IF  ( bullok )  THEN
C
C*  Convert the run date-time to integer.
C
	    rundt = rundt(1:9) // '00'
	    CALL TI_CTOI  ( rundt, irundt, ierctoi )
	    IF  ( ierctoi .ne. 0 )  THEN
	      CALL DC_WLOG  ( 2, 'TI', ierctoi, ' ', ierwlg )
	      bullok = .false.
	    END IF
	  END IF
C
	  IF ( bullok ) THEN
C
C*  Locate the next BUFR message within the bulletin,
C*  and store it within an equivalenced integer array.
C
	    ipt1 = INDEX ( bull ( ibptr : lenb ), 'BUFR' )
	    IF  ( ipt1 .eq. 0 )  THEN
	      bullok = .false.
	    ELSE
	      istart = ibptr + ipt1 - 1
	      cbull = bull ( istart : lenb )
	    END IF
	  END IF
C
          IF ( bullok ) THEN
C
C*  Determine which BUFR tables file this message
C*  corresponds to.  Retrieve the Section 3 descriptors
C*  from this TAMDAR BUFR message and compare it with
C*  the list of expected descriptors.
C
		ii = 1
		bullok = .false.
		DO WHILE  (  ( .not. bullok ) .and.
     +				( ii .le. NBUFRT )  )
		    CALL UT_CBS3  ( 4, ibull, cxdsc ( 1, ii ),
     +				    nxdsc ( ii ), iercs3 )
		    IF  ( iercs3 .eq. 0 )  THEN
			bullok = .true.
		    ELSE
			ii = ii + 1
		    END IF
		END DO
	  END IF
C
	  IF ( bullok ) THEN
		logmsg = 'message read using ' // bufrtf(ii)
		CALL DC_WLOG  ( 2, 'DC', 2, logmsg, ierwlg )
C
C*  Open the TAMDAR BUFR message.
C
		CALL READERME  ( ibull, iubfmf (ii),
     +				 bfstyp, ibfdt, ierrme )
		IF  ( ierrme .ne. 0 )  THEN
		    bullok = .false.
		ELSE
C
C*  Check update sequence number in Section 1 to see
C*  if correction is indicated.  If this bulletin is
C*  a correction, 'corn' will be set to a non-zero
C*  value so CORN may be set within TA_BUFR.
C
		    corn = FLOAT ( IUPBS01 ( ibull, 'USN' ) )
		    IF ( corn  .gt. 0.0 ) corn = 1.0
		    nrept = 0
		END IF
	  ELSE
		logmsg = 'message has unknown format'
		CALL DC_WLOG  ( 2, 'DC', 2, logmsg, ierwlg )
	  END IF
C
	  DO WHILE  ( bullok )
C
C*  Get the next report from this BUFR message.
C
	    CALL READSB  ( iubfmf(ii), ierrsb )
	    IF  ( ierrsb .ne. 0 )  THEN
C
C*  There are no more reports in this message.
C
	      bullok = .false.
C
C*  Print a count of the number of reports processed.
C
	      WRITE  ( UNIT = logmsg, FMT = '( A, I4, A )' )
     +		'contained ', nrept, ' reports'
	      CALL DC_WLOG  ( 2, 'DC', 2, logmsg, ierwlg )

C
C*  READSB - ierrsb = 0...
C
	    ELSE
	      nrept = nrept + 1
C
C*  Initialize the interface arrays.
C
	      CALL TA_IFIV  ( ier_ifiv )
C
C*  Decode the report into the interface arrays.
C
	      CALL TA_BFIF  ( iubfmf (ii) , ierbif )
C
C*  Write data for this report to the decoder log.
C
	      CALL TA_IFPT  ( rimnem, cimnem, ier_ifpt )
C
C*  Do not create output for reports that are more than
C*  NHOURS before or more than 3 hours after the run time.
C
	      IF  (  ( ERMISS ( rivals( iryear ) ) ) .or.
     +			 ( ERMISS ( rivals( irmnth ) ) ) .or.
     +			 ( ERMISS ( rivals( irdays ) ) ) .or.
     +			 ( ERMISS ( rivals( irhour ) ) ) .or.
     +			 ( ERMISS ( rivals( irminu ) ) )  ) THEN
		    iertmk = -1
	      ELSE
		irptdt(1) = INT ( rivals( iryear ) )
		irptdt(2) = INT ( rivals( irmnth ) )
		irptdt(3) = INT ( rivals( irdays ) )
		irptdt(4) = INT ( rivals( irhour ) )
		irptdt(5) = INT ( rivals( irminu ) )
		CALL DC_TMCK  ( 2, irundt, irptdt, nhours,
     +			    180, iertmk )
	      END IF
	      IF  ( iertmk .eq. 0 )  THEN
C
C*  Construct the aircraft id from aircraft number
C*  by appending flight phase character:
C*  - 'A' for ascending plane
C*  - 'D' for descending plane
C
		CALL ST_NUMB (civals ( icacrn ), icurra, iret)
		IF ( rivals (irpoaf) .eq. 5 ) THEN
		  cracrn = civals ( icacrn )(2:) // 'A'
		ELSE IF ( rivals (irpoaf) .eq. 6 ) THEN
		  cracrn = civals ( icacrn )(2:) // 'D'
		ELSE
		  cracrn = ' '
		ENDIF
		amatch = .false.
		ia    = 1
		iaircft = 0
C
C*  Start loop over the planes in the bulletin.
C
		DO WHILE ( .not. amatch )
		  IF ( ( cracrn .eq. acrn(ia) ) .and.
     +                 ( nmlev(ia) .ne. 0 ) ) THEN
C
C*  This is an existing plane - increment the number of levels.
C
		    amatch = .true.
		    nmlev(ia) = nmlev(ia) + 1
		    iaircft = ia
		  ELSE IF ( ( cracrn .ne. acrn(ia) ) .and.
     +                           ( rivals(irpoaf) .ne. RMISSD ) .and.
     +                           ( nmlev( ia ) .eq. 0 ) ) THEN
C
C*  This is a new plane - increment the total
C*  number of aircrafts and construct the aircraft number.
C
		    naircft = naircft + 1
		    acrn(naircft) = cracrn
		    IF ( cracrn(8:8) .eq. 'D') THEN
		      ipoaf = 6
		    ELSE IF (cracrn(8:8) .eq. 'A') THEN
		      ipoaf = 5
		    END IF
		    iacrn(naircft) = ipoaf*100000 + icurra
		    nmlev(naircft) = 1
		    iaircft = naircft
		    amatch = .true.
		  END IF
		  ia = ia + 1
		  IF ( ia .gt. NACRN ) amatch = .true.
		END DO
C
C*  Decode the records into the individual arrays.
C
		IF ( iaircft .gt. 0 ) THEN
		  il = nmlev(iaircft)
		  alat(iaircft, il) = rivals(irslat)
		  alon(iaircft, il) = rivals(irslon)
		  wspd(iaircft, il) = rivals(irwspd)
		  wdir(iaircft, il) = rivals(irwdir)
		  hght(iaircft, il) = rivals(irhmsl)
		  pres(iaircft, il) = rivals(irprlc)/100.
		  tmpk(iaircft, il) = rivals(irtmdb)
		  hmdt(iaircft, il) = rivals(irrehu)
		  CALL TI_ITOC ( irptdt, cdum, iret)
		  afdate(iaircft, il) = cdum
		END IF
	      END IF
C
		  logmsg = '-----------------------------------'
		  CALL DC_WLOG  ( 3, 'DC', 2, logmsg, ierwlg )
	    END IF

C
C*  END DO for DO WHILE ( bullok );
C*  looping over reports in a single BUFR message...
C
	  END DO

C
C======================================================================
C
C*  Write out the data for each aircraft for this BUFR message...
C
	  DO ia = 1, naircft
	    IF ( acrn(ia)(8:8) .eq. 'A' ) THEN
	      bot = 1
	      top = nmlev(ia)
	      inc = 1
	    ELSE
	      bot = nmlev(ia)
	      top = 1
	      inc = -1
	    END IF

C
C*  Create a "Nominal Time" of observation.
C
	    dattim = afdate(ia,bot)(1:9) // '00'

C
C*  Create a unique "observation ID" from time and aircraft ID.
C
	    acrn(ia) = dattim(8:9) // acrn(ia)(3:8)

C
C*  Make a file name from the template and the aircraft report time.
C
	    CALL FL_MNAM ( dattim(1:6), gemfil, filnam, ier )

C
C*  Check for the requested upper air data file.
C
	    CALL FL_INQR ( filnam, exist, newfil, ier )
	    IF  ( exist )  THEN

C
C*  Open the merged upper air data file.
C
	      CALL SN_OPNR ( filnam, isnfln, ifsc, nparm, parms, ivert,
     +				mrgdat, ieropnr )

C
C*  Check that the file was opened properly.
C
	      IF  ( ieropnr .ne. 0 )  THEN

C
C*  If not, write an error to the decoder log file.
C
	        CALL DC_WLOG ( 0, 'SN', ierr, filnam, iret )
	        RETURN
	      ELSE

C
C*  Check for existence of time/aircraft in file; if not present
C*  add the time/aircraft; if present, has data been added?
C
		CALL TA_TMST ( isnfln, dattim , acrn(ia),
     +				iacrn(ia), alat(ia,bot), alon(ia,bot),
     +				RMISSD, addstn, cirflg, datflg, ier_tmst)
		IF ( ier_tmst .eq. 0 ) THEN
		  IF ( datflg ) THEN

C
C*  Data already present in sounding file, read it in.
C
		    CALL SN_RDAT(isnfln, nlevi, datai, ihhmm, ier)
		    nlevo = nlevi + nmlev(ia)
		    IF ( acrn(ia)(8:8) .eq. 'A' ) THEN

C
C*  Ascending aircraft so can insert new data at the end of datai.
C
		      ii = nlevi * NPARMS
		      DO il = 1, nmlev(ia)
			datai(ii+1) = pres(ia,il)
			datai(ii+2) = tmpk(ia,il)
			datai(ii+3) = hmdt(ia,il)

C
C*  It appears that wdir equal to 360.0 is reserved for missing...
C
			IF ( wdir(ia,il) .eq. 360.0 )  THEN
			  datai(ii+4) = RMISSD
			  datai(ii+5) = RMISSD
			ELSE
			  datai(ii+4) = wdir(ia,il)
			  datai(ii+5) = wspd(ia,il)
			END IF
			datai(ii+6) = hght(ia,il)
			ii = ii + 6
		      END DO
		      CALL SN_WDAT(isnfln, ihhmm, nlevo, datai, ier)
		    ELSE

C
C*  Descending aircraft so must put new data at head of datai.
C*  Use datao for the output data.
C
		      ii = 0
		      DO il = bot, top, inc
			datao(ii+1) = pres(ia,il)
			datao(ii+2) = tmpk(ia,il)
			datao(ii+3) = hmdt(ia,il)

C
C*  It appears that wdir equal to 360.0 is reserved for missing...
C
			IF ( wdir(ia,il) .eq. 360.0 )  THEN
			  datao(ii+4) = RMISSD
			  datao(ii+5) = RMISSD
			ELSE
			  datao(ii+4) = wdir(ia,il)
			  datao(ii+5) = wspd(ia,il)
			END IF
			datao(ii+6) = hght(ia,il)
			ii = ii + 6
		      END DO
		      kk = 1
		      DO jj = nmlev(ia)*NPARMS+1, nlevo*NPARMS
			datao(jj) = datai(kk)
			kk = kk +1
		      END DO
		    END IF

C
C*  ELSE for "IF ( datflg ) THEN"
C
		  ELSE

C
C*  New aircraft for this time period, just add it...
C
		    CALL TI_CTOI(afdate(ia,bot), iarcft, ierctoi )
		    ihhmm = 100*iarcft(4) + iarcft(5)
		    ii = 0
		    DO il = bot, top, inc
		      datao(ii+1) = pres(ia,il)
		      datao(ii+2) = tmpk(ia,il)
		      datao(ii+3) = hmdt(ia,il)

C
C*  It appears that wdir equal to 360.0 is reserved for missing...
C
		      IF ( wdir(ia,il) .eq. 360.0 )  THEN
		        datao(ii+4) = RMISSD
		        datao(ii+5) = RMISSD
		      ELSE
		        datao(ii+4) = wdir(ia,il)
		        datao(ii+5) = wspd(ia,il)
		      END IF
		      datao(ii+6) = hght(ia,il)
		      ii = ii + 6
		    END DO
		    CALL SN_WDAT(isnfln, ihhmm, nmlev(ia), datao, ier)
		  END IF

C
C*  END IF for "IF ( iertmst .eq. 0 ) THEN"
C
		END IF
		CALL SN_CLOS(isnfln, ierclos)

C
C*  END IF for "IF  ( ieropnr .ne. 0 )  THEN"
C
	      END IF

C
C*  ELSE for "IF  ( exist )  THEN"
C
	    ELSE

C
C*  File does NOT exist, create the merged upper air data file.
C
	      CALL SN_CRFP ( filnam, prmfil, iflsrc, maxstn, maxtim,
     +               stmflg, isnfln, nparm, parms, pkflg, iercrfp )

C
C*  Check that the file was created properly.
C
              IF  ( iercrfp .ne. 0 )  THEN

C
C*  If not, write an error to the decoder log file.
C
                CALL DC_WLOG ( 0, 'SN', ierr, filnam, iret )
                RETURN
	      ELSE

C
C*  Add station to the data file.
C
		CALL TA_TMST ( isnfln, dattim, acrn(ia),
     +				iacrn(ia), alat(ia,bot), alon(ia,bot),
     +				RMISSD, addstn, cirflg, datflg, iertmst )
		IF ( iertmst .eq. 0 ) THEN
		  CALL TI_CTOI(afdate(ia,bot), iarcft, ierctoi )
		  ihhmm = 100*iarcft(4) + iarcft(5)
		  ii = 0
		  DO il = bot, top, inc
		    datao(ii+1) = pres(ia,il)
		    datao(ii+2) = tmpk(ia,il)
		    datao(ii+3) = hmdt(ia,il)

C
C*  It appears that wdir equal to 360.0 is reserved for missing...
C
		    IF ( wdir(ia,il) .eq. 360.0 )  THEN
		      datao(ii+4) = RMISSD
		      datao(ii+5) = RMISSD
		    ELSE
		      datao(ii+4) = wdir(ia,il)
		      datao(ii+5) = wspd(ia,il)
		    END IF
		    datao(ii+6) = hght(ia,il)
		    ii = ii + 6
		  END DO
		  CALL SN_WDAT (isnfln, ihhmm, nmlev(ia), datao, ierwdat)
		END IF
		CALL SN_CLOS(isnfln, ierclos)
	      END IF

C
C*  END IF for "IF  ( exist )  THEN"
C
	    END IF

C
C*  END DO for "DO ia = 1, naircft"
C
	  END DO

C
C*  END DO for "DO WHILE ( .true. )"; looping over BUFR bulltins...
C
	END DO
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
c> \file errclb.f

c> \brief This subroutine checks the validity of the input data.
c>
c> This subroutine checks the validity of the input data.
c>
c> @param n number of parameters
c> @param m history size of approximated Hessian
c> @param factr convergence criterion on function value
c> @param l lower bounds for parameters
c> @param u upper bounds for parameters
c> @param nbd indicates which bounds are present
c> @param task if an error occurs, contains a human-readable error message
c> @param info =0 on success; =-6 if nbd(k) was invalid; =-7 if both limits are given but l(k) > u(k)
c> @param k index of last errournous parameter
      subroutine errclb(n, m, factr, l, u, nbd, task, info, k)

      character*60     task
      integer          n, m, info, k, nbd(n)
      double precision factr, l(n), u(n)
c
c                           *  *  *
c
c     NEOS, November 1994. (Latest revision June 1996.)
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer          i
      double precision one,zero
      parameter        (one=1.0d0,zero=0.0d0)

c     Check the input arguments for errors.

      if (n .le. 0) task = 'ERROR: N .LE. 0'
      if (m .le. 0) task = 'ERROR: M .LE. 0'
      if (factr .lt. zero) task = 'ERROR: FACTR .LT. 0'

c     Check the validity of the arrays nbd(i), u(i), and l(i).

      do 10 i = 1, n
         if (nbd(i) .lt. 0 .or. nbd(i) .gt. 3) then
c                                                   return
            task = 'ERROR: INVALID NBD'
            info = -6
            k = i
         endif
         if (nbd(i) .eq. 2) then
            if (l(i) .gt. u(i)) then
c                                    return
               task = 'ERROR: NO FEASIBLE SOLUTION'
               info = -7
               k = i
            endif
         endif
  10  continue

      return

      end

max_stars_repo_name: jonathanschilling/L-BFGS-B

max_issues_repo_name: jonathanschilling/L-BFGS-B

max_forks_repo_name: jonathanschilling/L-BFGS-B


text:
c-------------------------------------------------------
	implicit none
	include 'netcdf.inc'
	include 'hbuf.inc'

	integer flag(HBUF_MAX_LENGTH)
c-------------------------------------------------------
	character caseid*80, filename*88, filename2*88
	character name*80
	integer npres, ntime,ntau
	real time(10000)
        real, allocatable ::  fq_isccp(:,:,:)  !accumulated fq_isccp
        real, allocatable ::  taulim(:), preslim(:)

        real, allocatable ::  totalcldarea(:),meantaucld(:)
        real, allocatable ::  lowcldarea(:),midcldarea(:),hghcldarea(:)
        real, allocatable ::  meanptop(:)

	integer i,j,k,l,m,nparms
	integer vdimids(4), ndimids
	integer ncid,err,tid,zid,timeid,varid
        character(80) long_name
        character(30) abbr_name
        character(12) units

	integer iargc
	external iargc

	time=0.
c-------------------------------------------------------
	m=COMMAND_ARGUMENT_COUNT()
	if(m.eq.0.) then
	 print*,'you forgot to specify the name of the file.'
	 stop
	endif
	call getarg(1,filename)
	print *,'open file: ',filename
	open(2,file=filename,status='old',form='unformatted')

	do i=1,83
	 if(filename(i:i+5).eq.'.isccp') then
	  filename2=filename(1:i-1)//'.isccp.nc'
	  print*,filename2
	 endif
	end do
c-------------------------------------------------------
c Count how many time points in a file:

        ntime = 0
        do while(.true.)
         read(2,end=111)  caseid
         ntime=ntime+1
         read (2) time(ntime)
         read (2) ntau,npres
         read (2)
         read (2)
         read (2)
         read (2)
         read (2)
	 print*,time(ntime)
        end do
 111    continue
        rewind(2)

	print*,caseid
	print*,'ntau=',ntau
	print*,'npres=',npres
	print*,'ntime=',ntime

	allocate (preslim(npres), taulim(ntau))
	allocate (fq_isccp(ntau,npres,ntime))
        allocate (totalcldarea(ntime))
        allocate (lowcldarea(ntime),midcldarea(ntime),hghcldarea(ntime))
        allocate (meantaucld(ntime),meanptop(ntime))

	do i=1,ntime
	 read(2)
	 read(2) time(i)
         read(2) ntau,npres
         read(2) taulim(1:ntau)
         read(2) preslim(1:npres)
         read(2) fq_isccp(1:ntau,1:npres,i)
         read(2) totalcldarea(i),
     &            lowcldarea(i),midcldarea(i),hghcldarea(i)
         read(2) meantaucld(i),meanptop(i)
        end do

	err = NF_CREATE(filename2, NF_CLOBBER, ncid)

	err = NF_REDEF(ncid)

	err = NF_PUT_ATT_TEXT(ncid,NF_GLOBAL,'model',7,'CSU CEM')
	err = NF_PUT_ATT_TEXT(ncid,NF_GLOBAL,'caseid',
     &                                  len_trim(caseid),caseid)
	err = NF_PUT_ATT_TEXT(ncid,NF_GLOBAL,'MissingValue',2,'-1')

	err = NF_DEF_DIM(ncid, 'time', ntime, timeid)
	err = NF_DEF_DIM(ncid, 'pres', npres, zid)
	err = NF_DEF_DIM(ncid, 'tau', ntau, tid)

	err = NF_DEF_VAR(ncid, 'pres', NF_FLOAT, 1, zid, varid)
	err = NF_PUT_ATT_TEXT(ncid,varid,'units',3,'hPa')
	err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',7,'pressure')

	err = NF_DEF_VAR(ncid, 'time', NF_FLOAT, 1, timeid, varid)
	err = NF_PUT_ATT_TEXT(ncid,varid,'units',3,'day')
	err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',4,'time')

cbloss(2016-02-11): Changed dimension of tau to optical depth.
c  This does not usually matter with the ISCCP simulator, where
c  ntau=npres, but seems desirable in general.
	err = NF_DEF_VAR(ncid, 'tau', NF_FLOAT, 1, tid, varid)
	err = NF_PUT_ATT_TEXT(ncid,varid,'units',1,' ')
	err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                                 17,'Optical thickness')

	err = NF_ENDDEF(ncid)

        err = NF_INQ_VARID(ncid,'time',varid)
	err = NF_PUT_VAR_REAL(ncid, varid, time)
        err = NF_INQ_VARID(ncid,'pres',varid)
	err = NF_PUT_VAR_REAL(ncid, varid, preslim)
        err = NF_INQ_VARID(ncid,'tau',varid)
	err = NF_PUT_VAR_REAL(ncid, varid, taulim)
c--------------------------------------------------------

        long_name = 'ISCCP Total Cloud Area tau > 0.3'
        abbr_name = 'CLDTOT'
        units = ' '
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  1, timeid,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                  len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                       len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid, totalcldarea(1:ntime))

c--------------------------------------------------------
        long_name = 'ISCCP Low Cloud Area tau > 0.3'
        abbr_name = 'CLDLOW'
        units = ' '
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  1, timeid,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                  len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                       len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid, lowcldarea(1:ntime))

c--------------------------------------------------------
        long_name = 'ISCCP Middle Cloud Area tau > 0.3'
        abbr_name = 'CLDMID'
        units = ' '
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  1, timeid,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                  len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                       len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid, midcldarea(1:ntime))

c--------------------------------------------------------
        long_name = 'ISCCP High Cloud Area tau > 0.3'
        abbr_name = 'CLDHGH'
        units = ' '
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  1, timeid,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                  len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                       len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid, hghcldarea(1:ntime))

c--------------------------------------------------------
        long_name = 'ISCCP mean optical depth for tau > 0.3'
	abbr_name = 'TAUMEAN'
        units = ' '
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  1, timeid,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                  len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                       len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid, meantaucld(1:ntime))

c--------------------------------------------------------
        long_name = 'ISCCP mean cloud top pressure for tau > 0.3'
	abbr_name = 'PTOPMEAN'
        units = 'hPa'
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  1, timeid,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                  len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                       len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid, meanptop(1:ntime))

c--------------------------------------------------------
        ndimids=3
        vdimids(1) = tid
        vdimids(2) = zid
        vdimids(3) = timeid

        long_name = 'ISCCP Table'
        abbr_name = 'isccp'
        units = ' '
        err = NF_REDEF(ncid)
        err = NF_DEF_VAR(ncid,trim(abbr_name),NF_FLOAT,
     &                                  ndimids,vdimids,varid)
        err = NF_PUT_ATT_TEXT(ncid,varid,'long_name',
     &                   len_trim(long_name),trim(long_name))
        err = NF_PUT_ATT_TEXT(ncid,varid,'units',
     &                   len_trim(units),trim(units))
        err = NF_ENDDEF(ncid)
        err = NF_PUT_VAR_REAL(ncid, varid,
     &                     fq_isccp(1:ntau,1:npres,1:ntime))

	err = NF_CLOSE(ncid)

	end

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


MODULE basic_type_methods

  USE basic_type_factory

  implicit none

  PRIVATE
  PUBLIC :: allocateGS
  PUBLIC :: allocateGSArray, allocateGSArray_v2
  PUBLIC :: allocateGBKS
  PUBLIC :: allocaterGBKS
  PUBLIC :: allocatecGBKS

#ifdef REAL_VER
  double precision,parameter :: zero = 0.d0
#else
  complex(kind(0d0)),parameter :: zero = (0.d0,0.d0)
#endif
  complex(kind(0d0)),parameter :: z0 = (0.d0,0.d0)

CONTAINS

  SUBROUTINE allocateGSArray( gs )
    implicit none
    type( GSArray ) :: gs
    integer :: m1,m2,n1,n2
    m1=gs%g_range%head
    m2=gs%g_range%tail
    n1=gs%s_range%head
    n2=gs%s_range%tail
    allocate( gs%val(m1:m2,n1:n2) ) ; gs%val=0.0d0
  END SUBROUTINE allocateGSArray

  SUBROUTINE allocateGSArray_v2( gs )
    implicit none
    type( GSArray_v2 ) :: gs
    integer :: m1,m2,n1,n2
    m1=gs%g_range%alloc%head
    m2=gs%g_range%alloc%tail
    n1=gs%s_range%alloc%head
    n2=gs%s_range%alloc%tail
    allocate( gs%val(m1:m2,n1:n2) ) ; gs%val=0.0d0
  END SUBROUTINE allocateGSArray_v2

  SUBROUTINE allocateGS( gs )
    implicit none
    type( GSArray    ) ::  gs
    allocate( gs%val(gs%g_prange%head:gs%g_prange%tail, &
                     gs%s_srange%head:gs%s_srange%tail) )
    gs%val = 0.d0
  END SUBROUTINE allocateGS

  SUBROUTINE allocateGBKS( gbks )
    implicit none
    type( GBKSArray  ) ::  gbks
    allocate( gbks%val(gbks%g_prange%head:gbks%g_prange%tail, &
                       gbks%b_prange%head:gbks%b_prange%tail, &
                       gbks%k_prange%head:gbks%k_prange%tail, &
                       gbks%s_prange%head:gbks%s_prange%tail) )
    gbks%val = zero
  END SUBROUTINE allocateGBKS

  SUBROUTINE allocaterGBKS( gbks )
    implicit none
    type( rGBKSArray ) :: gbks
    allocate( gbks%val(gbks%g_prange%head:gbks%g_prange%tail, &
                       gbks%b_prange%head:gbks%b_prange%tail, &
                       gbks%k_prange%head:gbks%k_prange%tail, &
                       gbks%s_prange%head:gbks%s_prange%tail) )
    gbks%val = 0.d0
  END SUBROUTINE allocaterGBKS

  SUBROUTINE allocatecGBKS( gbks )
    implicit none
    type( cGBKSArray ) :: gbks
    allocate( gbks%val(gbks%g_prange%head:gbks%g_prange%tail, &
                       gbks%b_prange%head:gbks%b_prange%tail, &
                       gbks%k_prange%head:gbks%k_prange%tail, &
                       gbks%s_prange%head:gbks%s_prange%tail) )
    gbks%val = z0
  END SUBROUTINE allocatecGBKS

END MODULE basic_type_methods


text:
c*********************************************************************
c
c   Display_Error - Function
c
c   VERSION:  1.0
c   AUTHOR: J. L. Pierce
c   CREATION DATE:  11-Feb-1984
c
c FUNCTION: Displays error message on bottom of screen.
c
c USEAGE:  iret = display_error(message,row,column,flag)
c
c ARGUMENTS:	message - character string descriptor
c			  error message to be displayed
c		row	- integer*4
c			  row to begin display on
c		column	- integer*4
c			  column to begin display on
c		flag	- integer*4
c			  attribute bit vector (see LIB$PUT_SCREEN doc.)
c ERRORS:  returns same error codes as LIB$PUT_SCREEN and LIB$ERASE_LINE
c
c NOTES:  Written in VAX/VMS FORTRAN.
c
c PROCESS: clear a line using LIB$ERASE_LINE as indicated by row and column
c	   call LIB$PUT_SCREEN to display message
c	   call wait for SYS_DEL_TIME - an external parameter
c	   call LIB$ERASE_LINE to clear error message
c	   return
c
c history:
c	  v1.1 09-april-1985 add str$trim
c*********************************************************************
c
	function display_error(error_message,row,column,flag)
c
	implicit none
c
	include 'c3mdef.for'			! master includes
c
	character*(*)	error_message		! text to be displayed
	integer*4	row			! row to begin on
	integer*4	column			! column to begin on
	integer*4	flag			! attribute flag
	integer*4	display_error		! for return code
c
	external	lib$put_screen
	external	lib$erase_line

	integer*4	col			! column to begin writing
	integer*4	olen			! length of the error message
	character	obuff*80
c
c******************************************************************************
c
c	Code
c
c******************************************************************************
c
	display_error = lib$erase_line(row,1)
	if(.not.display_error) return
c
c version 1.1 add str$trim
	call str$trim( obuff, error_message, olen )
	col = 40 - (olen/2)		! center the error message on screen
	display_error = lib$put_screen(obuff(1:olen),row,col,flag)
c version 1.1 end
c
	if(.not.display_error) return
c
c version 1.1 comment out wait and add read
c	call wait(SYS_DEL_TIME)
	read*
c
	call lib$erase_line(row,1)		! no error - has to work again
	return
	end

ext: for


text:
!!
!> @brief      Tools for quadrature
!!
!! Subroutine and functions to perform integrals
!!
!! @author     Rodrigo Navarro Perez
!!
module quadrature
use precisions, only: dp
implicit none

private
public :: booles_quadrature
contains

!!
!> @brief      Booles quadrature
!!
!! Given an array containing equidistant evaluations of a function and
!! the distance between evaluation points, uses the Booles rule to
!! calculate the integral of the evaluated function.
!!
!! Given the number of points necessary for the Booles rule, the
!! size of the array has to be a multiple of 4 minus 1
!!
!! @return     Integral of a function
!!
!! @author     Rodrigo Navarro Perez
!!
real(dp) function booles_quadrature(fx, delta_x) result(s)
    implicit none
    real(dp), intent(in), dimension(1:) :: fx !< Array of equidistant evaluations of a function
    real(dp), intent(in) :: delta_x !< distance between evaluation points

    integer :: fx_size, i

    fx_size = size(fx)

    if (mod(fx_size-1,4) /= 0) then
        print *, 'fx array size in booles_quadrature has to be a multiple of 4 plus 1'
        stop
    endif

    s = 0._dp

    do i=1,fx_size-1,4
        s = s + booles_rule(fx(i:i+4), delta_x)
    enddo
end function booles_quadrature

!!
!> @brief      Booles rule
!!
!! Given an array containing 5 equidistant evaluations of a function and
!! the distance between evaluation points, calculates Booles rule using those 5 points.
!!
!! @return     Booles rule for 5 points
!!
!! @author     Rodrigo Navarro Perez
!!
real(dp) function booles_rule(fx, delta_x) result(s)
    implicit none
    real(dp), intent(in), dimension(1:) :: fx !< Array of 5 equidistant evaluations of a function
    real(dp), intent(in) :: delta_x !< distance between evaluation points

    real(dp) :: fx0, fx1, fx2, fx3, fx4

    if( size(fx) /= 5 ) then
        print *, 'fx array size in booles_rule has to be 5'
        stop
    endif

    fx0 = fx(1)
    fx1 = fx(2)
    fx2 = fx(3)
    fx3 = fx(4)
    fx4 = fx(5)

    s = delta_x*2/45._dp*(7*fx4 + 32*fx3 + 12*fx2 + 32*fx1 + 7*fx0)
end function booles_rule

end module quadrature

max_stars_repo_name: rnavarroperez/nn-scattering-fit

max_issues_repo_name: rnavarroperez/nn-scattering-fit

max_forks_repo_name: rnavarroperez/nn-scattering-fit


text:
! Riemann solver for the LWR traffic model:

! q_t + (u_max q(1-q))_x = 0

! waves: 1
! equations: 1

! Conserved quantities:
!       1 q

! The speed limit umax should be in the common block cparam.

! See http://www.clawpack.org/riemann.html for a detailed explanation
! of the Riemann solver API.

subroutine rp1(maxm,num_eqn,num_waves,num_aux,num_ghost,num_cells, &
               ql,qr,auxl,auxr,wave,s,amdq,apdq)

    implicit none

    ! Inputs
    integer, intent(in) :: maxm, num_eqn, num_waves, num_aux, num_ghost, num_cells
    double precision, intent(in), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: ql,qr
    double precision, intent(in), dimension(num_aux, 1-num_ghost:maxm+num_ghost) :: auxl,auxr
    ! Outputs
    double precision, intent(out) :: s(num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out) :: wave(num_eqn,num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: amdq,apdq
    ! Locals
    double precision :: sim1, si, flux0, fluxim1, fluxi
    integer :: i

    common /cparam/ umax
    double precision :: umax

    do i=2-num_ghost,num_cells+num_ghost
        ! Compute the wave and speed

        wave(1,1,i) = ql(1,i) - qr(1,i-1)
        s(1,i) = umax * (1.d0 - (qr(1,i-1) + ql(1,i)))

        ! compute left-going and right-going flux differences:
        !------------------------------------------------------
        amdq(1,i) = dmin1(s(1,i), 0.d0) * wave(1,1,i)
        apdq(1,i) = dmax1(s(1,i), 0.d0) * wave(1,1,i)

        ! entropy fix for transonic rarefactions:
        sim1 = umax*(1.d0 - 2.d0*qr(1,i-1))
        si = umax*(1.d0 - 2.d0*ql(1,i))
        if (sim1.lt.0.d0 .and. si.gt.0.d0) then
            flux0 = 0.25d0*umax
            fluxim1 = qr(1,i-1)*umax*(1.d0 - qr(1,i-1))
            fluxi   = ql(1,i)*umax*(1.d0 - ql(1,i))
            amdq(1,i) = flux0 - fluxim1
            apdq(1,i) = fluxi - flux0
        endif
    enddo

    return
end subroutine rp1

max_stars_repo_name: ranocha/riemann

max_issues_repo_name: ranocha/riemann

max_forks_repo_name: ranocha/riemann


text:
	module foo
	integer a, b, c
	end module
	program bar
	use foo, my_a=>a, my_c=>c
	implicit none
	a = 1
	my_a = b
	end

max_stars_repo_name: kurtsansom/ftnchek

max_issues_repo_name: kurtsansom/ftnchek

max_forks_repo_name: kurtsansom/ftnchek


text:
	SUBROUTINE IS_BDLN ( clat1, clat2, clon1, clon2, rlat, rlon,
     +			     npt, iret )
C************************************************************************
C* IS_BDLN 								*
C*									*
C* This subroutine decodes a bounded area of the form lat1 TO lat2 FROM *
C* lon1 TO lon2.                                                        *
C*                                                                      *
C* IS_BDLN ( CLAT1, CLAT2, CLON1, CLON2, RLAT, RLON, NPT, IRET )        *
C*									*
C* Input parameters:							*
C*	CLAT1		CHAR*		First latitude value            *
C*	CLAT2		CHAR*		Second latitude value           *
C*	CLON1		CHAR*		First longitude value           *
C*	CLON2		CHAR*		Second longitude value          *
C*									*
C* Output parameters:							*
C*	RLAT (*)	REAL  		Signed latitude values          *
C*	RLON (*)	REAL   		Signed longitude values         *
C*	NPT		INTEGER		Number of values		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = invalid format            *
C*									*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	11/99	                                        *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	clat1, clat2, clon1, clon2
	REAL		rlat (*), rlon (*)
C*
	CHARACTER 	work*20
C------------------------------------------------------------------------
	iret = 0
	npt  = 4
C
	CALL ST_LSTR ( clat1, lens, ier )
	work = clat1 ( :lens ) // clon1
	CALL IS_LTLN ( work, rlat ( 1 ), rlon ( 1 ), ier1 )
	work = clat1 ( :lens ) // clon2
	CALL IS_LTLN ( work, rlat ( 2 ), rlon ( 2 ), ier2 )
C
	CALL ST_LSTR ( clat2, lens, ier )
	work = clat2 ( :lens ) // clon2
	CALL IS_LTLN ( work, rlat ( 3 ), rlon ( 3 ), ier3 )
	work = clat2 ( :lens ) // clon1
	CALL IS_LTLN ( work, rlat ( 4 ), rlon ( 4 ), ier4 )
C
	IF ( ( ier1 + ier2 + ier3 + ier4 ) .ne. 0 ) iret = -4
C*
	RETURN
	END


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:38 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE TSTIPA__genmod
          INTERFACE
            SUBROUTINE TSTIPA
            END SUBROUTINE TSTIPA
          END INTERFACE
        END MODULE TSTIPA__genmod


text:
SUBROUTINE rod_mm(mm,length)
!
! This subroutine forms the consistent mass matrix of a 1-d "rod" element.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),intent(in)::length
 REAL(iwp),intent(out)::mm(:,:)
 REAL(iwp)::one=1.0_iwp,d3=3.0_iwp,d6=6.0_iwp
 mm(1,1)=one/d3
 mm(1,2)=one/d6
 mm(2,1)=one/d6
 mm(2,2)=one/d3
 mm=mm*length
RETURN
END SUBROUTINE rod_mm

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM