text:
C$Procedure      CHANGU ( Change units )

      SUBROUTINE CHANGU ( ANGLE, LENGTH, TIME, MASS, CHARGE, IN,
     .                    OUT,   ERROR )

C$ Abstract
C
C     Determine units having the same dimensions of angle, length,
C     time, mass and charge as some set of input units, but with
C     respect to a "standard" set of units of the user's choosing.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     CONSTANTS
C     CONVERSION
C     PARSING
C     UNITS
C     UTILITY
C
C$ Declarations
      IMPLICIT NONE
      CHARACTER*(*)         ANGLE
      CHARACTER*(*)         LENGTH
      CHARACTER*(*)         TIME
      CHARACTER*(*)         MASS
      CHARACTER*(*)         CHARGE
      CHARACTER*(*)         IN
      CHARACTER*(*)         OUT
      CHARACTER*(*)         ERROR

      INTEGER               ROOM
      PARAMETER           ( ROOM   = 128 )

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     ANGLE      I   Default unit to use for angles (see OUNITS).
C     LENGTH     I   Default unit to use for lengths (see OUNITS).
C     TIME       I   Default unit to use for time (see OUNITS).
C     MASS       I   Default unit to use for mass (see OUNITS).
C     CHARGE     I   Default unit to use for charge (see OUNITS).
C     IN         I   Units to be transformed to the "standard".
C     OUT        O   Units that the input will be transformed to.
C     ERROR      O   Contains a description of a problem if one occurs.
C     ROOM       P   Maximum number of components in a compound unit.
C
C$ Detailed_Input
C
C     See individual entry points
C
C     ANGLE      is a string indicating which angle unit should be
C                used for outputs.
C
C     LENGTH     is a string indicating which distance unit should
C                be used for outputs.
C
C     TIME       is a string indicating which time unit should be
C                used for outputs.
C
C     MASS       is a string indicating which mass unit should be
C                used for outputs.
C
C     CHARGE     is a string indicating which charge unit should be
C                used for outputs.
C
C     IN         is the set of units associated with some measurment.
C                The dimensionally equivalent "standard" units are
C                returned in OUT.
C
C$ Detailed_Output
C
C     See individual entry points.
C
C     OUT        is the set of "standard" units that are dimensionally
C                equivalent to the input units given by IN.
C
C     ERROR      Contains a descriptive error message if the
C                subroutine call can not be executed successfully.
C
C$ Parameters
C
C     ROOM       This routine uses internal storage to construct
C                the output for TRANSU.  ROOM is the parameter that
C                describes the maximum number of components that
C                are expected for any compound unit.  The components
C                of a compound unit are
C
C                   Left parenthesis  --- '('
C                   Right parenthesis --- ')'
C                   Exponentiation    --- '**'
C                   Multiplication    --- '*'
C                   Division          --- '/'
C                   Numbers
C                   Reconized units of angle, distance, time, mass or
C                   charge.
C
C                Thus  ((10**12*KG)*(10**9*KM)**3)/((2/3)*SEC**2)
C                      ^^ ^ ^ ^^ ^^^^ ^ ^^^ ^^ ^^^^^^^^^^^  ^ ^^^
C
C                Has 31 components. (Each '^' points to the end of a
C                component).
C
C                At the time this routine was written, it was assumed
C                that compound units would have fewer than 128
C                components.
C
C$ Exceptions
C
C     1) The units used as the "standard" set must be recognized.
C        If they are not the error 'SPICE(UNKNOWNUNITS)' is signalled
C        by the entry point OUNITS.
C
C     2) If the input string IN can not be parsed as a unit, the error
C        'SPICE(INVALIDUNITS)' is signalled by the entry point TRANSU.
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine (and its entry points) are utilities that work
C     in conjunction with the general units conversion routine
C     CONVRT_2.
C
C     Here's why it is needed.
C
C        For many applications it is convenient to have command
C        driven programs.  Such commands might look like
C
C           SET FIELD OF VIEW  <number> [units]
C
C        Where "<number>" is some number that represents the size of
C        the field of view and must be supplied for the command to
C        mean anything.  The field "[units]" is an optional argument
C        that specifies the units to associate with the numeric
C        part of the command.  For example you might type any of the
C        following:
C
C           SET FIELD OF VIEW 12 DEGREES
C
C           SET FIELD OF VIEW 5  10E-3*RADIANS
C
C           SET FIELD OF VIEW 12 NANORADIANS
C
C           SET FIELD OF VIEW 6 ARCSECONDS
C
C        Allowing this kind of flexibility for inputs, gives user's
C        a friendlier interface to the program.  Instead of spending
C        time converting to some standard set of inputs, the program
C        "understands" many different units.
C
C        Ultimately, the measurements written in these expressions
C        must be converted to a set of units that the program
C        "understands."  If the above command were the only one
C        recognized by the program, the problem of converting to
C        internal units would be relatively simple.  You could just
C        list the collection of recognized units and translate them.
C        For this command such a would probably not contain more than
C        30 different units.  However, when compound units are
C        allowed such as:
C
C           KM/SEC**2,  MILES/HOUR/DAY,  AU/(100*DAYS)**2, etc.
C
C        it is no longer practical to simply list all of the possible
C        compound expressions.  Instead it is much simpler to select a
C        set of primitive units in which all compound units will be
C        expressed and used internally.  For example you might decide
C        that the fundamental units best suited to your application are:
C
C           For angles   ---  Degrees
C           For distance ---  Astronomical Units (AU)
C           For time     ---  DAYS
C           For mass     ---  KG
C           For Charge   ---  ELECTRON_CHARGES
C
C        When a measurment is encountered, your program would convert
C        it to this set of standard units automatically.  For example
C        If an input had the form
C
C           3 KM/SEC
C
C        the program would automatically convert it to the appropriate
C        number of
C
C             AU/DAYS.
C
C        In terms of the primitive units of angle, length, time, mass
C        and charge.  These two quantities are dimensionally equivalent.
C
C
C        This routine serves as the umbrella for two functions:
C
C        1) Establishing what units to use as "standard" for the
C           fundamental quanities of angle, distance, time, mass and
C           charge.  (OUNITS)
C
C        2) Computing the standard units that are dimensionally
C           equivalent to any given input units.
C
C        With the dimensionally equivalent standard units in hand,
C        it is an easy matter (as the example below illustrates)
C        to convert inputs measurments to the standard units your
C        program needs.
C
C$ Examples
C
C        To set up your default units as above:
C
C           IF ( FIRST ) THEN
C
C              CALL OUNITS ( 'DEGREES',         'AU', 'DAYS', 'KG',
C          .                 'ELECTRON_CHARGES'                    )
C
C              FIRST = .FALSE.
C
C           END IF
C
C        To translate a measurement X UNITS to the default units.
C
C           CALL TRANSU   (    UNITS, MINE       )
C           CALL CONVRT_2 ( X, UNITS, MINE, MY_X )
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber     (JPL)
C
C$ Version
C
C-    Beta Version 1.0.0, 29-MAY-1991 (WLT)
C
C-&

C
C     SPICELIB functions.
C
      INTEGER               BSRCHC
      INTEGER               LASTNB
      EXTERNAL              SCAN
C
C     Local Parameters
C
      INTEGER               NMARKS
      PARAMETER           ( NMARKS =   6 )

      INTEGER               STRSIZ
      PARAMETER           ( STRSIZ = 2*ROOM )
C
C     Local Variables
C
      CHARACTER*(2)         OP     ( NMARKS )

      CHARACTER*(8)         TCLASS (   5 )

      CHARACTER*(32)        TYPE   ( 0:5 )

      CHARACTER*(256)       O
      CHARACTER*(256)       STRING

      DOUBLE PRECISION      VALUE

      INTEGER               B
      INTEGER               BEG    ( ROOM )
      INTEGER               BLANK
      INTEGER               CLASS
      INTEGER               DIV
      INTEGER               E
      INTEGER               END    ( ROOM )
      INTEGER               EXP
      INTEGER               F
      INTEGER               I
      INTEGER               IDENT  ( ROOM )
      INTEGER               SIZE   ( 0:5 )
      INTEGER               LPAREN
      INTEGER               MULT
      INTEGER               NEST
      INTEGER               NOP
      INTEGER               NTOKNS
      INTEGER               OPLEN ( NMARKS )

C
C       Here is the range of       Character      ASCII code
C       initial characters that    ---------      ----------
C       will be used by the        ' '             32
C       "known" marks.             '('             40
C                                  ')'             41
C                                  '*'             42
C                                  '/'             47
C
C     So the required number of pointers is 47 - 32 + 5 = 20.
C
      INTEGER               OPPTR ( 20 )
      INTEGER               PASS
      INTEGER               RPAREN
      INTEGER               S
      INTEGER               START

      LOGICAL               FIRST
      LOGICAL               FOUND

C
C     Saved variables
C

      SAVE

C
C     Initial Values
C
      DATA                  TCLASS  / 'ANGLE', 'LENGTH', 'TIME',
     .                                'MASS',  'CHARGE'          /

      DATA                  FIRST   / .TRUE.  /
      DATA                  NOP     /  NMARKS /
      DATA                  OP      / ' ', '(', ')', '*', '**', '/' /

      RETURN

      ENTRY OUNITS ( ANGLE, LENGTH, TIME, MASS, CHARGE, ERROR )

C
C        On the first pass through this routine, set up the stuff
C        required for scanning the input string.
C
         IF ( FIRST ) THEN

            FIRST = .FALSE.

            CALL SCANPR ( NOP, OP, OPLEN, OPPTR )

            BLANK  = BSRCHC ( ' ',  NOP, OP )
            LPAREN = BSRCHC ( '(',  NOP, OP )
            RPAREN = BSRCHC ( ')',  NOP, OP )
            MULT   = BSRCHC ( '*',  NOP, OP )
            EXP    = BSRCHC ( '**', NOP, OP )
            DIV    = BSRCHC ( '/',  NOP, OP )

         END IF

         TYPE(0) = '1'
         TYPE(1) = ANGLE
         TYPE(2) = LENGTH
         TYPE(3) = TIME
         TYPE(4) = MASS
         TYPE(5) = CHARGE

         I     = 1
         ERROR = ' '

         DO WHILE ( I .LE. 5 )

            CALL FNDUCV ( TYPE(I), FOUND, CLASS, VALUE )

            IF ( .NOT. FOUND ) THEN
               ERROR = 'Unrecognized unit: ' // TYPE(I)
            ELSE IF ( CLASS .NE. I ) THEN

               CALL SUFFIX ( 'The',                          1, ERROR )
               CALL SUFFIX (  TCLASS(I),                     1, ERROR )
               CALL SUFFIX ( 'argument is ''',               1, ERROR )
               CALL SUFFIX (  TYPE(I),                       1, ERROR )
               CALL SUFFIX ( '''. This is not a unit ',      0, ERROR )
               CALL SUFFIX ( 'of type',                      1, ERROR )
               CALL SUFFIX (  TCLASS(I),                     1, ERROR )
               CALL SUFFIX ( '.',                            0, ERROR )

            END IF

            I = I + 1

         END DO

         DO I = 0, 5
            SIZE(I) = LASTNB( TYPE(I) )
         END DO

      RETURN

C
C     Construct the units having the same dimensions as the input
C     but that have fundamentals (angle, length, time, ... ) in the
C     form that are expected by the calling program.
C
      ENTRY TRANSU ( IN, OUT )

         STRING = IN
         O      = ' '
         NEST   = 0
         START  = 1
         F      = 0

         CALL  SCAN   ( STRING,
     .                  OP,     OPLEN, OPPTR, ROOM,   START,
     .                  NTOKNS, IDENT, BEG,   END            )

         I = 1

         DO WHILE ( I .LE. NTOKNS )

            B = BEG(I)
            E = END(I)

            IF ( IDENT(I) .EQ. BLANK ) THEN
C
C              Don't do anything....
C

            ELSE IF ( IDENT(I) .NE. 0 ) THEN

               S      = F+1
               F      = S+E-B
               O(S:F) = STRING(B:E)
C
C              We have to excercise a bit of caution.  If this
C              is an exponentiation operation, we need to just copy
C              the exponent to the output string.
C
               IF ( IDENT(I) .EQ.  EXP ) THEN

                  NEST  = 0
                  PASS  = 0

                  DO WHILE (      ( PASS .LT. 1 )
     .                       .OR. ( NEST .GT. 0 ) )

                     I       = I + 1
                     PASS    = PASS + 1

                     B       = BEG(I)
                     E       = END(I)
                     S       = F + 1
                     F       = S + B - E

                     O(S:F)  = STRING(B:E)

                     IF ( IDENT(I) .EQ. RPAREN ) THEN

                        NEST = NEST - 1

                     ELSE IF ( IDENT(I) .EQ. LPAREN ) THEN

                        NEST = NEST + 1

                     END IF

                  END DO

               END IF

            ELSE

C
C              If you get to this point, just copy the units
C              associated with the class of this token.
C
               CALL FNDUCV ( STRING(B:E), FOUND, CLASS, VALUE )

               S      = F           + 1
               F      = SIZE(CLASS) - 1 + S
               O(S:F) = TYPE(CLASS)

            END IF

            I = I + 1

         END DO

         OUT = O

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

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

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

  ifault = 0
  nullty = 0

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

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

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

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

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

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

      l = l + 1

      if ( irow == icol ) then
        exit
      end if

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

        u(k) = w / u(l)

      else

        u(k) = 0.0D+00

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

      end if

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

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

    else

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

      u(k) = sqrt ( w )

    end if

    j = j + icol

  end do

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

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

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

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

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

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

  ifault = 2
  j = 1
  k = 0

  do icol = 1, n

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

    do irow = 1, icol

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

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

      l = l + 1

      if ( irow == icol ) then
        exit
      end if

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

        u(k) = w / u(l)

      else

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

        u(k) = 0.0D+00

      end if

    end do

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

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

      u(k) = sqrt ( w )

    else

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

    end if

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

  end do

  return
end subroutine subchl
end module cholf90

max_stars_repo_name: BingzhangChen/HOT

max_issues_repo_name: BingzhangChen/HOT

max_forks_repo_name: BingzhangChen/HOT


text:
C$Procedure      UNORM ( Unit vector and norm, 3 dimensional )

      SUBROUTINE UNORM ( V1, VOUT, VMAG )

C$ Abstract
C
C     Normalize a double precision 3-vector and return its magnitude.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     VECTOR
C
C$ Declarations

      DOUBLE PRECISION  V1   ( 3 )
      DOUBLE PRECISION  VOUT ( 3 )
      DOUBLE PRECISION  VMAG

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     V1         I   Vector to be normalized.
C     VOUT       O   Unit vector V1 / |V1|.
C                    If V1 is the zero vector, then VOUT will also
C                    be zero.
C     VMAG       O   Magnitude of V1, i.e. |V1|.
C
C$ Detailed_Input
C
C     V1      This variable may contain any 3-vector, including the
C             zero vector.
C
C$ Detailed_Output
C
C     VOUT    This variable contains the unit vector in the direction
C             of V1.  If V1 is the zero vector, then VOUT will also be
C             the zero vector.
C
C     VMAG    This is the magnitude of V1.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     UNORM references a function called VNORM (which itself is
C     numerically stable) to calculate the norm of the input vector V1.
C     If the norm is equal to zero, then each component of the output
C     vector VOUT is set to zero.  Otherwise, VOUT is calculated by
C     dividing V1 by the norm.
C
C$ Examples
C
C     The following table shows how selected V1 implies VOUT and MAG.
C
C        V1                    VOUT                   MAG
C        ------------------    ------------------     ----
C        (5, 12, 0)            (5/13, 12/13, 0)       13
C        (1D-7, 2D-7, 2D-7)    (1/3, 2/3, 2/3)        3D-7
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.M. Owen       (JPL)
C     W.L. Taber      (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.2, 23-APR-2010 (NJB)
C
C        Header correction: assertions that the output
C        can overwrite the input have been removed.
C
C-    SPICELIB Version 1.0.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.0.0, 31-JAN-1990 (WMO)
C
C-&

C$ Index_Entries
C
C     3-dimensional unit vector and norm
C
C-&

C$ Revisions
C
C-    Beta Version 1.0.1, 10-JAN-1989 (WLT)
C
C     Error free specification added.
C
C-&
C

C
C     SPICELIB functions
C
      DOUBLE PRECISION   VNORM

C
C  Obtain the magnitude of V1
C
      VMAG = VNORM(V1)
C
C   If VMAG is nonzero, then normalize.  Note that this process is
C   numerically stable: overflow could only happen if VMAG were small,
C   but this could only happen if each component of V1 were small.
C   In fact, the magnitude of any vector is never less than the
C   magnitude of any component.
C
      IF ( VMAG .GT. 0.D0 ) THEN
         VOUT(1) = V1(1) / VMAG
         VOUT(2) = V1(2) / VMAG
         VOUT(3) = V1(3) / VMAG
      ELSE
         VOUT(1) = 0.D0
         VOUT(2) = 0.D0
         VOUT(3) = 0.D0
      END IF

      RETURN
      END

ext: f


text:
!====================================================================
!  Particle occupation report support routines
!  Update by KSM, based on code by Calvin in blanczoslib1.f90
!
!  Original code was not totaling to all the nucleons
!
!  versions for 'new' parallelization scheme -- FALL 2011
!===========================================================================
module pocc_mod
   ! single particle occupation table
   real(kind=4), allocatable :: spoccx(:,:,:) ! spoccx(species, eigenstate, orbit)

   integer :: aerr
contains

subroutine pocc_cleanup()
   implicit none

   if(allocated(spoccx)) deallocate(spoccx)
end subroutine pocc_cleanup

subroutine pocc_init_spoccx()
   use sporbit
   use io
   use lanczos_info
   use flagger
   implicit none
   integer :: maxorb

   call pocc_cleanup()
   if(spoccflag .and. .not. densityflag) then
      maxorb = MAX(numorb(1), numorb(2))
      allocate(spoccx(2, nkeep, maxorb), stat=aerr)
      if(aerr /= 0) call memerror("pocc_init_spoccx 1")
   else
      allocate(spoccx(1,1,1), stat=aerr)  ! always init
      if(aerr /= 0) call memerror("pocc_init_spoccx 2")
   end if
   spoccx = 0.0
end subroutine pocc_init_spoccx

!=============================================================
! Print out an occupation vector.
!   fn is the file number to write to
!   spoccx(pn, vec, orbit) has the occupation data
!   pn selects proton-1, neutron-2
!   lbl is either "   p occ:  ", or "  n occ:  "
subroutine pocc_write_occvec(fn, spoccx, i, pn, lbl)
   use sporbit
   use nodeinfo
   implicit none
   integer, intent(in) :: fn
   real(kind=4), intent(in)  :: spoccx(:,:,:)  ! single-particle occupations
   integer, intent(in) :: i, pn
   character(len=*), intent(in) :: lbl
   integer :: iorb
   if(iproc/=0)return
   do iorb = 1, numorb(1)
      if(MOD(iorb-1, 10) == 0) then
         if(iorb == 1) then
            write(fn, "(X,A10)", advance='no') lbl
         else
            write(fn, "(X)") ! newline
            write(fn, "(11X)", advance='no')
         end if
      else if(iorb /= 1) then
!         write(fn, "(A)", advance='no') ", "
         write(fn, "(A)", advance='no') "  "    ! putting in comma makes it more difficult to read with fortran
      end if
      write(fn, "(f9.5)", advance='no') spoccx(pn, i, iorb)
   end do
   write(fn, "(X)") ! cause advance
end subroutine pocc_write_occvec

! write out orbits in a form suitable for loading into Mathematica
subroutine pocc_write_orbits_m(mfn)
   use sporbit
   implicit none
   integer :: mfn
   integer :: i
   character(len=2) :: endb = "} ", endbc="},", endv

   write(mfn, '(A)') 'densSetOrbits[{ (*  {#, n, j*2, L } *)'
   do i=1,numorb(1)
      endv = endbc
      if(i == numorb(1)) endv = endb
      write(mfn, '(A,I4,A,I4,A,I4,A,I4,A)') "   {", i, ",", orbqn(1,i)%nr, ",", orbqn(1,i)%j, ",", orbqn(1,i)%l, endv
   end do
   write(mfn, '(A)') '}];'
end subroutine pocc_write_orbits_m

!=============================================================
! write out orbit quantum numbers with label so that
! occupation report is understandable
subroutine pocc_write_orbits_sub(fn)
   use sporbit
   implicit none
   integer :: fn ! file to write to
   integer :: i
   write(fn,*) ' Single particle state quantum numbers'
   write(fn,"(A)", advance='NO') 'ORBIT : '
   do i=1,numorb(1)
      write(fn, 100, advance='NO') i
   end do
   write(fn, *) " "  ! newline at end
   write(fn, "(A)", advance='NO') '   N : '
   do i=1,numorb(1)
      write(fn, 100, advance='NO') orbqn(1,i)%nr
   end do
   write(fn, *) " "  ! newline at end
   write(fn, "(A)", advance='NO') '   J : '
   do i=1,numorb(1)
      write(fn, 100, advance='NO') orbqn(1,i)%j
   end do
   write(fn, *) " "  ! newline at end
   write(fn, "(A)", advance='NO') '   L : '
   do i=1,numorb(1)
      write(fn, 100, advance='NO') orbqn(1,i)%l
   end do
   write(fn, *) " "  ! newline at end
   write(fn, *) " "  ! blank line
100 format(i6)
end subroutine pocc_write_orbits_sub

!=============================================================
! write out orbit quantum numbers with label so that
! density matrix report is understandable
! written in alternate form
subroutine pocc_write_orbits_alt(fn)
   use sporbit
   implicit none
   integer :: fn ! file to write to
   integer :: i
   write(fn,*) ' Single particle state quantum numbers'
   write(fn,"(A)", advance='yes') 'ORBIT      N     L   2 x J  '
   do i=1,numorb(1)
      write(fn, 101, advance='yes') i,orbqn(1,i)%nr,orbqn(1,i)%l, orbqn(1,i)%j
   end do

   write(fn, *) " "  ! newline at end
   write(fn, *) " "  ! blank line
101 format(4i6)
end subroutine pocc_write_orbits_alt

!=============================================================
! write out orbit quantum numbers with label so that
! occupation report is understandable
! write to both the log file and the .res file
subroutine pocc_write_orbits()
   use io
   use nodeinfo
   use flagger
   implicit none
   if(iproc == 0) then
      call pocc_write_orbits_sub(6) ! to log file
      if(writeout)  call pocc_write_orbits_sub(resultfile)
   end if
end subroutine pocc_write_orbits

subroutine pocc_write_nkeep_msg_sub(fn)
   use nodeinfo
   use lanczos_info
   implicit none
   integer :: fn

   if(iproc /= 0) return
   write(fn, '(A,i8)') 'Number of kept vectors=', nkeep
   write(fn, '(A)') ' '
end subroutine pocc_write_nkeep_msg_sub

subroutine pocc_write_nkeep_msg()
   use io
   implicit none

   call pocc_write_nkeep_msg_sub(6)
   if(writeout) call pocc_write_nkeep_msg_sub(resultfile)
end subroutine pocc_write_nkeep_msg

subroutine pocc_write_table_header()
   use nodeinfo
   use io
   implicit none
   character (len=*), parameter :: msg1 = ' State      E        Ex         J       T '
   character (len=*), parameter :: msg2 = ' State      E        Ex         '     ! added in 7.6.8

   if(iproc /= 0) return
   if(get_JT)then    !added in 7.6.8
       write(6,*) msg1
       if(writeout) write(resultfile, *) msg1
       if(densityflag) write(denresultfile, *) msg1

   else
       write(6,*) msg2
       if(writeout) write(resultfile, *) msg2
       if(densityflag) write(denresultfile, *) msg2

   endif
   return
end subroutine pocc_write_table_header

! write out table entry i
!
subroutine pocc_write_ejt(i, e, xj, xt)
   use nodeinfo
   use io
   implicit none
   integer :: i
   real(kind=egv_prec) :: e(*)
   real    :: xj, xt
   if(iproc == 0) then
      write(6, 100) i, e(i), e(i)-e(1), xj, xt
      if(writeout) write(resultfile, 100) i, e(i), e(i)-e(1), xj, xt
	  if(densityflag)write(denresultfile, 100) i, e(i), e(i)-e(1), xj, xt
   end if
100 format(i5,3x,2f10.5,2x,2f8.3)
end subroutine pocc_write_ejt
!============================================================
!
! CALLED BY:
!   lanczos_output
!   exactdiag_p
!   particle_occupation_p
!
! SUBROUTINES CALLED:
!   makespeme
!   applyspoccbundled
!   setup4obsspe
subroutine pocc_compute_spocc(i, restorejt)
   use sporbit
   use io
   use nodeinfo
   use flagger
   use obs
   use coupledmatrixelements !interaction
   use system_parameters
   use btbme_mod
   use diagh
   use apply_obs_mod
   implicit none
   logical, intent(in) :: restorejt
   integer, intent(in) :: i  ! which eigenvector
   integer :: iorb
   integer :: it

   if(spoccflag .and. .not.densityflag)then
      do iorb = 1,MAX(numorb(1),numorb(2) )
         if(numorb(1) > 0 .and. iorb <= numorb(1))then
            ! set up interaction
            pspe = 0.0
            if(numorb(2)>1)nspe = 0.0
            pspe(iorb) = 1.0
            call makespeme(1,'J')
            call makespeme(2,'J')
         end if
         if(numorb(2) > 0 .and. iorb <= numorb(2))then
            ! set up interaction
            nspe = 0.0
            if(numorb(1)>0)pspe = 0.0
            nspe(iorb) = 1.0
            call makespeme(1,'T')

            call makespeme(2,'T')
         end if
         xj2 = 0.0e0_obs_prec
         xt2 = 0.0e0_obs_prec
         call applyspoccbundled(1)
         if(np(1) > 0 .and. iorb <= numorb(1)) spoccx(1,i,iorb)= xj2   ! NOTE here xj2,xt2 are just dummy variables
         if(np(2) > 0 .and. iorb <= numorb(2)) spoccx(2,i,iorb)= xt2   ! they are NOT J2 and T2

!.......... ADDED IN 7.6.0...... FOR PARTICLE-HOLE CONJUGATION..............
         do it =1,2
            if(phconj(it))then
				spoccx(it,i,iorb)= orbqn(it,iorb)%j+1-spoccx(it,i,iorb)    ! subtract off from filled
			end if
		end do  ! it
      end do
      if(restorejt) then
   !------------- RESTORE SINGLE PARTICLE ENERGIES FOR J2, T2 OPERATORS--------
         call setup4obsspe('J')
         call makespeme(1,'J')
         call makespeme(2,'J')
		 call setup4obsmaster('T')
      end if
   end if
end subroutine pocc_compute_spocc
!..==============================================================
!
!  CALLED BY: main routine
!
!  SUBROUTINES CALLED:
!
!   setup_localvectors
subroutine particle_occupation_p
   use io
   use localvectors
   use obs
   use sporbit
   use lanczos_info
   use bvectorlib_mod
   use wfn_mod
   use system_parameters
   use coupledmatrixelements !interaction
	implicit none
   integer :: i,ikeep  ! which eigenvector
   real (kind=egv_prec), allocatable :: evec(:)
   real :: e, xj, xt
   character (len=*), parameter :: hdr_msg = "Generating new occupation report from saved eigenvectors (*.wfn)"
   if(iproc == 0) then
      write(6,*) hdr_msg
      if(writeout) write(resultfile,*) hdr_msg
   end if
   call setup_localvectors

   if(iproc == 0) then
      print *, "nfragments=", nfragments
   end if
   ! print *, "iproc=", iproc, ", v1s=", v1s, ", v1e=", v1e, ", v2s=", v2s, ", v2e=", v2e
   call wfn_read_nkeep(oldwfnfile, nkeep)
   call pocc_write_nkeep_msg()
   call pocc_write_orbits()
   call pocc_write_table_header()
   allocate(evec(nkeep), stat=aerr)
   if(aerr /= 0) call memerror("particle_occupation_p 1")
   if(allocated(pspe))deallocate(pspe,nspe)
   allocate(pspe(numorb(1)), nspe(numorb(2)), stat=aerr)
   if(aerr /= 0) call memerror("particle_occupation_p 2")
   call pocc_init_spoccx()
   twoobsflag = .true.  ! compute both two observables at same time

   ! outchoice = 'd'
   ! if(np(1)>0 .and. np(2) > 0)outchoice='b'
   ! if(np(1)>0 .and. np(2) == 0)outchoice='p'
   ! if(np(1)==0 .and. np(2) > 0)outchoice='n'
   do ikeep=1,nkeep
      ! note: write and read actually store xt^2, have to convert
      ! new interface, we say which vec to read and it checks
      i = ikeep
      call wfn_readeigenvec(oldwfnfile, frag1, fcomm1, vec1, i, e, xj, xt) ! KSM: updated
      evec(ikeep) = e
      xt2 = xt
      xt = real(-0.5 + sqrt(xt2 + 0.25))
      call pocc_compute_spocc(i, .false.)
      if(iproc /= 0) cycle
      ! now we have to figure out what to print
      call pocc_write_ejt(i, evec, xj, xt)
      if(npeff(1) > 0) call pocc_write_occvec(6, spoccx, i, 1, "   p_occ:")
      if(npeff(2) > 0) call pocc_write_occvec(6, spoccx, i, 2, "   n_occ:")
      if(writeout) then
         if(npeff(1) > 0) call pocc_write_occvec(resultfile, spoccx, i, 1, "   p_occ:")
         if(npeff(2) > 0) call pocc_write_occvec(resultfile, spoccx, i, 2, "   n_occ:")
      end if
   end do
   ! add blank line for formatting
   if(iproc == 0) then
      write(6,*) " "
      if(writeout) write(resultfile, *) " "
   end if
   deallocate(evec)
   deallocate(pspe)
   deallocate(nspe)
   call pocc_cleanup()
end subroutine particle_occupation_p

!=============================================================
!  subroutine  particle_occupation_p_orig
!  based upon subroutine expectator_p
!       master subroutine to control calculation of particle occupations
!
!  NOTE: There will have to be modifications when sections of the lanczos vector are broken into pieces
!
! routines called:
!  setup_localvectors
!  pocc_write_orbits
!   wfn_read_nkeep
!  pocc_write_nkeep_msg()
!   pocc_write_table_header
!

subroutine particle_occupation_p_orig

  use sporbit
  use precisions
  use basis
  use io
  use localvectors
  use obs
  use fragments
  use nodeinfo
  use coupledmatrixelements ! interaction
  use wfn_mod
  use bvectorlib_mod
  use diagh

  implicit none
!  include 'binterfaces.inc'

  integer nkeep
  integer ikeep
  integer i,j,n
  real :: e, xj, xt
  real (kind=egv_prec), allocatable :: evec(:)
  real (kind = 8) :: dnorm
  logical zeroflag
  integer ierr
  real (kind = lanc_prec), allocatable :: exph(:)
   integer it
   integer iorb
  character (len=*), parameter :: hdr_msg = "Generating occupation report from saved eigenvectors (*.wfn)"

  if(iproc == 0) then
	  write(6,*) hdr_msg
	  if(writeout) write(resultfile,*) hdr_msg
  end if

!........... TEMP...... WILL NEED TO FIX LATER....
  basestart(1) = 1
  basestop(1) = dimbasis
  frag1 = 1
  frag2 = 1
  call setup_localvectors

  call pocc_write_orbits()

  call wfn_read_nkeep(oldwfnfile, nkeep)

  allocate(evec(nkeep), stat=aerr)
  if(aerr /= 0) call memerror("particle_occupation_p_orig 1")

  call pocc_write_nkeep_msg()
  call pocc_write_table_header()

  twoobsflag = .false.  ! may change
  !print *, "iproc=", iproc, ", Reading eigenvectors, nkeep=", nkeep
  do ikeep = 1,nkeep
     ! new interface, we say which vec to read and it checks
     i = ikeep
     call wfn_readeigenvec(oldwfnfile, frag1, fcomm1, vec1, i, e, xj, xt) ! KSM: updated
     evec(ikeep) = e
     xt2 = xt
     xt = real(-0.5 + sqrt(xt2 + 0.25))
     ! print *, "iproc=", iproc, ", read ev, xj=", xj, ", xt=", xt
     call dnormvec_p('n','i',dnorm,zeroflag)
     if(allocated(pspe))deallocate(pspe,nspe)
     allocate(pspe(numorb(1)), nspe(numorb(2)), stat=aerr)
     if(aerr /= 0) call memerror("particle_occupation_p_orig 2")
     do it = 1,2

        allocate(exph( numorb(it)), stat=aerr)
        if(aerr /= 0) then
           call memerror("particle_occupation_p_orig 3")
           stop 5 ! for -Wuninitialized
        end if
        exph = 0.0
        do iorb = 1,numorb(it)
!--------------------- SET UP ARRAY
        pspe = 0.0
        nspe = 0.0
        if(it == 1)pspe(iorb) = 1.0
        if(it == 2)nspe(iorb) = 1.0
        call makespeme(1,'H')
        call makespeme(2,'H')

!---------------------- COMPUTE OCCUPATION
        call applyspes_expect(vec1,exph(iorb))
     enddo  ! iorb

!------------------ WRITE OUT
     if(iproc==0)then

       if(it == 1)then
          call pocc_write_ejt(i, evec, xj, xt)

          write(6,101)(exph(iorb),iorb = 1,numorb(it))
          if(writeout)write(resultfile,101)(exph(iorb),iorb = 1,numorb(it))
       else
          write(6,102)(exph(iorb),iorb = 1,numorb(it))
          if(writeout)write(resultfile,102)(exph(iorb),iorb = 1,numorb(it))
       endif

    end if
101 format(' p : ',10f9.5)
102 format(' n : ',10f9.5)
!1001  format(i5,2x,f11.4,2x,f5.1,2x,f5.1)

     deallocate(exph)
   enddo  ! it
  end do ! ikeep
  deallocate(evec)

  return
end subroutine particle_occupation_p_orig

!==================================================================

!
!  subroutine applyspes_obs
!
!  applies single-particle energies -- purely diagonal
!

subroutine applyspes_expect(vecin,exph)

  use basis
  use sectors
  use diagh
  use precisions
  use nodeinfo
  use localvectors
  use fragments
  use bmpi_mod
  implicit none

  real (kind = lanc_prec) :: vecin(dimbasis)
  real (kind = lanc_prec) :: v
  real (kind = lanc_prec) :: exph
  integer :: is,isc,jsc
  integer(kind=basis_prec) :: ip, in, ibasis
  real pspe,nspe
  integer:: ierr

!---------- LOOP OVER PROTON SECTORS-----------

  do is = 1,nsectors(1)
!................ LOOP OVER PROTON SDs in that sector

     do ip = xsd(1)%sector(is)%xsdstart,xsd(1)%sector(is)%xsdend
        pspe = pspe_h(ip)   ! the proton contribution to the single-particle energies
!-------- LOOP OVER NEUTRON SECTORS conjugate to is
        do isc = 1,xsd(1)%sector(is)%ncsectors
           jsc = xsd(1)%sector(is)%csector(isc)
!............... LOOP OVER NEUTRON SDS in sector jsc
           do in = xsd(2)%sector(jsc)%xsdstart,xsd(2)%sector(jsc)%xsdend
              nspe = nspe_h(in)  ! neutron contributions to s.p.e.s
              ibasis = pstart(ip)+nstart(in)  ! find basis state
			  if(ibasis >= v1s .and. ibasis <= v1e) then
				  v = vecin(ibasis)
				  exph = exph + v*v*(pspe+nspe)
			  end if
           enddo
        enddo !isc
     enddo  !ip

  enddo  !is

  if(.not. isfragroot) exph = 0.0
  call BMPI_ALLREDUCE(exph, 1,  MPI_SUM, icomm, ierr) ! in place

  return
end subroutine applyspes_expect
!======================================================================
end module pocc_mod

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick


text:
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
subroutine output_frame()
  use amr_commons
  use pm_commons
  use hydro_commons
  implicit none
#ifndef WITHOUTMPI
  include "mpif.h"
#endif

  integer::dummy_io,info
  integer,parameter::tag=100

  character(len=5) :: istep_str
  character(len=100) :: moviedir, moviecmd, moviefile
  character(len=100) :: moviefile1,moviefile2,moviefile3

  integer::icell,ncache,iskip,ngrid,nlevelmax_frame
  integer::ilun,nx_loc,ipout,npout,npart_out,ind,ix,iy,iz
  integer::imin,imax,jmin,jmax,ii,jj
  character(LEN=80)::fileloc
  character(LEN=5)::nchar
  real(dp)::scale,scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v
  real(dp)::xcen,ycen,zcen,delx,dely,delz
  real(dp)::xleft_frame,xright_frame,yleft_frame,yright_frame,zleft_frame,zright_frame
  real(dp)::xleft,xright,yleft,yright,zleft,zright
  real(dp)::xxleft,xxright,yyleft,yyright,zzleft,zzright
  real(dp)::dx_frame,dy_frame,dx,dx_loc
  real(dp)::dx_cell,dy_cell,dz_cell,dvol
  real(kind=8)::cell_value
  integer ,dimension(1:nvector)::ind_grid,ind_cell
  logical,dimension(1:nvector)::ok
  real(dp),dimension(1:3)::skip_loc
  real(dp),dimension(1:twotondim,1:3)::xc
  real(dp),dimension(1:nvector,1:ndim)::xx
  real(kind=8),dimension(:,:,:),allocatable::data_frame,data_frame_all
  real(kind=4),dimension(:,:),allocatable::data_single
  real(kind=8) :: z1,z2,om0in,omLin,hubin,Lbox
  real(kind=8) :: observer(3),thetay,thetaz,theta,phi,temp,ekk
  integer::igrid,jgrid,ipart,jpart,idim,icpu,ilevel
  integer::i,ig,ip,npart1
  integer::nalloc1,nalloc2

  integer,dimension(1:nvector),save::ind_part
  logical::opened
  opened=.false.

#if NDIM > 1

  ! Update counter
  imov=imov+1
  if(imov>imovout)return

  ! Determine the filename, dir, etc
  if(myid==1)write(*,*)'Computing and dumping movie frame'

  ! Determine the filename, dir, etc
  if(myid==1)write(*,*)'Computing and dumping movie frame'

  call title(imov, istep_str)
  moviedir = 'movie/'
  moviecmd = 'mkdir -p '//trim(moviedir)
  if(myid==1) write(*,*) "Writing frame ", istep_str
#ifdef NOSYSTEM
  if(myid==1)call PXFMKDIR(TRIM(moviedir),LEN(TRIM(moviedir)),O'755',info)
#else
  if(myid==1)call system(moviecmd)
#endif

  moviefile = trim(moviedir)//'info_'//trim(istep_str)//'.txt'
  if(myid==1)call output_info(moviefile)

  moviefile1 = trim(moviedir)//'dens_'//trim(istep_str)//'.map'
  moviefile2 = trim(moviedir)//'temp_'//trim(istep_str)//'.map'
  moviefile3 = trim(moviedir)//'metal_'//trim(istep_str)//'.map'

  if(levelmax_frame==0)then
     nlevelmax_frame=nlevelmax
  else
     nlevelmax_frame=levelmax_frame
  endif

  ! Conversion factor from user units to cgs units
  call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2)

  ! Local constants
  nx_loc=(icoarse_max-icoarse_min+1)
  skip_loc=(/0.0d0,0.0d0,0.0d0/)
  if(ndim>0)skip_loc(1)=dble(icoarse_min)
  if(ndim>1)skip_loc(2)=dble(jcoarse_min)
  if(ndim>2)skip_loc(3)=dble(kcoarse_min)
  scale=boxlen/dble(nx_loc)

  ! Compute frame boundaries
  xcen=xcentre_frame(1)+xcentre_frame(2)*aexp+xcentre_frame(3)*aexp**2+xcentre_frame(4)*aexp**3
  ycen=ycentre_frame(1)+ycentre_frame(2)*aexp+ycentre_frame(3)*aexp**2+ycentre_frame(4)*aexp**3
  zcen=zcentre_frame(1)+zcentre_frame(2)*aexp+zcentre_frame(3)*aexp**2+zcentre_frame(4)*aexp**3
  delx=deltax_frame(1)+deltax_frame(2)/aexp !+deltax_frame(3)*aexp**2+deltax_frame(4)*aexp**3  !Essentially comoving or physical
  dely=deltay_frame(1)+deltay_frame(2)/aexp !+deltay_frame(3)*aexp**2+deltay_frame(4)*aexp**3
  delz=deltaz_frame(1)+deltaz_frame(2)/aexp !+deltaz_frame(3)*aexp**2+deltaz_frame(4)*aexp**3
  xleft_frame=xcen-delx/2.
  xright_frame=xcen+delx/2.
  yleft_frame=ycen-dely/2.
  yright_frame=ycen+dely/2.
  zleft_frame=zcen-delz/2.
  zright_frame=zcen+delz/2.

  ! Allocate image
  allocate(data_frame(1:nx_frame,1:ny_frame,1:4))
  data_frame=0d0
  dx_frame=delx/dble(nx_frame)
  dy_frame=dely/dble(ny_frame)

  ! Loop over levels
  do ilevel=levelmin,nlevelmax_frame

     ! Mesh size at level ilevel in coarse cell units
     dx=0.5D0**ilevel

     ! Set position of cell centres relative to grid centre
     do ind=1,twotondim
        iz=(ind-1)/4
        iy=(ind-1-4*iz)/2
        ix=(ind-1-2*iy-4*iz)
        if(ndim>0)xc(ind,1)=(dble(ix)-0.5D0)*dx
        if(ndim>1)xc(ind,2)=(dble(iy)-0.5D0)*dx
        if(ndim>2)xc(ind,3)=(dble(iz)-0.5D0)*dx
     end do

     dx_loc=dx*scale
     ncache=active(ilevel)%ngrid

     ! Loop over grids by vector sweeps
     do igrid=1,ncache,nvector
        ngrid=MIN(nvector,ncache-igrid+1)
        do i=1,ngrid
           ind_grid(i)=active(ilevel)%igrid(igrid+i-1)
        end do
        ! Loop over cells
        do ind=1,twotondim
           ! Gather cell indices
           iskip=ncoarse+(ind-1)*ngridmax
           do i=1,ngrid
              ind_cell(i)=iskip+ind_grid(i)
           end do
           ! Gather cell centre positions
           do idim=1,ndim
              do i=1,ngrid
                 xx(i,idim)=xg(ind_grid(i),idim)+xc(ind,idim)
              end do
           end do
           ! Rescale position from code units to user units
           do idim=1,ndim
              do i=1,ngrid
                 xx(i,idim)=(xx(i,idim)-skip_loc(idim))*scale
              end do
           end do

           ! Check if cell is to be considered
           do i=1,ngrid
              ok(i)=son(ind_cell(i))==0.or.ilevel==nlevelmax_frame
           end do

           do i=1,ngrid
              if(ok(i))then
                 ! Check if the cell intersect the domain
                 xleft=xx(i,1)-dx_loc/2.
                 xright=xx(i,1)+dx_loc/2.
                 yleft=xx(i,2)-dx_loc/2.
                 yright=xx(i,2)+dx_loc/2.
#if NDIM>2
                 zleft=xx(i,3)-dx_loc/2.
                 zright=xx(i,3)+dx_loc/2.
                 if(    xright.lt.xleft_frame.or.xleft.ge.xright_frame.or.&
                      & yright.lt.yleft_frame.or.yleft.ge.yright_frame.or.&
                      & zright.lt.zleft_frame.or.zleft.ge.zright_frame)cycle
#else
                 if(    xright.lt.xleft_frame.or.xleft.ge.xright_frame.or.&
                      & yright.lt.yleft_frame.or.yleft.ge.yright_frame)cycle
#endif
                 ! Compute map indices for the cell
                 if(xleft>xleft_frame)then
                    imin=min(int((xleft-xleft_frame)/dx_frame)+1,nx_frame)
                 else
                    imin=1
                 endif
                 imax=min(int((xright-xleft_frame)/dx_frame)+1,nx_frame)
                 if(yleft>yleft_frame)then
                    jmin=min(int((yleft-yleft_frame)/dy_frame)+1,ny_frame)
                 else
                    jmin=1
                 endif
                 jmax=min(int((yright-yleft_frame)/dy_frame)+1,ny_frame)

                 ! Fill up map with projected mass
#if NDIM>2
                 dz_cell=min(zright_frame,zright)-max(zleft_frame,zleft)
#endif
                 do ii=imin,imax
                    xxleft=xleft_frame+dble(ii-1)*dx_frame
                    xxright=xxleft+dx_frame
                    dx_cell=min(xxright,xright)-max(xxleft,xleft)
                    do jj=jmin,jmax
                       yyleft=yleft_frame+dble(jj-1)*dy_frame
                       yyright=yyleft+dy_frame
                       dy_cell=min(yyright,yright)-max(yyleft,yleft)
                       ! Intersection volume
                       dvol=dx_cell*dy_cell
#if NDIM>2
                       dvol=dvol*dz_cell
#endif
                       data_frame(ii,jj,1)=data_frame(ii,jj,1)+dvol*uold(ind_cell(i),1)
                       data_frame(ii,jj,2)=data_frame(ii,jj,2)+dvol*uold(ind_cell(i),1)**2

                       !Get temperature
                       ekk=0.0d0
                       do idim=1,3
                          ekk=ekk+0.5*uold(ind_cell(i),idim+1)**2/uold(ind_cell(i),1)
                       enddo
                       temp=(gamma-1.0)*(uold(ind_cell(i),5)-ekk) !pressure
                       temp=temp/uold(ind_cell(i),1)*scale_T2 !temperature in K

                       data_frame(ii,jj,3)=data_frame(ii,jj,3)+dvol*uold(ind_cell(i),1)*temp !mass weighted temperature

                       if(metal)then
                       data_frame(ii,jj,4)=data_frame(ii,jj,4)+dvol*uold(ind_cell(i),6)
                       endif
                    end do
                 end do
              end if
           end do

        end do
        ! End loop over cells

     end do
     ! End loop over grids

  end do
  ! End loop over levels

  ! Convert into mass weighted
!  do ii=1,nx_frame
!     do jj=1,ny_frame
!        data_frame(ii,jj,2)=data_frame(ii,jj,2)/data_frame(ii,jj,1)
!        data_frame(ii,jj,3)=data_frame(ii,jj,3)/data_frame(ii,jj,1)
!        if(metal)then
!        data_frame(ii,jj,4)=data_frame(ii,jj,4)/data_frame(ii,jj,1)
!        endif
!     end do
!  end do
#ifndef WITHOUTMPI
  allocate(data_frame_all(1:nx_frame,1:ny_frame,1:4))
  call MPI_ALLREDUCE(data_frame,data_frame_all,nx_frame*ny_frame*4,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info)
  data_frame=data_frame_all
  deallocate(data_frame_all)
#endif
  ! Convert into mass weighted
  do ii=1,nx_frame
     do jj=1,ny_frame
        data_frame(ii,jj,2)=data_frame(ii,jj,2)/data_frame(ii,jj,1)
        data_frame(ii,jj,3)=data_frame(ii,jj,3)/data_frame(ii,jj,1)
        if(metal)then
        data_frame(ii,jj,4)=data_frame(ii,jj,4)/data_frame(ii,jj,1)
        endif
     end do
  end do

!     write(*,*) 'testing1', data_frame(100,100,1),data_frame(100,100,2)

  if(myid==1)then
     ilun=10
     allocate(data_single(1:nx_frame,1:ny_frame))
     ! Output mass weighted density
     open(ilun,file=TRIM(moviefile1),form='unformatted')
     data_single=data_frame(:,:,2)
     rewind(ilun)
     if(tendmov>0)then
        write(ilun)t,delx,dely,delz
     else
        write(ilun)aexp,delx,dely,delz
     endif
     write(ilun)nx_frame,ny_frame
     write(ilun)data_single
     close(ilun)
     ! Output mass weighted temperature
     open(ilun,file=TRIM(moviefile2),form='unformatted')
     data_single=data_frame(:,:,3)
!     write(*,*) 'testing', data_single(100,100)
     rewind(ilun)
     if(tendmov>0)then
        write(ilun)t,delx,dely,delz
     else
        write(ilun)aexp,delx,dely,delz
     endif
     write(ilun)nx_frame,ny_frame
     write(ilun)data_single
     close(ilun)
     ! Output mass weighted metal fraction
     if(metal)then
        open(ilun,file=TRIM(moviefile3),form='unformatted')
        data_single=data_frame(:,:,4)
        rewind(ilun)
        if(tendmov>0)then
           write(ilun)t,delx,dely,delz
        else
           write(ilun)aexp,delx,dely,delz
        endif
        write(ilun)nx_frame,ny_frame
        write(ilun)data_single
        close(ilun)
     endif
     deallocate(data_single)
  endif

  deallocate(data_frame)
#endif

end subroutine output_frame

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


text:
program helloworld
  print *, "helloworld"
end program helloworld

max_stars_repo_name: qzpm/Hello-language

max_issues_repo_name: qzpm/Hello-language

max_forks_repo_name: qzpm/Hello-language


text:
C    @(#)fnd_ptio.f	20.1 11/11/97
C****************************************************************
C
C     File: fnt_ptio.f
C
C     Purpose: Routine to find the PTI onwer hash index from the
C              owner number
C
c     Return code:  n > 0 : PTI owner hash number
c                   n < 0 : Not found
c
C     Note:  Arrays pti_onam() and pti_onum() are parallel, but are
C            constructed such that either name or number hashing points
c            to the same location.
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      integer function fnd_ptio (numowner)
      integer numowner
      character ownercode*3

      include 'ipfinc/parametr.inc'
      include 'ipfinc/pti_data.inc'

      integer        h, p

      if (numowner .eq. 0) then
        fnd_ptio = 0
      else
c
c       Hash owner number
c
        h = mod (numowner, MAXOWN-1) + 1
        p = htable_o(h)
        do while (p .gt. 0)
          if (numowner .eq. pti_onum(p)) then
            p = -p
          else
            p = nextptr_o(p)
          endif
        enddo
        fnd_ptio = -p
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
The wiki:wikipedia:Lolita style is popular in Japan as a street fashion, a trend in couture and a genre of pop culture.  Cosplay is related to Lolita, but the difference that it is NOT a costume, but a style of dress, and a lifestyle to some.  Its like dressing wiki:wikipedia:Punk or wiki:wikipedia:Hipster.  Lolita fashion is influenced by older styles of dress, such as the Roccoco and Victorian Eras.  There is an emphasis on modesty and innocence.

wiki:wikipedia:Maid_cafe maid cafés are influenced by the Victorian and sometimes Edwardian working class.  The style is very different from the costumey and showy French Maid costumes.  The employees of Crepe House Uni in Davis wore maidstyle or butler style but stopped because of inacceptance and misunderstanding of the uniform.

There are several subgenres of Lolita style, which are also prone to crossover:
  Sweet
  Gothic
  Country
  Kuro/Shiro
  Oji

A popular place to dress lolita is at Ciocolat

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c Copyright 2005-2021 Free Software Foundation, Inc.

c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 3 of the License, or
c (at your option) any later version.
c
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program.  If not, see <http://www.gnu.org/licenses/>.

c Ihis file is the Fortran source file for subarray.exp.  It was written
c by Wu Zhou. (woodzltc@cn.ibm.com)

        PROGRAM subarray

        character *7 str
        integer array(7)

c Initialize character array "str" and integer array "array".
        str = 'abcdefg'
        do i = 1, 7
          array(i) = i
        end do

        write (*, *) str(2:4)
        write (*, *) str(:3)
        write (*, *) str(5:)
        write (*, *) str(:)

        END PROGRAM


text:
	SUBROUTINE G2T_GSZONE ( nz, ns, iret )
C************************************************************************
C* G2T_GSZONE								*
C*									*
C* This subroutine gets the number of subzones for the NZth zone area.	*
C*									*
C* G2T_GSZONE ( NZ, NS, IRET )						*
C*									*
C* Input parameters:							*
C*	NZ		INTEGER		Zone number			*
C*									*
C* Output parameters:							*
C*	NS		INTEGER		Number of subzones in zone area	*
C*	IRET		INTEGER		Return code			*
C*					 0 = normal return		*
C**									*
C* Log:									*
C* T. Lee/SAIC		06/07	Created					*
C************************************************************************
	INCLUDE		'goftxt.cmn'
C------------------------------------------------------------------------
	iret = 0
C
C*	Retrieve zone info from the common.
C
	ns = nsubzn ( nz )
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! RUN: not %flang_fc1 %s 2>&1 | FileCheck %s
! Check errors found in folding
! TODO: test others emitted from flang/lib/Evaluate
module m
 contains
  subroutine s1(a,b)
    real :: a(*), b(:)
    !CHECK: error: DIM=1 dimension is out of range for rank-1 assumed-size array
    integer :: ub1(ubound(a,1))
    !CHECK-NOT: error: DIM=1 dimension is out of range for rank-1 assumed-size array
    integer :: lb1(lbound(a,1))
    !CHECK: error: DIM=0 dimension is out of range for rank-1 array
    integer :: ub2(ubound(a,0))
    !CHECK: error: DIM=2 dimension is out of range for rank-1 array
    integer :: ub3(ubound(a,2))
    !CHECK: error: DIM=0 dimension is out of range for rank-1 array
    integer :: lb2(lbound(b,0))
    !CHECK: error: DIM=2 dimension is out of range for rank-1 array
    integer :: lb3(lbound(b,2))
  end subroutine
  subroutine s2
    integer, parameter :: array(2,3) = reshape([(j, j=1, 6)], shape(array))
    integer :: x(2, 3)
    !CHECK: error: Invalid 'dim=' argument (0) in CSHIFT
    x = cshift(array, [1, 2], dim=0)
    !CHECK: error: Invalid 'shift=' argument in CSHIFT: extent on dimension 1 is 2 but must be 3
    x = cshift(array, [1, 2], dim=1)
  end subroutine
  subroutine s3
    integer, parameter :: array(2,3) = reshape([(j, j=1, 6)], shape(array))
    integer :: x(2, 3)
    !CHECK: error: Invalid 'dim=' argument (0) in EOSHIFT
    x = eoshift(array, [1, 2], dim=0)
    !CHECK: error: Invalid 'shift=' argument in EOSHIFT: extent on dimension 1 is 2 but must be 3
    x = eoshift(array, [1, 2], dim=1)
    !CHECK: error: Invalid 'boundary=' argument in EOSHIFT: extent on dimension 1 is 3 but must be 2
    x = eoshift(array, 1, [0, 0, 0], 2)
  end subroutine
  subroutine s4
    integer, parameter :: array(2,3) = reshape([(j, j=1, 6)], shape(array))
    logical, parameter :: mask(*,*) = reshape([(.true., j=1,3),(.false., j=1,3)], shape(array))
    integer :: x(3)
    !CHECK: error: Invalid 'vector=' argument in PACK: the 'mask=' argument has 3 true elements, but the vector has only 2 elements
    x = pack(array, mask, [0,0])
  end subroutine
  subroutine s5
    logical, parameter :: mask(2,3) = reshape([.false., .true., .true., .false., .false., .true.], shape(mask))
    integer, parameter :: field(3,2) = reshape([(-j,j=1,6)], shape(field))
    integer :: x(2,3)
    !CHECK: error: Invalid 'vector=' argument in UNPACK: the 'mask=' argument has 3 true elements, but the vector has only 2 elements
    x = unpack([1,2], mask, 0)
  end subroutine
  subroutine s6
    !CHECK: error: POS=-1 out of range for BTEST
    logical, parameter :: bad1 = btest(0, -1)
    !CHECK: error: POS=32 out of range for BTEST
    logical, parameter :: bad2 = btest(0, 32)
    !CHECK-NOT: error: POS=33 out of range for BTEST
    logical, parameter :: ok1 = btest(0_8, 33)
    !CHECK: error: POS=64 out of range for BTEST
    logical, parameter :: bad4 = btest(0_8, 64)
  end subroutine
  subroutine s7
    !CHECK: error: SHIFT=-33 count for ishft is less than -32
    integer, parameter :: bad1 = ishft(1, -33)
    integer, parameter :: ok1 = ishft(1, -32)
    integer, parameter :: ok2 = ishft(1, 32)
    !CHECK: error: SHIFT=33 count for ishft is greater than 32
    integer, parameter :: bad2 = ishft(1, 33)
    !CHECK: error: SHIFT=-65 count for ishft is less than -64
    integer(8), parameter :: bad3 = ishft(1_8, -65)
    integer(8), parameter :: ok3 = ishft(1_8, -64)
    integer(8), parameter :: ok4 = ishft(1_8, 64)
    !CHECK: error: SHIFT=65 count for ishft is greater than 64
    integer(8), parameter :: bad4 = ishft(1_8, 65)
  end subroutine
  subroutine s8
    !CHECK: error: SHIFT=-33 count for shiftl is negative
    integer, parameter :: bad1 = shiftl(1, -33)
    !CHECK: error: SHIFT=-32 count for shiftl is negative
    integer, parameter :: bad2 = shiftl(1, -32)
    integer, parameter :: ok1 = shiftl(1, 32)
    !CHECK: error: SHIFT=33 count for shiftl is greater than 32
    integer, parameter :: bad3 = shiftl(1, 33)
    !CHECK: error: SHIFT=-65 count for shiftl is negative
    integer(8), parameter :: bad4 = shiftl(1_8, -65)
    !CHECK: error: SHIFT=-64 count for shiftl is negative
    integer(8), parameter :: bad5 = shiftl(1_8, -64)
    integer(8), parameter :: ok2 = shiftl(1_8, 64)
    !CHECK: error: SHIFT=65 count for shiftl is greater than 64
    integer(8), parameter :: bad6 = shiftl(1_8, 65)
  end subroutine
  subroutine s9
    integer, parameter :: rank15(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) = 1
    !CHECK: error: SOURCE= argument to SPREAD has rank 15 but must have rank less than 15
    integer, parameter :: bad1 = spread(rank15, 1, 1)
    integer, parameter :: matrix(2, 2) = reshape([1, 2, 3, 4], [2, 2])
    !CHECK: error: DIM=0 argument to SPREAD must be between 1 and 3
    integer, parameter :: bad2 = spread(matrix, 0, 1)
    !CHECK: error: DIM=4 argument to SPREAD must be between 1 and 3
    integer, parameter :: bad3 = spread(matrix, 4, 1)
  end subroutine
  subroutine warnings
    real, parameter :: ok1 = scale(0.0, 99999) ! 0.0
    real, parameter :: ok2 = scale(1.0, -99999) ! 0.0
    !CHECK: SCALE intrinsic folding overflow
    real, parameter :: bad1 = scale(1.0, 99999)
  end subroutine
end module


text:
      SUBROUTINE getbfld(amat_i, amat_p, data_array, r12sqr,
     1  gsqrt, orsq, gobsr1, gobsz1, r12, z12, kcbfld)
      USE vmec_main
      USE vmec_params, ONLY: signgs
      USE vsvd
      USE realspace, ONLY: wint
      USE vspline, ONLY: hthom, hstark
      USE mgrid_mod, ONLY: nbfldn, nbsets, nbcoils, needbfld, abcoil,
     1                     plbfld
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER kcbfld
      REAL(rprec), DIMENSION(isnodes,*) :: amat_i
      REAL(rprec), DIMENSION(ipnodes,*) :: amat_p
      REAL(rprec), DIMENSION(*) :: data_array
      REAL(rprec), DIMENSION(nrzt) :: r12sqr,
     1  gsqrt, orsq, gobsr1, gobsz1, r12, z12
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: js, n1, m1, iobs, isym, iobsmax, iloop,
     1   msym, nsym, indexbfld, l, lk
      REAL(rprec), DIMENSION(:), ALLOCATABLE :: gobsz2, gobsr2
      REAL(rprec):: tpisg, sumir, sumiz, wscaleb, sumpr,
     1    sumpz, deltab, plasbfld, coscoil, sincoil, t2, t1
C----------------------------------------------

!     IRESIDUE > 0, OTHERWISE FLUX_INIT NOT CALLED YET
      kcbfld = 0
      IF (iresidue.le.0 .or. (nbfldn.eq.0 .and. iequi.eq.0))return

      ALLOCATE( gobsz2(nrzt), gobsr2(nrzt), stat=l)
      IF (l .ne. 0) STOP 'allocation problem in getbfld'
!
!       COMPUTE "GREEN'S FUNCTION" KERNEL ONLY FOR NEEDED OBSERVATION POINTS
!       (OR FOR FINAL OUTPUT IF IEQUI=1)
!       GOBSR1 = BR, GOBSZ1 = BZ
!       IF A (BR,BZ) OR (BRHO,BTHETA) PAIR, JUST CALL GRNBFLD ONCE
!
      IF (ANY(rm2(2:ns) .eq. zero)) STOP 'rm2 = 0'
      tpisg = twopi * signgs                !Positive volume integral
      DO n1 = 1, nbsets
         DO m1 = 1, nbcoils(n1)
            isym = needbfld(m1,n1)
            IF (isym.eq.needit .or. iequi.eq.1) THEN
               CALL grnbfld(r12sqr,r12,z12,gobsr1,gobsz1,nrzt,m1,n1)
               DO l = 2,nrzt
                 gobsr2(l) = gobsr1(l)*orsq(l)
                 gobsr1(l) = gobsr1(l)*gsqrt(l)
                 gobsz2(l) = gobsz1(l)*orsq(l)
                 gobsz1(l) = gobsz1(l)*gsqrt(l)
               END DO
!
!       DO INTEGRAL OVER ANGLES (ALL INTEGRALS ARE FROM THETA=0,TWOPI)
!
               DO js = 2, ns
                  sumir = zero
                  sumiz = zero
                  sumpr = zero
                  sumpz = zero
                  DO lk = js,nrzt,ns
                    sumir = sumir + gobsr2(lk)*wint(lk)
                    sumpr = sumpr + gobsr1(lk)*wint(lk)
                    sumiz = sumiz + gobsz2(lk)*wint(lk)
                    sumpz = sumpz + gobsz1(lk)*wint(lk)
                  ENDDO
                  imb(js,m1,n1,1) = tpisg*sumir
                  pmb(js,m1,n1,1) = (-tpisg*sumpr) + imb(js,m1,n1,1)
     1               /rm2(js)
                  imb(js,m1,n1,2) = tpisg*sumiz
                  pmb(js,m1,n1,2) = (-tpisg*sumpz) + imb(js,m1,n1,2)
     1               /rm2(js)
               END DO

            ELSE IF (isym .eq. ISAMECOIL) THEN       !Same coil position as previous coil
              DO js = 2,ns
                imb(js,m1,n1,1) = imb(js,m1-1,n1,1)
                pmb(js,m1,n1,1) = pmb(js,m1-1,n1,1)
                imb(js,m1,n1,2) = imb(js,m1-1,n1,2)
                pmb(js,m1,n1,2) = pmb(js,m1-1,n1,2)
              ENDDO
            ENDIF
          ENDDO   !m1
        ENDDO   !n1

!
!       CHECK FOR SYMMETRIC COIL (MAY BE IN DIFFERENT COIL SET,
!       SO HAD TO MOVE OUT OF M1,N1 LOOP ABOVE)
!
        DO n1 = 1,nbsets
          DO m1 = 1,nbcoils(n1)
            isym = needbfld(m1,n1)
            IF (isym .ge. ISYMCOIL) THEN
              msym = 1 + (isym-1)/nbsets
              nsym = 1 + MOD(isym-1,nbsets)
              DO js = 2,ns                     !BR(-Z) = -BR(Z), BZ(-Z) = BZ(Z)
                imb(js,m1,n1,1) =-imb(js,msym,nsym,1)
                pmb(js,m1,n1,1) =-pmb(js,msym,nsym,1)
                imb(js,m1,n1,2) = imb(js,msym,nsym,2)
                pmb(js,m1,n1,2) = pmb(js,msym,nsym,2)
              ENDDO
            ENDIF
          ENDDO
        ENDDO

!
!       COMPUTE SPLINE MATRIX ELEMENTS BY INTEGRATING OVER RADIUS
!
        DO 2000 iloop = 0,iequi                !iequi = 0 normally, = 1 at END
          DO n1 = 1, nbsets
            iobsmax = nbfld(n1)
            IF (iloop .eq. 1) iobsmax = nbcoils(n1)
            IF (iobsmax .gt. 0) THEN
              DO 1000 iobs = 1, iobsmax
                indexbfld = indxbfld(iobs,n1)
                IF (iloop .eq. 1) indexbfld = iobs
                IF (indexbfld .le. 0) GOTO 1000
                coscoil = COS( abcoil(indexbfld,n1) )
                sincoil = SIN( abcoil(indexbfld,n1) )
                DO js = 2,ns
                  pmb(js,0,n1,1) = ochip(js) *
     >            (pmb(js,indexbfld,n1,1)*coscoil +
     >             pmb(js,indexbfld,n1,2)*sincoil)
                  imb(js,0,n1,1) = ovrm2(js) *
     >            (imb(js,indexbfld,n1,1)*coscoil +
     >             imb(js,indexbfld,n1,2)*sincoil)
                END DO

                     IF (iloop .eq. 0) THEN
                        deltab = plbfld(indexbfld,n1)
                        kcbfld = kcbfld + 1

                        CALL splinint (imb(1,0,n1,1), current,
     1                     amat_i(1,kcbfld), hstark, u_ib, u1_ib,
     2                     w_ib, w1_ib, nk_ib, isnodes, intder, ns)

                        CALL splinint (pmb(1,0,n1,1), presint,
     1                     amat_p(1,kcbfld), hthom, u_pb, u1_pb, w_pb,
     2                     w1_pb, nk_pb, ipnodes, intder, ns)

                        wscaleb = one/sigma_b(indexbfld,n1)
                        data_array(kcbfld) = wscaleb*deltab
                        t2 = mu0*pthommax      !!*pfac moved to getthom

                        amat_i(:,kcbfld) = wscaleb*amat_i(:,kcbfld)
                        wscaleb = wscaleb*t2
                        amat_p(:,kcbfld) = wscaleb*amat_p(:,kcbfld)

                     ELSE      !Store plasma fluxes in EXTFLX for output
                        plasbfld = zero
                        DO js = 2, ns
                           t1 = current(js)*iotaf(js) - current(js-1)*
     1                        iotaf(js - 1)
                           plasbfld = plasbfld + pmb(js,0,n1,1)*
     1                        presph(js) + imb(js,0,n1,1)*t1
                        END DO
                        plbfld(iobs,n1) = plasbfld
                     ENDIF
 1000         CONTINUE
            ENDIF
          ENDDO       !n1
 2000   CONTINUE

      DEALLOCATE( gobsz2, gobsr2, stat=l)

      END SUBROUTINE getbfld


text:
SUBROUTINE prm_orbext (PRMfname, found)

! ----------------------------------------------------------------------
! SUBROUTINE: prm_orbext.f90
! ----------------------------------------------------------------------
! Purpose:
!  External Orbit settings
!  Precise orbit data (sp3) interpolation or computation of Keplerian orbits
! is set here.
!
!  The formulated orbit may be used for the following purposes:
!  1. External orbit comparison
!  2. Pseudo-observations within the dynamic orbit estimation procedure
! ----------------------------------------------------------------------
! Input arguments:
! - PRMfname:		Configuration file name for the orbit parameterization
!
! Output arguments:
! 	External Orbit is stored in the allocatable arrays orbext_ICRF, orbext_ITRF, orbext_kepler
! 	that are set as global variables in the module mdl_param.f03.
!
! - orbext_ICRF: 	Orbit array (Nx8) in inertial frame (ICRF) including the state vector per epoch
! 					  Collumns desciption per epoch:
!               	- Modified Julian Day number (including the fraction of the day)
!					- Seconds since 00h
!					- Position vector (m)
!					- Velocity vector (m/sec)
! - orbext_ITRF:	Orbit array (Nx8) in terrestrial frame (ITRF)
! 					  Collumns desciption per epoch: similar to the orbext_ICRF
! - orbext_kepler:	Orbit array (Nx8) in inertial frame including the Keplerian elements per epoch
! 					  Collumns desciption per epoch:
!               	- Modified Julian Day number (including the fraction of the day)
!					- Seconds since 00h
!          			- semi-major axis  (m)
!         			- eccentricity
!    	     		- inclination (degrees)
!     				- right ascension of the ascending node (degrees)
!     				- argument of perigee (degrees)
!         			- eccentric anomaly (degrees)
! ----------------------------------------------------------------------
! Author :	Dr. Thomas Papanikolaou, Cooperative Research Centre for Spatial Information, Australia
! Created:	20 September 2017
! ----------------------------------------------------------------------

      USE mdl_precision
      USE mdl_num
      USE mdl_param
      !USE mdl_arr
      USE m_keplerorb
      USE m_rso
      USE m_interporb
      USE m_orbT2C
      USE m_orbC2T
      USE m_sp3
      use pod_yaml
      use mdl_config
      IMPLICIT NONE

! ----------------------------------------------------------------------
! Dummy arguments declaration
! ----------------------------------------------------------------------
! IN
      CHARACTER (LEN=100), INTENT(IN) :: PRMfname
! OUT
      LOGICAL, INTENT (OUT) :: found
!
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Local Variables declaration
! ----------------------------------------------------------------------
!      REAL (KIND = prec_d), DIMENSION(:,:), ALLOCATABLE :: orbext_ICRF, orbext_ITRF, orbext_kepler
      INTEGER (KIND = prec_int2) :: data_opt
      CHARACTER (LEN=10) :: FRAME_EXT
      CHARACTER (LEN=300) :: fname_orb
      CHARACTER (LEN=300) :: fname_orb_0, fname_orb_1, fname_orb_2
      CHARACTER (LEN=300) :: fname_orbint
      CHARACTER (LEN=300) :: fname_write
      INTEGER (KIND = prec_int8) :: NPint
      INTEGER (KIND = prec_int8) :: interpstep
      INTEGER (KIND = prec_int8) :: sz1, sz2
      INTEGER (KIND = prec_int8) :: Ndays
      INTEGER (KIND = prec_int4) :: Zo_el
	  REAL (KIND = prec_q) :: GMearth
	  REAL (KIND = prec_d) :: Zo(6), Sec0, MJDo
      INTEGER (KIND = prec_int2) :: AllocateStatus, DeAllocateStatus
! ----------------------------------------------------------------------
      INTEGER (KIND = prec_int8) :: i, read_i
      INTEGER (KIND = prec_int2) :: UNIT_IN, ios
      INTEGER (KIND = prec_int2) :: ios_line, ios_key, ios_data
      INTEGER (KIND = prec_int2) :: space_i
      CHARACTER (LEN=7) :: Format1, Format2, Format3
      CHARACTER (LEN=500) :: line_ith
      CHARACTER (LEN=150) :: word1_ln, word_i, t0
! ----------------------------------------------------------------------
      CHARACTER (LEN=30) :: fmt_line
      CHARACTER (LEN=1) :: GNSSid
      INTEGER (KIND = prec_int4) :: PRN_no
! ----------------------------------------------------------------------
      CHARACTER (LEN=3) :: time
	  REAL (KIND = prec_d) :: mjd , mjd_TT, mjd_GPS, mjd_TAI, mjd_UTC
! ----------------------------------------------------------------------
      REAL (KIND = prec_q), DIMENSION(:,:), ALLOCATABLE :: clock_matrix
      INTEGER(KIND=prec_int2) time_system

! ----------------------------------------------------------------------
! Orbit parameterization INPUT file read:
! ----------------------------------------------------------------------
found = .false.

! ----------------------------------------------------------------------
      UNIT_IN = 9
      Format1 = '(A)'
      Format2 = '(F)'
      Format3 = '(I100)'
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Read input file
i = 0
fname_orb = yml_ext_orbit_filename
interpstep = yml_ext_orbit_steps
NPint = yml_ext_orbit_points
data_opt = yml_ext_orbit_opt

if (.not. yaml_found) then
! ----------------------------------------------------------------------
! Open .in file
      OPEN (UNIT = UNIT_IN, FILE = TRIM (PRMfname), IOSTAT = ios)
      IF (ios /= 0) THEN
         PRINT *, "Error in opening file:", PRMfname
         PRINT *, "OPEN IOSTAT=", ios
      END IF

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

DO

READ (UNIT=UNIT_IN,FMT=Format1,IOSTAT=ios_line) line_ith
i = i + 1
! PRINT *, "READ Line (i,ios):", i, ios_line

! ----------------------------------------------------------------------
! End of file
         IF (ios_line < 0) THEN
!            PRINT *, "End of file, i=", i
            EXIT
         END IF
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! 1st Word of Line ith
READ (line_ith, * , IOSTAT=ios_data) word1_ln  ! 1st word
!READ (line_ith, * , IOSTAT=ios_data) word1_ln, charN
! ----------------------------------------------------------------------
!PRINT *, "word1_ln: ", word1_ln

! ----------------------------------------------------------------------
! Parameters Keywords read
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! External Orbit
! ----------------------------------------------------------------------
! 1. Orbit data in sp3 format (including position and velocity vectors) ! Old option: RSO orbit data by GFZ
! 2. Interpolated orbit based on Lagrange interpolation of daily sp3 data
! 3: Keplerian orbit
! 4. Interpolated orbit (3 days arc) based on Lagrange interpolation of 3 boundary sp3 data files
! 5. Orbit: Position vector from sp3 data; Velocity vector approximated through sp3 position differences without applying interpolation
IF (word1_ln == "orbit_external_opt") THEN
   READ ( line_ith, FMT = * , IOSTAT=ios_key ) word_i, data_opt
END IF
yml_ext_orbit_opt = data_opt
! ----------------------------------------------------------------------

! Reference Frame of external orbit
! 1. ITRF
! 2. ICRF
IF (word1_ln == "orbit_ext_frame") THEN
   READ ( line_ith, FMT = * , IOSTAT=ios_key ) word_i, FRAME_EXT
   if (FRAME_EXT == 'ITRF') yml_ext_orbit_frame = ITRF
   if (FRAME_EXT == 'ICRF') yml_ext_orbit_frame = ICRF
END IF

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

! ----------------------------------------------------------------------
! GNSS orbit data (sp3) file name
IF (word1_ln == "orbit_filename" .and. .not. ooption_on_command_line) THEN
   READ ( line_ith, FMT = * , IOSTAT=ios_key ) word_i, fname_orb
END IF
! Case 4:
!fname_orb_0 = 'wum19283.sp3'
!fname_orb_1 = 'wum19284.sp3'
!fname_orb_2 = 'wum19285.sp3'
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Numerical Interpolation or Kepler orbit interval (sec)
IF (word1_ln == "orbit_interp_step") THEN
   READ ( line_ith, FMT = * , IOSTAT=ios_key ) word_i, interpstep
END IF
!interpstep = integstep ! Module mdl_param.f03
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Interpolated Orbit:	Cases: 2,4
! ----------------------------------------------------------------------
! Number of data points used in Lagrange interpolation
IF (word1_ln == "orbit_interp_points") THEN
   READ ( line_ith, FMT = * , IOSTAT=ios_key ) word_i, NPint
END IF
!NPint = 12

! Output file name for writing interpolated orbit (Cases: 2, 4)
!IF (word1_ln == "orbit_interp_fname") THEN
!   READ ( line_ith, FMT = * , IOSTAT=ios_key ) word_i, fname_orbint
!END IF
!fname_orbint = 'orb_interp.out'
! ----------------------------------------------------------------------

END DO
CLOSE (UNIT=UNIT_IN)
! Close of input parameterization file
! ----------------------------------------------------------------------
end if

data_opt = yml_ext_orbit_opt
!if (yml_interpolate_start /= 0.d0) data_opt = TYPE_INTERP

! ----------------------------------------------------------------------
! Compute the interpolated orbit arc (based on GNSS sp3 orbits) or the Keplerian orbit arc
! Form the orbit arrays (satellite state vector per epoch) to be used for external orbit comparison
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Case 1: Orbit SP3 data with position and velocity vectors (Old: RSO data by GFZ)
! ----------------------------------------------------------------------
If (data_opt == TYPE_SP3) then

! PRN: GNSS constellation ID letter + Satellite number
fmt_line = '(A1,I2.2)'
READ (PRN, fmt_line , IOSTAT=ios) GNSSid, PRN_no
! Read RSO data and produce the orbit array (orbext_ITRF) for the input satellite PRN
!CALL rso (fname_orb, PRN_no, orbext_ITRF)

! Read sp3 orbit data that include position and velocity vectors
If (yml_ext_orbit_frame == ICRF) Then
! ICRF
Call sp3 (fname_orb, PRN, orbext_ICRF, yml_interpolate_start, clock_matrix, time_system, found)
! Orbit transformation ICRF to ITRF
if (found) Call orbC2T (orbext_ICRF, time_system, orbext_ITRF)
Else If (yml_ext_orbit_frame == ITRF) Then
! ITRF
Call sp3 (fname_orb, PRN, orbext_ITRF, yml_interpolate_start, clock_matrix, time_system, found)
! Orbit transformation ITRF to ICRF
if (found) Call orbT2C (orbext_ITRF, time_system, orbext_ICRF)
End IF
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Case 2: Orbit based on Lagrange interpolation of sp3 data e.g. IGS final/rapid (15 min); MGEX (5 min) IGS rapid orbits; interval 15 min
! ----------------------------------------------------------------------
Else if (data_opt == TYPE_INTERP) then

! Interpolated Orbit: Read sp3 orbit data and apply Lagrange interpolation
If (yml_ext_orbit_frame == ICRF) Then
CALL interp_orb (fname_orb, PRN, interpstep, NPint, yml_interpolate_start, orbext_ICRF, time_system, found)
! Orbit transformation ICRF to ITRF
if (found) Call orbC2T (orbext_ICRF, time_system, orbext_ITRF)

Else If (yml_ext_orbit_frame == ITRF) Then

CALL interp_orb (fname_orb, PRN, interpstep, NPint, yml_interpolate_start, orbext_ITRF, time_system, found)
! Orbit transformation ITRF to ICRF
if (found) Call orbT2C (orbext_ITRF, time_system, orbext_ICRF)

END IF
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Case 3:	Keplerian orbit
! ----------------------------------------------------------------------
Else if (data_opt == TYPE_KEPLER) then

! Keplerian orbit arc length in number of days
Ndays = INT(orbarc / (24.0D0 * 3600.0D0)) + 1 ! orbarc: global in Module mdl_param.f03

! ----------------------------------------------------------------------
! Initial Epoch
MJDo = MJD_to ! Module mdl_param.f03
Sec0 = SEC_to ! Module mdl_param.f03

! Time scale transformation
mjd = MJD_to
CALL time_TT (mjd , mjd_TT, mjd_GPS, mjd_TAI, mjd_UTC)

!If (time == 'TT') Then
If (yml_time_scale == GPS_time) then
    MJDo = mjd_GPS
Else if (yml_time_scale == UTC_time) then
    MJDo = mjd_UTC
Else if (yml_time_scale == TAI_time) then
    MJDo = mjd_TAI
End If
Sec0 = (MJDo - INT(MJDo)) * (86400.D0)
! ----------------------------------------------------------------------

! Initial State Vector in ICRF
Zo = SVEC_Zo ! Module mdl_param.f03

! Initial state vector type: SVEC_Zo is computed in ICRF >> Set Zo_el to 2
! 1. Keplerian Elements
! 2. State Vector (Position & Velocity) cartesian coordinates (in ICRF)
Zo_el = 2

! Earth Gravity constant:  Module mdl_num.f90
GMearth = GM_global

! Computation of the Kepler orbit arc
CALL keplerorb (GMearth, MJDo, Sec0, Zo, Zo_el, Ndays, interpstep, orbext_kepler, orbext_ICRF, found)

! ----------------------------------------------------------------------
! Temporary
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! Orbit Transformation to ITRF is not applied
! The orbit comparison based on Kepler orbits is perfromed only in the inertial frame
if (found) then
sz1 = size(orbext_ICRF, DIM = 1)
sz2 = size(orbext_ICRF, DIM = 2)
ALLOCATE (orbext_ITRF(sz1,sz2), STAT = AllocateStatus)
orbext_ITRF = orbext_ICRF
end if
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Deactivated cases
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Case 4: Interpolated orbit 3 days arc based on Lagrange interpolation of 3 daily sp3 data files
! ----------------------------------------------------------------------
else if (data_opt == TYPE_3LINTERP) then
	! Performing interpolation individually for each sp3 orbit file
	! Extrapolation is applied in order to avoid the gaps at the end of the daily file (sp3 last epoch is 23h 45min 00sec)
	!CALL orb3sp3(fname_orb_0, fname_orb_1, fname_orb_2, PRN, interpstep, NPint, fname_orbint)
! or
	! Performing interpolation after forming the orbit array (position vector) based on the mutliple sp3 files
	!CALL orb3sp3_2(fname_orb_0, fname_orb_1, fname_orb_2, PRN, interpstep, NPint, fname_orbint)
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Case 5: Orbit: Position from sp3 data; Velocity approximated through sp3 position differences without applying interpolation
! ----------------------------------------------------------------------
Else if (data_opt == TYPE_POSSP3) then
! fname_write = 'orb_sp3.out'
! CALL orb_sp3 (fname_orb, PRN, fname_write)
! ----------------------------------------------------------------------

End if

if (.false.) then
! ----------------------------------------------------------------------
PRINT *,"--------------------- INPUT ----------------------------"
PRINT *, "External orbit case:              ", data_opt
PRINT *, "External orbit file name:         ", TRIM(fname_orb)
PRINT *, "External orbit interval           ", interpstep
PRINT *, "External orbit inerpolation points", NPint
PRINT *, "External orbit output file name   ", TRIM(fname_orbint)
PRINT *,"--------------------------------------------------------"
PRINT *," "
! ----------------------------------------------------------------------
end if

END

max_stars_repo_name: HiTMonitor/ginan

max_issues_repo_name: HiTMonitor/ginan

max_forks_repo_name: HiTMonitor/ginan


text:
C    %W% %G%
      subroutine svsss1

C     THIS SUBROUTINE DECODES AND ERROR CHECKS TYPE W SUPPLEMENTARY
C     SIGNAL DATA CARD.  IT IS CALLED BY INPUT3.

      include 'tspinc/param.inc'
      include 'tspinc/svs.inc'
      include 'tspinc/inp3.inc'
      include 'tspinc/inp3a.inc'
      include 'tspinc/tzro.inc'
      include 'tspinc/prt.inc'
      dimension temp(16)
      character*8 name1, name2
      character*1 id1, id2

C     -     begin     begin     begin     begin     begin     begin
      read (work80(icard), 10000) (temp(i),i = 1, 4),(temp(i),i=7,10)
     &                    ,ierr, name1, bkv1, id1, name2, bkv2, id2
10000 format (bz, 16x, f4.0, 3f4.4, f6.2, 2f4.2, f4.3, i1, 2x, 2(a8,
     & f4.0, a1))

      svserr(ksv) = 0.0
      if (ierr .ne. 0) svserr(ksv) = 1.

      if (temp(4) .le. 0.0) then
        write (errbuf(1), 10040) nbname, bkv, nid
10040   format ('0', 2x, a8, 2x, f5.1, 2x, a1, 5x,
     &   'TS9 MUST BE GREATER THAN ZERO')
        call prterr('E', 1)
        imchn = 1
      endif

      if (temp(7) .le. 0.0) then
        write (errbuf(1), 10050) nbname, bkv, nid
10050   format ('0', 2x, a8, 2x, f5.1, 2x, a1, 5x,
     &   'TS13 MUST BE GREATER THAN ZERO')
        call prterr('E', 1)
        imchn = 1
      endif

      if (temp(8) .le. 0.0) then
        write (errbuf(1), 10060) nbname, bkv, nid
10060   format ('0', 2x, a8, 2x, f5.1, 2x, a1, 5x,
     &   'TS14 MUST BE GREATER THAN ZERO')
        call prterr('E', 1)
        imchn = 1
      endif

      if (temp(9) .eq. 0.0) then
        write (errbuf(1), 10070) nbname, bkv, nid
10070   format ('0', 2x, a8, 2x, f5.1, 2x, a1, 5x,
     &   'KS3 MUST NOT BE ZERO')
        call prterr('E', 1)
        imchn = 1
      endif

      if (temp(10) .le. 0.0) then
        write (errbuf(1), 10080) nbname, bkv, nid
10080   format ('0', 2x, a8, 2x, f5.1, 2x, a1, 5x,
     &   'VSCSMAX MUST BE GREATER THAN ZERO')
        call prterr('E', 1)
        imchn = 1
      endif

      if (imchn .ne. 1) then
        if (subtyp .eq. 'A') then
          isupt(ksv) = 1
          if (name1 .eq. ' ') then
            jbus = 0
          else
            kb = nambas(bkv1)
            jbus = inam(name1, kb)
          endif
          irbus(ksv) = jbus
          irmote(ksv) = jbus
          igenno(ksv) = jbus
          iparsv(ksv) = id1
        elseif (subtyp .eq. 'B') then
          isupt(ksv) = 2
          kb1 = nambas(bkv1)
          kb2 = nambas(bkv2)
          irmote(ksv) = inam(name1, kb1)
          jrmote(ksv) = inam(name2, kb2)
          iparsv(ksv) = id2
        else
          isupt(ksv) = 3
          kb1 = nambas(bkv1)
          irbus(ksv) = inam(name1, kb1)
          irmote(ksv) = irbus(ksv)
        endif

        cks1(ksv) = temp(1)
        cks3(ksv) = temp(9)
        vscsmx(ksv) = temp(10)
        xcon = 2.0*frqbse/dtsvs
        as7(ksv) = xcon*temp(2) + 1.
        as8(ksv) = xcon*temp(3) + 1.
        as9(ksv) = xcon*temp(4) + 1.
        as13(ksv) = xcon*temp(7)
        as14(ksv) = xcon*temp(8) + 1.
        ispcde(ksv) = 1
      endif
      return
      end


text:
c
c -------------------------------------------------------------
c
      subroutine ginit(msave, first, nvar, quad, gradThreshold)
c
      implicit double precision (a-h,o-z)
      logical first, quad
      logical            graf
      parameter  (maxgr = 192, maxlv=12)
      common /nodal/ rnode(12,maxgr),node(17,maxgr),lstart(maxlv),
     *  newstl(maxlv),
     *  listsp(maxlv),intrat(maxlv),tol,bzonex,bzoney,mstart,ndfree,
     *  lfine,iorder,mxnest,kcheck,lwidth,
     *  maxcl, graf, lhead
      include "calloc.i"
c
c ginit = init's all grids at 'level'  by calling to finit
c         if first = true, (first call to init), then allocate a
c         storage area, else use one already allocated.
c
      mptr  = msave
 10   if (mptr .eq. 0) go to 99
          maxi    = node(5,mptr)
          maxj    = node(6, mptr)
          maxip1  = maxi + 1
          maxjp1  = maxj + 1
          mitot   = maxi-1 + 2*lwidth
          mjtot   = maxj-1 + 2*lwidth
          corn1   = rnode(1,mptr)
          corn2   = rnode(2,mptr)
          hx      = rnode(9,mptr)
          hy      = rnode(10,mptr)
          if(.not. (first)) go to 20
              loc             = igetsp(maxip1*maxjp1*nvar)
              node(7,mptr)    = loc
              rnode(12, mptr) = 0.0
c             get space for irr array, and for shifted values for quad.recon.
              !! get space for ncount and numHoods too
              node(14, mptr) = igetsp(3*mitot*mjtot)
              call setirr(alloc(node(14,mptr)),mitot,mjtot,mptr,quad,
     .                    gradThreshold)
              go to 30
 20       continue
c
c  if 2nd time through, put initial values in 8 so finer grids
c  can be advanced with interpolation of their boundary values.
c  new time soln should still be in location 7.
c
          loc     = node(8,mptr)
c
   30     continue
c         # added irr,mitot,mjtot,lstgrd to finit call
c         # is this lstgrd value right???
          lstgrd  = node(17,mptr)
          call finit(alloc(loc),nvar,maxip1,maxjp1,
     1               corn1,corn2,hx,hy,alloc(node(14,mptr)),
     2               mitot,mjtot,lstgrd)
c
          mptr  = node(10, mptr)
      go to 10
c
c
 99   return
      end


text:
! The MIT License (MIT)
!
! Copyright (c) 2016 odlomax
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

! module containing compile-time parameters
module m_kind_parameters

   use iso_c_binding

   implicit none

   ! variable types
   integer,parameter :: rel_kind=c_double               ! c double/float kind
   integer,parameter :: int_kind=c_int                  ! c int kind
   integer,parameter :: cpx_kind=c_double_complex       ! c double complex kind
   integer,parameter :: chr_kind=c_char                 ! c character kind
   integer,parameter :: log_kind=c_bool                 ! c boolean kind
   integer,parameter :: string_length=2048              ! default character length

   integer(kind=int_kind),parameter :: n_dim=3          ! number of dimensions

end module

module m_sph_parameters

   use m_kind_parameters

   implicit none

   ! sph parameters
   integer(kind=int_kind),parameter :: n_leaf=2**n_dim          ! max particles per leaf cell
   integer(kind=int_kind),parameter :: h_iterations=10          ! maximum number of h iterations
   integer(kind=int_kind),parameter :: n_lookup=128             ! number of entries in kernel lookup table
   real(kind=rel_kind),parameter :: delta_h=0.01_rel_kind       ! smoothing length convergence parameter

   ! ray parameters
   integer(kind=int_kind),parameter :: n_ray=128                ! initial size of ray_particle array (can be resized)
   real(kind=rel_kind),parameter :: ray_eta=1.2_rel_kind        ! ray length over estimation factor
   real(kind=rel_kind),parameter :: ray_conv=0.01_rel_kind      ! ray convergence parameter
   real(kind=rel_kind),parameter :: mrw_max_len=3._rel_kind     ! maximum length of MRW step in units of h

end module

module m_constants_parameters

   use m_kind_parameters

   implicit none

   ! numbers
   real(kind=rel_kind),parameter :: pi=4._rel_kind*atan(1._rel_kind)    ! 3.1412...
   real(kind=rel_kind),parameter :: one_over_four_pi=0.25_rel_kind/pi   ! 1/ 4*pi
   real(kind=rel_kind),parameter :: one_third=1._rel_kind/3._rel_kind   ! 1/3
   real(kind=rel_kind),parameter :: root_two=sqrt(2._rel_kind)          ! sqrt of two
   real(kind=rel_kind),parameter :: root_three=sqrt(3._rel_kind)        ! sqrt of three
   real(kind=rel_kind),parameter :: arcsec_rad=pi/6.48e+5_rel_kind      ! 1 arcsec in radians
   real(kind=rel_kind),parameter :: fwhm_sigma=2._rel_kind*sqrt(2._rel_kind*log(2._rel_kind))   ! FWHM in units of sigma

   ! physical constants
   real(kind=rel_kind),parameter :: k_boltz=1.3806488e-23_rel_kind      ! Boltzmann constant [J/K]
   real(kind=rel_kind),parameter :: c_light=2.99792458e+8_rel_kind      ! speed of light [m/s]
   real(kind=rel_kind),parameter :: h_planck=6.6260696e-34_rel_kind     ! Planck constant [J/Hz]
   real(kind=rel_kind),parameter :: sigma_sb=5.6703668e-08_rel_kind     ! Stefan-Boltzmann constant [W m^-2 K^-4]

end module

max_stars_repo_name: odlomax/spamcart-dev

max_issues_repo_name: odlomax/spamcart-dev

max_forks_repo_name: odlomax/spamcart-dev


text:
Describe Users/Tammydockter here.

Hi Tammy, Welcome to the Wiki! Please dont delete comments.  Myra

Tammyplease stop deleting other peoples comments.  Users/PeteB