text:
C$Procedure      M2EUL ( Matrix to Euler angles )

      SUBROUTINE M2EUL (  R,  AXIS3,   AXIS2,   AXIS1,
     .                        ANGLE3,  ANGLE2,  ANGLE1  )

C$ Abstract
C
C     Factor a rotation matrix as a product of three rotations about
C     specified coordinate axes.
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     ROTATION
C
C$ Keywords
C
C     ANGLE
C     MATRIX
C     ROTATION
C     TRANSFORMATION
C
C$ Declarations

      DOUBLE PRECISION      R      ( 3, 3 )

      INTEGER               AXIS3
      INTEGER               AXIS2
      INTEGER               AXIS1

      DOUBLE PRECISION      ANGLE3
      DOUBLE PRECISION      ANGLE2
      DOUBLE PRECISION      ANGLE1

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     R          I   A rotation matrix to be factored.
C     AXIS3,
C     AXIS2,
C     AXIS1      I   Numbers of third, second, and first rotation axes.
C     ANGLE3,
C     ANGLE2,
C     ANGLE1     O   Third, second, and first Euler angles, in radians.
C
C$ Detailed_Input
C
C     R              is a 3x3 rotation matrix that is to be factored as
C                    a product of three rotations about a specified
C                    coordinate axes.  The angles of these rotations are
C                    called `Euler angles'.
C
C     AXIS3,
C     AXIS2,
C     AXIS1          are the indices of the rotation axes of the
C                    `factor' rotations, whose product is R.  R is
C                    factored as
C
C                       R = [ ANGLE3 ]     [ ANGLE2 ]     [ ANGLE1 ]   .
C                                    AXIS3          AXIS2         AXIS1
C
C                    The axis numbers must belong to the set {1, 2, 3}.
C                    The second axis number MUST differ from the first
C                    and third axis numbers.
C
C                    See the $ Particulars section below for details
C                    concerning this notation.
C
C$ Detailed_Output
C
C     ANGLE3,
C     ANGLE2,
C     ANGLE1         are the Euler angles corresponding to the matrix
C                    R and the axes specified by AXIS3, AXIS2, and
C                    AXIS1.  These angles satisfy the equality
C
C                       R = [ ANGLE3 ]     [ ANGLE2 ]     [ ANGLE1 ]
C                                   AXIS3          AXIS2          AXIS1
C
C
C                    See the $ Particulars section below for details
C                    concerning this notation.
C
C                    The range of ANGLE3 and ANGLE1 is (-pi, pi].
C
C                    The range of ANGLE2 depends on the exact set of
C                    axes used for the factorization.  For
C                    factorizations in which the first and third axes
C                    are the same,
C
C                       R = [r]  [s]  [t] ,
C                              a    b    a
C
C                    the range of ANGLE2 is [0, pi].
C
C
C                    For factorizations in which the first and third
C                    axes are different,
C
C                       R = [r]  [s]  [t] ,
C                              a    b    c
C
C                    the range of ANGLE2 is [-pi/2, pi/2].
C
C                    For rotations such that ANGLE3 and ANGLE1 are not
C                    uniquely determined, ANGLE3 will always be set to
C                    zero; ANGLE1 is then uniquely determined.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)   If any of AXIS3, AXIS2, or AXIS1 do not have values in
C
C             { 1, 2, 3 },
C
C          then the error SPICE(BADAXISNUMBERS) is signaled.
C
C     2)   An arbitrary rotation matrix cannot be expressed using
C          a sequence of Euler angles unless the second rotation axis
C          differs from the other two.  If AXIS2 is equal to AXIS3 or
C          AXIS1, then then error SPICE(BADAXISNUMBERS) is signaled.
C
C     3)   If the input matrix R is not a rotation matrix, the error
C          SPICE(NOTAROTATION) is signaled.
C
C     4)   If ANGLE3 and ANGLE1 are not uniquely determined, ANGLE3
C          is set to zero, and ANGLE1 is determined.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     A word about notation:  the symbol
C
C        [ x ]
C             i
C
C     indicates a coordinate system rotation of x radians about the
C     ith coordinate axis.  To be specific, the symbol
C
C        [ x ]
C             1
C
C     indicates a coordinate system rotation of x radians about the
C     first, or x-, axis; the corresponding matrix is
C
C        +-                    -+
C        |  1      0       0    |
C        |                      |
C        |  0    cos(x)  sin(x) |.
C        |                      |
C        |  0   -sin(x)  cos(x) |
C        +-                    -+
C
C     Remember, this is a COORDINATE SYSTEM rotation by x radians; this
C     matrix, when applied to a vector, rotates the vector by -x
C     radians, not x radians.  Applying the matrix to a vector yields
C     the vector's representation relative to the rotated coordinate
C     system.
C
C     The analogous rotation about the second, or y-, axis is
C     represented by
C
C        [ x ]
C             2
C
C     which symbolizes the matrix
C
C        +-                    -+
C        | cos(x)   0   -sin(x) |
C        |                      |
C        |  0       1      0    |,
C        |                      |
C        | sin(x)   0    cos(x) |
C        +-                    -+
C
C     and the analogous rotation about the third, or z-, axis is
C     represented by
C
C        [ x ]
C             3
C
C     which symbolizes the matrix
C
C        +-                    -+
C        |  cos(x)  sin(x)   0  |
C        |                      |
C        | -sin(x)  cos(x)   0  |.
C        |                      |
C        |  0        0       1  |
C        +-                    -+
C
C
C     The input matrix is assumed to be the product of three
C     rotation matrices, each one of the form
C
C        +-                    -+
C        |  1      0       0    |
C        |                      |
C        |  0    cos(r)  sin(r) |     (rotation of r radians about the
C        |                      |      x-axis),
C        |  0   -sin(r)  cos(r) |
C        +-                    -+
C
C
C        +-                    -+
C        | cos(s)   0   -sin(s) |
C        |                      |
C        |  0       1      0    |     (rotation of s radians about the
C        |                      |      y-axis),
C        | sin(s)   0    cos(s) |
C        +-                    -+
C
C     or
C
C        +-                    -+
C        |  cos(t)  sin(t)   0  |
C        |                      |
C        | -sin(t)  cos(t)   0  |     (rotation of t radians about the
C        |                      |      z-axis),
C        |  0        0       1  |
C        +-                    -+
C
C     where the second rotation axis is not equal to the first or
C     third.  Any rotation matrix can be factored as a sequence of
C     three such rotations, provided that this last criterion is met.
C
C     This routine is related to the SPICELIB routine EUL2M, which
C     produces a rotation matrix, given a sequence of Euler angles.
C     This routine is a `right inverse' of EUL2M:  the sequence of
C     calls
C
C        CALL M2EUL ( R,  AXIS3,   AXIS2,   AXIS1,
C       .                 ANGLE3,  ANGLE2,  ANGLE1     )
C
C        CALL EUL2M (     ANGLE3,  ANGLE2,  ANGLE1,
C       .                 AXIS3,   AXIS2,   AXIS1,   R )
C
C     preserves R, except for round-off error.
C
C
C     On the other hand, the sequence of calls
C
C        CALL EUL2M (     ANGLE3,  ANGLE2,  ANGLE1,
C       .                 AXIS3,   AXIS2,   AXIS1,   R )
C
C        CALL M2EUL ( R,  AXIS3,   AXIS2,   AXIS1,
C       .                 ANGLE3,  ANGLE2,  ANGLE1     )
C
C
C     preserve ANGLE3, ANGLE2, and ANGLE1 only if these angles start
C     out in the ranges that M2EUL's outputs are restricted to.
C
C$ Examples
C
C     1)  Conversion of instrument pointing from a matrix representation
C         to Euler angles:
C
C         Suppose we want to find camera pointing in alpha, delta, and
C         kappa, given the inertial-to-camera coordinate transformation
C
C
C    +-                                                               -+
C    |  0.49127379678135830  0.50872620321864170  0.70699908539882417  |
C    |                                                                 |
C    | -0.50872620321864193 -0.49127379678135802  0.70699908539882428  |
C    |                                                                 |
C    |  0.70699908539882406 -0.70699908539882439  0.01745240643728360  |
C    +-                                                               -+
C
C
C         We want to find angles alpha, delta, kappa such that
C
C            TICAM  =  [ kappa ]  [ pi/2 - delta ]  [ pi/2 + alpha ] .
C                               3                 1                 3
C
C         We can use the following small program to do this computation:
C
C
C            PROGRAM EX1
C            IMPLICIT NONE
C
C            DOUBLE PRECISION      DPR
C            DOUBLE PRECISION      HALFPI
C            DOUBLE PRECISION      TWOPI
C
C            DOUBLE PRECISION      ALPHA
C            DOUBLE PRECISION      ANG1
C            DOUBLE PRECISION      ANG2
C            DOUBLE PRECISION      DELTA
C            DOUBLE PRECISION      KAPPA
C            DOUBLE PRECISION      TICAM  ( 3, 3 )
C
C
C            DATA TICAM /  0.49127379678135830D0,
C           .             -0.50872620321864193D0,
C           .              0.70699908539882406D0,
C           .              0.50872620321864170D0,
C           .             -0.49127379678135802D0,
C           .             -0.70699908539882439D0,
C           .              0.70699908539882417D0,
C           .              0.70699908539882428D0,
C           .              0.01745240643728360D0  /
C
C
C            CALL M2EUL ( TICAM, 3, 1, 3, KAPPA, ANG2, ANG1 )
C
C            DELTA = HALFPI() - ANG2
C            ALPHA = ANG1     - HALFPI()
C
C            IF ( KAPPA .LT. 0.D0 ) THEN
C               KAPPA = KAPPA + TWOPI()
C            END IF
C
C            IF ( ALPHA .LT. 0.D0 ) THEN
C               ALPHA = ALPHA + TWOPI()
C            END IF
C
C            WRITE (*,'(1X,A,F24.14)') 'Alpha (deg): ', DPR() * ALPHA
C            WRITE (*,'(1X,A,F24.14)') 'Delta (deg): ', DPR() * DELTA
C            WRITE (*,'(1X,A,F24.14)') 'Kappa (deg): ', DPR() * KAPPA
C
C            END
C
C
C         The program's output should be something like
C
C            Alpha (deg):       315.00000000000000
C            Delta (deg):         1.00000000000000
C            Kappa (deg):        45.00000000000000
C
C         possibly formatted a little differently, or degraded slightly
C         by round-off.
C
C
C     2)  Conversion of instrument pointing angles from a non-J2000,
C         not necessarily inertial frame to J2000-relative RA, Dec,
C         and Twist.
C
C         Suppose that we have pointing for some instrument expressed as
C
C            [ gamma ]  [ beta ]  [ alpha ]
C                     3         2          3
C
C         with respect to some coordinate system S.  For example, S
C         could be a spacecraft-fixed system.
C
C         We will suppose that the transformation from J2000
C         coordinates to system S coordinates is given by the rotation
C         matrix J2S.
C
C         The rows of J2S are the unit basis vectors of system S, given
C         in J2000 coordinates.
C
C         We want to express the pointing with respect to the J2000
C         system as the sequence of rotations
C
C            [ kappa ]  [ pi/2 - delta ]  [ pi/2 + alpha ] .
C                     3                 1                 3
C
C         First, we use subroutine EUL2M to form the transformation
C         from system S to instrument coordinates S2INST.
C
C            CALL EUL2M ( GAMMA, BETA, ALPHA, 3, 2, 3, S2INST )
C
C         Next, we form the transformation from J2000 to instrument
C         coordinates J2INST.
C
C            CALL MXM ( S2INST, J2S, J2INST )
C
C         Finally, we express J2INST using the desired Euler angles, as
C         in the first example:
C
C            CALL M2EUL ( J2INST, 3, 1, 3, TWIST, ANG2, ANG3 )
C
C            RA   =  ANG3 - HALFPI()
C            DEC  =  HALFPI() - ANG2
C
C         If we wish to make sure that RA, DEC, and TWIST are in
C         the ranges [0, 2pi), [-pi/2, pi/2], and [0, 2pi)
C         respectively, we may add the code
C
C            IF ( RA .LT. 0.D0 ) THEN
C               RA = RA + TWOPI()
C            END IF
C
C            IF ( TWIST .LT. 0.D0 ) THEN
C               TWIST = TWIST + TWOPI()
C            END IF
C
C         Note that DEC is already in the correct range, since ANG2
C         is in the range [0, pi] when the first and third input axes
C         are equal.
C
C         Now RA, DEC, and TWIST express the instrument pointing
C         as RA, Dec, and Twist, relative to the J2000 system.
C
C         A warning note:  more than one definition of RA, Dec, and
C         Twist is extant.  Before using this example in an application,
C         check that the definition given here is consistent with that
C         used in your application.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.2.1, 21-DEC-2006 (NJB)
C
C        Error corrected in header example:  input matrix
C        previously did not match shown outputs.  Offending
C        example now includes complete program.
C
C-    SPICELIB Version 1.2.0, 15-OCT-2005 (NJB)
C
C        Updated to remove non-standard use of duplicate arguments
C        in MXM and MTXM calls.  A short error message cited in
C        the Exceptions section of the header failed to match
C        the actual short message used; this has been corrected.
C
C-    SPICELIB Version 1.1.2, 13-OCT-2004 (NJB)
C
C        Fixed header typo.
C
C-    SPICELIB Version 1.1.1, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.1.0, 02-NOV-1990 (NJB)
C
C        Header upgraded to describe notation in more detail.  Argument
C        names were changed to describe the use of the arguments more
C        accurately.  No change in functionality was made; the operation
C        of the routine is identical to that of the previous version.
C
C-    SPICELIB Version 1.0.0, 03-SEP-1990 (NJB)
C
C-&

C$ Index_Entries
C
C     matrix to euler angles
C
C-&

C$ Revisions
C
C-    SPICELIB Version 1.2.0, 26-AUG-2005 (NJB)
C
C        Updated to remove non-standard use of duplicate arguments
C        in MXM and MTXM calls.  A short error message cited in
C        the Exceptions section of the header  failed to match
C        the actual short message used; this has been corrected.
C
C-    SPICELIB Version 1.1.0, 02-NOV-1990 (NJB)
C
C        Argument names were changed to describe the use of the
C        arguments more accurately.  The axis and angle numbers
C        now decrease, rather than increase, from left to right.
C        The current names reflect the order of operator application
C        when the Euler angle rotations are applied to a vector:  the
C        rightmost matrix
C
C           [ ANGLE1 ]
C                     AXIS1
C
C        is applied to the vector first, followed by
C
C           [ ANGLE2 ]
C                     AXIS2
C
C        and then
C
C           [ ANGLE3 ]
C                     AXIS3
C
C        Previously, the names reflected the order in which the Euler
C        angle matrices appear on the page, from left to right.  This
C        naming convention was found to be a bit obtuse by a various
C        users.
C
C        No change in functionality was made; the operation of the
C        routine is identical to that of the previous version.
C
C        Also, the header was upgraded to describe the notation in more
C        detail.  The symbol
C
C           [ x ]
C                i
C
C        is explained at mind-numbing length.  An example was added
C        that shows a specific set of inputs and the resulting output
C        matrix.
C
C        The angle sequence notation was changed to be consistent with
C        Rotations required reading.
C
C          1-2-3  and  a-b-c
C
C        have been changed to
C
C          3-2-1  and  c-b-a.
C
C       Also, one `)' was changed to a `}'.
C
C       The phrase `first and third' was changed to `first or third'
C       in the $ Particulars section, where the criterion for the
C       existence of an Euler angle factorization is stated.
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               ISROT
      LOGICAL               RETURN

C
C     Local parameters
C

C
C     NTOL and DETOL are used to determine whether R is a rotation
C     matrix.
C
C     NTOL is the tolerance for the norms of the columns of R.
C
C     DTOL is the tolerance for the determinant of a matrix whose
C     columns are the unitized columns of R.
C
C
      DOUBLE PRECISION      NTOL
      PARAMETER           ( NTOL = 0.1D0 )

      DOUBLE PRECISION      DTOL
      PARAMETER           ( DTOL = 0.1D0 )

C
C     Local variables
C
      DOUBLE PRECISION      CHANGE ( 3, 3 )
      DOUBLE PRECISION      SIGN
      DOUBLE PRECISION      TMPMAT ( 3, 3 )
      DOUBLE PRECISION      TMPROT ( 3, 3 )

      INTEGER               C
      INTEGER               I
      INTEGER               NEXT   ( 3 )

      LOGICAL               DEGEN

C
C     Saved variables
C
      SAVE                  NEXT

C
C     Initial values
C
      DATA  NEXT / 2, 3, 1 /

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

C
C     The first order of business is to screen out the goofy cases.
C
C     Make sure the axis numbers are all right:  They must belong to
C     the set {1, 2, 3}...
C
      IF (      (  ( AXIS3 .LT. 1 ) .OR. ( AXIS3 .GT. 3 )  )
     .     .OR. (  ( AXIS2 .LT. 1 ) .OR. ( AXIS2 .GT. 3 )  )
     .     .OR. (  ( AXIS1 .LT. 1 ) .OR. ( AXIS1 .GT. 3 )  )   ) THEN

         CALL SETMSG ( 'Axis numbers are #,  #,  #. ' )
         CALL ERRINT ( '#',  AXIS3                    )
         CALL ERRINT ( '#',  AXIS2                    )
         CALL ERRINT ( '#',  AXIS1                    )
         CALL SIGERR ( 'SPICE(BADAXISNUMBERS)'        )
         CALL CHKOUT ( 'M2EUL'                        )
         RETURN

C
C     ...and the second axis number must differ from its neighbors.
C
      ELSE IF ( ( AXIS3 .EQ. AXIS2 ) .OR. ( AXIS1 .EQ. AXIS2 ) ) THEN

         CALL SETMSG ( 'Middle axis matches neighbor: # # #.' )
         CALL ERRINT ( '#',  AXIS3                            )
         CALL ERRINT ( '#',  AXIS2                            )
         CALL ERRINT ( '#',  AXIS1                            )
         CALL SIGERR ( 'SPICE(BADAXISNUMBERS)'                )
         CALL CHKOUT ( 'M2EUL'                                )
         RETURN

C
C     R must be a rotation matrix, or we may as well forget it.
C
      ELSE IF (  .NOT. ISROT ( R, NTOL, DTOL )  ) THEN

         CALL SETMSG ( 'Input matrix is not a rotation.' )
         CALL SIGERR ( 'SPICE(NOTAROTATION)'             )
         CALL CHKOUT ( 'M2EUL'                           )
         RETURN

      END IF

C
C     AXIS3, AXIS2, AXIS1 and R have passed their tests at this
C     point.  We take the liberty of working with TMPROT, a version
C     of R that has unitized columns.
C
      DO I = 1, 3
         CALL VHAT ( R(1,I), TMPROT(1,I) )
      END DO

C
C     We now proceed to recover the promised Euler angles from
C     TMPROT.
C
C        The ideas behind our method are explained in excruciating
C        detail in the ROTATION required reading, so we'll be terse.
C        Nonetheless, a word of explanation is in order.
C
C        The sequence of rotation axes used for the factorization
C        belongs to one of two categories:  a-b-a or c-b-a.  We
C        wish to handle each of these cases in one shot, rather than
C        using different formulas for each sequence to recover the
C        Euler angles.
C
C        What we're going to do is use the Euler angle formula for the
C        3-1-3 factorization for all of the a-b-a sequences, and the
C        formula for the 3-2-1 factorization for all of the c-b-a
C        sequences.
C
C        How can we get away with this?  The Euler angle formulas for
C        each factorization are different!
C
C        Our trick is to apply a change-of-basis transformation to the
C        input matrix R.  For the a-b-a factorizations, we choose a new
C        basis such that a rotation of ANGLE3 radians about the basis
C        vector indexed by AXIS3 becomes a rotation of ANGLE3 radians
C        about the third coordinate axis, and such that a rotation of
C        ANGLE2 radians about the basis vector indexed by AXIS2 becomes
C        a rotation of ANGLE2 radians about the first coordinate axis.
C        So R can be factored as a 3-1-3 rotation relative to the new
C        basis, and the Euler angles we obtain are the exact ones we
C        require.
C
C        The c-b-a factorizations can be handled in an analogous
C        fashion.  We transform R to a basis where the original axis
C        sequence becomes a 3-2-1 sequence.  In some cases, the angles
C        we obtain will be the negatives of the angles we require.  This
C        will happen if and only if the ordered basis (here the e's are
C        the standard basis vectors)
C
C            { e        e        e      }
C               AXIS3    AXIS2    AXIS1
C
C        is not right-handed.  An easy test for this condition is that
C        AXIS2 is not the successor of AXIS3, where the ordering is
C        cyclic.
C

      IF ( AXIS3 .EQ. AXIS1 ) THEN

C
C        The axis order is a-b-a.  We're going to find a matrix CHANGE
C        such that
C
C                 T
C           CHANGE  R  CHANGE
C
C        gives us R in the a useful basis, that is, a basis in which
C        our original a-b-a rotation is a 3-1-3 rotation, but where the
C        rotation angles are unchanged. To achieve this pleasant
C        simplification, we set column 3 of CHANGE to to e(AXIS3),
C        column 1 of CHANGE to e(AXIS2), and column 2 of CHANGE to
C
C          (+/-) e(C),
C
C        (C is the remaining index) depending on whether
C        AXIS3-AXIS2-C is a right-handed sequence of axes:  if it
C        is, the sign is positive.  (Here e(1), e(2), e(3) are the
C        standard basis vectors.)
C
C        Determine the sign of our third basis vector, so that we can
C        ensure that our new basis is right-handed.  The variable NEXT
C        is just a little mapping that takes 1 to 2, 2 to 3, and 3 to
C        1.
C
         IF ( AXIS2 .EQ. NEXT(AXIS3) ) THEN
            SIGN =  1.D0
         ELSE
            SIGN = -1.D0
         END IF

C
C        Since the axis indices sum to 6,
C
         C = 6 - AXIS3 - AXIS2

C
C        Set up the entries of CHANGE:
C
         CALL CLEARD ( 9, CHANGE )

         CHANGE ( AXIS3, 3 ) =         1.D0
         CHANGE ( AXIS2, 1 ) =         1.D0
         CHANGE ( C,     2 ) =  SIGN * 1.D0

C
C        Transform TMPROT.
C
         CALL MXM  ( TMPROT,  CHANGE,  TMPMAT )
         CALL MTXM ( CHANGE,  TMPMAT,  TMPROT )

C
C        Now we're ready to find the Euler angles, using a
C        3-1-3 factorization.  In general, the matrix product
C
C           [ a1 ]   [ a2 ]   [ a3 ]
C                 3        1        3
C
C        has the form
C
C     +-                                                              -+
C     |         cos(a1)cos(a3)          cos(a1)sin(a3)  sin(a1)sin(a2) |
C     | -sin(a1)cos(a2)sin(a3)  +sin(a1)cos(a2)cos(a3)                 |
C     |                                                                |
C     |        -sin(a1)cos(a3)         -sin(a1)sin(a3)  cos(a1)sin(a2) |
C     | -cos(a1)cos(a2)sin(a3)  +cos(a1)cos(a2)cos(a3)                 |
C     |                                                                |
C     |         sin(a2)sin(a3)         -sin(a2)cos(a3)         cos(a2) |
C     +-                                                              -+
C
C
C        but if a2 is 0 or pi, the product matrix reduces to
C
C
C     +-                                                              -+
C     |         cos(a1)cos(a3)          cos(a1)sin(a3)               0 |
C     | -sin(a1)cos(a2)sin(a3)  +sin(a1)cos(a2)cos(a3)                 |
C     |                                                                |
C     |        -sin(a1)cos(a3)         -sin(a1)sin(a3)               0 |
C     | -cos(a1)cos(a2)sin(a3)  +cos(a1)cos(a2)cos(a3)                 |
C     |                                                                |
C     |                      0                       0         cos(a2) |
C     +-                                                              -+
C
C
C        In this case, a1 and a3 are not uniquely determined.  If we
C        arbitrarily set a1 to zero, we arrive at the matrix
C
C           +-                                         -+
C           |         cos(a3)         sin(a3)      0    |
C           | -cos(a2)sin(a3)  cos(a2)cos(a3)      0    |
C           |               0            0      cos(a2) |
C           +-                                         -+
C
C        We take care of this case first.  We test three conditions
C        that are mathematically equivalent, but may not be satisfied
C        simultaneously because of round-off:
C
C
         DEGEN =      (        (       TMPROT(1,3)    .EQ. 0.D0  )
     .                   .AND. (       TMPROT(2,3)    .EQ. 0.D0  )  )
     .           .OR. (        (       TMPROT(3,1)    .EQ. 0.D0  )
     .                   .AND. (       TMPROT(3,2)    .EQ. 0.D0  )  )
     .           .OR. (           ABS( TMPROT(3,3) )  .EQ. 1.D0     )

C
C        In the following block of code, we make use of the fact that
C
C           SIN ( ANGLE2 )   >  0
C                            -
C        in choosing the signs of the ATAN2 arguments correctly.  Note
C        that ATAN2(x,y) = -ATAN2(-x,-y).
C
C
         IF ( DEGEN ) THEN

            ANGLE3 = 0.D0
            ANGLE2 = ACOS  (  TMPROT(3,3)                 )
            ANGLE1 = ATAN2 (  TMPROT(1,2),   TMPROT(1,1)  )

         ELSE
C
C           The normal case.
C
            ANGLE3 = ATAN2 (  TMPROT(1,3),   TMPROT(2,3)  )
            ANGLE2 = ACOS  (  TMPROT(3,3)                 )
            ANGLE1 = ATAN2 (  TMPROT(3,1),  -TMPROT(3,2)  )

         END IF

      ELSE

C
C        The axis order is c-b-a.  We're going to find a matrix CHANGE
C        such that
C
C                 T
C           CHANGE  R  CHANGE
C
C        gives us R in the a useful basis, that is, a basis in which
C        our original c-b-a rotation is a 3-2-1 rotation, but where the
C        rotation angles are unchanged, or at worst negated.  To
C        achieve this pleasant simplification, we set column 1 of
C        CHANGE to to e(AXIS3), column 2 of CHANGE to e(AXIS2), and
C        column 3 of CHANGE to
C
C          (+/-) e(AXIS1),
C
C        depending on whether AXIS3-AXIS2-AXIS1 is a right-handed
C        sequence of axes:  if it is, the sign is positive.  (Here
C        e(1), e(2), e(3) are the standard basis vectors.)
C
C        We must be cautious here, because if AXIS3-AXIS2-AXIS1 is a
C        right-handed sequence of axes, all of the rotation angles will
C        be the same in our new basis, but if it's a left-handed
C        sequence, the third angle will be negated.  Let's get this
C        straightened out right now.  The variable NEXT is just a
C        little mapping that takes 1 to 2, 2 to 3, and 3 to 1.
C
         IF ( AXIS2 .EQ. NEXT(AXIS3) ) THEN
            SIGN =  1.D0
         ELSE
            SIGN = -1.D0
         END IF

C
C        Set up the entries of CHANGE:
C
         CALL CLEARD ( 9, CHANGE )

         CHANGE ( AXIS3, 1 ) =         1.D0
         CHANGE ( AXIS2, 2 ) =         1.D0
         CHANGE ( AXIS1, 3 ) =  SIGN * 1.D0

C
C        Transform TMPROT.
C
         CALL MXM  ( TMPROT,  CHANGE,  TMPMAT )
         CALL MTXM ( CHANGE,  TMPMAT,  TMPROT )

C
C        Now we're ready to find the Euler angles, using a
C        3-2-1 factorization.  In general, the matrix product
C
C           [ a1 ]   [ a2 ]   [ a3 ]
C                 1        2        3
C
C        has the form
C
C
C     +-                                                              -+
C     |         cos(a2)cos(a3)          cos(a2)sin(a3)        -sin(a2) |
C     |                                                                |
C     |        -cos(a1)sin(a3)          cos(a1)cos(a3)  sin(a1)cos(a2) |
C     | +sin(a1)sin(a2)cos(a3)  +sin(a1)sin(a2)sin(a3)                 |
C     |                                                                |
C     |         sin(a1)sin(a3)         -sin(a1)cos(a3)  cos(a1)cos(a2) |
C     | +cos(a1)sin(a2)cos(a3)  +cos(a1)sin(a2)sin(a3)                 |
C     +-                                                              -+
C
C
C        but if a2 is -pi/2 or pi/2, the product matrix reduces to
C
C
C     +-                                                              -+
C     |                      0                       0        -sin(a2) |
C     |                                                                |
C     |        -cos(a1)sin(a3)          cos(a1)cos(a3)               0 |
C     | +sin(a1)sin(a2)cos(a3)  +sin(a1)sin(a2)sin(a3)                 |
C     |                                                                |
C     |         sin(a1)sin(a3)         -sin(a1)cos(a3)               0 |
C     | +cos(a1)sin(a2)cos(a3)  +cos(a1)sin(a2)sin(a3)                 |
C     +-                                                              -+
C
C
C        In this case, a1 and a3 are not uniquely determined.  If we
C        arbitrarily set a1 to zero, we arrive at the matrix
C
C           +-                                              -+
C           |               0                 0    -sin(a2)  |
C           |        -sin(a3)           cos(a3)          0   |,
C           |  sin(a2)cos(a3)    sin(a2)sin(a3)          0   |
C           +-                                              -+
C
C
C        We take care of this case first.  We test three conditions
C        that are mathematically equivalent, but may not be satisfied
C        simultaneously because of round-off:
C
C
         DEGEN =      (        (       TMPROT(1,1)    .EQ. 0.D0  )
     .                   .AND. (       TMPROT(1,2)    .EQ. 0.D0  )  )
     .           .OR. (        (       TMPROT(2,3)    .EQ. 0.D0  )
     .                   .AND. (       TMPROT(3,3)    .EQ. 0.D0  )  )
     .           .OR. (           ABS( TMPROT(1,3) )  .EQ. 1.D0     )

C
C        In the following block of code, we make use of the fact that
C
C           COS ( ANGLE2 )   >  0
C                            -
C        in choosing the signs of the ATAN2 arguments correctly.  Note
C        that ATAN2(x,y) = -ATAN2(-x,-y).
C
C
         IF ( DEGEN ) THEN

            ANGLE3 =          0.D0
            ANGLE2 =          ASIN  ( -TMPROT(1,3)                 )
            ANGLE1 = SIGN  *  ATAN2 ( -TMPROT(2,1),   TMPROT(2,2)  )

         ELSE
C
C           The normal case.
C
            ANGLE3 =          ATAN2 (  TMPROT(2,3),   TMPROT(3,3)  )
            ANGLE2 =          ASIN  ( -TMPROT(1,3)                 )
            ANGLE1 = SIGN  *  ATAN2 (  TMPROT(1,2),   TMPROT(1,1)  )

         END IF

      END IF

      CALL CHKOUT ( 'M2EUL' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
*DECK S1MERG
      SUBROUTINE S1MERG (TCOS, I1, M1, I2, M2, I3)
C***BEGIN PROLOGUE  S1MERG
C***SUBSIDIARY
C***PURPOSE  Merge two strings of ascending real numbers.
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (S1MERG-S, D1MERG-D, C1MERG-C, I1MERG-I)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   This subroutine merges two ascending strings of numbers in the
C   array TCOS.  The first string is of length M1 and starts at
C   TCOS(I1+1).  The second string is of length M2 and starts at
C   TCOS(I2+1).  The merged string goes into TCOS(I3+1).
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  SCOPY
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   901120  Modified to use IF-THEN-ELSE.  Previous spaghetti code did
C           not compile correctly with optimization on the IBM RS6000.
C           (RWC)
C   920130  Code name changed from MERGE to S1MERG.  (WRB)
C***END PROLOGUE  S1MERG
      INTEGER I1, I2, I3, M1, M2
      REAL TCOS(*)
C
      INTEGER J1, J2, J3
C
C***FIRST EXECUTABLE STATEMENT  S1MERG
      IF (M1.EQ.0 .AND. M2.EQ.0) RETURN
C
      IF (M1.EQ.0 .AND. M2.NE.0) THEN
         CALL SCOPY (M2, TCOS(I2+1), 1, TCOS(I3+1), 1)
         RETURN
      ENDIF
C
      IF (M1.NE.0 .AND. M2.EQ.0) THEN
         CALL SCOPY (M1, TCOS(I1+1), 1, TCOS(I3+1), 1)
         RETURN
      ENDIF
C
      J1 = 1
      J2 = 1
      J3 = 1
C
   10 IF (TCOS(I1+J1) .LE. TCOS(I2+J2)) THEN
         TCOS(I3+J3) = TCOS(I1+J1)
         J1 = J1+1
         IF (J1 .GT. M1) THEN
            CALL SCOPY (M2-J2+1, TCOS(I2+J2), 1, TCOS(I3+J3+1), 1)
            RETURN
         ENDIF
      ELSE
         TCOS(I3+J3) = TCOS(I2+J2)
         J2 = J2+1
         IF (J2 .GT. M2) THEN
            CALL SCOPY (M1-J1+1, TCOS(I1+J1), 1, TCOS(I3+J3+1), 1)
            RETURN
         ENDIF
      ENDIF
      J3 = J3+1
      GO TO 10
      END

ext: f


text:
      subroutine ivimp(i1,i2,pas,iv)
c     generate iv=i1:pas:i2
c     Copyright INRIA
      integer i1,i2,pas,iv(*)
      k=0
      do 10 i=i1,i2,pas
         k=k+1
         iv(k)=i
 10   continue
      end

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
	SUBROUTINE AF_DUPE  ( lunf, rpid, dattim, xlat, xlon, xelv,
     +				iret )
C************************************************************************
C* AF_DUPE								*
C*									*
C* This routine will check for a duplicate station entry in the ship	*
C* format data file. A duplicate station entry is a report at the same	*
C* lat/lon location and flight level (station elevation). 		*
C*									*
C* AF_DUPE  ( LUNF, RPID, DATTIM, XLAT, XLON, XELV, IRET )		*
C*									*
C* Input parameters:							*
C*	LUNF		INTEGER		Surface file number		*
C*	RPID		CHAR*8		Station identifier		*
C*	DATTIM		CHAR*		GEMPAK time			*
C*	XLAT		REAL		Station latitude		*
C*	XLON		REAL		Station longitude		*
C*	XELV		REAL		Station elevation		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					   0 = normal return		*
C*					   1 = duplicate found		*
C**									*
C* Log:									*
C* S. Chiswell/Unidata	 5/01	Initial development			*
C* S. Danz/AWC	 	 4/04	Changed epsilon to 0.011 to catch round *
C*				off in values coming back from the db   *
C************************************************************************
	CHARACTER*(*)	rpid, dattim
	REAL		xlat, xlon, xelv
	INTEGER		lunf, iret
C
	INTEGER		isdupe
	CHARACTER	stid*9, stn*8, area*12
C------------------------------------------------------------------------
	iret = 0
	isdupe = 0
C
C*	See if this time already exists in the output file
C
	CALL SF_STIM ( lunf, dattim, ier )
	IF ( ier .ne. 0 ) RETURN
C
	stid = '@' // rpid
	area = ' '
	CALL SF_UARE ( lunf, stid, .true., area, stn, ier )

	DO WHILE ( ( ier .eq. 0 ) .and. ( isdupe .eq. 0 ) )

	   CALL SF_SNXT ( lunf, stn, istnm, slat, slon,
     +			selv, ispri, ier )
	   IF ( ier .eq. 0) THEN
	      dlat = abs(slat - xlat)
	      dlon = abs(slon - xlon)
	      delv = abs(selv - xelv)
	      IF(( dlat .lt. 0.011 ) .and. ( dlon .lt. 0.011 ) .and.
     +		 ( delv .lt. 1 ) ) THEN
		 isdupe = 1
	      END IF
	   END IF
	END DO

	IF ( isdupe .eq. 1 ) THEN
	   iret = isdupe
	END IF

	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
module bc_state_HP_wall
    use mod_kinds,          only: rk,ik
    use mod_constants,      only: ZERO, ONE
    use type_bc_state,      only: bc_state_t
    use type_chidg_worker,  only: chidg_worker_t
    use type_properties,    only: properties_t
    use mpi_f08,            only: mpi_comm
    use DNAD_D
    implicit none

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   11/09/2018
    !!
    !-----------------------------------------------------------------------------------------
    type, public, extends(bc_state_t) :: HP_wall_t

    contains

        procedure   :: init
        procedure   :: compute_bc_state    !> bc implementation

    end type HP_wall_t
    !******************************************************************************************

contains

    !------------------------------------------------------------------------------------------
    subroutine init(self)
        class(HP_wall_t),  intent(inout)   :: self

        !
        ! Set name
        !
        call self%set_name('HP Wall')
        call self%set_family('Scalar')

    end subroutine init
    !******************************************************************************************

    !> Specialized compute routine for Extrapolation Boundary Condition
    !!
    !!  @author Nathan A. Wukie
    !!
    !!  @param[in]      eqnset  Equation Set type governing the current domain
    !!  @param[in]      mesh    Mesh data containing elements and faces for the domain
    !!  @param[inout]   sdata   Solver data containing solution vector, rhs, linearization, etc.
    !!  @param[in]      ielem   Index of the element being computed
    !!  @param[in]      iface   Index of the face being computed
    !!  @param[in]      iblk    Index of the linearization block being computed
    !---------------------------------------------------------------------------------------------
    subroutine compute_bc_state(self,worker,prop,bc_COMM)
        class(HP_wall_t),   intent(inout)   :: self
        type(chidg_worker_t),                   intent(inout)   :: worker
        class(properties_t),                    intent(inout)   :: prop
        type(mpi_comm),                         intent(in)      :: bc_COMM

        type(AD_D), allocatable, dimension(:)   ::  &
            u_bc, p_bc, q_bc, r_bc,                 &
            grad1_u_bc, grad2_u_bc, grad3_u_bc,     &
            grad1_p_bc, grad2_p_bc, grad3_p_bc,     &
            grad1_q_bc, grad2_q_bc, grad3_q_bc,     &
            grad1_r_bc, grad2_r_bc, grad3_r_bc

        ! Get u_m from face interior to initialize derivatives
        u_bc       = worker%get_field('u','value', 'face interior')
        grad1_u_bc = worker%get_field('u','grad1', 'face interior')
        grad2_u_bc = worker%get_field('u','grad2', 'face interior')
        grad3_u_bc = worker%get_field('u','grad3', 'face interior')

        p_bc       = worker%get_field('p','value', 'face interior')
        grad1_p_bc = worker%get_field('p','grad1', 'face interior')
        grad2_p_bc = worker%get_field('p','grad2', 'face interior')
        grad3_p_bc = worker%get_field('p','grad3', 'face interior')

        q_bc       = worker%get_field('q','value', 'face interior')
        grad1_q_bc = worker%get_field('q','grad1', 'face interior')
        grad2_q_bc = worker%get_field('q','grad2', 'face interior')
        grad3_q_bc = worker%get_field('q','grad3', 'face interior')

        r_bc       = worker%get_field('r','value', 'face interior')
        grad1_r_bc = worker%get_field('r','grad1', 'face interior')
        grad2_r_bc = worker%get_field('r','grad2', 'face interior')
        grad3_r_bc = worker%get_field('r','grad3', 'face interior')

        ! ZERO dirichlet condition on walls
        u_bc = ZERO

        ! Store boundary condition state, Value
        call worker%store_bc_state('u', u_bc, 'value')
        call worker%store_bc_state('p', p_bc, 'value')
        call worker%store_bc_state('q', q_bc, 'value')
        call worker%store_bc_state('r', r_bc, 'value')

        ! Store boundary condition state, gradient
        call worker%store_bc_state('u', grad1_u_bc, 'grad1')
        call worker%store_bc_state('u', grad2_u_bc, 'grad2')
        call worker%store_bc_state('u', grad3_u_bc, 'grad3')

        call worker%store_bc_state('p', grad1_p_bc, 'grad1')
        call worker%store_bc_state('p', grad2_p_bc, 'grad2')
        call worker%store_bc_state('p', grad3_p_bc, 'grad3')

        call worker%store_bc_state('q', grad1_q_bc, 'grad1')
        call worker%store_bc_state('q', grad2_q_bc, 'grad2')
        call worker%store_bc_state('q', grad3_q_bc, 'grad3')

        call worker%store_bc_state('r', grad1_r_bc, 'grad1')
        call worker%store_bc_state('r', grad2_r_bc, 'grad2')
        call worker%store_bc_state('r', grad3_r_bc, 'grad3')

    end subroutine compute_bc_state
    !*********************************************************************************************

end module bc_state_HP_wall

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
!
! Copyright (c) 2014, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
! This tests replacement inside directives
!
#define FillInDirective FLUSH

program p
    logical :: res(1) = .false., expect(1) = .true.
    !$omp FillInDirective
    print *, "Testing a single quote ' yep! That's a quote alright"
    res(1) = .true.
    call check(res, expect, 1)
end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


! chk_allocstring --
!     Check: does the compiler support character strings with allocatable length?
!
program chk_allocstring
    character(len=:), allocatable :: string

    allocate( character(len=10) :: string )

    write( *, '(a,i0)' ) 'Allocated length: ', len(string)

    string = '1'

    write( *, '(a,i0)' ) 'Length after assignment: ', len(string)
    if ( len(string) /= 10 ) then
        write( *, '(a)' ) '    Be aware that the length may be influenced by assignment!'
        write( *, '(a)' ) '    (automatic reallocation is taking place)'
    endif
end program chk_allocstring


text:
	SUBROUTINE IM_CALIB( iarray, ilen, ioff, iret)
C************************************************************************
C* IM_CALIB                                                             *
C*                                                                      *
C* This subroutine reads the calibration block of an AREA file, and	*
C* stores the colorbar levels in the appropriate common block variables.*
C* If the calibration block is a PROD, then the pixel and data values	*
C* are retrieved for the colorbar levels.				*
C*                                                                      *
C* IM_CALIB  ( IARRAY, ILEN, IOFF, IRET )				*
C*                                                                      *
C* Input parameters:                                                    *
C*      IARRAY          INTEGER*        Header array			*
C*	ILEN		INTEGER		Header array length		*
C*	IOFF		INTEGER		Byte offset for calibration	*
C*									*
C* Output parameters:                                                   *
C*      IRET            INTEGER         Return code                     *
C*                                        0 = normal return             *
C**                                                                     *
C* Log:                                                                 *
C* Chiz/Unidata		07/00		Created for CIMSS products	*
C*					Initially, only PROD is handled.*
C* S. Chiswell/Unidata	 2/04		Added imcmn.cmn for imcalbar	*
C************************************************************************

	INTEGER		iarray(*), ilen, ioff, iret

	CHARACTER	tval*4

	INCLUDE		'IMGDEF.CMN'
	INCLUDE		'imcmn.cmn'

	iret = 0

	istart = (ioff / 4) + 1
C
C*	See if this is a known Calibration type
C
	CALL ST_ITOC ( iarray (istart), 1, tval, ier )

C	IF ( ier .eq. 0 ) THEN
        if((tval .eq. 'PROD').or.(tval .eq. 'N0R').or.
     +		(tval .eq. 'FIRE' ) ) then

c          cmcalb = tval
           cmcalb = 'PROD'
C
C*	   Read in min and max pixel and corresponding data
C*	   values. Scale value is for data points.
C
	   iminval = iarray( istart + 1 )
	   imaxval = iarray( istart + 2 )
	   iminpix = iarray( istart + 3 )
	   imaxpix = iarray( istart + 4 )

C
C*	   Un-Swap calibration text values if necessary
C
	   iloc = istart + 5
	   if (imbswp .eq. 1)
     +		ier = MV_SWP4 ( 1, iarray (iloc), iarray (iloc) )
	   CALL ST_ITOC ( iarray (iloc), 1, tval, ier )

           if(ier .eq. 0) then
              cmbunt = tval
           else
              cmbunt = 'Unk'
           endif

	   iscaleval = iarray( istart + 6 )
	   if(iscaleval .eq. 0) iscaleval = 1

C
C*	   Determine colorbar levels and store commonblock values.
C
	   ratio = float(imaxval - iminval)/float(imaxpix - iminpix)

	   immnpx = iminpix
	   immxpx = imaxpix
           iminpix = iminpix + 1
           imaxpix = imaxpix + 1
           do i=1,256
              if(i.eq.iminpix) then
                 CALL ST_INCH ( iminval/iscaleval, cmblev (i), ier )
              else if(i.eq.imaxpix) then
                 CALL ST_INCH ( imaxval/iscaleval, cmblev (i), ier )
              else if((i.gt.iminpix).and.(i.lt.imaxpix)) then
                    IF ( mod ( i-1,16) .eq. 0) THEN
                        level = nint( (i-iminpix) * ratio) + iminval
                        CALL ST_INCH ( level/iscaleval,
     +                                 cmblev (i), ier )
                    ELSE
                        cmblev (i) = ' '
                    END IF
              else
                 cmblev(i) = ' '
              endif
           end do
        endif

	imcalbar = 1

	RETURN
	END


text:
! { dg-do run }
! { dg-options "-finit-local-zero -fbackslash" }

program init_flag_1
  call real_test
  call logical_test
  call int_test
  call complex_test
  call char_test
end program init_flag_1

! Test some initializations for both implicitly and
! explicitly declared local variables.
subroutine real_test
  real r1
  real r2(10)
  dimension r3(10,10)
  if (r1 /= 0.0) call abort
  if (r2(2) /= 0.0) call abort
  if (r3(5,5) /= 0.0) call abort
  if (r4 /= 0.0) call abort
end subroutine real_test

subroutine logical_test
  logical l1
  logical l2(2)
  if (l1 .neqv. .false.) call abort
  if (l2(2) .neqv. .false.) call abort
end subroutine logical_test

subroutine int_test
  integer i1
  integer i2(10)
  dimension i3(10,10)
  if (i1 /= 0) call abort
  if (i2(2) /= 0) call abort
  if (i3(5,5) /= 0) call abort
  if (i4 /= 0) call abort
end subroutine int_test

subroutine complex_test
  complex c1
  complex c2(20,20)
  if (c1 /= (0.0,0.0)) call abort
  if (c2(1,1) /= (0.0,0.0)) call abort
end subroutine complex_test

subroutine char_test
  character*1 c1
  character*8 c2, c3(5)
  character c4(10)
  if (c1 /= '\0') call abort
  if (c2 /= '\0\0\0\0\0\0\0\0') call abort
  if (c3(1) /= '\0\0\0\0\0\0\0\0') call abort
  if (c3(5) /= '\0\0\0\0\0\0\0\0') call abort
  if (c4(5) /= '\0') call abort
end subroutine char_test

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Beverages and More!, better known as BevMo, is your source for any and all things alcohol.  Large beer selection, larger wine selection, and a good selection of other alcoholic beverages.  They also have brewing kits if you want to make your own.  Throw in a bunch of cheeses, chips, and other snacks and accessories, and you have what the name of the store describes.  If you like cigars, they have a few of those as well.

There is another location in Arden, but the Natomas location (linked above) is closest to Davis, just off the Truxel exit of I80.

20080716 00:52:17 nbsp Disney Land for Alcoholics Users/DavisResident23

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE DELFIL (FILE_NAME,NOT_EX_ERR)
      IMPLICIT NONE

*     FORMAL_PARAMETERS:
      CHARACTER*(*) FILE_NAME
      LOGICAL NOT_EX_ERR

**    Local variables
      LOGICAL THERE
      INTEGER UN,GETUN,IL,IOS
      CHARACTER*132 LFNAME
      SAVE

      IL = LEN_TRIM (FILE_NAME)
      IF (IL.EQ.0) CALL FATALERR ('DELFIL','empty file name')

      LFNAME = FILE_NAME
      CALL FLNAME (LFNAME)
      INQUIRE (FILE=LFNAME(1:IL),EXIST=THERE)

      IF (THERE) THEN

         UN = GETUN (10,99)
         OPEN (UN,FILE=LFNAME(1:IL),STATUS='OLD',
     &         ACCESS='SEQUENTIAL',FORM='UNFORMATTED',IOSTAT=IOS)
         IF (IOS.NE.0) CALL FATALERR ('DELFIL','cannot open file')

         CLOSE (UN,STATUS='DELETE',IOSTAT=IOS)
         IF (IOS.NE.0) CALL FATALERR ('DELFIL','cannot open file')

      ELSE IF (.NOT.THERE.AND.NOT_EX_ERR) THEN

         CALL FATALERR ('DELFIL','file to be deleted does not exist')

      END IF

      RETURN
      END

ext: for

max_issues_repo_name: zhumiansanren/FIRI-WHS


text:
      program equiv01

C     Check that unusable information is not used

      equivalence (i,j), (k,x)

      read*, i, j, k, x

      if(i.gt.n) then
         i = i - 1
         print *, i
      endif

      if(j.gt.n) then
         j = j - 1
         print *, j
      endif

      if(k.gt.n) then
         k = k - 1
         print *, k
      endif

      print *, i, j, k, x

      end


text:
MODULE mo_ts_procedures
  USE grid, ONLY : nzp,nxp,nyp,level
  USE mo_aux_state, ONLY : dzt,zt
  USE mo_stats_finder
  IMPLICIT NONE

  PRIVATE

  PUBLIC :: globalMax, globalMean3d, globalAvgWaterPaths,    &
            globalAvgCloudBoundaries, globalAvgCloudFraction

  CONTAINS

    ! Check if variable exists in the FieldArray instances, then check if the data is stored
    ! or if it is associated with an onDemand procedure. Finally calculate the statistics.

    !
    ! Global max values
    !
    SUBROUTINE globalMax(name,output)
      USE util, ONLY : get_gmax
      CHARACTER(len=*), INTENT(in) :: name
      REAL, INTENT(out) :: output
      REAL :: fvar(nzp,nxp,nyp)
      fvar = 0.
      output = 0.
      CALL stats_get3d(name,fvar)
      CALL get_gmax(nzp,nxp,nyp,fvar,output)
    END SUBROUTINE globalMax

    !
    ! Global mean from 3d
    !
    SUBROUTINE globalMean3d(name,output)
      USE util, ONLY : get_avg4_root
      CHARACTER(len=*), INTENT(in) :: name
      REAL, INTENT(out) :: output
      REAL :: fvar(nzp,nxp,nyp)
      fvar = 0.
      output = 0.
      CALL stats_get3d(name,fvar)
      CALL get_avg4_root(nzp,nxp,nyp,fvar,output,dzt%d)
    END SUBROUTINE globalMean3d

    !
    ! Global average water paths
    !
    SUBROUTINE globalAvgWaterPaths(name,output)
      USE util, ONLY : get_avg3_root
      USE mo_diag_state, ONLY : a_rc, a_srp, a_ri, a_riri, a_dn
      USE mo_progn_state, ONLY : a_rpp
      CHARACTER(len=*), INTENT(in) :: name
      REAL, INTENT(out) :: output
      REAL :: field(nzp,nxp,nyp)
      REAL :: prof(nzp)
      prof = 0.
      field = 0.
      output = 0.
      SELECT CASE(name)
      CASE('lwp_bar')
         field = a_rc%d*a_dn%d
      CASE('rwp_bar')
         IF (level < 4) THEN
            field = a_rpp%d*a_dn%d
         ELSE
            field = a_srp%d*a_dn%d
         END IF
      CASE('iwp_bar')
         field = (a_riri%d+a_ri%d)*a_dn%d
      END SELECT
      CALL get_avg3_root(nzp,nxp,nyp,field,prof)
      output = SUM(prof/dzt%d)
    END SUBROUTINE globalAvgWaterPaths

    !
    ! Global average cloud boundaries
    ! Take the global mean profile of cloud water and diagnose from there
    !
    SUBROUTINE globalAvgCloudBoundaries(name,output)
      USE util, ONLY : get_avg3_root
      USE mo_diag_state, ONLY : a_rc,a_ri,a_riri
      CHARACTER(len=*), INTENT(in) :: name
      REAL, INTENT(out) :: output
      REAL :: field(nzp,nxp,nyp)
      REAL :: prof(nzp)
      LOGICAL :: cldmask(nzp)
      REAL, PARAMETER :: TH = 1.e-6
      INTEGER :: k
      field = 0.
      prof = 0.
      output = 0.
      IF (level < 5) THEN
         field = a_rc%d
      ELSE
         field = a_rc%d + a_ri%d + a_riri%d
      END IF
      CALL get_avg3_root(nzp,nxp,nyp,field,prof)
      cldmask = ( prof > TH  )
      IF (ANY(cldmask)) THEN ! any clouds found?
         SELECT CASE(name)
         CASE('ctop')  ! Catch the first cloud top starting from model top
            DO k = nzp,1,-1
               IF (cldmask(k)) THEN
                  output = zt%d(k)
                  EXIT
               END IF
            END DO
         CASE('cbase') ! Catch the first cloud base starting from surface
            DO k = 2,nzp
               IF (cldmask(k)) THEN
                  output = zt%d(k)
                  EXIT
               END IF
            END DO
         END SELECT
      END IF
    END SUBROUTINE globalAvgCloudBoundaries

    !
    ! Global average cloud fraction
    ! Mask individual columns where there is any cloud as 1, the rest as 0 and take the global mean
    !
    SUBROUTINE globalAvgCloudFraction(name,output)
      USE util, ONLY : get_avg2_root
      USE mo_diag_state, ONLY : a_rc,a_ri,a_riri
      CHARACTER(len=*), INTENT(in) :: name
      REAL, INTENT(out) :: output
      REAL :: field(nzp,nxp,nyp)
      REAL :: lateral(nxp,nyp)
      REAL, PARAMETER :: TH = 1.e-6
      field = 0.
      output = 0.
      lateral = 0.
      field = a_rc%d
      IF (level == 5) field = field + a_riri%d + a_riri%d
      lateral(:,:) = MERGE(1.,0., ANY( field(:,:,:) > TH, DIM=1 ))
      CALL get_avg2_root(nxp,nyp,lateral,output)
    END SUBROUTINE globalAvgCloudFraction

END MODULE mo_ts_procedures

max_stars_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_issues_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_forks_repo_name: UCLALES-SALSA/UCLALES-SALSA


text:
***********************************************************************
c                    ORBEL_FGET.F
***********************************************************************
*     PURPOSE:  Solves Kepler's eqn. for hyperbola using hybrid approach.
*
*             Input:
*                           e ==> eccentricity anomaly. (real scalar)
*                        capn ==> hyperbola mean anomaly. (real scalar)
*             Returns:
*                  orbel_fget ==>  eccentric anomaly. (real scalar)
*
*     ALGORITHM: Based on pp. 70-72 of Fitzpatrick's book "Principles of
*           Cel. Mech. ".  Quartic convergence from Danby's book.
*     REMARKS:
*     AUTHOR: M. Duncan
*     DATE WRITTEN: May 11, 1992.
*     REVISIONS: 2/26/93 hfl
***********************************************************************

	real*8 function orbel_fget(e,capn)

      include '../swift.inc'

c...  Inputs Only:
	real*8 e,capn

c...  Internals:
	integer i,IMAX
	real*8 tmp,x,shx,chx
	real*8 esh,ech,f,fp,fpp,fppp,dx
	PARAMETER (IMAX = 10)

c----
c...  Executable code

c Function to solve "Kepler's eqn" for F (here called
c x) for given e and CAPN.

c  begin with a guess proposed by Danby
	if( capn .lt. 0.d0) then
	   tmp = -2.d0*capn/e + 1.8d0
	   x = -log(tmp)
	else
	   tmp = +2.d0*capn/e + 1.8d0
	   x = log( tmp)
	endif

	orbel_fget = x

	do i = 1,IMAX
	  call orbel_schget(x,shx,chx)
	  esh = e*shx
	  ech = e*chx
	  f = esh - x - capn
c	  write(6,*) 'i,x,f : ',i,x,f
	  fp = ech - 1.d0
	  fpp = esh
	  fppp = ech
	  dx = -f/fp
	  dx = -f/(fp + dx*fpp/2.d0)
	  dx = -f/(fp + dx*fpp/2.d0 + dx*dx*fppp/6.d0)
	  orbel_fget = x + dx
c   If we have converged here there's no point in going on
	  if(abs(dx) .le. TINY) RETURN
	  x = orbel_fget
	enddo

	write(6,*) 'FGET : RETURNING WITHOUT COMPLETE CONVERGENCE'
	return
	end   ! orbel_fget
c------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
    subroutine carbo(Patm,temperature,runoff,xmass_carbonate,j)
!   ======================================================
    implicit none

    include 'combine_foam.inc'

    Pmax=1.+0.302483*runoff**0.8
    !Patm=var(12,nbasin)/0.0508d+18 ! obsolete, Patm is now an input argument
    pco2soil=Patm+Pmax/(1+dexp(1.315-0.119*(temperature)))
    SO4soil=SO4_conc(j)

    temper=temperature+273.15

    x1=1.d-7
    x2=1.d-1
    xacc=1.d-10
!    if (value_keep1(j).ne.0) then
!        xn=value_keep1(j)
!    else
    xn=1.d-4
!    endif

!    call rtbis(x1,x2,xacc,out,temper)
    call newton(xn,fx,dfx,xacc,temper)

    out=xn
!    value_keep1(j)=xn

    h2co3soil=betasoil*pco2soil
    hco3soil=ak1soil*h2co3soil/out
    co3soil=hco3soil*ak2soil/out
    Calcium=0.5*(hco3soil+2*co3soil+2.*SO4soil)  !mol/m3
    xmass_carbonate=Calcium !mol/m3

    return
    end


text:
!> @file
!! Define the maxdiff operation via mpi wrappers
!! @author
!!    Copyright (C) 2012-2018 BigDFT group
!!    This file is distributed under the terms of the
!!    GNU General Public License, see ~/COPYING file
!!    or http://www.gnu.org/copyleft/gpl.txt .
!!    For the list of contributors, see ~/AUTHORS
!!$module f_maxdiff
!!$  use f_precisions
!!$  use f_allreduce
!!$  implicit none
!!$
!!$  private
!!$
!!$
!!$
!!$contains

  !> Detect the maximum difference between arrays all over a given communicator
  function mpimaxdiff_i0(n,array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    integer, intent(in) :: n !<number of elements to be controlled
    integer(f_integer) :: array !< starting point of the array
    integer(f_integer), dimension(:,:), allocatable :: array_glob
    integer(f_integer) :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = n
    maxdiff=0

    include 'maxdiff-inc.f90'
  end function mpimaxdiff_i0

  !> Detect the maximum difference between arrays all over a given communicator
  function mpimaxdiff_li0(n,array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    integer(f_long), intent(in) :: n !<number of elements to be controlled
    integer(f_long) :: array !< starting point of the array
    integer(f_long), dimension(:,:), allocatable :: array_glob
    integer(f_long) :: maxdiff
    include 'maxdiff-decl-inc.f90'
    ndims = int(n,kind=4)
    maxdiff=int(0,f_long)
    include 'maxdiff-inc.f90'
  end function mpimaxdiff_li0

  function mpimaxdiff_c1(array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    !> array to be checked
    character, dimension(:), intent(inout) :: array
    integer, dimension(:,:), allocatable :: array_glob
    integer :: maxdiff
    include 'maxdiff-decl-inc.f90'
    ndims = size(array)
    maxdiff=0
    include 'maxdiff-arr-inc.f90'
  end function mpimaxdiff_c1

  function mpimaxdiff_d0(n,array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    integer, intent(in) :: n !<number of elements to be controlled
    double precision :: array !< starting point of the array
    double precision, dimension(:,:), allocatable :: array_glob
    double precision :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = n
    maxdiff=0.d0

    include 'maxdiff-inc.f90'
  end function mpimaxdiff_d0

  function mpimaxdiff_d1(array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    !> array to be checked
    double precision, dimension(:), intent(inout) :: array
    double precision, dimension(:,:), allocatable :: array_glob
    double precision :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = size(array)

    maxdiff=0.d0

    include 'maxdiff-arr-inc.f90'
  end function mpimaxdiff_d1

  function mpimaxdiff_i1(array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    !> array to be checked
    integer(f_integer), dimension(:), intent(inout) :: array
    integer(f_integer), dimension(:,:), allocatable :: array_glob
    integer(f_integer) :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = size(array)

    maxdiff=0

    include 'maxdiff-arr-inc.f90'
  end function mpimaxdiff_i1

  function mpimaxdiff_li1(array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    !> array to be checked
    integer(f_long), dimension(:), intent(inout) :: array
    integer(f_long), dimension(:,:), allocatable :: array_glob
    integer(f_long) :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = size(array)

    maxdiff=0

    include 'maxdiff-arr-inc.f90'
  end function mpimaxdiff_li1

  function mpimaxdiff_i2(array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    !> array to be checked
    integer(f_integer), dimension(:,:), intent(inout) :: array
    integer(f_integer), dimension(:,:), allocatable :: array_glob
    integer :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = size(array)

    maxdiff=0

    include 'maxdiff-arr-inc.f90'
  end function mpimaxdiff_i2

  function mpimaxdiff_d2(array,root,source,comm,bcast) result(maxdiff)
    use dynamic_memory
    implicit none
    !> array to be checked
    double precision, dimension(:,:), intent(inout) :: array
    double precision, dimension(:,:), allocatable :: array_glob
    double precision :: maxdiff
    include 'maxdiff-decl-inc.f90'

    ndims = size(array)

    maxdiff=0.d0

    include 'maxdiff-arr-inc.f90'
  end function mpimaxdiff_d2

!!$end module f_maxdiff

max_stars_repo_name: adegomme/aiida-bigdft-plugin


text:
      function xnorm(sigma)
c   xnorm  - function which generates successive random samples from a
c             normal population at each call
c        xnorm  uses random number generator ranf and a transformation
c             algorithm to generate the random samples
c) the seed for xnorm can be changed by calling ranset(x)
c)   with x as the new desired seed (can use either a
c)   call statement or dummy function statement
c        the maximum value of xnorm  is limited in the program to
c             approximately three times sigma
c        usage
c           x = xnorm(sigma)
c             x is the variable into which the next sample is stored
c             sigma is standard deviation of desired sequence
	integer seed
	integer LARGE
      data c0,c1,c2,d1,d2,d3/2.515517,.802853,.010328,1.432788,.189269,.
     $001308/
c  Original code uses random number from 0 to 1
c100  u = ranf()-0.5
c     if (abs(u).lt.1.0e-03) go to 100
c
c replacement for Sun
c
  100 u = ((random(seed))/LARGE) - 0.5
      if (abs(u).lt.1.0e-03) go to 100
      sgn = u/abs(u)
      u = 1./(u*u)
      t = alog(u)
      t = sqrt(t)
      a = ((c2*t)+c1)*t+c0
      b = (((d3*t)+d2)*t+d1)*t+1.
      xnorm = (t-(a/b))*sgn*sigma
      return
      end

c $Id$


text:
      SUBROUTINE spline_akima(x,y,xx,yy,npts,iflag)
! 2012-07-27: Typo correction, J Geiger
      USE stel_kinds
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer, intent(in) :: npts        ! number of active points
      integer, intent(inout) :: iflag
      real(rprec), intent(in)  :: x
      real(rprec), intent(out) :: y
      real(rprec), dimension(npts), intent(in) :: xx, yy
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      real(rprec), dimension(-1:size(xx)+2) :: a, b, c, d
      real(rprec), dimension(-1:size(xx)+2) :: xloc, yloc
      real(rprec), dimension(-1:size(xx)+2) :: m, t, dm, p, q
      integer :: i, ix, iv, ivm1
      real(rprec) :: cl,bl,cr,br,dx
!-----------------------------------------------

!     print * , x, y
!     print * , npts
!     print * , xx
!     print * , yy
!     stop 'Test!'

      iflag = 0  !initialization
      ix = size(xx)
!     print *, size(xx)
!     print *, size(yy)
      if(npts > ix)
     >   stop'spline_akima: more active points requested than available'
      if(ix /= size(yy)) stop 'size mismatch of xx and yy!'
      iv=npts
      ivm1  = iv-1
! initialize local variables
      a = 0._DP ; b = 0._DP ; c = 0._DP ; d = 0._DP
      xloc = 0._DP ; yloc = 0._DP
      m = 0._DP ; t = 0._DP ; dm = 0._DP
      p = 0._DP ; q = 0._DP

      xloc(1:iv)=xx
      xloc(-1)= 2*xloc(1)-xloc(3)
      xloc( 0)= xloc(1)+xloc(2)-xloc(3)
      xloc(iv+2)= 2*xloc(iv)-xloc(iv-2)
      xloc(iv+1)= xloc(iv)+xloc(iv-1)-xloc(iv-2)
      yloc(1:iv)=yy
! calculate linear derivatives as far as existent
      m(1:iv-1) = (yloc(2:iv)-yloc(1:iv-1))/
     >            (xloc(2:iv)-xloc(1:iv-1))
! values for i=0, -1 and iv, iv+1 by quadratic extrapolation:
      cl = (m(2)-m(1))/(xloc(3)-xloc(1))
      bl = m(1) - cl*(xloc(2)-xloc(1))
      cr = (m(iv-2)-m(iv-1))/(xloc(iv)-xloc(iv-2))
      br = m(iv-2) - cr*(xloc(iv-1)-xloc(iv-2))
      yloc( 0)=yloc(1)+bl*(xloc( 0)-xloc(1))+
     >               cl*(xloc( 0)-xloc(1))**2
      yloc(-1)=yloc(1)+bl*(xloc(-1)-xloc(1))+
     >               cl*(xloc(-1)-xloc(1))**2
      yloc(iv+1)=yloc(iv)+br*(xloc(iv+1)-xloc(iv))+
     >               cr*(xloc(iv+1)-xloc(iv))**2
      yloc(iv+2)=yloc(iv)+br*(xloc(iv+2)-xloc(iv))+
     >               cr*(xloc(iv+2)-xloc(iv))**2
! rest of linear derivatives
      m(-1) = (yloc(0)-yloc(-1))/(xloc(0)-xloc(-1))
      m( 0) = (yloc(1)-yloc( 0))/(xloc(1)-xloc( 0))
      m(iv  ) = (yloc(iv+1)-yloc(iv  ))/(xloc(iv+1)-xloc(iv  ))
      m(iv+1) = (yloc(iv+2)-yloc(iv+1))/(xloc(iv+2)-xloc(iv+1))
! calculate weights for derivatives
      dm(-1:iv)= abs(m(0:iv+1)-m(-1:iv))
      where (dm /= 0._DP) !exclude division by zero
        p(1:iv) = dm(1:iv)/(dm(1:iv)+dm(-1:iv-2))
      end where
      where (dm /= 0._DP) !exclude division by zero
        q(1:iv) = dm(-1:iv-2)/(dm(1:iv)+dm(-1:iv-2))
      end where
      t(1:iv) = p(1:iv)*m(0:iv-1)+q(1:iv)*m(1:iv)
      where ( p(1:iv)+q(1:iv) < TINY(1._DP)) ! in case of two zeros give equal weight
        t(1:iv) = 0.5_DP*m(0:iv-1)+0.5_DP*m(1:iv)
      end where
! fix coefficients
      a = yloc
      b = t
      c(1:iv-1) = (3*m(1:iv-1)-t(2:iv)-2*t(1:iv-1))/
     >               (xloc(2:iv)-xloc(1:iv-1))
      d(1:iv-1) = (t(2:iv)+t(1:iv-1)-2*m(1:iv-1))/
     >               (xloc(2:iv)-xloc(1:iv-1))**2

! calculation
      if(x<xloc(1) .or. x>xloc(iv))then
        y=0.0
        iflag=-1
        return
      endif
      if(x == xloc(iv)) then
        y = yy(iv)
        iflag = 0
        return
      endif
      do i=1,iv-1
        if(x >= xloc(i) .and. x < xloc(i+1))then
          dx=x-xloc(i)
          y=a(i)+dx*(b(i)+dx*(c(i)+d(i)*dx))
          iflag = 0
          return
        endif
      enddo

      END SUBROUTINE spline_akima

max_stars_repo_name: ORNL-Fusion/libstell

max_issues_repo_name: ORNL-Fusion/libstell

max_forks_repo_name: ORNL-Fusion/LIBSTELL


text:
module control

  use shared_data

  implicit none

  contains

  subroutine user_control

    x_min = 0.0_num
    x_max = 1.0_num
    y_min = x_min
    y_max = x_max

    nx = 32
    ny = nx

    CFL = 0.8_num
    t_end = 10.0_num

    nsteps = 10 !set to <0 to run till t_end
    dumpfreq = 1
    use_minmod = .false.

    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = no_slip
    bc_ymax = no_slip

    grav_y  = -1.0_num

    ! uncomment one of these to overwrite user_control
    ! and users initial condition with one of the test cases

!    minion_test = .true.
!    shear_test = .true.
!    rti1_test = .true.
    circular_drop = .true.
!    blob1_test = .true.
!    vardens_adv_test = .true.

    ! DONT set more than one of the above true

  end subroutine user_control

  subroutine shear_test_control
    x_min = 0.0_num
    x_max = 1.0_num
    y_min = x_min
    y_max = x_max
    nx = 32
    ny = nx !but do force nx=ny
    CFL = 0.8_num
    t_end = 1.0_num
    dt_snapshot = 0.1_num
    nsteps = -1
    use_minmod = .false.
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = periodic
    bc_ymax = periodic
    dumpfreq = -1
    grav_x = 0.0_num
    grav_y = 0.0_num
  end subroutine shear_test_control

  subroutine minion_test_control
    ! periodic every unit length in each direction - can
    ! use to test Lx /= Ly scenarios also (dx=dy must be maintained though)
    x_min = 0.0_num
    x_max = 1.0_num
    y_min = x_min
    y_max = x_max
    nx = 32
    ny = nx !but do force nx=ny
    CFL = 0.5_num
    t_end = 0.5_num
    nsteps = -1
    use_minmod = .false.
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = periodic
    bc_ymax = periodic
    dumpfreq = 10000
    grav_x = 0.0_num
    grav_y = 0.0_num
  end subroutine minion_test_control

  subroutine vardens_adv_test_control
    x_min = 0.0_num
    x_max = 1.0_num
    y_min = x_min
    y_max = x_max
    !nx = 32   allow it to use whatever was set for easy overwriting in
    ! user control. Dont change it here in practice.
    ny = nx !but do force nx=ny
    CFL = 1.0_num
    t_end = 1000.0_num
    dt_snapshot = 10.0_num
    nsteps = -1
    use_minmod = .false.
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = periodic
    bc_ymax = periodic
    dumpfreq = -1
  end subroutine vardens_adv_test_control

  subroutine rti1_control
    x_min = -0.5_num
    x_max = 0.5_num
    y_min = -2.0_num
    y_max = 2.0_num
    nx = 16
    ny = 4*nx ! ensure nx=ny !!
    CFL = 0.5_num
    t_end = 10.0_num
    nsteps = -1
    use_minmod = .false.
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = slip_hse
    bc_ymax = outflow_hse
    dumpfreq = -1
    dt_snapshot = 1.0_num
    grav_x = 0.0_num
    grav_y = -1.0_num
  end subroutine rti1_control

  subroutine circular_drop_control
    x_min = 0.0_num
    x_max = 1.0_num
    y_min = x_min
    y_max = x_max
    nx = 256
    ny = nx !but do force nx=ny
    CFL = 0.5_num
    t_end = 5.0_num
    nsteps = -1
    use_minmod = .true.
    bc_xmin = slip
    bc_xmax = slip
    bc_ymin = slip_hse
    bc_ymax = slip_hse !outflow_hse
    dumpfreq = -1
    dt_snapshot = 0.1_num
  end subroutine circular_drop_control

  subroutine blob1_control
    x_min = 0.0_num
    x_max = 1.0_num
    y_min = x_min
    y_max = x_max
    nx = 512
    ny = nx !but do force nx=ny
    CFL = 0.8_num
    t_end = 3.0_num
    nsteps = -1
    use_minmod = .false.
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = no_slip
    bc_ymax = outflow !no_slip
    dumpfreq = -1
    dt_snapshot = 0.1_num
  end subroutine blob1_control

end module control

max_stars_repo_name: JOThurgood/SimpleCFD

max_issues_repo_name: JOThurgood/SimpleCFD

max_forks_repo_name: JOThurgood/SimpleCFD


text:
C-----CALCULATION OF THE BOUND-FREE ABSORPTION COEFFICIENT
      REAL*4 FUNCTION AHYDBF(DIVTE,WLAM,Z,EMI)

      REAL*8 WLAM,WLAMZ
      REAL*4 GAUNT(7),GHBF1(5,3),GHBF2(5),DIVTE,Z,EMI
      DATA GAUNT/7*1./
      DATA GHBF1/
     *     1.105,1.101,1.101,1.102,1.099,
     *     -7.922E-05,-3.290E-05,-1.923E-05,-1.304E-05,-9.020E-06,
     *     4.536E-09,1.152E-09,5.110E-10,2.638E-10,1.367E-10/

      DATA GHBF2/2.02848E5,1.31854E4,2.06536E-1,9.36244E4,2.21980E3/

C-----DETERMINATION OF THE LOWER AND UPPER SUMMATION INDICES N1, N2 FOR THE
C-----THE FUNCTION F(N1,N2,1,WLAM,TE) AND THE GAUNT-FACTORS FOR LEVELS N = 1 - 15
      ZZ = Z * Z
      WLAMZ = ZZ * (WLAM * 1.E8)             ! cm -> A
      N1 = 1. + SQRT( WLAMZ / 911.76 )
      N2 = MAX(N1,10)
      IF ( N1 .LE. 6 ) THEN
         DO 100 J = N1, 6
            IF ( J .EQ. 1 ) THEN
               GAUNT(J) = ( GHBF2(1) + WLAMZ * ( GHBF2(2) - WLAMZ ) ) *
     *                      GHBF2(3) / ( GHBF2(4)
     *                                  + WLAMZ * ( GHBF2(5) + WLAMZ ) )
            ELSE
               I = J - 1
               GAUNT(J) = GHBF1(I,1) + WLAMZ * ( GHBF1(I,2)
     *                                            + WLAMZ * GHBF1(I,3) )
            ENDIF
 100     CONTINUE
      ENDIF
C
C-----DETERMINATION OF THE FUNCTION F(N1,N2,1,WLAM,TE)
C
      X1 = 1.57804E+05 * ZZ * DIVTE
      DIVX1 = 1. / X1
      IF ( N1 .LT. 10 ) THEN
         N3 = N2 + 1
         G = 1. / ( N3 * N3 )
         B0 = X1 * G
         EB0 = EXP(B0)
         F = .5 * ( EB0 - 1. ) * DIVX1
         DO 200  K = N1, N2
            KINDEX = MIN(K,7)
            K2 = K * K
            K3 = K2 * K
            FK = GAUNT(KINDEX) * EXP( X1 / K2 ) / K3
            F = F + FK
 200     CONTINUE
      ELSE
         DIV = 1. / EMI
         F = .5 * DIVX1 * ( DIV - 1. )
      ENDIF
      FAC = 2.08952E-18 * WLAM * WLAMZ * WLAMZ * EXP(-X1)
      AHYDBF = FAC * F
      RETURN
      END

max_stars_repo_name: andycasey/original-oracle

max_issues_repo_name: andycasey/original-oracle

max_forks_repo_name: andycasey/original-oracle


module netcdf_layer_base
  use netcdf
  implicit none
  include "mpif.h"

  type, abstract :: NetCDF_layer_
     procedure (nf90_open), pointer, nopass :: open_file => nf90_open
     procedure (nf90_def_dim), pointer, nopass :: def_dim => nf90_def_dim
     procedure (nf90_inq_varid), pointer, nopass :: inq_varid => nf90_inq_varid
     procedure (nf90_close), pointer, nopass :: close_file => nf90_close

     procedure (integer), pointer, nopass :: put_var => nf_put_var
     procedure (integer), pointer, nopass :: get_var => nf_get_var
     procedure (integer), pointer, nopass :: put_att => nf_put_att
     procedure (integer), pointer, nopass :: def_var => nf_def_var
   contains
     procedure (create_file_signature), pass(object), deferred :: create_file
  end type NetCDF_layer_

  integer, external :: nf_put_att
  integer, external :: nf_def_var
  integer, external :: nf_put_var
  integer, external :: nf_get_var

  abstract interface

     function create_file_signature(object, path, cmode, initialsize, chunksize, ncid) result(res)
       import NetCDF_layer_
       class(NetCDF_layer_), intent(in   ) :: object
       character (len = *), intent(in   ) :: path
       integer,             intent(in   ) :: cmode
       integer, optional,   intent(in   ) :: initialsize
       integer, optional,   intent(inout) :: chunksize
       integer,             intent(  out) :: ncid
       integer                            :: res
     end function create_file_signature

  end interface

  type, extends(NetCDF_layer_) :: NetCDF_serial_
   contains
     procedure, pass(object) :: create_file => create_file_serial
  end type NetCDF_serial_

  type, extends(NetCDF_layer_) :: NetCDF_parallel_
     integer :: MPI_communicator
     integer :: default_info = MPI_INFO_NULL
   contains
     procedure, pass(object) :: create_file => create_file_parallel
  end type NetCDF_parallel_

contains

  function create_file_serial (object, path, cmode, initialsize, chunksize, ncid) result(res)
    class(NetCDF_serial_),  intent(in) :: object
    character (len = *), intent(in   ) :: path
    integer,             intent(in   ) :: cmode
    integer, optional,   intent(in   ) :: initialsize
    integer, optional,   intent(inout) :: chunksize
    integer,             intent(  out) :: ncid
    integer                            :: res

    res = nf90_create(path = path, cmode = cmode, ncid = ncid)

  end function create_file_serial

  function create_file_parallel(object, path, cmode, initialsize, chunksize, ncid) result(res)
    class(NetCDF_parallel_),intent(in) :: object
    character (len = *), intent(in   ) :: path
    integer,             intent(in   ) :: cmode
    integer, optional,   intent(in   ) :: initialsize
    integer, optional,   intent(inout) :: chunksize
    integer,             intent(  out) :: ncid
    integer                            :: res

    res = nf90_create(path = path, cmode = cmode, ncid = ncid, &
         &  comm = object%mpi_communicator, info = object%default_info)

  end function create_file_parallel

end module netcdf_layer_base


text:
INTEGER(kind = INT8), DIMENSION(0:5), PARAMETER                                 :: const_factorial_INT8 = (/ &
                                                                                1_INT8, &
                                                                                1_INT8, &
                                                                                2_INT8, &
                                                                                6_INT8, &
                                                                                24_INT8, &
                                                                                120_INT8  &
                                                                                /)

max_stars_repo_name: Guymer/fortranlib

max_issues_repo_name: Guymer/fortranlib

max_forks_repo_name: Guymer/fortranlib


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE FNDDIS (NAMECO, NAMENV, DEFOK, IXDEF,
     &                   IYDEF, IZDEF, DEFFAC, NAMLEN)
C=======================================================================

C   --*** FNDDIS *** (MESH) Find displacement variables
C   --   Written by Amy Gilkey - revised 01/25/88
C   --
C   --FNDDIS finds the displacement variables.  The first two/three nodal
C   --variables are displacement variables if and only if they begin with
C   --'D' and end with the last character of the corresponding coordinate
C   --name.
C     NAMING CONVENTION FOR DISPLACEMENT VARIABLES
C     FIRST OPTION
C     1) The nodal variables must be the first two or three variables
C        in the nodal variable names list.
C     2) Nodal variables must begin with a 'D'
C     3) The last character of the nodal variable must be equal to
C        the last character in the coordinate name.
C     4) The names of the displacement variables must be the same
C        except for the last character
C     SECOND OPTION
C     1) The nodal variables must appear together in the nodal
C        variable names list.
C     2) The nodal variables must begin with "DISP"
C     3) The last character of the nodal variable must be equal to
C        the last character in the coordinate name.
C     4) The names of the displacement variables must be the same
C        except for the last character
C     Sample of valid displacement variable names
C     Coordinate Name: X, Y, Z
C     Nodal Variable Names: DISPLX, DISPLY, DISPLZ
C     Note: Nodal variable names should be written with capital letters
C   --
C   --Parameters:
C   --   NAMECO - IN  - the coordinate names
C   --   NAMENV - IN  - the nodal variable names
C   --   DEFOK  - OUT - true iff the displacement variables were found
C   --   IXDEF  - OUT - index of x displacement variable
C   --   IYDEF  - OUT - index of y displacement variable
C   --   IZDEF  - OUT - index of z displacement variable
C   --   DEFFAC - OUT - the default displacement magnification
C   --
C   --Common Variables:
C   --   Uses NDIM, NVARNP, NSTEPW of /DBNUMS/
C   --   Uses IS3DIM of /D3NUMS/

      include 'params.blk'
      include 'dbnums.blk'
      include 'd3nums.blk'

      CHARACTER*(MXSTLN) NAMECO(*)
      CHARACTER*(NAMLEN) NAMENV(*)
      LOGICAL DEFOK

      CHARACTER CDUM
      CHARACTER*1024 STRING

      IF ((NSTEPW .GT. 0) .AND. (NVARNP .GE. NDIM)) THEN

C      --Locate displacement variables

         DEFOK = .TRUE.
         NDEF0 = 0
         LN = MAX (LENSTR (NAMENV(1)), 2)
         DO 100 I = 1, NDIM
            LC = LENSTR (NAMECO(I))
C           Check if first character is a 'D'
C           Check if all displacement names are the same other
C           than last character
C           Check if the last character in the displacement nodal variable
C           name is equal to the last character in the coordinate name.
            IF ((NAMENV(I)(1:1) .NE. 'D')
     &         .OR. (NAMENV(I)(1:LN-1) .NE. NAMENV(1)(1:LN-1))
     &         .OR. (NAMENV(I)(LN:LN) .NE. NAMECO(I)(LC:LC)))
     &         DEFOK = .FALSE.
  100    CONTINUE

         IF (.NOT. DEFOK) THEN
           DO 300 J = 2, NVARNP
             IF (LENSTR(NAMENV(J)) .GE. 5) THEN
               IF ((NAMENV(J)(1:4) .EQ. "DISP")
     &           .AND. (J+NDIM-1 .LE. NVARNP)) THEN
                 DEFOK = .TRUE.
                 NDEF0 = J-1
                 LN = MAX (LENSTR (NAMENV(J)), 2)
                 DO 200 I = 1, NDIM
                   LC = LENSTR (NAMECO(I))
                   IF ((NAMENV(I+J-1)(1:4) .NE. "DISP")
     &               .OR. (NAMENV(I+J-1)(1:LN-1) .NE.
     &                     NAMENV(J)(1:LN-1))
     &               .OR. (NAMENV(I+J-1)(LN:LN) .NE.
     &                     NAMECO(I)(LC:LC))) THEN
                     DEFOK = .FALSE.
                   ENDIF
 200             CONTINUE
                 IF (DEFOK) THEN
                    EXIT
                 ENDIF
               ENDIF
             ENDIF
 300       CONTINUE
         ENDIF

         IF (DEFOK) THEN
C           find the index for the x displacement variable
            CALL DBVIX_BL ('N', NDEF0+1, IXDEF)
C           find the index for the y displacement variable
            CALL DBVIX_BL ('N', NDEF0+2, IYDEF)
            IF (IS3DIM) THEN
C              find the index for the z displacement variable
               CALL DBVIX_BL ('N', NDEF0+3, IZDEF)
            ELSE
               IZDEF = 0
            END IF
         END IF

      ELSE
C      --If no time steps, no displacement variables
         DEFOK = .FALSE.
      END IF

      IF (.NOT. DEFOK) THEN
         WRITE (*, 10000)
10000     FORMAT (/,
     &      ' A valid set of displacement functions cannot be found.', /
     &      ' Your mesh plots will not be deformed.')

         IXDEF = 0
         IYDEF = 0
         IZDEF = 0

         DEFFAC = 0.0

      ELSE

         WRITE (*, *)
         CALL DBVTYP_BL (IXDEF, CDUM, NXDEF)
         CALL DBVTYP_BL (IYDEF, CDUM, NYDEF)
         IF (.NOT. IS3DIM) THEN
            WRITE (STRING, '(10A)') 'Displacement variables are ',
     &         NAMENV(NXDEF), ' ', NAMENV(NYDEF)
         ELSE
            CALL DBVTYP_BL (IZDEF, CDUM, NZDEF)
            WRITE (STRING, '(10A)') 'Displacement variables are ',
     &         NAMENV(NXDEF), ' ', NAMENV(NYDEF), ' ', NAMENV(NZDEF)
         END IF
         CALL SQZSTR (STRING, LSTR)
         WRITE (*, 10010) STRING(:LSTR)
10010     FORMAT (1X, 5A)

         DEFFAC = -999.0
      END IF

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      subroutine sumrun(wmelt,hascon)
c********************************************************************
c                                                                   *
c  This subroutine is called from SR CONTIN to generate a summary   *
c  of the number of runoff events and total runoff volume generated *
c  during the simulation period.                                    *
c                                                                   *
c********************************************************************
c
      real wmelt
      integer hascon
c
      include 'pmxelm.inc'
      include 'pmxpln.inc'
      include 'pmxprt.inc'
      include 'pmxhil.inc'
c
c********************************************************************
c                                                                   *
c   Common Blocks                                                   *
c                                                                   *
c********************************************************************
c
      include 'cavloss.inc'
c
      include 'cclim.inc'
      include 'cefflen.inc'
c      read: efflen(mxplan)
      include 'chydrol.inc'
c
      include 'cirriga.inc'
c
      include 'cslpopt.inc'
c
      include 'cstmflg.inc'
c
      include 'cstruc.inc'
      include 'csumirr.inc'
      include 'csumout.inc'
c
c********************************************************************
c                                                                   *
c   sumout variables updated                                        *
c     nrunot(mxplan),nrunoy(mxplan),nrunom(13,mxplan)               *
c     trunot(mxplan),trunoy(mxplan),trunom(13,mxplan)               *
c                                                                   *
c********************************************************************
cReza put IF here 3/4/94
c Change to include consideration of irrigation-induced runoff.
c dcf 8/25/94
c
      if(rain(iplane) .gt. 0.0 .or. irrund(iplane) .gt. 0.0) then
        nrunot(iplane) = nrunot(iplane) + 1
        nrunoy(iplane) = nrunoy(iplane) + 1
        nrunom(iplane) = nrunom(iplane) + 1
c
        if (hascon.ne.0) then
           trunot(iplane) = trunot(iplane) + (runoff(iplane)*1000.0)
           trunoy(iplane) = trunoy(iplane) + (runoff(iplane)*1000.0)
           trunom(iplane) = trunom(iplane) + (runoff(iplane)*1000.0)
        else
          trunot(iplane) = trunot(iplane) + (runoff(iplane)*1000.0) *
     1    efflen(iplane) / totlen(iplane)
          trunoy(iplane) = trunoy(iplane) + (runoff(iplane)*1000.0) *
     1    efflen(iplane) / totlen(iplane)
          trunom(iplane) = trunom(iplane) + (runoff(iplane)*1000.0) *
     1    efflen(iplane) / totlen(iplane)
        endif
c
c    monthly rainfall runoff totals 06-27-94 sjl
c
      else
c Reza 3/7/94.
c
c   XXX this appears to be incorrect.  Will calculate melt runoff if
c       irrigation has occured on a day without rainfall in the middle
c       of the summer.  06-22-94 04:49pm  sjl
c       however there is probably a flag for frozen soil
c       that i do not know about, i added following condition for my work
c
c
c    ...if 5 day avg temp less than freezing count runoff as melt runoff
c    otherwise it will be calculated as irrigation runoff06-27-94 10:49am  sjl
c
c XXX   The following causes problems with the water balance, because
c       melt runoff events on days in which average temp is greater
c       than 0 degrees Centigrade are not added in.  Change this so that
c       it checks melt values (wmelt(iplane)) -     dcf  7/6/94
c       if(tmnavg.lt.0)then
c
        if(wmelt .gt. 0.0)then
          nmunot(iplane) = nmunot(iplane) + 1
          nmunoy(iplane) = nmunoy(iplane) + 1
          nmunom(iplane) = nmunom(iplane) + 1
c
          if (hascon.ne.0) then
            tmunot(iplane) = tmunot(iplane) + (runoff(iplane)*1000.0)
            tmunoy(iplane) = tmunoy(iplane) + (runoff(iplane)*1000.0)
            tmunom(iplane) = tmunom(iplane) + (runoff(iplane)*1000.0)
          else
            tmunot(iplane) = tmunot(iplane) + (runoff(iplane)*1000.0) *
     1      efflen(iplane) / totlen(iplane)
            tmunoy(iplane) = tmunoy(iplane) + (runoff(iplane)*1000.0) *
     1      efflen(iplane) / totlen(iplane)
            tmunom(iplane) = tmunom(iplane) + (runoff(iplane)*1000.0) *
     1      efflen(iplane) / totlen(iplane)
          endif
        endif

      endif
c
      return
      end

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
C -----------------------------------------------------------------------
C  Program to run a least squares minimization to standard star
C  photometry to find transformation coefficients:
C  instrumental --> standard
C
C  ISYA 2004                               (c) Ignasi Ribas (IEEC, Spain)
C -----------------------------------------------------------------------
      PROGRAM FINDCOEFS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAX=1000,N=2)
      DIMENSION YINST(MAX),BINST(MAX)
      DIMENSION XMA(MAX,N),XDEP(MAX)
      DIMENSION XV(MAX),XCV(N,N),XCR(N,N),XR(N),XERR(MAX)
      DIMENSION YMA(MAX,N),YDEP(MAX)
      DIMENSION YV(MAX),YCV(N,N),YCR(N,N),YR(N),YERR(MAX)

      Open(1,file='standards.lst',status='old')
      Do 1 k=1,MAX
       Read(1,*,end=2) yinst(k),xerr(k),binst(k),sigb,ystd,ydep(k)
       xma(k,1)=1.d0
       yma(k,1)=1.d0
       xma(k,2)=binst(k)-yinst(k)
       yma(k,2)=binst(k)-yinst(k)
       yerr(k)=dsqrt(xerr(k)*xerr(k)+sigb*sigb)
       xdep(k)=ystd-yinst(k)
  1   Continue
  2   Close(1)
      ns=k-1

      CALL MINQ(XMA,XDEP,XERR,NS,XR,XCV,XCR,XV,XSIG,XR1)
      CALL MINQ(YMA,YDEP,YERR,NS,YR,YCV,YCR,YV,YSIG,YR1)

      Open(2,file='coefficients.dat')
      Write(2,50)
  50  Format(20x,'FIT RESULTS',/,20x,'-----------')
      Write(2,*)
      Write(2,100) XR(1),XR(2)
      Write(2,200) SQRT(XCV(1,1)),SQRT(XCV(2,2))
      Write(2,*)
      Write(2,400) XSIG,XR1
      Write(2,*)
      Write(2,*)
      Write(2,300) YR(1),YR(2)
      Write(2,200) SQRT(YCV(1,1)),SQRT(YCV(2,2))
      Write(2,*)
      Write(2,400) YSIG,YR1
 100  Format('y_std-y_inst =    ',f7.3,'  + ',f7.3,' * (b-y)_inst')
 200  Format('               +/-',f7.3,' +/-',f7.3)
 300  Format('   (b-y)_std =    ',f7.3,'  + ',f7.3,' * (b-y)_inst')
 400  Format(5x,'sigma = ',f7.3,5x,'correlation coeff = ',f7.3)
      Write(2,*)
      Write(2,500)
 500  Format('  y_inst','   O-C ','  sigma',' (b-y)_i','   O-C ',
     $ '  sigma')
      Do 40 i=1,ns
       calc1=xr(1)+xma(i,2)*xr(2)
       res1=xdep(i)-calc1
       calc2=yr(1)+yma(i,2)*yr(2)
       res2=ydep(i)-calc2
       Write(2,'(2(f8.3,2f7.3))') yinst(i),-res1,xerr(i),yma(i,2),res2,
     $  yerr(i)
 40   Continue
      Close(2)

C      WRITE(*,900)
C      Do 10 i=1,n
C        WRITE(*,800) (XCV(i,j),j=1,n)
C 10   Continue
C      WRITE(*,901)
C      Do 20 i=1,n
C        WRITE(*,800) (YCV(i,j),j=1,n)
C 20   Continue

 800  Format(8x,4f9.5)
 900  Format(/,/,9x,'Covariance matrix 1')
 901  Format(/,/,9x,'Covariance matrix 2')

      Stop
      End
*
*     -------------
*     LEAST SQUARES
*     -------------
*
      SUBROUTINE MINQ (A,ALTINCR,ERRY,NS,R,CV,CR,V,SIGMA,R2)
      PARAMETER (MAX=1000,N=2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(MAX,N),ALTINCR(MAX),ERRY(MAX)
      DIMENSION P(MAX,MAX),ATP(N,MAX),AT(N,MAX),C(N,N)
      DIMENSION V(MAX),CV(N,N),CR(N,N)
      DIMENSION D(N,MAX),E(MAX),F(MAX),G(MAX),R(N)
      DIMENSION LL(N),MM(N)
*
*     -------------------------------
*     INITIALIZATION OF WEIGHT MATRIX
*     -------------------------------
*
      PT=0.
      DO 5 I=1,NS
         DO 6 J=1,NS
            IF (I.EQ.J) THEN
               P(I,J)=1.D0/ERRY(I)
            ELSE
               P(I,J)=0.D0
            ENDIF
6        CONTINUE
         PT=PT+P(I,I)
5     CONTINUE
      DO 7 I=1,NS
7     P(I,I)=NS*(P(I,I)/PT)
*
*     ------------------------------------
*     CALCULATION OF THE TRANSPOSED MATRIX
*     -------------------------------------
*
      DO 10 I=1,NS
         DO 15 J=1,N
           AT(J,I)=A(I,J)
15       CONTINUE
10    CONTINUE
*
*     ----------------------------------
*     PRODUCT OF TRANSPOSED TIMES WEIGHT
*     ----------------------------------
*
      DO 20 I=1,N
         DO 25 J=1,NS
            ATP(I,J)=0.D0
           DO 30 K=1,NS
               ATP(I,J)=ATP(I,J)+AT(I,K)*P(K,J)
30         CONTINUE
25       CONTINUE
20    CONTINUE
*
*     --------------------------
*     PRODUCT WITH SYSTEM MATRIX
*     --------------------------
*
      DO 35 I=1,N
         DO 40 J=1,N
           C(I,J)=0.D0
            DO 45 K=1,NS
              C(I,J)=C(I,J)+ATP(I,K)*A(K,J)
45          CONTINUE
40       CONTINUE
35    CONTINUE
*
*     -------
*     INVERSE
*     -------
*
      CALL INVERT (C,N,DET,LL,MM)
*
*     -----------------------------------
*     PRODUCT OF INVERSE TIMES TRANSPOSED
*     -----------------------------------
*
      DO 50 I=1,N
         DO 55 J=1,NS
            D(I,J)=0.D0
            DO 60 K=1,N
              D(I,J)=D(I,J)+C(I,K)*AT(K,J)
60          CONTINUE
55       CONTINUE
50    CONTINUE
*
*     -----------------------------------
*     PRODUCT OF WEIGHT TIMES DATA MATRIX
*     -----------------------------------
*
      em1=0.
      em2=0.
      DO 65 I=1,NS
         E(I)=0.D0
         em1=em1+ALTINCR(I)
         em2=em2+A(I,2)
         DO 70 J=1,NS
            E(I)=E(I)+P(I,J)*ALTINCR(J)
70    CONTINUE
65    CONTINUE
      em1=em1/NS
      em2=em2/NS
      dyt=0.
      dxt=0.
      dxy=0.
      Do 77 I=1,NS
      DXY=DXY+(P(I,I)/NS)*(A(I,2)-em2)*(ALTINCR(I)-em1)
      DXT=DXT+(P(I,I)/NS)*(A(I,2)-em2)*(A(I,2)-em2)
77    DYT=DYT+(P(I,I)/NS)*(ALTINCR(I)-em1)*(ALTINCR(I)-em1)
      R2=DXY/(SQRT(DXT*DYT))
*
*     -------------------------
*     SOLUTION BY LEAST SQUARES
*     -------------------------
*
      DO 75 I=1,N
         R(I)=0.D0
         DO 80 J=1,NS
            R(I)=R(I)+D(I,J)*E(J)
80       CONTINUE
75    CONTINUE
*
*     ------------------------------------
*     CALCULATION OF THE COVARIANCE MATRIX
*     ------------------------------------
*
      DO 85 I=1,NS
         F(I)=0.D0
         DO 90 J=1,N
           F(I)=F(I)+A(I,J)*R(J)
90       CONTINUE
        V(I)=F(I)-ALTINCR(I)
85    CONTINUE
*
      XNUMSIG=0.D0
      DO 95 I=1,NS
         G(I)=0.D0
         DO 100 J=1,NS
            G(I)=G(I)+V(J)*P(J,I)
100      CONTINUE
         XNUMSIG=XNUMSIG+G(I)*V(I)
95    CONTINUE
      SIGMAQ=XNUMSIG/(NS-N)
      SIGMA=DSQRT(SIGMAQ)
*
      DO 105 I=1,N
         DO 110 J=1,N
            CV(I,J)=SIGMAQ*C(I,J)
110      CONTINUE
105   CONTINUE
*
*     ------------------
*     CORRELATION MATRIX
*     ------------------
*
      DO 115 I=1,N
         DO 120 J=1,N
            CR(I,J)=C(I,J)/DSQRT(C(I,I)*C(J,J))
120      CONTINUE
115   CONTINUE
*
      RETURN
      END
*
*     ----------------
*     MATRIX INVERSION
*     ----------------
*
      SUBROUTINE INVERT (A,N,D,LL,MM)
      DIMENSION A(2*N),LL(N),MM(N)
      DOUBLE PRECISION A,D,BIGA,HOLD,DABS
*     --------------------------
*     SEARCH FOR LARGEST ELEMENT
*     --------------------------
      D=1.0D0
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      LL(K)=K
      MM(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
10    IF (DABS(BIGA)-DABS(A(IJ))) 15,20,20
15    BIGA=A(IJ)
      LL(K)=I
      MM(K)=J
20    CONTINUE
*     ----------------
*     INTERCHANGE ROWS
*     ----------------
      J=LL(K)
      IF(J-K) 35,35,25
25    KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
30    A(JI) =HOLD
*     -------------------
*     INTERCHANGE COLUMNS
*     -------------------
35    I=MM(K)
      IF(I-K) 45,45,38
38    JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
40    A(JI) =HOLD
*     -------------------------------------------------------
*     DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS
*     CONTAINED IN BIGA)
*     -------------------------------------------------------
45    IF(BIGA) 48,46,48
46    D=0.0D0
      RETURN
48    DO 55 I=1,N
      IF(I-K) 50,55,50
50    IK=NK+I
      A(IK)=A(IK)/(-BIGA)
55    CONTINUE
*     -------------
*     REDUCE MATRIX
*     -------------
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I-K) 60,65,60
60    IF(J-K) 62,65,62
62    KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
65    CONTINUE
*     -------------------
*     DIVIDE ROW BY PIVOT
*     -------------------
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J-K) 70,75,70
70    A(KJ)=A(KJ)/BIGA
75    CONTINUE
*     -----------------
*     PRODUCT OF PIVOTS
*     -----------------
      D=D*BIGA
*     ---------------------------
*     REPLACE PIVOT BY RECIPROCAL
*     ---------------------------
      A(KK)=1.0D0/BIGA
80    CONTINUE
*     --------------------------------
*     FINAL ROW AND COLUMN INTERCHANGE
*     --------------------------------
      K=N
100   K=(K-1)
      IF(K) 150,150,105
105   I=LL(K)
      IF(I-K) 120,120,108
108   JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
110   A(JI) =HOLD
120   J=MM(K)
      IF(J-K) 100,100,125
125   KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
130   A(JI) =HOLD
      GO TO 100
150   RETURN
      END

max_stars_repo_name: mbadenas/Optical-Astronomy

max_issues_repo_name: mbadenas/Optical-Astronomy

max_forks_repo_name: mbadenas/Optical-Astronomy