text:
!************************************************************
!
!  This example shows a way to recursively traverse the file
!  using h5literate and h5literate_by_name_f.  The method shown
!  here guarantees that
!  the recursion will not enter an infinite loop, but does
!  not prevent objects from being visited more than once.
!  The program prints the directory structure of the file
!  specified in filename.  The default file used by this example
!  implements the structure described in the User's Guide,
!  chapter 4, figure 26.
!
!  This file is intended for use with HDF5 Library version 1.8
!  with --enable-fortran2003
! ************************************************************

MODULE g_traverse

  USE HDF5
  USE ISO_C_BINDING
  IMPLICIT NONE

  CHARACTER(LEN=18) :: filename = "h5ex_g_traverse.h5"

  !
  ! Define operator data structure type for H5Literate callback.
  ! During recursive iteration, these structures will form a
  ! linked list that can be searched for duplicate groups,
  ! preventing infinite recursion.
  !
  TYPE :: opdata
     INTEGER :: recurs              ! Recursion level.  0=root
     TYPE(opdata), POINTER :: prev  ! Pointer to previous opdata
     INTEGER(haddr_t) :: addr       ! Group address
  END TYPE opdata

CONTAINS

  !
  ! OPERATOR FUNCTION TO BE CALLED BY H5LITERATE_F
  !
  ! ************************************************************
  !
  !  Operator function.  This function prints the name and type
  !  of the object passed to it.  If the object is a group, it
  !  is first checked against other groups in its path using
  !  the group_check function, then if it is not a duplicate,
  !  H5Literate is called for that group.  This guarantees that
  !  the program will not enter infinite recursion due to a
  !  circular path in the file.
  !
  ! ************************************************************

  RECURSIVE INTEGER FUNCTION op_func(loc_id, name, info, operator_data) RESULT(ret_val) BIND(C)

    USE HDF5
    USE ISO_C_BINDING
    IMPLICIT NONE

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

    INTEGER :: status, return_val
    TYPE(h5o_info_t), TARGET :: infobuf
    TYPE(C_PTR) :: ptr
    CHARACTER(LEN=10) :: name_string
    INTEGER :: i
    TYPE(opdata), POINTER :: od
    TYPE(opdata), TARGET :: nextod
    INTEGER(HSIZE_T) :: idx

    TYPE(C_PTR) :: ptr2
    TYPE(C_FUNPTR) :: funptr

    CHARACTER(LEN=10) :: space
    INTEGER :: spaces ! Number of whitespaces to prepend to output
    INTEGER :: len

    ret_val = 0

    name_string(1:10) = " "
    len = 0
    DO
       len = len + 1
       IF(name(len)(1:1).EQ.C_NULL_CHAR) EXIT
       name_string(len:len) = name(len)(1:1)
    ENDDO
    len = len - 1 ! subtract NULL character

    space(1:10) = " "

    CALL C_F_POINTER(operator_data, od)
    !
    ! Get type of the object and display its name and type.
    ! The name of the object is passed to this function by
    ! the Library.
    !
    CALL H5Oget_info_by_name_f(loc_id, name_string, infobuf, status)

    spaces = 2*(od%recurs+1)

    WRITE(*,'(A)', ADVANCE='NO') space(1:spaces) !  Format output

    IF(infobuf%type.EQ.H5O_TYPE_GROUP_F)THEN

       WRITE(*,'("Group: ",A," {")') name_string(1:len)

!
!              Check group address against linked list of operator
!              data structures.  We will always run the check, as the
!              reference count cannot be relied upon if there are
!              symbolic links, and H5Oget_info_by_name always follows
!              symbolic links.  Alternatively we could use H5Lget_info
!              and never recurse on groups discovered by symbolic
!              links, however it could still fail if an object's
!              reference count was manually manipulated with
!              H5Odecr_refcount.
!

       i = group_check(od, infobuf%addr)

       IF(i.EQ.1)THEN
          WRITE(*,'(A)') space(1:spaces)//"  Warning: Loop detected!"
       ELSE

          nextod%recurs = od%recurs + 1
          nextod%prev => od
          nextod%addr = infobuf%addr
          idx = 0
          ptr2 = C_LOC(nextod)
          funptr = C_FUNLOC(op_func)
          CALL h5literate_by_name_f(loc_id, name_string, H5_INDEX_NAME_F, H5_ITER_NATIVE_F, idx, &
               funptr, ptr2, ret_val, status)

       ENDIF
       WRITE(*,'(A)') space(1:spaces)//"}"
       RETURN
    ELSE IF(infobuf%type.EQ.H5O_TYPE_DATASET_F)THEN
       WRITE(*,'("Dataset: ",A)') name_string(1:len)
    ELSE IF(infobuf%type.EQ.H5O_TYPE_NAMED_DATATYPE_F)THEN
       WRITE(*,'("Datatype: ",A)') name_string(1:len)
    ELSE
       WRITE(*,'("Unknown: ",A)') name_string(1:len)
    ENDIF

END FUNCTION op_func

!************************************************************
!
!  This function recursively searches the linked list of
!  opdata structures for one whose address matches
!  target_addr.  Returns 1 if a match is found, and 0
!  otherwise.
!
! ************************************************************/

  INTEGER RECURSIVE FUNCTION group_check(od, target_addr) result(g_c)

    IMPLICIT NONE
    TYPE(opdata), POINTER :: od
    INTEGER(haddr_t) :: target_addr

    IF (od%addr .EQ. target_addr)THEN
       g_c = 1       ! Addresses match
    ELSE IF (od%recurs.EQ.0)THEN
       g_c = 0       ! Root group reached with no matches
    ELSE
       ! Recursively examine the next node
       g_c = group_check(od%prev, target_addr)
    END IF
  END FUNCTION group_check

END MODULE g_traverse

PROGRAM main

  USE HDF5
  USE ISO_C_BINDING

  USE g_traverse

  IMPLICIT NONE

  INTEGER(hid_t) :: file ! Handle
  INTEGER :: status
  TYPE(h5o_info_t) :: infobuf
  TYPE(opdata), TARGET :: od
  TYPE(C_PTR) :: ptr
  INTEGER(hsize_t) :: idx
  INTEGER :: ret_value
  TYPE(C_FUNPTR) :: funptr
  !
  ! Open file and initialize the operator data structure.
  !
  CALL H5Fopen_f(filename, H5F_ACC_RDONLY_F, file, status)

  CALL h5oget_info_by_name_f(file, "/", infobuf, status)

  od%recurs = 0
  od%prev => NULL()
  od%addr = infobuf%addr
  !
  ! Print the root group and formatting, begin iteration.
  !
  idx = 0
  funptr = C_FUNLOC(op_func)
  ptr    = C_LOC(od)

  WRITE(*,'(A)') "/ {"
  CALL  H5Literate_f(file, H5_INDEX_NAME_F,  H5_ITER_NATIVE_F, idx, funptr, ptr, ret_value, status)
  WRITE(*,'(A)') "}"

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

END PROGRAM main

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
*> \brief \b SGBBRD
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SGBBRD + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE SGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
*                          LDQ, PT, LDPT, C, LDC, WORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          VECT
*       INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
*       ..
*       .. Array Arguments ..
*       REAL               AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
*      $                   PT( LDPT, * ), Q( LDQ, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SGBBRD reduces a real general m-by-n band matrix A to upper
*> bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
*>
*> The routine computes B, and optionally forms Q or P**T, or computes
*> Q**T*C for a given matrix C.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] VECT
*> \verbatim
*>          VECT is CHARACTER*1
*>          Specifies whether or not the matrices Q and P**T are to be
*>          formed.
*>          = 'N': do not form Q or P**T;
*>          = 'Q': form Q only;
*>          = 'P': form P**T only;
*>          = 'B': form both.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] NCC
*> \verbatim
*>          NCC is INTEGER
*>          The number of columns of the matrix C.  NCC >= 0.
*> \endverbatim
*>
*> \param[in] KL
*> \verbatim
*>          KL is INTEGER
*>          The number of subdiagonals of the matrix A. KL >= 0.
*> \endverbatim
*>
*> \param[in] KU
*> \verbatim
*>          KU is INTEGER
*>          The number of superdiagonals of the matrix A. KU >= 0.
*> \endverbatim
*>
*> \param[in,out] AB
*> \verbatim
*>          AB is REAL array, dimension (LDAB,N)
*>          On entry, the m-by-n band matrix A, stored in rows 1 to
*>          KL+KU+1. The j-th column of A is stored in the j-th column of
*>          the array AB as follows:
*>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
*>          On exit, A is overwritten by values generated during the
*>          reduction.
*> \endverbatim
*>
*> \param[in] LDAB
*> \verbatim
*>          LDAB is INTEGER
*>          The leading dimension of the array A. LDAB >= KL+KU+1.
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is REAL array, dimension (min(M,N))
*>          The diagonal elements of the bidiagonal matrix B.
*> \endverbatim
*>
*> \param[out] E
*> \verbatim
*>          E is REAL array, dimension (min(M,N)-1)
*>          The superdiagonal elements of the bidiagonal matrix B.
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is REAL array, dimension (LDQ,M)
*>          If VECT = 'Q' or 'B', the m-by-m orthogonal matrix Q.
*>          If VECT = 'N' or 'P', the array Q is not referenced.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.
*>          LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
*> \endverbatim
*>
*> \param[out] PT
*> \verbatim
*>          PT is REAL array, dimension (LDPT,N)
*>          If VECT = 'P' or 'B', the n-by-n orthogonal matrix P'.
*>          If VECT = 'N' or 'Q', the array PT is not referenced.
*> \endverbatim
*>
*> \param[in] LDPT
*> \verbatim
*>          LDPT is INTEGER
*>          The leading dimension of the array PT.
*>          LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          C is REAL array, dimension (LDC,NCC)
*>          On entry, an m-by-ncc matrix C.
*>          On exit, C is overwritten by Q**T*C.
*>          C is not referenced if NCC = 0.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of the array C.
*>          LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (2*max(M,N))
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup realGBcomputational
*
*  =====================================================================
      SUBROUTINE SGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
     $                   LDQ, PT, LDPT, C, LDC, WORK, INFO )
*
*  -- LAPACK computational routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          VECT
      INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
*     ..
*     .. Array Arguments ..
      REAL               AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
     $                   PT( LDPT, * ), Q( LDQ, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            WANTB, WANTC, WANTPT, WANTQ
      INTEGER            I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
     $                   KUN, L, MINMN, ML, ML0, MN, MU, MU0, NR, NRT
      REAL               RA, RB, RC, RS
*     ..
*     .. External Subroutines ..
      EXTERNAL           SLARGV, SLARTG, SLARTV, SLASET, SROT, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      WANTB = LSAME( VECT, 'B' )
      WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB
      WANTPT = LSAME( VECT, 'P' ) .OR. WANTB
      WANTC = NCC.GT.0
      KLU1 = KL + KU + 1
      INFO = 0
      IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) )
     $     THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NCC.LT.0 ) THEN
         INFO = -4
      ELSE IF( KL.LT.0 ) THEN
         INFO = -5
      ELSE IF( KU.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDAB.LT.KLU1 ) THEN
         INFO = -8
      ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
         INFO = -12
      ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
         INFO = -14
      ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
         INFO = -16
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGBBRD', -INFO )
         RETURN
      END IF
*
*     Initialize Q and P**T to the unit matrix, if needed
*
      IF( WANTQ )
     $   CALL SLASET( 'Full', M, M, ZERO, ONE, Q, LDQ )
      IF( WANTPT )
     $   CALL SLASET( 'Full', N, N, ZERO, ONE, PT, LDPT )
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
*
      MINMN = MIN( M, N )
*
      IF( KL+KU.GT.1 ) THEN
*
*        Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
*        first to lower bidiagonal form and then transform to upper
*        bidiagonal
*
         IF( KU.GT.0 ) THEN
            ML0 = 1
            MU0 = 2
         ELSE
            ML0 = 2
            MU0 = 1
         END IF
*
*        Wherever possible, plane rotations are generated and applied in
*        vector operations of length NR over the index set J1:J2:KLU1.
*
*        The sines of the plane rotations are stored in WORK(1:max(m,n))
*        and the cosines in WORK(max(m,n)+1:2*max(m,n)).
*
         MN = MAX( M, N )
         KLM = MIN( M-1, KL )
         KUN = MIN( N-1, KU )
         KB = KLM + KUN
         KB1 = KB + 1
         INCA = KB1*LDAB
         NR = 0
         J1 = KLM + 2
         J2 = 1 - KUN
*
         DO 90 I = 1, MINMN
*
*           Reduce i-th column and i-th row of matrix to bidiagonal form
*
            ML = KLM + 1
            MU = KUN + 1
            DO 80 KK = 1, KB
               J1 = J1 + KB
               J2 = J2 + KB
*
*              generate plane rotations to annihilate nonzero elements
*              which have been created below the band
*
               IF( NR.GT.0 )
     $            CALL SLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA,
     $                         WORK( J1 ), KB1, WORK( MN+J1 ), KB1 )
*
*              apply plane rotations from the left
*
               DO 10 L = 1, KB
                  IF( J2-KLM+L-1.GT.N ) THEN
                     NRT = NR - 1
                  ELSE
                     NRT = NR
                  END IF
                  IF( NRT.GT.0 )
     $               CALL SLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
     $                            AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
     $                            WORK( MN+J1 ), WORK( J1 ), KB1 )
   10          CONTINUE
*
               IF( ML.GT.ML0 ) THEN
                  IF( ML.LE.M-I+1 ) THEN
*
*                    generate plane rotation to annihilate a(i+ml-1,i)
*                    within the band, and apply rotation from the left
*
                     CALL SLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ),
     $                            WORK( MN+I+ML-1 ), WORK( I+ML-1 ),
     $                            RA )
                     AB( KU+ML-1, I ) = RA
                     IF( I.LT.N )
     $                  CALL SROT( MIN( KU+ML-2, N-I ),
     $                             AB( KU+ML-2, I+1 ), LDAB-1,
     $                             AB( KU+ML-1, I+1 ), LDAB-1,
     $                             WORK( MN+I+ML-1 ), WORK( I+ML-1 ) )
                  END IF
                  NR = NR + 1
                  J1 = J1 - KB1
               END IF
*
               IF( WANTQ ) THEN
*
*                 accumulate product of plane rotations in Q
*
                  DO 20 J = J1, J2, KB1
                     CALL SROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
     $                          WORK( MN+J ), WORK( J ) )
   20             CONTINUE
               END IF
*
               IF( WANTC ) THEN
*
*                 apply plane rotations to C
*
                  DO 30 J = J1, J2, KB1
                     CALL SROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
     $                          WORK( MN+J ), WORK( J ) )
   30             CONTINUE
               END IF
*
               IF( J2+KUN.GT.N ) THEN
*
*                 adjust J2 to keep within the bounds of the matrix
*
                  NR = NR - 1
                  J2 = J2 - KB1
               END IF
*
               DO 40 J = J1, J2, KB1
*
*                 create nonzero element a(j-1,j+ku) above the band
*                 and store it in WORK(n+1:2*n)
*
                  WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
                  AB( 1, J+KUN ) = WORK( MN+J )*AB( 1, J+KUN )
   40          CONTINUE
*
*              generate plane rotations to annihilate nonzero elements
*              which have been generated above the band
*
               IF( NR.GT.0 )
     $            CALL SLARGV( NR, AB( 1, J1+KUN-1 ), INCA,
     $                         WORK( J1+KUN ), KB1, WORK( MN+J1+KUN ),
     $                         KB1 )
*
*              apply plane rotations from the right
*
               DO 50 L = 1, KB
                  IF( J2+L-1.GT.M ) THEN
                     NRT = NR - 1
                  ELSE
                     NRT = NR
                  END IF
                  IF( NRT.GT.0 )
     $               CALL SLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA,
     $                            AB( L, J1+KUN ), INCA,
     $                            WORK( MN+J1+KUN ), WORK( J1+KUN ),
     $                            KB1 )
   50          CONTINUE
*
               IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
                  IF( MU.LE.N-I+1 ) THEN
*
*                    generate plane rotation to annihilate a(i,i+mu-1)
*                    within the band, and apply rotation from the right
*
                     CALL SLARTG( AB( KU-MU+3, I+MU-2 ),
     $                            AB( KU-MU+2, I+MU-1 ),
     $                            WORK( MN+I+MU-1 ), WORK( I+MU-1 ),
     $                            RA )
                     AB( KU-MU+3, I+MU-2 ) = RA
                     CALL SROT( MIN( KL+MU-2, M-I ),
     $                          AB( KU-MU+4, I+MU-2 ), 1,
     $                          AB( KU-MU+3, I+MU-1 ), 1,
     $                          WORK( MN+I+MU-1 ), WORK( I+MU-1 ) )
                  END IF
                  NR = NR + 1
                  J1 = J1 - KB1
               END IF
*
               IF( WANTPT ) THEN
*
*                 accumulate product of plane rotations in P**T
*
                  DO 60 J = J1, J2, KB1
                     CALL SROT( N, PT( J+KUN-1, 1 ), LDPT,
     $                          PT( J+KUN, 1 ), LDPT, WORK( MN+J+KUN ),
     $                          WORK( J+KUN ) )
   60             CONTINUE
               END IF
*
               IF( J2+KB.GT.M ) THEN
*
*                 adjust J2 to keep within the bounds of the matrix
*
                  NR = NR - 1
                  J2 = J2 - KB1
               END IF
*
               DO 70 J = J1, J2, KB1
*
*                 create nonzero element a(j+kl+ku,j+ku-1) below the
*                 band and store it in WORK(1:n)
*
                  WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
                  AB( KLU1, J+KUN ) = WORK( MN+J+KUN )*AB( KLU1, J+KUN )
   70          CONTINUE
*
               IF( ML.GT.ML0 ) THEN
                  ML = ML - 1
               ELSE
                  MU = MU - 1
               END IF
   80       CONTINUE
   90    CONTINUE
      END IF
*
      IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
*
*        A has been reduced to lower bidiagonal form
*
*        Transform lower bidiagonal form to upper bidiagonal by applying
*        plane rotations from the left, storing diagonal elements in D
*        and off-diagonal elements in E
*
         DO 100 I = 1, MIN( M-1, N )
            CALL SLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA )
            D( I ) = RA
            IF( I.LT.N ) THEN
               E( I ) = RS*AB( 1, I+1 )
               AB( 1, I+1 ) = RC*AB( 1, I+1 )
            END IF
            IF( WANTQ )
     $         CALL SROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC, RS )
            IF( WANTC )
     $         CALL SROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
     $                    RS )
  100    CONTINUE
         IF( M.LE.N )
     $      D( M ) = AB( 1, M )
      ELSE IF( KU.GT.0 ) THEN
*
*        A has been reduced to upper bidiagonal form
*
         IF( M.LT.N ) THEN
*
*           Annihilate a(m,m+1) by applying plane rotations from the
*           right, storing diagonal elements in D and off-diagonal
*           elements in E
*
            RB = AB( KU, M+1 )
            DO 110 I = M, 1, -1
               CALL SLARTG( AB( KU+1, I ), RB, RC, RS, RA )
               D( I ) = RA
               IF( I.GT.1 ) THEN
                  RB = -RS*AB( KU, I )
                  E( I-1 ) = RC*AB( KU, I )
               END IF
               IF( WANTPT )
     $            CALL SROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
     $                       RC, RS )
  110       CONTINUE
         ELSE
*
*           Copy off-diagonal elements to E and diagonal elements to D
*
            DO 120 I = 1, MINMN - 1
               E( I ) = AB( KU, I+1 )
  120       CONTINUE
            DO 130 I = 1, MINMN
               D( I ) = AB( KU+1, I )
  130       CONTINUE
         END IF
      ELSE
*
*        A is diagonal. Set elements of E to zero and copy diagonal
*        elements to D.
*
         DO 140 I = 1, MINMN - 1
            E( I ) = ZERO
  140    CONTINUE
         DO 150 I = 1, MINMN
            D( I ) = AB( 1, I )
  150    CONTINUE
      END IF
      RETURN
*
*     End of SGBBRD
*
      END

ext: f

max_issues_repo_name: WatsonZhouAnda/Cheetah-Software

max_forks_repo_name: WatsonZhouAnda/Cheetah-Software


!##############################################################################
! PROGRAM creditriskmanagement_MC
!
! ## Credit risk management with Monte Carlo Simulations
!
! This code is published under the GNU General Public License v3
!                         (https://www.gnu.org/licenses/gpl-3.0.en.html)
!
! Authors: Hans Fehr and Fabian Kindermann
!          contact@ce-fortran.com
!
! #VC# VERSION: 1.0  (23 January 2018)
!
!##############################################################################
program creditriskmanagement_MC

    use toolbox

    implicit none
    integer, parameter :: KK = 8               ! number of ratings
    integer, parameter :: JJ = 5               ! number of industries
    integer, parameter :: BB = 16              ! number of bonds in portfolio
    integer, parameter :: NN = 100000          ! number of MC simulations
    real*8, parameter  :: sig_RR = 0.25d0      ! volatility of recovery rate
    real*8, parameter  :: r_f = 0.03d0         ! riskfree rate
    logical, parameter :: correlated = .false. ! should shocks be correlated

    real*8 :: RM_data(KK-1, KK), cs_data(KK-1), RR_data(JJ)
    real*8 :: pi_D, NE, mu_RR, cs, cr, pi_k(KK)
    real*8 :: EL_d(BB+1), UL_d(BB+1), EL_m(BB+1), UL_m(BB+1)
    real*8 :: VaR_d, CVaR_d, VaR_m, CVaR_m
    real*8 :: p_0, p_prime, sum_NE
    real*8 :: u(NN), L_d(BB+1, NN), L_m(BB+1, NN), L_temp(NN)
    real*8 :: psi, alpha, beta
    integer :: ij, ik, ik_p, it, TT, ib, ii, ii_VaR

    ! define a bond type
    type :: bond_type
        character(len=1) :: industry
        character(len=3) :: rating
        real*8 :: couponrate
        integer :: maturity
        real*8 :: exposure
    end type

    ! declare the bond variable
    type(bond_type) :: portfolio(BB)

    ! initialize bond data
    call init_data

    ! open file for output
    open(12, file = 'output.out')
    write(12,'(a)')'RISK ANALYSIS OF CORPORATE BOND PORTFOLIO'
    write(12, *)
    write(12,'(a)')'I  R        cr   T         p_0          NE        EL_D        UL_D        EL_M        UL_M'

    ! set up the portfolio
    portfolio(1)  = bond_type('F', 'BBB', 5.00d0,  4, 1000000d0)
    portfolio(2)  = bond_type('M', 'A'  , 4.63d0,  3, 3000000d0)
    portfolio(3)  = bond_type('M', 'BBB', 3.13d0,  6, 1000000d0)
    portfolio(4)  = bond_type('S', 'BBB', 2.88d0, 10, 1000000d0)
    portfolio(5)  = bond_type('T', 'BB' , 6.50d0,  2, 2000000d0)
    portfolio(6)  = bond_type('M', 'A'  , 5.00d0,  6, 5000000d0)
    portfolio(7)  = bond_type('M', 'A'  , 2.63d0,  5, 4000000d0)
    portfolio(8)  = bond_type('M', 'A'  , 1.88d0,  7,  500000d0)
    portfolio(9)  = bond_type('F', 'AA' , 6.50d0, 11, 8000000d0)
    portfolio(10) = bond_type('U', 'A'  , 6.38d0,  3, 4500000d0)
    portfolio(11) = bond_type('F', 'A'  , 4.00d0,  5, 1000000d0)
    portfolio(12) = bond_type('M', 'A'  , 2.00d0,  9,  600000d0)
    portfolio(13) = bond_type('M', 'A'  , 4.25d0,  5,  400000d0)
    portfolio(14) = bond_type('U', 'BB' , 5.13d0,  4, 2000000d0)
    portfolio(15) = bond_type('M', 'B'  , 4.38d0,  3, 1000000d0)
    portfolio(16) = bond_type('M', 'AA' , 5.63d0,  4,10000000d0)

    sum_NE = 0d0

    ! draw one common series of shocks when correlated
    if(correlated)then
        call simulate_uniform(u(1:NN/2), 0d0, 1d0)
        u(NN/2+1:NN) = 1d0 - u(1:NN/2)
    endif

    ! iterate over all bonds in portfolio
    do ib = 1, BB

        ! draw new series of random shocks when uncorrelated shocks
        if(.not. correlated)then
            call simulate_uniform(u(1:NN/2), 0d0, 1d0)
            u(NN/2+1:NN) = 1d0 - u(1:NN/2)
        endif

        ! retrieve bond's properties
        cr    = portfolio(ib)%couponrate/100d0
        TT    = portfolio(ib)%maturity
        NE    = portfolio(ib)%exposure
        ij    = ij_code(portfolio(ib)%industry)
        ik    = ik_code(portfolio(ib)%rating)
        pi_D  = RM_data(ik, KK)
        pi_k  = RM_data(ik, :)
        cs    = cs_data(ik)
        mu_RR = RR_data(ij)

        ! compute price under current rating
        p_0 = 0d0
        do it = 1, TT-1
            p_0 = p_0 + cr/(1d0 + r_f + cs)**it
        enddo
        p_0 = p_0 + (1d0 + cr)/(1d0 + r_f + cs)**TT

        ! define beta distribution parameters and sum up NE
        alpha = mu_RR**2*(1.0d0 - mu_RR)/sig_RR**2 - mu_RR
        beta = alpha/mu_RR - alpha
        sum_NE = sum_NE + NE

        ! start MC simulation
        do ii = 1, NN

            ! check whether bond is in default
            if(u(ii) <= 1d0 - pi_D)then

                ! where does the bond migrate
                do ik_p = 1, KK-1
                    if(u(ii) <= sum(pi_k(1:ik_p)))then
                        cs = cs_data(ik_p)
                        exit
                    endif
                enddo

                ! calculate price change if credit rating changes
                p_prime = 0d0
                do it = 1, TT-1
                    p_prime = p_prime + cr/(1d0 + r_f + cs)**it
                enddo
                p_prime = p_prime + (1d0 + cr)/(1d0 + r_f + cs)**TT

                L_d(ib, ii) = 0d0
                L_m(ib, ii) = NE*(p_0 - p_prime)
            else

                ! simulate recovery rate under default
                call simulate_beta(psi, alpha, beta)
                L_d(ib, ii) = NE*(p_0 - psi)
                L_m(ib, ii) = NE*(p_0 - psi)
            endif
        enddo
    enddo

    ! get loss of total portfolio
    L_d(BB+1, :) = sum(L_d(1:BB, :), 1)
    L_m(BB+1, :) = sum(L_m(1:BB, :), 1)

    ! plot migration mode loss histogram
    call plot_hist(L_m(BB+1, :)/sum_NE*10000d0, 99, -200d0, 600d0)
    call execplot(xlim=(/-200d0, 600d0/), xlabel='Portfolio loss in bp', ylabel='Relative frequency')

    ! plot tail distribution
    call plot_hist(L_m(BB+1, :)/sum_NE*10000d0, 99, 200d0, 600d0)
    call execplot(xlim=(/200d0, 600d0/), xlabel='Portfolio loss in bp', ylabel='Relative frequency')

    open(928,file='fig49_data.m')
    write(928, '(a)')'L_m = ['
    do ii = 1, NN
        write(928,'(f40.10)')L_m(BB+1, ii)
    enddo
    write(928, '(a)')'];'
    close(928)

    ! calculate EL and UL
    do ib = 1, BB+1
        EL_d(ib) = sum(L_d(ib, :))/dble(NN)
        UL_d(ib) = sqrt( sum((L_d(ib, :)-EL_d(ib))**2)/dble(NN-1) )
        EL_m(ib) = sum(L_m(ib, :))/dble(NN)
        UL_m(ib) = sqrt( sum((L_m(ib, :)-EL_m(ib))**2)/dble(NN-1) )
    enddo

    ! calculate value-at-risk
    call tic

    ii_VaR = ceiling(0.9d0*NN)

    ! sort portfolio losses (default mode)
    L_temp = L_d(BB+1, :)
    call sort(L_temp)
    VaR_d  = L_temp(ii_VaR)
    CVaR_d = sum(L_temp(ii_VaR:NN))/dble(NN-ii_VaR+1)

    ! sort portfolio losses (migration mode)
    L_temp = L_m(BB+1, :)
    call sort(L_temp)
    VaR_m  = L_temp(ii_VaR)
    CVaR_m = sum(L_temp(ii_VaR:NN))/dble(NN-ii_VaR+1)

    call toc

    ! print results
    do ib = 1, BB

        cr    = portfolio(ib)%couponrate/100d0
        TT    = portfolio(ib)%maturity
        NE    = portfolio(ib)%exposure
        ik    = ik_code(portfolio(ib)%rating)
        cs    = cs_data(ik)

        p_0 = 0d0
        do it = 1, TT-1
            p_0 = p_0 + cr/(1d0 + r_f + cs)**it
        enddo
        p_0 = p_0 + (1d0 + cr)/(1d0 + r_f + cs)**TT

        write(12, '(a1,a5,f8.2,i4,f12.6,f12.0,4f12.2)')portfolio(ib)%industry, portfolio(ib)%rating, &
                        cr*100d0, TT, p_0, NE, EL_d(ib), UL_d(ib), EL_m(ib), UL_m(ib)
    enddo
    write(12,'(30x,f12.0,4f12.2)')sum_NE, EL_d(BB+1), UL_d(BB+1), EL_m(BB+1), UL_m(BB+1)
    write(12,'(44x,4(f10.2,a))')EL_d(BB+1)/sum_NE*10000d0, 'bp', UL_d(BB+1)/sum_NE*10000d0, 'bp', &
                                EL_m(BB+1)/sum_NE*10000d0, 'bp', UL_m(BB+1)/sum_NE*10000d0, 'bp'
    close(12)

    write(*,'(a)')'                          EL          UL         VaR        CVaR'
    write(*,'(a,4f12.2)')'Default mode  : ', EL_d(BB+1), UL_d(BB+1), VaR_d, CVaR_d
    write(*,'(a,4f12.2)')'        in bp : ', (/EL_d(BB+1), UL_d(BB+1), VaR_d, CVaR_d/)/sum_NE*10000d0
    write(*,*)
    write(*,'(a,4f12.2)')'Migration mode: ', EL_m(BB+1), UL_m(BB+1), VaR_m, CVaR_m
    write(*,'(a,4f12.2)')'        in bp : ', (/EL_m(BB+1), UL_m(BB+1), VaR_m, CVaR_m/)/sum_NE*10000d0

contains

    ! initialize the probability parameters
    subroutine init_data()

        implicit none

        RM_data(1, :) = (/ 90.24d0,  8.99d0,  0.56d0,  0.05d0,  0.08d0,  0.03d0,  0.05d0,  0.00d0 /)
        RM_data(2, :) = (/  0.58d0, 90.03d0,  8.65d0,  0.56d0,  0.06d0,  0.08d0,  0.02d0,  0.02d0 /)
        RM_data(3, :) = (/  0.04d0,  2.00d0, 91.58d0,  5.71d0,  0.40d0,  0.17d0,  0.02d0,  0.08d0 /)
        RM_data(4, :) = (/  0.01d0,  0.13d0,  3.89d0, 90.69d0,  4.18d0,  0.68d0,  0.16d0,  0.26d0 /)
        RM_data(5, :) = (/  0.02d0,  0.04d0,  0.18d0,  5.81d0, 84.15d0,  7.97d0,  0.83d0,  1.00d0 /)
        RM_data(6, :) = (/  0.00d0,  0.05d0,  0.15d0,  0.25d0,  6.31d0, 83.16d0,  5.01d0,  5.07d0 /)
        RM_data(7, :) = (/  0.00d0,  0.00d0,  0.20d0,  0.30d0,  0.91d0, 15.96d0, 51.31d0, 31.32d0 /)
        RM_data = RM_data/100d0

        CS_data = (/0.00d0, 0.05d0, 1.00d0, 2.00d0, 5.00d0, 10.00d0, 20.00d0/)/100d0

        RR_data = (/0.437d0, 0.246d0, 0.575d0, 0.327d0, 0.493d0/)

    end subroutine init_data

    ! return rating code
    function ik_code(rating)

        implicit none
        character(len=3), intent(in) :: rating
        integer :: ik_code

        select case (rating)
            case('AAA')
                ik_code = 1
            case('AA')
                ik_code = 2
            case('A')
                ik_code = 3
            case('BBB')
                ik_code = 4
            case('BB')
                ik_code = 5
            case('B')
                ik_code = 6
            case('CCC')
                ik_code = 7
            case default
                stop 'Credit rating does not exists'
        end select

    end function

    ! return industry code
    function ij_code(industry)

        implicit none
        character(len=1), intent(in) :: industry
        integer :: ij_code

        select case (industry)
            case('M')
                ij_code = 1
            case('F')
                ij_code = 2
            case('U')
                ij_code = 3
            case('T')
                ij_code = 4
            case('S')
                ij_code = 5
            case default
                stop 'Industry does not exists'
        end select

    end function

end program


        !COMPILER-GENERATED INTERFACE MODULE: Wed Jul 21 10:02:25 2021
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE WARNINGS__genmod
          INTERFACE
            SUBROUTINE WARNINGS
            END SUBROUTINE WARNINGS
          END INTERFACE
        END MODULE WARNINGS__genmod


text:
      subroutine geneq2_wtsi(i)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C**********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 To generate equations isothermal air-water solution at each node.
CD1
C**********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Initial implementation: 22-May-02, Programmer: George Zyvoloski
CD2
CD2 $Log:   /pvcs.config/fehm90/src/geneq2_wtsi.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:06   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:06:00   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
CD2
C***********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3  2.3.2 Heat- and mass-transfer equations
CD3  2.3.3 Noncondensible gas flow equations
CD3
C***********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C**********************************************************************
c
c generate equations for 3-d schemes'isothermal air-water ,
c full derivatives'
c
      use comflow
      use davidi
      use comji
      use comfi
      use comgi
      use comei
      use comdi
      use comci
      use combi
      use comii
      use comdti
      use comai
      use comwt
      use comriv
      implicit none

      integer i,imm
      integer icd
      integer ii1
      integer ii2
      integer idg
      integer iq
      integer jmi
      integer jml
      integer jmia
      integer jm
      integer neqp1
      integer ij
      integer ij1
      integer ij2
      integer iz
      integer kb,iunsat,kb_unsat
      integer neighc
      integer iau
      integer ial
      integer kz
      integer nmatavw,dry_nodes
      real*8 sx1d
      real*8 axi
      real*8 ayi
      real*8 azi
      real*8 alxi
      real*8 alyi
      real*8 alzi
      real*8 avxi
      real*8 avyi
      real*8 avzi
      real*8 pvii
      real*8 phii
      real*8 swi
      real*8 dlpi
      real*8 dlpkb
      real*8 dvpi
      real*8 dvpkb
      real*8 dlei,df
      real*8 dvei
      real*8 dlekb
      real*8 dvekb
      real*8 dlapi
      real*8 dvapi
      real*8 dlapkb
      real*8 dvapkb
      real*8 dili
      real*8 dilkb
      real*8 divi,dz
      real*8 divkb
      real*8 axkb
      real*8 aykb
      real*8 azkb
      real*8 alxkb
      real*8 alykb
      real*8 alzkb
      real*8 sx2c
      real*8 sx4d
c     real*8 sxzc
      real*8 pvikb
      real*8 phikb
      real*8 radi
      real*8 radkb
      real*8 fid
      real*8 fid1
      real*8 axyd,axyds
      real*8 axy
      real*8 axyf
      real*8 heatc
      real*8 pvxy
      real*8 pxy
      real*8 pxyh
      real*8 pxyi
      real*8 sx4h
      real*8 vxyd
      real*8 vxy
      real*8 vxyf
      real*8 dpvti
      real*8 dilpi
      real*8 dilei
      real*8 divpi
      real*8 divei
      real*8 dilpkb
      real*8 divpkb
      real*8 divekb
      real*8 dilekb
      real*8 sx2t
      real*8 sx3t
      real*8 sxzt
      real*8 dlaei
      real*8 dlaekb
      real*8 dvaei
      real*8 dvaekb
      real*8 ti
      real*8 dis2,dis_tol,sx_min,dxy2
      real*8 delx2
      real*8 dely2
      real*8 delz2
      real*8 reduction_factor
      real*8 grav_air,s_wt_fac,ds_fac_wti,ds_fac_wtkb
      real*8 rlpfree,dfid,dfid1,phi_0,cap_fac
      real*8 rlxyd, drlxyi, drlxykb, dfi, dfkb, tiny
      parameter(dis_tol=1.d-12, cap_fac = 0.0d00, tiny = 1.d-20)

      logical bit, test_bit
c     integer isl,isw_term
      integer isl
      integer iz4m1
      integer imd,iwd

c
c additional terms
c
      real*8 tcap1(nn),tcap2(nn),dtcap1(nn),dtcap2(nn)
      real*8 pcai, axycap, pcaij,pcakb
c
      phi_0 = crl(4,1) + phi_inc
c     water column height correction
      pcai = cap_fac*dzrg(i)*(rlxyf(i)-0.5)

c changed by avw -- entered here by seh
      neqp1=neq+1
      if(i.gt.neq) then
         nmatavw=ldna
      else
         nmatavw=0
      endif
      if(icons.le.abs(maxit)) then
       grav_air=0.0
      else
       grav_air=grav
      endif

c
c storage for upwind
c

      sx1d=sx1(i)
      axi=pnx(i)
      ayi=pny(i)
      azi=pnz(i)
      alxi=axi
      avxi=axi
      alyi=ayi
      avyi=ayi
      alzi=azi
      avzi=azi
      pvii=phi(i)
      dili=dil(i)
      dilpi=dilp(i)
      if(irdof.ne.13) then
        phii=pvii-pcp(i)
        dpvti=dpcef(i)
        divi=div(i)
        divpi=divp(i)
        divei=dive(i)
        dilei=dile(i)
      else
        phii=pvii
      endif
      ti=t(i)
      swi=s(i)
c
c form constants for i>neq
c
      if(i.gt.neq.and.idualp.eq.0) then
         icd=neq
      else
         icd=0
      endif
      iz=i-icd
c
      iz4m1 = 4*(iz-1)+1
c
      ii1=nelm(i-icd)+1
      ii2=nelm(i-icd+1)
      idg=nelmdg(i-icd)-ii1+1
      neqp1=neq+1
      iq=0
      jmi=nelmdg(i-icd)
      jmia=jmi-neqp1

c Take care of source/sink term by putting it into empty slot (from
c same node to same node).
c If this is an isothermal air-water simulation
      a_axy(jmia+nmatavw)=sk(i)
c Take care of source/sink term
c If this is an isothermal air-water simulation
      if (irdof .ne. 13) then
         a_vxy(jmia+nmatavw)=qh(i)
      end if

      do 58 jm=jmi+1,ii2
         kb=nelm(jm)+icd
         if(iriver.eq.2.and.kb.gt.neq_primary) go to 58
         iq=iq+1
         it8(iq)=kb
         it9(iq)=jm-ii1+1
         it10(iq)=istrw(jm-neqp1)
c gaz 11-18-2001
c         if(imdnode.ne.0) then
c           imd = mdnodes(i) + mdnodes(kb)
c           if(imd.lt.2) it10(iq) = -abs(it10(iq))
c         endif
         it11(iq)=jm-neqp1
         ij1=nelm(nelm(jm))+1
         ij2=nelmdg(nelm(jm))-1
         do 68 ij=ij1,ij2
            if(nelm(ij).eq.iz) then
               it12(iq)=ij-neqp1
            endif
 68      continue
 58   continue
c
c 3-d geometry
c
      if(icnl.eq.0) then
         do 59 jm=1,iq
            kb=it8(jm)
            kz=kb-icd
            neighc=it9(jm)
            iwd=it10(jm)
            iw =abs(iwd)
            axkb=pnx(kb)
            aykb=pny(kb)
            azkb=pnz(kb)
            alxkb=axkb
            alykb=aykb
            alzkb=azkb
            reduction_factor = red_factor(istrw_itfc(it11(jm)))
            perml(1)=2.*alxkb*alxi/(alxkb+alxi)
            perml(2)=2.*alykb*alyi/(alykb+alyi)
            perml(3)=2.*alzkb*alzi/(alzkb+alzi)
            sx2c=sx(iw,isox)+sx(iw,isoy)+sx(iw,isoz)

            pvikb=phi(kb)
            if (irdof .ne. 13) then
               phikb=pvikb-pcp(kb)
            else
               phikb=pvikb
            end if
            delx2=(cord(kz,1)-cord(iz,1))**2
            dely2=(cord(kz,2)-cord(iz,2))**2
            delz2=(cord(kz,3)-cord(iz,3))**2
            dis2=delx2+dely2+delz2
            if(dis2.gt.dis_tol.and.iwd.gt.0) then
               pxy=sx2c*dis2/(delx2/perml(1)+
     &              dely2/perml(2)+delz2/perml(3))
            else
               pxy=sx2c*sx_mult*max(perml(1),perml(2),perml(3))
            endif
	    pxy = pxy*reduction_factor
            pxyi=pxy*(phikb-phii)
            pxyh=pxy*(pvikb-pvii)
            t1(neighc)=pxyi
            t2(neighc)=pxyh
            t3(neighc)=pxy
            t4(neighc)=pxy
            t6(neighc)=-grav*t3(neighc)
            t7(neighc)=-grav_air*t4(neighc)
 59      continue
c
c
c 2-d geometry
c
      elseif(icnl.ne.0) then
         radi=cord(iz,3)
         do 69 jm=1,iq
            kb=it8(jm)
            kz=kb-icd
            iwd=it10(jm)
            iw = abs(iwd)
            neighc=it9(jm)
            axkb=pnx(kb)
            aykb=pny(kb)
            alxkb=axkb
            alykb=aykb
            reduction_factor = red_factor(istrw_itfc(it11(jm)))
            perml(1)=2.*alxkb*alxi/(alxkb+alxi)
            perml(2)=2.*alykb*alyi/(alykb+alyi)
            radkb=0.5*(radi+cord(kz,3))
            sx2c=radkb*(sx(iw,isox)+sx(iw,isoy))
            pvikb=phi(kb)
            if (irdof .ne. 13) then
               phikb=pvikb-pcp(kb)
            else
               phikb=pvikb
            end if
            delx2=(cord(kz,1)-cord(iz,1))**2
            dely2=(cord(kz,igrav)-cord(iz,igrav))**2
            dis2=delx2+dely2
            if(dis2.gt.dis_tol.and.iwd.gt.0) then
              pxy=sx2c*dis2/(delx2/perml(1)+
     &           dely2/perml(2))
            else
              pxy=sx2c*sx_mult*max(perml(1),perml(2))
            endif
	    pxy = pxy*reduction_factor
            pxyi=pxy*(phikb-phii)
            pxyh=pxy*(pvikb-pvii)
            t1(neighc)=pxyi
            t2(neighc)=pxyh
            t3(neighc)=pxy
            t4(neighc)=pxy
            t6(neighc)=-grav*t3(neighc)
            t7(neighc)=-grav_air*t4(neighc)
 69      continue
      endif
c
c     liquid phase calculations

         do 70 jm=1,iq
            kb=it8(jm)
            kz=kb-icd
            neighc=it9(jm)
            pxyi=t1(neighc)
		  pxy=t3(neighc)
            sx4d=t6(neighc)
            pcakb = cap_fac*dzrg(kb)*(rlxyf(kb)-0.5)
            axyd=pxyi+0.5*sx4d
     &           *(rolf(i)+rolf(kb))
     &           *((cord(kz,igrav)-cord(iz,igrav))
     &           +(pcakb-pcai))
            t8(neighc)=axyd

 70      continue
c
c     determine upwind nodes and if liquid phase exists
c     wtsi_isot = 1 isotropic; wtsi_isot = 0 anisotropic
         isl=1
         do 71 jm=1,iq
            kb=it8(jm)
            kz=kb-icd
            neighc=it9(jm)
            if(wtsi_isot.ne.0) then
c calculate information needed for direction cosines
               delx2=(cord(kz,1)-cord(iz,1))**2
               if(icnl.eq.0) then
c 3-d problem
                  dely2=(cord(kz,2)-cord(iz,2))**2
               else
c 2-d problem
                  dely2=0.
               endif
               delz2=(cord(kz,igrav)-cord(iz,igrav))**2
               dxy2=delx2+dely2
			   dis2=delx2+dely2+delz2+dis_tol
            endif
            if(iad.le.iad_up) then
               axyd=t8(neighc)
               if(axyd.gt.0.0) then
                  imm=iz
               else
                  imm=kz
               endif

               if(wtsi_isot.eq.2) then
c
c calculate harmonic weighted relative permeability
c

                  rlxyd = 2.*rlxyf(iz)*rlxyf(kz)/(rlxyf(iz)+rlxyf(kz)
     &                 +tiny)
                  drlxyi = 2.*drlxyf(iz)*rlxyf(kz)/(rlxyf(iz)+rlxyf(kz)
     &                 +tiny) - 2.*rlxyf(iz)*rlxyf(kz)/(rlxyf(iz)
     &                 +rlxyf(kz)+tiny)**2*drlxyf(iz)
                  drlxykb= 2.*drlxyf(kz)*rlxyf(iz)/(rlxyf(iz)+rlxyf(kz)
     &                 +tiny) - 2.*rlxyf(iz)*rlxyf(kz)/(rlxyf(iz)
     &                 +rlxyf(kz)+tiny)**2*drlxyf(kz)
c not correct
                  rlpfree=rlxyd*(dxy2/dis2)+rlzf(imm)*(delz2/dis2)
                  dfi=drlxyi*(dxy2/dis2)+drlzf(imm)*(delz2/dis2)
                  dfkb=drlxykb*(dxy2/dis2)+drlzf(imm)*(delz2/dis2)
               else if(wtsi_isot.eq.1) then

                  rlpfree=rlxyf(imm)*(dxy2/dis2)+rlzf(imm)*(delz2/dis2)
                  df=drlxyf(imm)*(dxy2/dis2)+drlzf(imm)*(delz2/dis2)

               else
                  rlpfree=rlxyf(imm)
                  df = drlxyf(imm)
               endif

               if(axyd.gt.0.0) then
c water is flowing out of node i
                  if(wtsi_isot.eq.2) then
                     t91(neighc)=rlpfree
                     t9(neighc)=0.0
                     dfid1f(neighc)=dfi
                     dfidf(neighc)=dfkb
                  else
                     t91(neighc)=rlpfree
                     t9(neighc)=0.0
                     dfid1f(neighc)=df
                     dfidf(neighc)=0.0
                  end if
               else
c water is flowing into node i
                  if(wtsi_isot.eq.2) then
                     t9(neighc)=rlpfree
                     t91(neighc)=0.0
                     dfid1f(neighc)=dfi
                     dfidf(neighc)=dfkb
                  else
                     t9(neighc)=rlpfree
                     t91(neighc)=0.0
                     dfid1f(neighc)=0.0
                     dfidf(neighc)=df
                  end if
               end if

               call setbit(nbits,neighc,upwind_v(iz4m1),fid)
            else
c placeholder if_block:  we should never go here
               test_bit = bit(nbits,neighc,upwind_v(iz4m1))
               if(.not.test_bit) then
c placeholder if_block:  we should never go here
               endif
            endif
 71      continue
c
c     form equations
c
         if(isl.ne.0) then
            do 72 jm=1,iq
               kb=it8(jm)
               kz=kb-icd
               neighc=it9(jm)
               iau=it11(jm)
               ial=it12(jm)
               jml=nelmdg(kz)-neqp1
               heatc=t5(neighc)
               axyd=t8(neighc)

               fid=t9(neighc)
               fid1=t91(neighc)
               dfid=dfidf(neighc)
               dfid1=dfid1f(neighc)
               pxyi=t1(neighc)
               pxy=t3(neighc)
               sx4d=t6(neighc)
               dilkb=dil(kb)
               dilpkb=dilp(kb)
               pcakb = cap_fac*dzrg(kb)*(rlxyf(kb)-0.5)
	         axyd = axyd
	         dlpi=-pxy+0.5*sx4d*(dglp(i))
     &              *((cord(kz,igrav)-cord(iz,igrav))+(pcakb-pcai))
     &              -0.5*sx4d*(rolf(i)+rolf(kb))
     &              *cap_fac*dzrg(i)*drlxyf(i)

               dlpkb=pxy+0.5*sx4d*(dglp(kb))
     &              *((cord(kz,igrav)-cord(iz,igrav))+(pcakb-pcai))
     &              +0.5*sx4d*(rolf(i)+rolf(kb))
     &              *cap_fac*dzrg(kb)*drlxyf(kb)

c
c put derivative wrt free surface in
c

               axyf=(fid*dilkb+fid1*dili)
               axy=axyd*axyf
c
               dlapi=(dlpi*axyf+axyd*fid1*dilpi
     &         + axyd*dili*dfid1)

               dlapkb=(dlpkb*axyf+axyd*fid*dilpkb
     &         + axyd*dilkb*dfid)
               a_axy(iau+nmatavw)=axy
               a_axy(ial+nmatavw)=-axy

               bp(iz+nrhs(1))=bp(iz+nrhs(1))+axy
               bp(kz+nrhs(1))=bp(kz+nrhs(1))-axy
               a(jmia+nmat(1))=a(jmia+nmat(1))+dlapi
               a(ial+nmat(1))=a(ial+nmat(1))-dlapi
               a(iau+nmat(1))=a(iau+nmat(1))+dlapkb
               a(jml+nmat(1))=a(jml+nmat(1))-dlapkb
               if(kb.gt.neq_primary) then
                write(*,*) i, kb, jmia, ial, iau,jml, axy, dlapi,dlapkb
                pause
               endif
 72         continue
	endif
      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
!
! LBLRTM_r1p3_Module
!
! Module containing procedures for LBLRTM Record r1p3.
!
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, 24-Dec-2012
!                     paul.vandelst@noaa.gov
!

MODULE LBLRTM_r1p3_Module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE File_Utility   , ONLY: File_Open
  USE Message_Handler, ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  ! Line-by-line model parameters
  USE LBL_Parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Visibility
  ! ----------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: LBLRTM_r1p3_type
  ! Procedures
  PUBLIC :: LBLRTM_r1p3_Write

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Message string length
  INTEGER, PARAMETER :: ML = 256
  ! The record I/O format
  CHARACTER(*), PARAMETER :: LBLRTM_R1P3_FMT = '(8(es10.3),4x,i1,5x,es10.3,3x,i2)'

  ! -------------
  ! Derived types
  ! -------------
  TYPE :: LBLRTM_r1p3_type
    REAL(fp) :: v1        = 0.0_fp   !  Beginning wavenumber value for the calculation
    REAL(fp) :: v2        = 0.0_fp   !  Ending wavenumber value for the calculation
    REAL(fp) :: sample    = 4.0_fp   !  number of sample points per mean halfwidth (between 1 and 4)
    REAL(fp) :: dvset     = 0.0_fp   !  Selected DV for the final monochromatic calculation
    REAL(fp) :: alfal0    = 0.04_fp  !  Average collision broadened halfwidth
    REAL(fp) :: avmass    = 36.0_fp  !  Average molecular mass (amu) for Doppler halfwidth
    REAL(fp) :: dptmin    = 0.0_fp   !  Minimum molecular optical depth below which lines will be rejected
    REAL(fp) :: dptfac    = 0.0_fp   !  Factor multiplying molecular continuum optical depth for rejection check
    INTEGER  :: ilnflg    = 0        !  Flag for binary record of line rejection information
    REAL(fp) :: dvout     = 0.0_fp   !  Selected DV grid for the optical depth 'monochromatic' output spacing
    INTEGER  :: nmol_scal = 0        !  Enables the scaling of the atmospheric profile for selected species
  END TYPE LBLRTM_r1p3_type

CONTAINS

  FUNCTION LBLRTM_r1p3_Write(r1p3,fid) RESULT(err_stat)

    ! Arguments
    TYPE(LBLRTM_r1p3_type), INTENT(IN) :: r1p3
    INTEGER               , INTENT(IN) :: fid
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_r1p3_Write'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Setup
    err_stat = SUCCESS
    ! ...Check unit is open
    IF ( .NOT. File_Open(fid) ) THEN
      msg = 'File unit is not connected'
      CALL Cleanup(); RETURN
    END IF

    ! Write the record
    WRITE( fid,FMT=LBLRTM_R1P3_FMT,IOSTAT=io_stat,IOMSG=io_msg) r1p3
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing record - '//TRIM(io_msg)
      CALL Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE CleanUp()
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE CleanUp

  END FUNCTION LBLRTM_r1p3_Write

END MODULE LBLRTM_r1p3_Module

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


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

C     -*- Mode: fortran -*-
C=======================================================================
      SUBROUTINE SPOFF (XN, YN, ZN, NSPL, ZS, XS, XS2, YS, YS2,
     &     SLTOP, SLBOT, NUMNP, NDIM )
C=======================================================================

C     --*** SPLXYZ *** (GEN3D) Calculate 3D coordinates for Spline projection
C     --   Written by Amy Gilkey - revised 05/09/88
C     --   Modified by Greg Sjaardema - 02/06/89
C     --
C     --WRPXYZ calculates the coordinate array for the 3D warp translations.
C     --
C     --Parameters:
C     --   XN, YN, ZN - IN/OUT - the coordinates, modified
C     --   NSPL - Number of points on spline curve
C     --   ZS, XS, YS - IN - spline point coordinates
C     --   XS2, YS2 - spline slopes
C     --   SLTOP, SLBOT - slopes at top and bottom of spline
C     --   NUMNP - number of nodal points
C     --   NDIM - spatial dimension (2 or 3)
C     --

      REAL XN(*), YN(*), ZN(*)
      REAL ZS(NSPL), XS(NSPL), XS2(NSPL), YS(NSPL), YS2(NSPL)
      REAL SLTOP(2), SLBOT(2)

C$$$      CALL RDSPLN (.TRUE., NDB, SLTOP, SLBOT, NSPL, ZS, XS, YS)

      CALL SPLINE (ZS, XS, NSPL, SLTOP(1), SLBOT(1), XS2)
      CALL SPLINE (ZS, YS, NSPL, SLTOP(2), SLBOT(2), YS2)

      KLO = 1

      DO 80 INP = 1, NUMNP
         Z = ZN(INP)
         CALL HUNT (ZS, NSPL, Z, KLO)
         IF (KLO .EQ. NSPL) KLO = KLO - 1
         IF (KLO .EQ.    0) KLO = KLO + 1

         H  = ZS(KLO+1) - ZS(KLO)
         A  = (ZS(KLO+1)-Z)/H
         B  = (Z-ZS(KLO))/H

         XOFF = A * XS(KLO) + B * XS(KLO+1) +
     *        ((A**3-A) * XS2(KLO)+(B**3-B)*XS2(KLO+1)) * (H**2) / 6.

         YOFF = A * YS(KLO) + B * YS(KLO+1) +
     *        ((A**3-A) * YS2(KLO)+(B**3-B)*YS2(KLO+1)) * (H**2) / 6.

         XN(INP) = XN(INP) + XOFF
         YN(INP) = YN(INP) + YOFF
   80 CONTINUE

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!DESCRIPTION:                                                               !
!'quadsolveplot' is a program that solves and graphs quadratic functions    !
!depending on user entered coefficients a, b, and c of the general form     !
!being ax^2 + bx + c = 0. The program will work with the function even if   !
!coefficients equalling zero make the function linear. However, the program !
!will NOT work with unsolvable forms or ones where all real numbers         !
!are the solutions to x.                                                    !
!                                                                           !
!ORIGIN INFO:                                                               !
!Code written in gedit (pseudocode) and emacs (hard code)                   !
!Created Monday, September 13, 2013                                         !
!Author: Michael Conroy                                                     !
!Contact: sietekk@gmail.com                                                 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program quadsolvplot
use quadsolvplotmodule !Graphing module
implicit none

!DECLARE VARIABLES, ARRAYS, AND CONSTANTS
!The integers hold single bit values corresponding to choices made by queries
integer :: varcheck !Holds integer values from 1 to 4 according to
integer :: discrim !Holds value for discriminant
integer :: redoentry !If unsolvable, redoentry will be 1 = new vars, 0 = quit program
integer :: datfilecreate !For creating a graphable data file, either yes = 1 or no = 0
integer :: graphoption !For graphing query, either yes = 1 or no = 0
integer, parameter :: n=101 !For incremented values in x and y data/graphing arrays; 51 => extra place for 0.0 in arrays

!Variables to hold the entered coefficients, real roots, linear roots, and complex roots, and the discriminant (external function)
!a, b, c: user coefficiets; root1, root2: real roots; complex1, complex2: complex roots; linroot (a=0), linroot2 (a=c=0): linear roots
real :: a, b, c, root1, root2, linroot1, linroot2
complex :: comproot1, comproot2
real, external :: discrimcalcfunc !Computes the determinant

!Declare arrays for data file creation and graphing
real, dimension(n) :: x, y !2D arrays that hold x- and y-axis data

!!Subroutines (Not declared like functions => listed here only for my own benefit)
!!varentrysub(a,b,c)
!!realrootsub(a,b,discrim,root1,root2)
!!complexrootsub(a,b,discrim,complex1,complex2)
!!datafilesub(a,b,c,varcheck,n,x,y)

!DATA ENTRY
!Subroutine varentrysub is called on to prompt user to enter coefficients
10 call varentrysub(a,b,c) !10 is the goto location for the COEFFICIENT CHECK section

!COEFFICIENT CHECK
!1) Depending on coefficient values entered and the combinations of those values, a value of 1-7 is
!   given to 'varcheck', which will tell the program how to solve the function
!2) Detailed description of 'varcheck' values
!   If 1, then a=0 (linear); or 2, then b=0 (quadratic, double root); or 3, then c=0 (quadratic);
!   or 4, then a=b=0 (Unsolvable); or 5, then a=c=0 (linear); or 6, then b=c=0 (quadratic, double root);
!   or 7, then a=b=c=0 (all real numbers are roots)
if (a==0) then !Linear, one root at x=(-c)/b
	varcheck = 1
else if (b==0) then !Quadratic, double root at x=(+/-)sqrt((-c)/a)
	varcheck = 2
else if (c==0) then !Quadratic
	varcheck = 3
else if (a==0.and.b==o) then !Unsolvable, ASK TO RE-ENTER
	varcheck = 4
else if (a==0.and.c==0) then !Linear, one root at x=(-b)
	varcheck = 5
else if (b==0.and.c==0) then !Quadratic, double root at x=0
	varcheck = 6
else !For condition: a==0.and.b==0.and.c==0 = true, where a=b=c=0 (all are roots)
	varcheck = 7
end if

!CALL TO SOLUTION STRUCTURES AND PRINT ROOT(S)
!A case statement calls on a function or subroutine to solve for the discriminant, linear roots,
!real roots, and/or the complex roots. The quadratic function entered according to the
!coefficients is displayed along with the number of and the values of the root(s).
select case (varcheck)
	case (1) !a=0 (linear)
		print *, 'Equation is linear if a=0!'
		print *, 'Solving as a linear equation...'
		linroot1 = (-c)/b
		print *, 'You entered coefficients for the linear function ', b, 'x', ' + ', c, ' = 0,'
		print *, 'and this function has one real root at x = ', linroot1, '.'
11 	case (4) !a=b=0 (unsolvable)
		print *, 'Equation is unsolvable if a=b=0!'
		print *, 'Do you want to re-enter coefficients or quit? ('1'=re-entry, '0'=quit program):'
		read *, redoentry
		if (redoentry==0) then
			print *, 'Quitting now...'
			stop
		else if (redoentry==1) then
			goto 10 !Returns program back to varentrysub() subroutine to re-enter variables
		else !For condition: redoentry/=1.or.redoentry/=0
			print *, 'WARNING: INVALID INPUT! Please try again...'
			goto 11
		end if
	case (5) !a=c=0 (linear)
		print *, 'Equation is linear if a=c=0!'
		print *, 'Solving as a linear equation...'
		linroot1 = (-b)
		print *, 'You entered coefficients for the linear function ', b, 'x', ' + ', c, ' = 0,'
		print *, 'and this function has one real root at x = ', linroot2, '.'
	case (7) !a=b=c=0 (all are roots)
		print *, 'All real numbers are solutions!'
	case default !All quadratic options: cases (2, 3, and 6)
		print *, 'Equation is quadratic!'
		print *, 'Solving as a quadratic equation...'
		discrim = discrimcalcfunc(a,b,c)
		if (discrim>0) then !Real roots
			call realrootsub(a,b,discrim,root1,root2)
			print *, 'You entered coefficients for the quadratic function ', a, 'x^2 + ', b, 'x', ' + ', c, ' = 0,'
			print *, 'and this function has two real roots at x = ', root1, ' and x = ', root2, '.'
		else if (discrim==0) !Works for both b=0 and b=c=0
			call realrootsub(a,b,discrim,root1,root2)
			print *, 'You entered coefficients for the quadratic function ', a, 'x^2 + ', b, 'x', ' + ', c, ' = 0,'
			print *, 'and this function has a double real root at x = ', root1, ' and x = ', root2, '.'
		else !For condition: discrim<0
			call complexrootsub(a,b,discrim,complex1,complex2)
			print *, 'You entered coefficients for the linear function ', a, 'x^2 + ', b, 'x', ' + ', c, ' = 0,'
			print *, 'and this function has two complex roots at x = ', complex1, ' and x = ', complex2, '.'
end select

!DATA FILE CREATION
!Data file creation query
if (varcheck==1.or.varcheck==5) then !Linear
	print *, 'Your function is linear! Do you wish to save graphable data for your function? ('1'=re-entry, '0'=quit program):'
	read *, datfilecreate
12 	if (datfilecreate==0) then
		print *, 'Exiting data file creation...'
		stop
	else if (datfilecreate==1) then
		call datafilesub(a,b,c,varcheck,n,x,y)
	else !For condition: datfilecreate/=1.or.datfile/=0
		print *, 'WARNING: INVALID INPUT! Please try again...'
		goto 12
	end if
else if (varcheck==2.or.varcheck==3.or.varcheck==6) then
	print *, 'Your function is quadratic! Do you wish to save graphable data for your function? ('1'=re-entry, '0'=quit program):'
	read *, datfilecreate
13 	if (datfilecreate==0) then
		print *, 'Exiting data file creation...'
		stop
	else if (datfilecreate==1) then
		print *, 'Creating data file...'
		call datafilesub(a,b,c,varcheck,n,x,y)
	else !For condition: datfilecreate/=1.or.datfile/=0
		print *, 'WARNING: INVALID INPUT! Please try again...'
		goto 13
	end if
else !For condition: varcheck==4.or.varcheck==7
	print *, 'Sorry, your function is neither quadratic nor linear, so data cannot be generated for you.'
	print *, 'Exiting data file creation...'
end if

!GRAPHING THE FUNCTION
14 print *, 'Graph the function with the new data file?  ('1'= yes, '0'= no):'
read *, graphoption
if (graphoption==0) then
	print *, 'Exiting the graphing module...'
	stop
else if (graphoption==1) then
	print *, 'Graphing data...'
	call graph(x,y,n)
else !For condition: datfilecreate/=1.or.datfile/=0
	print *, 'WARNING: INVALID INPUT! Please try again...'
	goto 14
end if

!END OF THE MAIN PROGRAM
end program quadsolveplot

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

!SUBROUTINES
subroutine varentrysub(a,b,c)
	implicit none

	!DEFINE VARIANLES
	real,intent(inout) :: a,b,c

	!BODY OF SUBROUTINE
	print *, 'This program solves for the roots of a quadratic equation of the form ax^2 + bx + c = 0,'
	print *, 'and will solve for real and complex roots. It will also identify if the equation is'
	print *, 'unsolvable or the roots are all real numbers. The option to save a file with graphable'
	print *, 'data will be presented along with the option to graph the functioned entered.'
	print *, '...'
	print *, 'Please enter the values of the coefficients a, b, and c separated by commas or spaces:'
	read *, a,b,c
	print *, 'You entered the values a = ', a, ', b = ', b, ', c = ', c, '.' !Not sure if this prints variable values inside text
end subroutine varentrysub

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

subroutine datafilesub(a,b,c,varcheck,n,x,y)
	implicit none

	!DEFINE VARIANLES
	real, intent(in) :: a, b, c
	integer, intent(in) :: varcheck, n
	real, dimension(n) intent(inout) :: x,y
	integer :: bitcheck !For internal select case statement
	integer :: i !x,y populator do loop counters

	!BODY OF SUBROUTINE
	if (varcheck==1.or.varcheck==5) checkbit=1
	if (varcheck==2.or.varcheck==3.or.varcheck==6) checkbit=2
	print *, 'Initializing data file creation...'
	do i=-50,50 !Do loop to populate rest of x array
		x(i) = i
	end do
	y(i)=0.0 !Define zero for y values in y array population
	select case (bitcheck)
		case (1) !Linear, data will be constructed accordingly
			do i=-50,50
				y(i)=b*i+c
			end do

		case (2) !Quadratic, data will be constructed accordingly
			do i=-50,50
				y(i)=a(i-root1)(i-root2)
			end do
	end select
	open(unit=99,file='quadsolvplot.data',status='replace')
	do i=-50,50
		write(99,*) x(i), ' ', y(i)
	end do
	close(unit=99)
	print *, 'Data arrays generated and saved!'
end subroutine datafilesub

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

max_stars_repo_name: sietekk/quadsolveplot

max_issues_repo_name: sietekk/quadsolveplot

max_forks_repo_name: sietekk/quadsolveplot


text:
module m_std
    use iso_c_binding;
end module m_std

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text: #redirect J. Alex Tomaich, DDS MD

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine pebi3 (px,py,pz, cpb,vol)
!***********************************************************************
!  Copyright, 1993, 2004,  The  Regents of the University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1  PURPOSE
CD1
CD1  To calculate internodal volumes using perpendicular bisectors.
CD1
C***********************************************************************
CD2
CD2  REVISION HISTORY
CD2
CD2 $Log:   /pvcs.config/fehm90/src/pebi3.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:36   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:11:50   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:11:30   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:25:46   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:05:50   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:44:12 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.3   Thu Jan 18 10:06:16 1996   zvd
CD2 Modifications to prolog
CD2
CD2    Rev 1.2   Wed Jan 10 15:26:58 1996   hend
CD2 Added Prolog
CD2
CD2    Rev 1.1   03/18/94 15:57:52   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:26:12   pvcs
CD2 original version in process of being certified
CD2
C***********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3  2.2 Finite-Element Coefficient Generation
CD3
C***********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C***********************************************************************
CD5
CD5 INTERFACES
CD5
CD5 Formal Calling Parameters
CD5
CD5   Identifier      Type     Use  Description
CD5
CD5   cpb             REAL*8   O    Perpindicular bisector
CD5   px              REAL*8   I    X nodal coordinates
CD5   py              REAL*8   I    Y nodal coordinates
CD5   pz              REAL*8   I    Z nodal coordinates
CD5   vol             REAL*8   O    Internodal volume
CD5
CD5 Interface Tables
CD5
CD5   None
CD5
CD5 Files
CD5
CD5   None
CD5
C***********************************************************************
CD6
CD6 GLOBAL OBJECTS
CD6
CD6 Global Constants
CD6
CD6   None
CD6
CD6 Global Types
CD6
CD6   None
CD6
CD6 Global Variables
CD6
CD6   None
CD6
CD6 Global Subprograms
CD6
CD6   Identifier      Type     Description
CD6
CD6   lubksub0                 Solves the set of n linear equations A*x = b.
CD6   ludcmp0                  Given an n by n matrix A, with physical
CD6                              dimension np, this routine replaces it
CD6                              by the LU decompostion of a row wise
CD6                              permutation of itself.
CD6
C***********************************************************************
CD7
CD7 LOCAL IDENTIFIERS
CD7
CD7 Local Constants
CD7
CD7   None
CD7
CD7 Local Types
CD7
CD7   None
CD7
CD7 Local variables
CD7
CD7   Identifier      Type     Description
CD7
CD7   amat            REAL*8   Array
CD7   b               REAL*8   Vector, right-hand side
CD7   cx              REAL*8
CD7   cy              REAL*8
CD7   cz              REAL*8
CD7   d               REAL*8   Flag denoting even or odd number of
CD7                              row interchanges
CD7   i               INT      Loop index
CD7   i1              INT      Array index
CD7   i2              INT      Array index
CD7   i3              INT      Array index
CD7   i4              INT      Array index
CD7   iar             INT
CD7   indx            INT      Permutation vector
CD7   iperm           INT
CD7   j               INT      Loop index
CD7   k               INT      Array index
CD7   l               INT      Array index
CD7   n3              INT      Array dimension
CD7   qx              REAL*8
CD7   qy              REAL*8
CD7   qz              REAL*8
CD7   vol1            REAL*8
CD7   vol2            REAL*8
CD7   vol3            REAL*8
CD7   vol4            REAL*8
CD7   vol5            REAL*8
CD7   x1              REAL*8
CD7   x2              REAL*8
CD7   x3              REAL*8
CD7   x4              REAL*8
CD7   xm              REAL*8
CD7   y1              REAL*8
CD7   y2              REAL*8
CD7   y3              REAL*8
CD7   y4              REAL*8
CD7   ym              REAL*8
CD7   z1              REAL*8
CD7   z2              REAL*8
CD7   z3              REAL*8
CD7   z4              REAL*8
CD7   zm              REAL*8
CD7
CD7 Local Subprograms
CD7
CD7   None
CD7
C***********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN pebi3
CPS
CPS   FOR each coordinate position
CPS       FOR each incremental coordinate position
CPS           calculate coordinate midpoints
CPS       END FOR
CPS   END FOR
CPS
CPS   FOR each coordinate position
CPS       set indices and calculate internodal lengths, etc.
CPS       call ludcmp0 and lubksb0 to solve for bisector coordinates
CPS   END FOR
CPS
CPS   calculate internodal lengths, etc. from first coordinate to each . . .
CPS   . . . other coordinate
CPS   call ludcmp0 and lubksb0 to solve for coordinates
CPS
CPS   FOR each coordinate position
CPS       FOR each incremental coordinate position
CPS           IF x internode length is greater than y and z internode length
CPS              use y and z bisector coordinates to calculate . . .
CPS              . . . perpindicular bisector
CPS           ELSE IF y internode length is greater than z internode length
CPS              use x and z bisector coordinates to calculate . . .
CPS              . . . perpindicular bisector
CPS           ELSE
CPS              use x and y bisector coordinates to calculate . . .
CPS              . . . perpindicular bisector
CPS           END IF
CPS           set perpindicular bisector
CPS       END FOR
CPS   END FOR
CPS
CPS   FOR each coordinate position
CPS       set indices
CPS       calculate internodal volumes
CPS   END FOR
CPS
CPS END pebi3
CPS
C***********************************************************************

      implicit none

      integer n3, i, i1, i2, i3, i4, j, k, l
      parameter (n3 = 3)
      integer indx(n3), iperm(4, 4), iar(4, 4)
      real*8 cpb(4, 4), vol(4), px(4), py(4), pz(4)
      real*8 cx(4), cy(4), cz(4), xm(4, 4), ym(4, 4), zm(4, 4)
      real*8 amat(n3, n3), b(n3), d
      real*8 arx, ary, arz, qx, qy, qz
      real*8 vol1, vol2, vol3, vol4, vol5
      real*8 x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4

      data iperm/1,2,3,4, 2,1,4,3, 3,1,2,4, 4,1,3,2/
      data iar /0,4,2,3, 3,0,4,1, 4,1,0,2, 2,3,1,0/

      do i = 1, 3
         do j = i+1, 4
	    xm(i,j) = 0.5 * (px(i)+px(j))
	    ym(i,j) = 0.5 * (py(i)+py(j))
	    zm(i,j) = 0.5 * (pz(i)+pz(j))
	    xm(j,i) = xm(i,j)
	    ym(j,i) = ym(i,j)
	    zm(j,i) = zm(i,j)
         enddo
      enddo

      do i = 1, 4
         i1 = iperm(1,i)
         i2 = iperm(2,i)
         i3 = iperm(3,i)
         i4 = iperm(4,i)

         amat(1,1) = px(i2)-px(i3)
         amat(1,2) = py(i2)-py(i3)
         amat(1,3) = pz(i2)-pz(i3)
         amat(2,1) = px(i2)-px(i4)
         amat(2,2) = py(i2)-py(i4)
         amat(2,3) = pz(i2)-pz(i4)
         amat(3,1) = (py(i2)-py(i3))*(pz(i2)-pz(i4)) -
     &        (pz(i2)-pz(i3))*(py(i2)-py(i4))
         amat(3,2) = (pz(i2)-pz(i3))*(px(i2)-px(i4)) -
     &        (px(i2)-px(i3))*(pz(i2)-pz(i4))
         amat(3,3) = (px(i2)-px(i3))*(py(i2)-py(i4)) -
     &        (py(i2)-py(i3))*(px(i2)-px(i4))

         b(1) = xm(i2,i3)*amat(1,1)+ym(i2,i3)*amat(1,2)+zm(i2,i3)*
     &        amat(1,3)
         b(2) = xm(i2,i4)*amat(2,1)+ym(i2,i4)*amat(2,2)+zm(i2,i4)*
     &        amat(2,3)
         b(3) = px(i2)*amat(3,1)+py(i2)*amat(3,2)+pz(i2)*amat(3,3)

         call ludcmp0 (amat,n3,n3,indx,d)
         call lubksb0 (amat,n3,n3,indx,b)
         cx(i1) = b(1)
         cy(i1) = b(2)
         cz(i1) = b(3)
      enddo

      amat(1,1) = px(1)-px(2)
      amat(1,2) = py(1)-py(2)
      amat(1,3) = pz(1)-pz(2)
      amat(2,1) = px(1)-px(3)
      amat(2,2) = py(1)-py(3)
      amat(2,3) = pz(1)-pz(3)
      amat(3,1) = px(1)-px(4)
      amat(3,2) = py(1)-py(4)
      amat(3,3) = pz(1)-pz(4)

      b(1) = xm(1,2)*amat(1,1) + ym(1,2)*amat(1,2) + zm(1,2)*amat(1,3)
      b(2) = xm(1,3)*amat(2,1) + ym(1,3)*amat(2,2) + zm(1,3)*amat(2,3)
      b(3) = xm(1,4)*amat(3,1) + ym(1,4)*amat(3,2) + zm(1,4)*amat(3,3)

      call ludcmp0 (amat,n3,n3,indx,d)
      call lubksb0 (amat,n3,n3,indx,b)
      qx = b(1)
      qy = b(2)
      qz = b(3)

      do i = 1, 3
         do j = i+1, 4
	    k = iar(i,j)
	    l = iar(j,i)
	    if (abs(px(i)-px(j)) .ge. abs(py(i)-py(j)) .and.
     &           abs(px(i)-px(j)) .ge. abs(pz(i)-pz(j))) then
               arx = (qy-ym(i,j))*(cz(k)-cz(l)) -
     &              (qz-zm(i,j))*(cy(k)-cy(l))
               cpb(i,j) = 0.5 * arx / (px(i)-px(j))
	    else if (abs(py(i)-py(j)) .ge. abs(pz(i)-pz(j))) then
               ary = (qz-zm(i,j))*(cx(k)-cx(l)) -
     &              (qx-xm(i,j))*(cz(k)-cz(l))
               cpb(i,j) = 0.5 * ary / (py(i)-py(j))
	    else
               arz = (qx-xm(i,j))*(cy(k)-cy(l)) -
     &              (qy-ym(i,j))*(cx(k)-cx(l))
               cpb(i,j) = 0.5 * arz / (pz(i)-pz(j))
	    end if
	    cpb(j,i) = cpb(i,j)
         enddo
      enddo

      do i = 1, 4
         i1 = iperm(1,i)
         i2 = iperm(2,i)
         i3 = iperm(3,i)
         i4 = iperm(4,i)
         x1 = qx
         y1 = qy
         z1 = qz
         x2 = cx(i4)
         y2 = cy(i4)
         z2 = cz(i4)
         x3 = cx(i3)
         y3 = cy(i3)
         z3 = cz(i3)
         x4 = cx(i2)
         y4 = cy(i2)
         z4 = cz(i2)
         vol1 = (y1*z2-y2*z1)*(x3-x4) + (y1*z3-y3*z1)*(x4-x2) +
     &        (y1*z4-y4*z1)*(x2-x3) + (y2*z3-y3*z2)*(x1-x4) +
     &        (y2*z4-y4*z2)*(x3-x1) + (y3*z4-y4*z3)*(x1-x2)

         x1 = px(i1)
         y1 = py(i1)
         z1 = pz(i1)
         x3 = cx(i2)
         y3 = cy(i2)
         z3 = cz(i2)
         x4 = cx(i3)
         y4 = cy(i3)
         z4 = cz(i3)
         vol2 = (y1*z2-y2*z1)*(x3-x4) + (y1*z3-y3*z1)*(x4-x2) +
     &        (y1*z4-y4*z1)*(x2-x3) + (y2*z3-y3*z2)*(x1-x4) +
     &        (y2*z4-y4*z2)*(x3-x1) + (y3*z4-y4*z3)*(x1-x2)

         x2 = xm(i1,i3)
         y2 = ym(i1,i3)
         z2 = zm(i1,i3)
         x3 = cx(i2)
         y3 = cy(i2)
         z3 = cz(i2)
         x4 = cx(i4)
         y4 = cy(i4)
         z4 = cz(i4)
         vol3 = (y1*z2-y2*z1)*(x3-x4) + (y1*z3-y3*z1)*(x4-x2) +
     &        (y1*z4-y4*z1)*(x2-x3) + (y2*z3-y3*z2)*(x1-x4) +
     &        (y2*z4-y4*z2)*(x3-x1) + (y3*z4-y4*z3)*(x1-x2)

         x2 = xm(i1,i4)
         y2 = ym(i1,i4)
         z2 = zm(i1,i4)
         x3 = cx(i3)
         y3 = cy(i3)
         z3 = cz(i3)
         x4 = cx(i2)
         y4 = cy(i2)
         z4 = cz(i2)
         vol4 = (y1*z2-y2*z1)*(x3-x4) + (y1*z3-y3*z1)*(x4-x2) +
     &        (y1*z4-y4*z1)*(x2-x3) + (y2*z3-y3*z2)*(x1-x4) +
     &        (y2*z4-y4*z2)*(x3-x1) + (y3*z4-y4*z3)*(x1-x2)

         x2 = xm(i1,i2)
         y2 = ym(i1,i2)
         z2 = zm(i1,i2)
         x3 = cx(i4)
         y3 = cy(i4)
         z3 = cz(i4)
         x4 = cx(i3)
         y4 = cy(i3)
         z4 = cz(i3)
         vol5 = (y1*z2-y2*z1)*(x3-x4) + (y1*z3-y3*z1)*(x4-x2) +
     &        (y1*z4-y4*z1)*(x2-x3) + (y2*z3-y3*z2)*(x1-x4) +
     &        (y2*z4-y4*z2)*(x3-x1) + (y3*z4-y4*z3)*(x1-x2)

         vol(i1) = (vol1+vol2+vol3+vol4+vol5) / 6.
      enddo

      return
      end


text:
*----------------------------------------------------------------------*
      subroutine wr_formula(ffform,form)
*----------------------------------------------------------------------*
*     write next formula record to ffform
*----------------------------------------------------------------------*

      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
      include 'def_operator.h'
      include 'def_filinf.h'
      include 'def_formula_item.h'

      type(filinf), intent(in) ::
     &     ffform
      type(formula_item), intent(in) ::
     &     form
      integer ::
     &     lu

      lu = ffform%unit

      ! write command record
      write(lu) form%command, form%target

      select case(form%command)
      case(command_add_contribution)
        call rw_contr_kernel(-1,lu,form%contr)
      case(command_del_intermediate)
        write(lu) form%label
      case(command_new_intermediate)
        call rw_opdef_kernel(-1,lu,form%interm,form%incore,
     &                       form%parent1,form%parent2,
     &                       form%tra,form%tra1,form%tra2)
      case(command_reorder)
        call rw_reo_kernel(-1,lu,form%reo)
      case(command_add_bc_reo,command_bc_reo,command_add_reo)
        call rw_bcontr_kernel(-1,lu,form%bcontr)
        call rw_reo_kernel(-1,lu,form%reo)
      case(command_add_intm,command_cp_intm,command_add_bc,command_bc)
        call rw_bcontr_kernel(-1,lu,form%bcontr)
      end select

      return

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!
!===============================================================================
!-- time_marching Module
!> Module time_marching structures and code
!===============================================================================
!/****/h* modules|time_marching/time_marching_mod
!* SYNOPSIS
MODULE time_marching_mod
!* PURPOSE
!*   Module for global variables, included library, and time_marching code
!* INCLUDES
!*   Name                  Purpose
!* CONTAINS
!*   Subroutine              Purpose
!*
!* AUTHOR
!*   Pieter Boom
!* MODIFICATION HISTORY
!*   2019/08/20: Created (PB)
!* COPYRIGHT
!*   (c) University of Manchester
!******/
!> author: Pieter Boom
!> date: 2019/08/20
!> Module time_marching structures and code
!===============================================================================

  !-- use statements and implicit none --
  USE common_mod
  USE mpi_mod

  IMPLICIT NONE

#include "time_marching.inc"

!===============================================================================
CONTAINS
!===============================================================================
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!/****/s* time_marching_mod/initialise_time_marching
!* SYNOPSIS
  SUBROUTINE initialise_time_marching()
!* PURPOSE
!*   Deallocate global structures/variables
!* SIDE EFFECTS
!*   Global structure (lcl_complex) deallocated
!* AUTHOR
!*   Pieter Boom
!* MODIFICATION HISTORY
!*   2019/08/21: Created (PB)
!* COPYRIGHT
!*   (c) University of Manchester
!******/
!> author: Pieter Boom
!> date: 2019/08/21
!> Deallocate global structures/variables
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IMPLICIT NONE

    !-- local variables --
    INTEGER         :: i

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! MAIN EXECUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !-- deallocate primary elements --

    IF (ALLOCATED(time_marching%A)) DEALLOCATE(time_marching%A)
    IF (ALLOCATED(time_marching%B)) DEALLOCATE(time_marching%B)
    IF (ALLOCATED(time_marching%U)) DEALLOCATE(time_marching%U)
    IF (ALLOCATED(time_marching%V)) DEALLOCATE(time_marching%V)
    IF (ALLOCATED(time_marching%c)) DEALLOCATE(time_marching%c)
    IF (ALLOCATED(time_marching%dc)) DEALLOCATE(time_marching%dc)
    IF (ALLOCATED(time_marching%LHS_tvec)) DEALLOCATE(time_marching%LHS_tvec)
    IF (ALLOCATED(time_marching%q_step_indx)) DEALLOCATE(time_marching%q_step_indx)

    SELECTCASE(time_marching_method)

    CASE('explicit_euler','ee')
      time_marching%implicit = 0
      time_marching%p = 1;    time_marching%q = 1
      time_marching%r = 1;    time_marching%s = 1;    time_marching%s_initial = 1

      ALLOCATE(time_marching%A(1,1),time_marching%U(1,1),time_marching%B(1,1),&
        time_marching%V(1,1),time_marching%c(1),time_marching%q_step_indx(1))

      time_marching%A = transpose(reshape((/ 0.d0 /),(/ time_marching%s,time_marching%s /)))
      time_marching%B = transpose(reshape((/ 1.d0 /),(/ time_marching%s,time_marching%r /)))
      time_marching%U = transpose(reshape((/ 1.d0 /),(/ time_marching%r,time_marching%s /)))
      time_marching%V = transpose(reshape((/ 1.d0 /),(/ time_marching%r,time_marching%r /)))
      time_marching%c = (/ 0.d0 /)

    CASE('implicit_euler','ie')
      time_marching%implicit = 1
      time_marching%stifflyAcc = .TRUE.
      time_marching%p = 1;    time_marching%q = 1
      time_marching%r = 1;    time_marching%s = 1;    time_marching%s_initial = 1

      ALLOCATE(time_marching%A(1,1),time_marching%U(1,1),time_marching%B(1,1),&
        time_marching%V(1,1),time_marching%c(1),time_marching%q_step_indx(1))

      time_marching%A = transpose(reshape((/ 1.d0 /),(/ time_marching%s,time_marching%s /)))
      time_marching%B = transpose(reshape((/ 1.d0 /),(/ time_marching%s,time_marching%r /)))
      time_marching%U = transpose(reshape((/ 1.d0 /),(/ time_marching%r,time_marching%s /)))
      time_marching%V = transpose(reshape((/ 1.d0 /),(/ time_marching%r,time_marching%r /)))
      time_marching%c = (/ 1.d0 /)

    CASE('bdf2')
      time_marching%implicit = 1
      time_marching%stifflyAcc = .TRUE.
      time_marching%p = 2;  time_marching%q = 2
      time_marching%r = 2;  time_marching%s = 1;  time_marching%s_initial = 1

      !call AllocateUnsteady(unsteady)

      ALLOCATE(time_marching%A(1,1),time_marching%U(1,2),time_marching%B(2,1),&
        time_marching%V(2,2),time_marching%c(1),time_marching%q_step_indx(2))

      time_marching% A(1,1) = 2.d0/3.d0
      time_marching% B = transpose(reshape((/ &
           2.d0/3.d0, &
           0.d0 &
           /),(/ time_marching% s, time_marching% r /)))
      time_marching% U = transpose(reshape((/ &
           4.d0/3.d0, -1.d0/3.d0 &
           /),(/ time_marching% r, time_marching% s /)))
      time_marching% V = transpose(reshape((/ &
           4.d0/3.d0, -1.d0/3.d0, &
           1.d0, 0.d0 &
           /),(/ time_marching% r, time_marching% r /)))
           time_marching%c = (/ 1.d0 /)

       CASE('sdirk3')
          time_marching%implicit = 1
          time_marching%stifflyAcc = .TRUE.
          time_marching%p = 3;  time_marching%q = 1
          time_marching%r = 1;  time_marching%s = 4;  time_marching%s_initial = 1

          !call AllocateUnsteady(unsteady)

          ALLOCATE(time_marching%A(4,4),time_marching%U(4,1),time_marching%B(1,4),&
            time_marching%V(1,1),time_marching%c(4),time_marching%q_step_indx(1))

          time_marching%A = transpose(reshape( (/ &
               5.7281606245015908d-1,0.0000000000000000d0,0.0000000000000000d0,0.0000000000000000d0, &
               -4.2974075486898089d-2,5.7281606245015908d-1,0.0000000000000000d0,0.0000000000000000d0, &
               -3.8405397546217372d0,3.6079433958911875d0,5.7281606245015908d-1,0.0000000000000000d0, &
               -9.9166439219916533d0,1.1013593846890206d1,-6.6976598734871517d-1,5.7281606245015908d-1 &
               /),(/ time_marching%s, time_marching%s /)))
          time_marching%B = transpose(reshape( (/ &
               -9.9166439219916533d0,1.1013593846890206d1,-6.6976598734871517d-1,5.7281606245015908d-1 &
               /),(/ time_marching%s, time_marching%r /)))
          time_marching%U = transpose(reshape( (/ &
               1.0000000000000000d0, &
               1.0000000000000000d0, &
               1.0000000000000000d0, &
               1.0000000000000000d0 &
               /),(/ time_marching%r, time_marching%s /)))
          time_marching%V = transpose(reshape( (/ &
               1.0000000000000000d0 &
               /),(/ time_marching%r, time_marching%r /)))
           time_marching% c = SUM(time_marching% A,DIM=2)

       CASE('sdirk4')
          time_marching%implicit = 1
          time_marching%stifflyAcc = .FALSE.
          time_marching%p = 4;  time_marching%q = 1
          time_marching%r = 1;  time_marching%s = 3;  time_marching%s_initial = 1

          !call AllocateUnsteady(unsteady)

          ALLOCATE(time_marching%A(3,3),time_marching%U(3,1),time_marching%B(1,3),&
            time_marching%V(1,1),time_marching%c(3),time_marching%q_step_indx(1))

          time_marching%A = transpose(reshape( (/ &
               1.0685790213016289d0,0.0000000000000000d0,0.0000000000000000d0, &
               -5.6857902130162898d-1,1.0685790213016289d0,0.0000000000000000d0, &
               2.1371580426032577d0,-3.2743160852065154d0,1.0685790213016289d0 &
               /),(/ time_marching%s, time_marching%s /)))
          time_marching%B = transpose(reshape( (/ &
               1.2888640051572038d-1,7.4222719896855924d-1,1.2888640051572038d-1 &
               /),(/ time_marching%s, time_marching%r /)))
          time_marching%U = transpose(reshape( (/ &
               1.0000000000000000d0, &
               1.0000000000000000d0, &
               1.0000000000000000d0 &
               /),(/ time_marching%r, time_marching%s /)))
          time_marching%V = transpose(reshape( (/ &
               1.0000000000000000d0 &
               /),(/ time_marching%r, time_marching%r /)))
           time_marching% c = SUM(time_marching% A,DIM=2)

       CASE('trap')
         time_marching%implicit = 1
         time_marching%stifflyAcc = .TRUE.
         time_marching%p = 2;  time_marching%q = 2
         time_marching%r = 1;  time_marching%s = 2;  time_marching%s_initial = 2

         !call AllocateUnsteady(unsteady)

         ALLOCATE(time_marching%A(2,2),time_marching%U(2,1),time_marching%B(1,2),&
           time_marching%V(1,1),time_marching%c(2),time_marching%q_step_indx(1))

         time_marching%A = transpose(reshape( (/ &
              0.0000000000000000d0,0.0000000000000000d0, &
              5.0000000000000000d-1,5.0000000000000000d-1 &
              /),(/ time_marching%s, time_marching%s /)))
         time_marching%B = transpose(reshape( (/ &
              5.0000000000000000d-1,5.0000000000000000d-1 &
              /),(/ time_marching%s, time_marching%r /)))
         time_marching%U = transpose(reshape( (/ &
              1.0000000000000000d0, &
              1.0000000000000000d0 &
              /),(/ time_marching%r, time_marching%s /)))
         time_marching%V = transpose(reshape( (/ &
              1.0000000000000000d0 &
              /),(/ time_marching%r, time_marching%r /)))
         time_marching% c = SUM(time_marching% A,DIM=2)

       CASE('esdirk4')
         time_marching%implicit = 1
         time_marching%stifflyAcc = .TRUE.
         time_marching%p = 4;  time_marching%q = 2
         time_marching%r = 1;  time_marching%s = 6;  time_marching%s_initial = 2

         !call AllocateUnsteady(unsteady)

         ALLOCATE(time_marching%A(6,6),time_marching%U(6,1),time_marching%B(1,6),&
           time_marching%V(1,1),time_marching%c(6),time_marching%q_step_indx(1))

         time_marching% A = transpose(reshape( (/ &
              0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
              1.d0/4.d0, 1.d0/4.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
              8611.d0/62500.d0, -1743.d0/31250.d0, 1.d0/4.d0, 0.d0, 0.d0, 0.d0, &
              5012029.d0/34652500.d0, -654441.d0/2922500.d0, 174375.d0/388108.d0, &
                1.d0/4.d0, 0.d0, 0.d0, &
              15267082809.d0/155376265600.d0, -71443401.d0/120774400.d0, &
                730878875.d0/902184768.d0, 2285395.d0/8070912.d0, 1.d0/4.d0, 0.d0, &
              82889.d0/524892.d0, 0.d0, 15625.d0/83664.d0, 69875.d0/102672.d0, &
                -2260.d0/8211.d0, 1.d0/4.d0 &
              /),(/ time_marching% s, time_marching% s /)))
         time_marching% B = transpose(reshape((/ &
              82889.d0/524892.d0, 0.d0, 15625.d0/83664.d0, 69875.d0/102672.d0, &
                -2260.d0/8211.d0, 1.d0/4.d0 &
              /),(/ time_marching% s, time_marching% r /)))
         time_marching% U = transpose(reshape((/ &
              1.d0, &
              1.d0, &
              1.d0, &
              1.d0, &
              1.d0, &
              1.d0 &
              /),(/ time_marching% r, time_marching% s /)))
         time_marching% V(1,1) = 1.d0
         time_marching% c = SUM(time_marching% A,DIM=2)

    CASE DEFAULT
      CALL end_mpi()
    END SELECT

    DO i=1,time_marching%r;   time_marching%q_step_indx(i) = i;    END DO

    IF (time_marching%implicit == 0) THEN;  time_marching%A(i,:) = dt*time_marching%A(i,:)
    ELSEIF (time_marching%implicit == 1) THEN
      ALLOCATE(time_marching%LHS_tvec(time_marching%s))
      DO i = time_marching%s_initial,time_marching%s
        time_marching%LHS_tvec(i) = 1.d0/(time_marching%A(i,i) * dt)
        time_marching%U(i,:) = time_marching%U(i,:) * time_marching%LHS_tvec(i)
        time_marching%A(i,:) = time_marching%A(i,:) / time_marching%A(i,i)
      END DO
      ALLOCATE(time_marching%dc(time_marching%s+1))
      time_marching%dc(1) = time_marching%c(1)
      DO i = 2,time_marching%s
        time_marching%dc(i) = time_marching%c(i) - time_marching%c(i-1)
      END DO
      time_marching%dc(i) = 1.d0 - time_marching%c(i-1)
    END IF
    time_marching%B = dt*time_marching%B

    RETURN

  END SUBROUTINE
! initialise_time_marching
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

END MODULE
! time_marching_mod
!===============================================================================
!

max_stars_repo_name: ParaFEM/ParaGEMS

max_issues_repo_name: ParaFEM/ParaGEMS

max_forks_repo_name: ParaFEM/ParaGEMS


max_stars_repo_name: umarcor/svunit

max_issues_repo_name: umarcor/svunit

max_forks_repo_name: umarcor/svunit


text:
module legendre_mod
    implicit none

contains
    ! Computes the associated Legendre polynomials P_n^m,
    ! with m = m and n = 0, ..., n_max
    ! They include the Condon-Shortley phase
    subroutine associated_legendre(x, n_max, m, p_nm, p_nm_d)
        real(8), intent(in) :: x
        integer, intent(in) :: n_max
        integer, intent(in) :: m
        real(8), allocatable, intent(out) :: p_nm(:)
        real(8), allocatable, intent(out) :: p_nm_d(:)

        integer :: i
        real(8) :: theta
        real(8) :: parity_factor(0:n_max)
        real(8) :: dummy_array(n_max+1)
        integer :: n_array(0:n_max)

        allocate(p_nm(0:n_max))
        allocate(p_nm_d(0:n_max))

        if (x .le. -1d0 .or. x .ge. 1d0) then
            print *, "Error: associated_legendre argument (x) out of range"
            p_nm = 0d0
            p_nm_d = 0d0
            return
        end if

        if (x .lt. 0d0) then
            parity_factor = [((-1d0) ** (i + m), i=0,n_max)]
        else
            parity_factor = 1d0
        end if

        ! compute associated Legendre polynomial
        theta = acos(abs(x))
        call XDLEGF(0d0, n_max, m, m, theta, 3, p_nm, dummy_array)
        p_nm = p_nm * parity_factor

        ! compute first derivate using recurrence
        ! (x^2 - 1)(d/dx)p_nm(x) = n * x * p_nm(x) - (n + m) * p_{n-1, m}(x)
        n_array = [(i, i=0,n_max)]
        p_nm_d(1:n_max) = n_array(1:n_max) * x * p_nm(1:n_max)  &
                          - (n_array(1:n_max) + m) * p_nm(0:n_max-1)
        p_nm_d(0) = 0d0
        p_nm_d = p_nm_d / (x ** 2 - 1)

    end subroutine

end module

max_stars_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_issues_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_forks_repo_name: Riarrieta/WindowedGreenFunctionMethod


text:
MU East Conference Room is located across from the Information Booth in the east side of the Memorial Union. The room is rectangular and has very loud air conditioning which makes it a bad place for meetings. It is often used by ASUCD to announce the results of ASUCD Elections when there is bad weather. It is also used for storage by the UCD Bookstore during the bookbuy rush at the beginning of each quarter.

There are rumors that the UCD administration is considering bringing a bank to fill in the space.