text:
Module m_uniinv
Integer, Parameter :: kdp = selected_real_kind(15)
public :: uniinv
private :: kdp
private :: R_uniinv, I_uniinv, D_uniinv
private :: R_nearless, I_nearless, D_nearless, nearless
interface uniinv
  module procedure d_uniinv, r_uniinv, i_uniinv
end interface uniinv
interface nearless
  module procedure D_nearless, R_nearless, I_nearless
end interface nearless
contains

Subroutine D_uniinv (XDONT, IGOEST)
! __________________________________________________________
!   UNIINV = Merge-sort inverse ranking of an array, with removal of
!   duplicate entries.
!   The routine is similar to pure merge-sort ranking, but on
!   the last pass, it sets indices in IGOEST to the rank
!   of the value in the ordered set with duplicates removed.
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! __________________________________________________________
      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
      Integer, Dimension (:), Intent (Out) :: IGOEST
! __________________________________________________________
      Real (kind=kdp) :: XTST, XDONA, XDONB
!
! __________________________________________________________
      Integer, Dimension (SIZE(IGOEST)) :: JWRKT, IRNGT
      Integer :: LMTNA, LMTNC, IRNG, IRNG1, IRNG2, NUNI
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      NVAL = Min (SIZE(XDONT), SIZE(IGOEST))
!
      Select Case (NVAL)
      Case (:0)
         Return
      Case (1)
         IGOEST (1) = 1
         Return
      Case Default
         Continue
      End Select
!
!  Fill-in the index array, creating ordered couples
!
      Do IIND = 2, NVAL, 2
         If (XDONT(IIND-1) < XDONT(IIND)) Then
            IRNGT (IIND-1) = IIND - 1
            IRNGT (IIND) = IIND
         Else
            IRNGT (IIND-1) = IIND
            IRNGT (IIND) = IIND - 1
         End If
      End Do
      If (Modulo (NVAL, 2) /= 0) Then
         IRNGT (NVAL) = NVAL
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      LMTNA = 2
      LMTNC = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (NVAL <= 4) Exit
!
!   Loop on merges of A and B into C
!
         Do IWRKD = 0, NVAL - 1, 4
            If ((IWRKD+4) > NVAL) Then
               If ((IWRKD+2) >= NVAL) Exit
!
!   1 2 3
!
               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
!
!   1 3 2
!
               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
                  IRNG2 = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNG2
!
!   3 1 2
!
               Else
                  IRNG1 = IRNGT (IWRKD+1)
                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNG1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
!
!   1 3 x x
!
            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   1 3 2 4
                  IRNGT (IWRKD+3) = IRNG2
               Else
!   1 3 4 2
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+4) = IRNG2
               End If
!
!   3 x x x
!
            Else
               IRNG1 = IRNGT (IWRKD+1)
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
                  IRNGT (IWRKD+2) = IRNG1
                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   3 1 2 4
                     IRNGT (IWRKD+3) = IRNG2
                  Else
!   3 1 4 2
                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                     IRNGT (IWRKD+4) = IRNG2
                  End If
               Else
!   3 4 1 2
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+3) = IRNG1
                  IRNGT (IWRKD+4) = IRNG2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (2*LMTNA >= NVAL) Exit
         IWRKF = 0
         LMTNC = 2 * LMTNC
!
!   Loop on merges of A and B into C
!
         Do
            IWRK = IWRKF
            IWRKD = IWRKF + 1
            JINDA = IWRKF + LMTNA
            IWRKF = IWRKF + LMTNC
            If (IWRKF >= NVAL) Then
               If (JINDA >= NVAL) Exit
               IWRKF = NVAL
            End If
            IINDA = 1
            IINDB = JINDA + 1
!
!  One steps in the C subset, that we create in the final rank array
!
!  Make a copy of the rank array for the iteration
!
            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
            XDONA = XDONT (JWRKT(IINDA))
            XDONB = XDONT (IRNGT(IINDB))
!
            Do
               IWRK = IWRK + 1
!
!  We still have unprocessed values in both A and B
!
               If (XDONA > XDONB) Then
                  IRNGT (IWRK) = IRNGT (IINDB)
                  IINDB = IINDB + 1
                  If (IINDB > IWRKF) Then
!  Only A still with unprocessed values
                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
                     Exit
                  End If
                  XDONB = XDONT (IRNGT(IINDB))
               Else
                  IRNGT (IWRK) = JWRKT (IINDA)
                  IINDA = IINDA + 1
                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
                  XDONA = XDONT (JWRKT(IINDA))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 2 * LMTNA
      End Do
!
!   Last merge of A and B into C, with removal of duplicates.
!
      IINDA = 1
      IINDB = LMTNA + 1
      NUNI = 0
!
!  One steps in the C subset, that we create in the final rank array
!
      JWRKT (1:LMTNA) = IRNGT (1:LMTNA)
      If (IINDB <= NVAL) Then
        XTST = NEARLESS (Min(XDONT(JWRKT(1)), XDONT(IRNGT(IINDB))))
      Else
        XTST = NEARLESS (XDONT(JWRKT(1)))
      Endif
      Do IWRK = 1, NVAL
!
!  We still have unprocessed values in both A and B
!
         If (IINDA <= LMTNA) Then
            If (IINDB <= NVAL) Then
               If (XDONT(JWRKT(IINDA)) > XDONT(IRNGT(IINDB))) Then
                  IRNG = IRNGT (IINDB)
                  IINDB = IINDB + 1
               Else
                  IRNG = JWRKT (IINDA)
                  IINDA = IINDA + 1
               End If
            Else
!
!  Only A still with unprocessed values
!
               IRNG = JWRKT (IINDA)
               IINDA = IINDA + 1
            End If
         Else
!
!  Only B still with unprocessed values
!
            IRNG = IRNGT (IWRK)
         End If
         If (XDONT(IRNG) > XTST) Then
            XTST = XDONT (IRNG)
            NUNI = NUNI + 1
         End If
         IGOEST (IRNG) = NUNI
!
      End Do
!
      Return
!
End Subroutine D_uniinv

Subroutine R_uniinv (XDONT, IGOEST)
! __________________________________________________________
!   UNIINV = Merge-sort inverse ranking of an array, with removal of
!   duplicate entries.
!   The routine is similar to pure merge-sort ranking, but on
!   the last pass, it sets indices in IGOEST to the rank
!   of the value in the ordered set with duplicates removed.
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! _________________________________________________________
      Real, Dimension (:), Intent (In) :: XDONT
      Integer, Dimension (:), Intent (Out) :: IGOEST
! __________________________________________________________
      Real    :: XTST, XDONA, XDONB
!
! __________________________________________________________
      Integer, Dimension (SIZE(IGOEST)) :: JWRKT, IRNGT
      Integer :: LMTNA, LMTNC, IRNG, IRNG1, IRNG2, NUNI
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      NVAL = Min (SIZE(XDONT), SIZE(IGOEST))
!
      Select Case (NVAL)
      Case (:0)
         Return
      Case (1)
         IGOEST (1) = 1
         Return
      Case Default
         Continue
      End Select
!
!  Fill-in the index array, creating ordered couples
!
      Do IIND = 2, NVAL, 2
         If (XDONT(IIND-1) < XDONT(IIND)) Then
            IRNGT (IIND-1) = IIND - 1
            IRNGT (IIND) = IIND
         Else
            IRNGT (IIND-1) = IIND
            IRNGT (IIND) = IIND - 1
         End If
      End Do
      If (Modulo (NVAL, 2) /= 0) Then
         IRNGT (NVAL) = NVAL
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      LMTNA = 2
      LMTNC = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (NVAL <= 4) Exit
!
!   Loop on merges of A and B into C
!
         Do IWRKD = 0, NVAL - 1, 4
            If ((IWRKD+4) > NVAL) Then
               If ((IWRKD+2) >= NVAL) Exit
!
!   1 2 3
!
               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
!
!   1 3 2
!
               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
                  IRNG2 = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNG2
!
!   3 1 2
!
               Else
                  IRNG1 = IRNGT (IWRKD+1)
                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNG1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
!
!   1 3 x x
!
            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   1 3 2 4
                  IRNGT (IWRKD+3) = IRNG2
               Else
!   1 3 4 2
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+4) = IRNG2
               End If
!
!   3 x x x
!
            Else
               IRNG1 = IRNGT (IWRKD+1)
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
                  IRNGT (IWRKD+2) = IRNG1
                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   3 1 2 4
                     IRNGT (IWRKD+3) = IRNG2
                  Else
!   3 1 4 2
                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                     IRNGT (IWRKD+4) = IRNG2
                  End If
               Else
!   3 4 1 2
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+3) = IRNG1
                  IRNGT (IWRKD+4) = IRNG2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (2*LMTNA >= NVAL) Exit
         IWRKF = 0
         LMTNC = 2 * LMTNC
!
!   Loop on merges of A and B into C
!
         Do
            IWRK = IWRKF
            IWRKD = IWRKF + 1
            JINDA = IWRKF + LMTNA
            IWRKF = IWRKF + LMTNC
            If (IWRKF >= NVAL) Then
               If (JINDA >= NVAL) Exit
               IWRKF = NVAL
            End If
            IINDA = 1
            IINDB = JINDA + 1
!
!  One steps in the C subset, that we create in the final rank array
!
!  Make a copy of the rank array for the iteration
!
            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
            XDONA = XDONT (JWRKT(IINDA))
            XDONB = XDONT (IRNGT(IINDB))
!
            Do
               IWRK = IWRK + 1
!
!  We still have unprocessed values in both A and B
!
               If (XDONA > XDONB) Then
                  IRNGT (IWRK) = IRNGT (IINDB)
                  IINDB = IINDB + 1
                  If (IINDB > IWRKF) Then
!  Only A still with unprocessed values
                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
                     Exit
                  End If
                  XDONB = XDONT (IRNGT(IINDB))
               Else
                  IRNGT (IWRK) = JWRKT (IINDA)
                  IINDA = IINDA + 1
                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
                  XDONA = XDONT (JWRKT(IINDA))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 2 * LMTNA
      End Do
!
!   Last merge of A and B into C, with removal of duplicates.
!
      IINDA = 1
      IINDB = LMTNA + 1
      NUNI = 0
!
!  One steps in the C subset, that we create in the final rank array
!
      JWRKT (1:LMTNA) = IRNGT (1:LMTNA)
      If (IINDB <= NVAL) Then
        XTST = NEARLESS (Min(XDONT(JWRKT(1)), XDONT(IRNGT(IINDB))))
      Else
        XTST = NEARLESS (XDONT(JWRKT(1)))
      Endif
      Do IWRK = 1, NVAL
!
!  We still have unprocessed values in both A and B
!
         If (IINDA <= LMTNA) Then
            If (IINDB <= NVAL) Then
               If (XDONT(JWRKT(IINDA)) > XDONT(IRNGT(IINDB))) Then
                  IRNG = IRNGT (IINDB)
                  IINDB = IINDB + 1
               Else
                  IRNG = JWRKT (IINDA)
                  IINDA = IINDA + 1
               End If
            Else
!
!  Only A still with unprocessed values
!
               IRNG = JWRKT (IINDA)
               IINDA = IINDA + 1
            End If
         Else
!
!  Only B still with unprocessed values
!
            IRNG = IRNGT (IWRK)
         End If
         If (XDONT(IRNG) > XTST) Then
            XTST = XDONT (IRNG)
            NUNI = NUNI + 1
         End If
         IGOEST (IRNG) = NUNI
!
      End Do
!
      Return
!
End Subroutine R_uniinv
Subroutine I_uniinv (XDONT, IGOEST)
! __________________________________________________________
!   UNIINV = Merge-sort inverse ranking of an array, with removal of
!   duplicate entries.
!   The routine is similar to pure merge-sort ranking, but on
!   the last pass, it sets indices in IGOEST to the rank
!   of the value in the ordered set with duplicates removed.
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! __________________________________________________________
      Integer, Dimension (:), Intent (In)  :: XDONT
      Integer, Dimension (:), Intent (Out) :: IGOEST
! __________________________________________________________
      Integer :: XTST, XDONA, XDONB
!
! __________________________________________________________
      Integer, Dimension (SIZE(IGOEST)) :: JWRKT, IRNGT
      Integer :: LMTNA, LMTNC, IRNG, IRNG1, IRNG2, NUNI
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      NVAL = Min (SIZE(XDONT), SIZE(IGOEST))
!
      Select Case (NVAL)
      Case (:0)
         Return
      Case (1)
         IGOEST (1) = 1
         Return
      Case Default
         Continue
      End Select
!
!  Fill-in the index array, creating ordered couples
!
      Do IIND = 2, NVAL, 2
         If (XDONT(IIND-1) < XDONT(IIND)) Then
            IRNGT (IIND-1) = IIND - 1
            IRNGT (IIND) = IIND
         Else
            IRNGT (IIND-1) = IIND
            IRNGT (IIND) = IIND - 1
         End If
      End Do
      If (Modulo (NVAL, 2) /= 0) Then
         IRNGT (NVAL) = NVAL
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      LMTNA = 2
      LMTNC = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (NVAL <= 4) Exit
!
!   Loop on merges of A and B into C
!
         Do IWRKD = 0, NVAL - 1, 4
            If ((IWRKD+4) > NVAL) Then
               If ((IWRKD+2) >= NVAL) Exit
!
!   1 2 3
!
               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
!
!   1 3 2
!
               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
                  IRNG2 = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNG2
!
!   3 1 2
!
               Else
                  IRNG1 = IRNGT (IWRKD+1)
                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
                  IRNGT (IWRKD+2) = IRNG1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
!
!   1 3 x x
!
            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   1 3 2 4
                  IRNGT (IWRKD+3) = IRNG2
               Else
!   1 3 4 2
                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+4) = IRNG2
               End If
!
!   3 x x x
!
            Else
               IRNG1 = IRNGT (IWRKD+1)
               IRNG2 = IRNGT (IWRKD+2)
               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
                  IRNGT (IWRKD+2) = IRNG1
                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
!   3 1 2 4
                     IRNGT (IWRKD+3) = IRNG2
                  Else
!   3 1 4 2
                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
                     IRNGT (IWRKD+4) = IRNG2
                  End If
               Else
!   3 4 1 2
                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
                  IRNGT (IWRKD+3) = IRNG1
                  IRNGT (IWRKD+4) = IRNG2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (2*LMTNA >= NVAL) Exit
         IWRKF = 0
         LMTNC = 2 * LMTNC
!
!   Loop on merges of A and B into C
!
         Do
            IWRK = IWRKF
            IWRKD = IWRKF + 1
            JINDA = IWRKF + LMTNA
            IWRKF = IWRKF + LMTNC
            If (IWRKF >= NVAL) Then
               If (JINDA >= NVAL) Exit
               IWRKF = NVAL
            End If
            IINDA = 1
            IINDB = JINDA + 1
!
!  One steps in the C subset, that we create in the final rank array
!
!  Make a copy of the rank array for the iteration
!
            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
            XDONA = XDONT (JWRKT(IINDA))
            XDONB = XDONT (IRNGT(IINDB))
!
            Do
               IWRK = IWRK + 1
!
!  We still have unprocessed values in both A and B
!
               If (XDONA > XDONB) Then
                  IRNGT (IWRK) = IRNGT (IINDB)
                  IINDB = IINDB + 1
                  If (IINDB > IWRKF) Then
!  Only A still with unprocessed values
                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
                     Exit
                  End If
                  XDONB = XDONT (IRNGT(IINDB))
               Else
                  IRNGT (IWRK) = JWRKT (IINDA)
                  IINDA = IINDA + 1
                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
                  XDONA = XDONT (JWRKT(IINDA))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 2 * LMTNA
      End Do
!
!   Last merge of A and B into C, with removal of duplicates.
!
      IINDA = 1
      IINDB = LMTNA + 1
      NUNI = 0
!
!  One steps in the C subset, that we create in the final rank array
!
      JWRKT (1:LMTNA) = IRNGT (1:LMTNA)
      If (IINDB <= NVAL) Then
        XTST = NEARLESS (Min(XDONT(JWRKT(1)), XDONT(IRNGT(IINDB))))
      Else
        XTST = NEARLESS (XDONT(JWRKT(1)))
      Endif
      Do IWRK = 1, NVAL
!
!  We still have unprocessed values in both A and B
!
         If (IINDA <= LMTNA) Then
            If (IINDB <= NVAL) Then
               If (XDONT(JWRKT(IINDA)) > XDONT(IRNGT(IINDB))) Then
                  IRNG = IRNGT (IINDB)
                  IINDB = IINDB + 1
               Else
                  IRNG = JWRKT (IINDA)
                  IINDA = IINDA + 1
               End If
            Else
!
!  Only A still with unprocessed values
!
               IRNG = JWRKT (IINDA)
               IINDA = IINDA + 1
            End If
         Else
!
!  Only B still with unprocessed values
!
            IRNG = IRNGT (IWRK)
         End If
         If (XDONT(IRNG) > XTST) Then
            XTST = XDONT (IRNG)
            NUNI = NUNI + 1
         End If
         IGOEST (IRNG) = NUNI
!
      End Do
!
      Return
!
End Subroutine I_uniinv

Function D_nearless (XVAL) result (D_nl)
!  Nearest value less than given value
! __________________________________________________________
      Real (kind=kdp), Intent (In) :: XVAL
      Real (kind=kdp) :: D_nl
! __________________________________________________________
      D_nl = nearest (XVAL, -1.0_kdp)
      return
!
End Function D_nearless
Function R_nearless (XVAL) result (R_nl)
!  Nearest value less than given value
! __________________________________________________________
      Real, Intent (In) :: XVAL
      Real :: R_nl
! __________________________________________________________
      R_nl = nearest (XVAL, -1.0)
      return
!
End Function R_nearless
Function I_nearless (XVAL) result (I_nl)
!  Nearest value less than given value
! __________________________________________________________
      Integer, Intent (In) :: XVAL
      Integer :: I_nl
! __________________________________________________________
      I_nl = XVAL - 1
      return
!
End Function I_nearless

end module m_uniinv

max_stars_repo_name: ZoliQua/orthologs


!
!     --------------------------------------------------------------
!     R E C O U P L S 4
!     --------------------------------------------------------------
!                                                                  *
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville            February 1994   *
!
      SUBROUTINE RECOUPLS4(K, JA1, JA2, JA3, JA4, K1, K2, K3, K4, KA, IRE, IAT&
         , REC)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE CONSTS_C
      USE MEDEFN_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  11:00:59  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE dlsa3_I
      USE dlsa2_I
      USE dlsa4_I
      USE dlsa1_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: K
      INTEGER  :: JA1
      INTEGER  :: JA2
      INTEGER  :: JA3
      INTEGER  :: JA4
      INTEGER  :: K1
      INTEGER  :: K2
      INTEGER  :: K3
      INTEGER  :: K4
      INTEGER  :: KA
      INTEGER  :: IRE
      INTEGER  :: IAT
      REAL(DOUBLE) , INTENT(OUT) :: REC
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IA4, IB4, ISKR
      REAL(DOUBLE) :: S1, S2, S3, S4, S, RE
!-----------------------------------------------
      S1 = DBLE(J1QN1(JA1,K))
      S2 = DBLE(J1QN1(JA2,K))
      S3 = DBLE(J1QN1(JA3,K))
      S4 = DBLE(J1QN1(JA4,K))
      S = S1*S2*S3*S4
      REC = ONE/DSQRT(S)
      IA4 = J1QN1(JA4,K) - 1
      IB4 = J1QN2(JA4,K) - 1
      REC = REC*DSQRT(DBLE(IA4 + 1))/DSQRT(DBLE((K4 + 1)*(IB4 + 1)))
!
      ISKR = JA3 - JA2
      IF (ISKR > 1) THEN
         IAT = 0
         CALL DLSA3 (K, JA2, JA3, KA, IRE, IAT, RE)
         IF (IAT == 0) RETURN
         REC = RE*REC
      ENDIF
!
      ISKR = JA4 - JA3
      IF (ISKR > 1) THEN
         IAT = 0
         CALL DLSA3 (K, JA3, JA4, K4, IRE, IAT, RE)
         IF (IAT == 0) RETURN
         REC = RE*REC
      ENDIF
!
      IAT = 0
      CALL DLSA2 (K, JA1, JA4, K4, IRE, IAT, RE)
      IF (IAT == 0) RETURN
      REC = RE*REC
!
      IAT = 0
      CALL DLSA4 (K, JA1, JA2, K1, K2, KA, IRE, IAT, RE)
      IF (IAT == 0) RETURN
      REC = RE*REC
!
      IAT = 0
      CALL DLSA4 (K, JA2, JA3, KA, K3, K4, IRE, IAT, RE)
      IF (IAT == 0) RETURN
      REC = RE*REC
      IF (JA1==1 .AND. JA2==2) RETURN
!
      IAT = 0
      CALL DLSA1 (K, JA1, K1, IRE, IAT, RE)
      IF (IAT == 0) RETURN
      REC = RE*REC
!
      ISKR = JA2 - JA1
      IF (JA1 == 1) ISKR = JA2 - 1 - JA1
      IF (ISKR <= 1) RETURN
      IAT = 0
      CALL DLSA3 (K, JA1, JA2, K1, IRE, IAT, RE)
      REC = RE*REC
      RETURN
      END SUBROUTINE RECOUPLS4


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details
C=======================================================================
      SUBROUTINE RDSTEP (NDB, NCSTEP,
     &     NUMELB, IDELB, ISEVOK,
     &     TIME, VARGL, VARNP, VAREL,
     $     neldm)
C=======================================================================

C   --*** RDSTEP *** (EXPLORE) Read current database variables
C   --
C   --RDSTEP reads the variable data for the next time step from the
C   --database.
C   --
C   --Parameters:
C   --   NDB - IN - the database number
C   --   NCSTEP - IN/OUT - the current step number
C   --   NVARGL - IN - the number of global variables
C   --   NVARNP - IN - the number of nodal variables
C   --   NVAREL - IN - the number of element variables
C   --   NUMELB - IN - the number of elements per block
C   --   ISEVOK - IN - the element block variable truth table;
C   --      variable i of block j exists iff ISEVOK(i,j) is NOT 0
C   --   TIME - OUT - the time step time
C   --   VARGL - OUT - the global variables for the time step
C   --   VARNP - OUT - the nodal variables for the time step
C   --   VAREL - OUT - the element variables for the time step

      include 'exodusII.inc'
      include 'exp_dbnums.blk'

      INTEGER NUMELB(*)
      INTEGER IDELB(*)
      INTEGER ISEVOK(nvarel,*)
      REAL VARGL(*)
      REAL VARNP(numnp, *)
      REAL VAREL(numel, *)

      CHARACTER*80 ERRMSG

C ... Initialize variables to known value in case error on read
      TIME = -999.9
      CALL INIREA (NVARGL, -999.9, VARGL)
      CALL INIREA (NVARNP * NUMNP, -999.9, VARNP)
      CALL INIREA (NVAREL * NUMEL, -999.9, VAREL)

C ... Read the time for the specified step
      call exgtim(ndb, ncstep, time, ierr)
      if (ierr .ne. 0) go to 150

C ... Read the global variables (if any)
      if (nvargl .gt. 0) then
         call exggv (ndb, ncstep, nvargl, vargl, ierr)
         if (ierr .ne. 0) go to 170
      end if

C ... Read the nodal variables (if any)
      DO 120 I = 1, NVARNP
        call exgnv (ndb, ncstep, i, numnp, varnp(1,i), ierr)
        if (ierr .ne. 0) go to 180
 120  CONTINUE

C ... Read the element variables (if any)
      DO 140 I = 1, NVAREL
        IEL0 = 1
        DO 130 IELB = 1, NELBLK
          IF (ISEVOK(I,IELB) .NE. 0 .AND. numelb(ielb) .GT. 0) THEN
             call exgev (ndb, ncstep, i, idelb(ielb),
     &            numelb(ielb), varel(iel0,i), ierr)
             if (ierr .ne. 0) go to 190
          END IF
          IEL0 = IEL0 + NUMELB(IELB)
 130    CONTINUE
 140  CONTINUE

      RETURN

  150 CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     &   'TIME on step', NCSTEP
      GOTO 200
  170 CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     &   'GLOBAL VARIABLES on step', NCSTEP
      GOTO 200
  180 CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     &   'NODAL VARIABLE', I, ' on step', NCSTEP
      GOTO 200
  190 CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     &   'ELEMENT VARIABLE', I, ' in block', IELB, ' on step', NCSTEP
      GOTO 200
  200 CONTINUE
      CALL WDBERR (IERR, ERRMSG)

      RETURN

10000  FORMAT (5 (A, I12))
      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine hns(b,theta,rho,dg,sqrts0,xflow,h1,h2)
!
!     determine the flow depth h2 downstream of a hydraulic jump,
!     corresponding to a upstream flow depth of h1
!
      implicit none
!
      real*8 b,rho,dg,sqrts0,xflow,h1,h2,c2,f,df,dh2,hk,
     &  xflow2,tth,A1,yg1,A2,yg2,dA2dh2,dyg2dh2,theta
!
      call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
!
      h2=h1*(-1.d0+dsqrt(1.d0+8.d0*(hk/h1)**3))/2.d0
!
      if(dabs(theta).lt.1.d-10) return
!
!     hns for a trapezoid, non-rectangular cross section
!
      c2=rho*rho*dg*sqrts0
      xflow2=xflow*xflow
      tth=dtan(theta)
      A1=h1*(b+h1*tth)
      yg1=h1*(3.d0*b+2.d0*h1*tth)/(6.d0*(b+h1*tth))
!
!     Newton-Raphson iterations
!
      do
         A2=h2*(b+h2*tth)
         yg2=h2*(3.d0*b+2.d0*h2*tth)/(6.d0*(b+h2*tth))
         dA2dh2=b+2.d0*h2*tth
         dyg2dh2=((3.d0*b+4.d0*h2*tth)*(b+tth)
     &                    -tth*h2*(3.d0*b+2.d0*h2*tth))/
     &                    (6.d0*(b+h2*tth)**2)
         f=A2*xflow2+c2*(A1*A1*A2*yg1-A1*A2*A2*yg2)-A1*xflow2
         df=dA2dh2*xflow2+c2*(A1*A1*yg1*dA2dh2-2.d0*A1*A2*dA2dh2*yg2
     &                       -A1*A2*A2*dyg2dh2)
         dh2=f/df
         if(dabs(dh2)/h2.lt.1.d-3) exit
         h2=h2-dh2
      enddo
!
      write(*,*) 'hns ','h1= ',h1,'h2= ',h2,'hk= ',hk
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      subroutine oceaalbe(pws,paw,xsal,pcl,pwl,
     s           brdfalbe)
C INPUT:  pws=speed of wind (in m/s)
C         paw=azim. of sun - azim. of wind (in deg.)
C	  xsal=salinity (in ppt)
C	  pcl=pigment concentration (in mg.m-3)
C         pwl=wavelength of the computation (in micrometer)
C OUTPUT: brdfalbe=the spherical albedo of the ocean
C
      real Ref(39)
      real pwl,azw,pcl,wl,wspd,C,pws,brdfalbe,w,wlp,paw
      real ref_i,rwc,rw,rogalbe,a,rwb,xsal
      real nr,ni
      integer iwl

c effective reflectance of the whitecaps (Koepke, 1984)
      data Ref/
     &0.220,0.220,0.220,0.220,0.220,0.220,0.215,0.210,0.200,0.190,
     &0.175,0.155,0.130,0.080,0.100,0.105,0.100,0.080,0.045,0.055,
     &0.065,0.060,0.055,0.040,0.000,0.000,0.000,0.000,0.000,0.000,
     &0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/
C conversion of parameter
      C=pcl
      wspd=pws
      azw=paw
      wl=pwl

C COMPUTE WHITECAPS REFLECTANCE (LAMBERTIAN)
      W=2.95e-06*(wspd**3.52)
      iwl=1+int((wl-0.2)/0.1)
      wlp=0.5+(iwl-1)*0.1
      Ref_i=ref(iwl+1)+(wl-wlp)/0.1*(ref(iwl)-ref(iwl+1))
      Rwc=W*Ref_i
C COMPUTE WATER REFRACTION INDEX
      call indwat(wl,xsal,nr,ni)
C COMPUTE BACKSCATTERED REFLECTANCE FROM THE SEA WATER (LAMBERTIAN)
C  water reflectance below the sea surface
      call morcasiwat(wl,C,Rw)
C SUNGLINT spherical albedo
      call glitalbe(wspd,nr,ni,azw,rogalbe)
C  water reflectance above the sea surface, (albedo re=0.485)
C albedo is a=re is taken from table 2 of Austin,1974,The remote sensing
C of spectral radiance from below the ocean surface, in Optical
C Aspects of Oceanography (N.G. Jerlov and E. Steeman Nielsen,Eds),
C Academic,London,pp. 317-344
      a=0.485
      Rwb=(1.-Rogalbe)*(1.-a)*Rw/(1-a*Rw)
C SPHERICAL ALBEDO OF SEA WATER
      brdfalbe=Rwc+(1-W)*Rogalbe+(1-Rwc)*Rwb
      return
      end

max_stars_repo_name: xcorail/OTB

max_issues_repo_name: xcorail/OTB

max_forks_repo_name: xcorail/OTB


text:
      program prog
      character ld*1(3, 2)
      character ls*3(2)
      equivalence(ld, ls)
      data ls /'Yui', 'dFg'/
      write(6, *) ld
      end


text:
MODULE APF
  IMPLICIT NONE

  INTEGER, PRIVATE, PARAMETER :: APFMAXNUMARGS = 16, APFMAXARGLEN = 128
  INTEGER, PRIVATE, PARAMETER :: APFINT = 0, APFDBL = 1

  INTEGER, PRIVATE, SAVE :: N = 0

  TYPE INTPTR
     INTEGER, POINTER :: P
  ENDTYPE INTPTR

  TYPE DBLPTR
     DOUBLE PRECISION, POINTER :: P
  ENDTYPE DBLPTR

  INTEGER, PRIVATE, SAVE :: TYPES(APFMAXNUMARGS)
  TYPE(INTPTR), PRIVATE, SAVE :: INTPTRS (APFMAXNUMARGS)
  TYPE(DBLPTR), PRIVATE, SAVE :: DBLPTRS(APFMAXNUMARGS)
  CHARACTER(LEN=APFMAXARGLEN), PRIVATE, SAVE :: NAMES(APFMAXARGLEN)

CONTAINS

  SUBROUTINE APF_INT_ADD(IPTR, NAME)
    IMPLICIT NONE

    INTEGER, POINTER, INTENT(IN) :: IPTR
    CHARACTER(LEN=*), INTENT(IN) :: NAME

    N = N + 1
    TYPES(N) = APFINT
    INTPTRS(N)%P => IPTR
    NAMES(N) = NAME
  ENDSUBROUTINE APF_INT_ADD

  SUBROUTINE APF_DBL_ADD(DPTR, NAME)
    IMPLICIT NONE

    DOUBLE PRECISION, POINTER, INTENT(IN) :: DPTR
    CHARACTER(LEN=*), INTENT(IN) :: NAME

    N = N + 1
    TYPES(N) = APFDBL
    DBLPTRS(N)%P => DPTR
    NAMES(N) = NAME
  ENDSUBROUTINE APF_DBL_ADD

  SUBROUTINE APF_PARSE()
    IMPLICIT NONE

    INTEGER :: I, K, ARGC
    INTEGER :: IAUX, STAT
    DOUBLE PRECISION :: DAUX
    CHARACTER(LEN=APFMAXARGLEN) :: ARGV

    ARGC = COMMAND_ARGUMENT_COUNT()

    I = 1
    DO WHILE (I .LE. ARGC)
       CALL GET_COMMAND_ARGUMENT(I, ARGV)

       DO K = 1, N
          IF (ARGV .NE. NAMES(K)) CYCLE

          IF (TYPES(K) .EQ. APFINT) THEN
             IF (K+1 .GT. ARGC) THEN
                PRINT *, 'Expected integer after ', TRIM(ARGV)
                STOP
             ENDIF

             I = I + 1
             CALL GET_COMMAND_ARGUMENT(I, ARGV)
             READ (ARGV, *, IOSTAT=STAT) IAUX

             IF (STAT .NE. 0) THEN
                PRINT *, 'Invalid integer found'
                STOP
             ENDIF

             INTPTRS(K)%P = IAUX
          ELSEIF (TYPES(K) .EQ. APFDBL) THEN
             IF (K+1 .GT. ARGC) THEN
                PRINT *, 'Expected double after ', TRIM(ARGV)
                STOP
             ENDIF

             I = I + 1
             CALL GET_COMMAND_ARGUMENT(I, ARGV)

             READ (ARGV, *, IOSTAT=STAT) DAUX

             IF (STAT .NE. 0) THEN
                PRINT *, 'Invalid double found'
                STOP
             ENDIF

             DBLPTRS(K)%P = DAUX
          ENDIF
          EXIT
       ENDDO

       I = I + 1
    ENDDO
  ENDSUBROUTINE APF_PARSE

ENDMODULE APF

max_stars_repo_name: vsart/apf

max_issues_repo_name: vsart/apf

max_forks_repo_name: vsart/apf


!****h* root/fortran/test/tf.f90
!
! NAME
!  tf.f90
!
! FUNCTION
!  Contains subroutines which are needed in all the hdf5 fortran tests
!
! COPYRIGHT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the COPYING file, which can be found at the root of the source code       *
!   distribution tree, or in https://support.hdfgroup.org/ftp/HDF5/releases.  *
!   If you do not have access to either file, you may request a copy from     *
!   help@hdfgroup.org.                                                        *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! CONTAINS SUBROUTINES
!  write_test_status, check, verify, verifyLogical, verifyString, h5_fixname_f,
!  h5_cleanup_f, h5_exit_f, h5_env_nocleanup_f,dreal_eqv
!
!*****

#include "H5config_f.inc"

MODULE TH5_MISC

  USE, INTRINSIC :: ISO_C_BINDING

  IMPLICIT NONE

  INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(5)  ! This should map to REAL*4 on most modern processors
  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(10) ! This should map to REAL*8 on most modern processors

  ! generic compound datatype
  TYPE :: comp_datatype
    SEQUENCE
    REAL :: a
    INTEGER :: x
    DOUBLE PRECISION :: y
    CHARACTER(KIND=C_CHAR) :: z
  END TYPE comp_datatype

  PUBLIC :: H5_SIZEOF
  INTERFACE H5_SIZEOF
     MODULE PROCEDURE H5_SIZEOF_CMPD
     MODULE PROCEDURE H5_SIZEOF_CHR
     MODULE PROCEDURE H5_SIZEOF_I
     MODULE PROCEDURE H5_SIZEOF_SP,H5_SIZEOF_DP
  END INTERFACE

CONTAINS

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: write_test_status
!DEC$endif
  SUBROUTINE write_test_status( test_result, test_title, total_error)

    ! Writes the results of the tests

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: test_result  ! negative,  --skip --
                                        ! 0       ,   passed
                                        ! positive,   failed

    CHARACTER(LEN=*), INTENT(IN) :: test_title ! Short description of test
    INTEGER, INTENT(INOUT) :: total_error ! Accumulated error

! Controls the output style for reporting test results

    CHARACTER(LEN=8) :: error_string
    CHARACTER(LEN=8), PARAMETER :: success = ' PASSED '
    CHARACTER(LEN=8), PARAMETER :: failure = '*FAILED*'
    CHARACTER(LEN=8), PARAMETER :: skip    = '--SKIP--'

    error_string = failure
    IF (test_result ==  0) THEN
       error_string = success
    ELSE IF (test_result == -1) THEN
       error_string = skip
    ENDIF

    WRITE(*, fmt = '(A, T80, A)') test_title, error_string

    IF(test_result.GT.0) total_error = total_error + test_result

  END SUBROUTINE write_test_status

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: check
!DEC$endif
  SUBROUTINE check(string,error,total_error)
    CHARACTER(LEN=*) :: string
    INTEGER :: error, total_error
    IF (error .LT. 0) THEN
       total_error=total_error+1
       WRITE(*,*) string, " FAILED"
    ENDIF
    RETURN
  END SUBROUTINE check

!----------------------------------------------------------------------
! Name:		h5_fixname_f
!
! Purpose:	Create a file name from the a file base name.
!               It is a fortran counterpart for the h5_fixname in ../../test/h5test.c
!
! Inputs:
!		base_name	- base name of the file
!               fapl		- file access property list
! Outputs:
!		full_name	- full file name
!		hdferr:		- error code
!				 	Success:  0
!				 	Failure: -1
!
! Programmer:	Elena Pourmal
!		September 13, 2002
!
!
!----------------------------------------------------------------------
  SUBROUTINE h5_fixname_f(base_name, full_name, fapl, hdferr)
!
!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_fixname_f
!DEC$endif
    USE H5GLOBAL
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: base_name   ! base name
    CHARACTER(LEN=*), INTENT(IN) :: full_name   ! full name
    INTEGER, INTENT(OUT) :: hdferr         ! Error code
    INTEGER(HID_T), INTENT(IN) :: fapl ! file access property list

    INTEGER(SIZE_T) :: base_namelen ! Length of the base name character string
    INTEGER(SIZE_T) :: full_namelen ! Length of the full name character string
!            INTEGER(HID_T) :: fapl_default

    INTERFACE
       INTEGER FUNCTION h5_fixname_c(base_name, base_namelen, fapl, &
            full_name, full_namelen)
         USE H5GLOBAL
         !DEC$ IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ ATTRIBUTES C,reference,decorate,alias:'H5_FIXNAME_C':: h5_fixname_c
         !DEC$ ENDIF
         !DEC$ATTRIBUTES reference :: base_name
         !DEC$ATTRIBUTES reference :: full_name
         CHARACTER(LEN=*), INTENT(IN) :: base_name
         INTEGER(SIZE_T) :: base_namelen
         INTEGER(HID_T), INTENT(IN) :: fapl
         CHARACTER(LEN=*), INTENT(IN) :: full_name
         INTEGER(SIZE_T) :: full_namelen
       END FUNCTION h5_fixname_c
    END INTERFACE

    base_namelen = LEN(base_name)
    full_namelen = LEN(full_name)
    hdferr = h5_fixname_c(base_name, base_namelen, fapl, &
         full_name, full_namelen)

  END SUBROUTINE h5_fixname_f

!----------------------------------------------------------------------
! Name:		h5_cleanup_f
!
! Purpose:	Cleanups tests files
!               It is a fortran counterpart for the h5_cleanup in ../../test/h5test.c
!
! Inputs:
!		base_name	- base name of the file
!               fapl		- file access property list
! Outputs:
!		hdferr:		- error code
!				 	Success:  0
!				 	Failure: -1
!
! Programmer:	Elena Pourmal
!		September 19, 2002
!
!
!----------------------------------------------------------------------
  SUBROUTINE h5_cleanup_f(base_name, fapl, hdferr)
!
!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_cleanup_f
!DEC$endif
    USE H5GLOBAL
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: base_name   ! base name
    INTEGER, INTENT(OUT) :: hdferr         ! Error code
    INTEGER(HID_T), INTENT(IN) :: fapl ! file access property list

    INTEGER(SIZE_T) :: base_namelen ! Length of the base name character string

    INTERFACE
       INTEGER FUNCTION h5_cleanup_c(base_name, base_namelen, fapl)
         USE H5GLOBAL
         !DEC$ IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ ATTRIBUTES C,reference,decorate,alias:'H5_CLEANUP_C':: h5_cleanup_c
         !DEC$ ENDIF
         !DEC$ATTRIBUTES reference :: base_name
         CHARACTER(LEN=*), INTENT(IN) :: base_name
         INTEGER(SIZE_T) :: base_namelen
         INTEGER(HID_T), INTENT(IN) :: fapl
       END FUNCTION h5_cleanup_c
    END INTERFACE

    base_namelen = LEN(base_name)
    hdferr = h5_cleanup_c(base_name, base_namelen, fapl)

  END SUBROUTINE h5_cleanup_f

!----------------------------------------------------------------------
! Name:		h5_exit_f
!
! Purpose:	Exit application
!               It is a fortran counterpart for the standard C 'exit()' routine
! 		Be careful not to overflow the exit value range since
! 		UNIX supports a very small range such as 1 byte.
!		Therefore, exit(256) may end up as exit(0).
!
! Inputs:
!		status	- Status to return from application
!
! Outputs:
!		none
!
! Programmer:	Quincey Koziol
!		December 14, 2004
!
!
!----------------------------------------------------------------------
  SUBROUTINE h5_exit_f(status)
!
!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_exit_f
!DEC$endif
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: status         ! Return code

    INTERFACE
       SUBROUTINE h5_exit_c(status)
         !DEC$ IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ ATTRIBUTES C,reference,decorate,alias:'H5_EXIT_C':: h5_exit_c
         !DEC$ ENDIF
         INTEGER, INTENT(IN) :: status
       END SUBROUTINE h5_exit_c
    END INTERFACE

    CALL h5_exit_c(status)

  END SUBROUTINE h5_exit_f

!----------------------------------------------------------------------
! Name:		h5_env_nocleanup_f
!
! Purpose:	Uses the HDF5_NOCLEANUP environment variable in Fortran
!               tests to determine if the output files should be removed
!
! Inputs:
!
! Outputs:      HDF5_NOCLEANUP:  .true. - don't remove test files
!		                .false. - remove test files
!
! Programmer:	M.S. Breitenfeld
!               September 30, 2008
!
!----------------------------------------------------------------------
  SUBROUTINE h5_env_nocleanup_f(HDF5_NOCLEANUP)
!
!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_env_nocleanup_f
!DEC$endif
    IMPLICIT NONE
    LOGICAL, INTENT(OUT) :: HDF5_NOCLEANUP ! Return code
    INTEGER :: status

    INTERFACE
       SUBROUTINE h5_env_nocleanup_c(status)
         !DEC$ IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ ATTRIBUTES C,reference,decorate,alias:'H5_ENV_NOCLEANUP_C':: h5_env_nocleanup_c
         !DEC$ ENDIF
         INTEGER :: status
       END SUBROUTINE h5_env_nocleanup_c
    END INTERFACE

    CALL h5_env_nocleanup_c(status)

    HDF5_NOCLEANUP = .FALSE.
    IF(status.EQ.1) HDF5_NOCLEANUP = .TRUE.

  END SUBROUTINE h5_env_nocleanup_f

! ---------------------------------------------------------------------------------------------------
! H5_SIZEOF routines
!
! NOTES
!   (1) The Sun/Oracle compiler has the following restrictions on the SIZEOF intrinsic function:
!
!     "The SIZEOF intrinsic cannot be applied to arrays of an assumed size, characters of a
!      length that is passed, or subroutine calls or names. SIZEOF returns default INTEGER*4 data.
!      If compiling for a 64-bit environment, the compiler will issue a warning if the result overflows
!      the INTEGER*4 data range. To use SIZEOF in a 64-bit environment with arrays larger
!      than the INTEGER*4 limit (2 Gbytes), the SIZEOF function and
!      the variables receiving the result must be declared INTEGER*8."
!
!    Thus, we can not overload the H5_SIZEOF function to handle arrays (as used in tH5P_F03.f90), or
!    characters that do not have a set length (as used in tH5P_F03.f90), sigh...
!
!   (2) F08+TS29113 requires C interoperable variable as argument for C_SIZEOF.
!
!   (3) Unfortunately we need to wrap the C_SIZEOF/STORAGE_SIZE functions to handle different
!       data types from the various tests.
!
! ---------------------------------------------------------------------------------------------------

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_sizeof_cmpd
!DEC$endif
  INTEGER(C_SIZE_T) FUNCTION H5_SIZEOF_CMPD(a)
    IMPLICIT NONE
    TYPE(comp_datatype), INTENT(in) :: a

#ifdef H5_FORTRAN_HAVE_STORAGE_SIZE
    H5_SIZEOF_CMPD = storage_size(a, c_size_t)/storage_size(c_char_'a',c_size_t)
#else
    H5_SIZEOF_CMPD = SIZEOF(a)
#endif

  END FUNCTION H5_SIZEOF_CMPD

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_sizeof_chr
!DEC$endif
  INTEGER(C_SIZE_T) FUNCTION H5_SIZEOF_CHR(a)
    IMPLICIT NONE
    CHARACTER(LEN=1), INTENT(in) :: a

#ifdef H5_FORTRAN_HAVE_STORAGE_SIZE
    H5_SIZEOF_CHR = storage_size(a, c_size_t)/storage_size(c_char_'a',c_size_t)
#else
    H5_SIZEOF_CHR = SIZEOF(a)
#endif

  END FUNCTION H5_SIZEOF_CHR

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_sizeof_i
!DEC$endif
  INTEGER(C_SIZE_T) FUNCTION H5_SIZEOF_I(a)
    IMPLICIT NONE
    INTEGER, INTENT(in):: a

#ifdef H5_FORTRAN_HAVE_STORAGE_SIZE
    H5_SIZEOF_I = storage_size(a, c_size_t)/storage_size(c_char_'a',c_size_t)
#else
    H5_SIZEOF_I = SIZEOF(a)
#endif

  END FUNCTION H5_SIZEOF_I

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_sizeof_sp
!DEC$endif
  INTEGER(C_SIZE_T) FUNCTION H5_SIZEOF_SP(a)
    IMPLICIT NONE
    REAL(sp), INTENT(in):: a

#ifdef H5_FORTRAN_HAVE_STORAGE_SIZE
    H5_SIZEOF_SP = storage_size(a, c_size_t)/storage_size(c_char_'a',c_size_t)
#else
    H5_SIZEOF_SP = SIZEOF(a)
#endif

  END FUNCTION H5_SIZEOF_SP

!This definition is needed for Windows DLLs
!DEC$if defined(BUILD_HDF5_TEST_DLL)
!DEC$attributes dllexport :: h5_sizeof_dp
!DEC$endif
  INTEGER(C_SIZE_T) FUNCTION H5_SIZEOF_DP(a)
    IMPLICIT NONE
    REAL(dp), INTENT(in):: a

#ifdef H5_FORTRAN_HAVE_STORAGE_SIZE
    H5_SIZEOF_DP = storage_size(a, c_size_t)/storage_size(c_char_'a',c_size_t)
#else
    H5_SIZEOF_DP = SIZEOF(a)
#endif

  END FUNCTION H5_SIZEOF_DP

END MODULE TH5_MISC


text:
      SUBROUTINE PATHK
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
************************************************************************
*
*  		Written by Manyin Yi, Aug 1989.
*       Restartable reaction_path calulation.
*	The number of path_step and step value are read in through
*       keyword POINT and STEP.
*	The reaction profile is archived.
*
************************************************************************
      COMMON /GEOM  / GEO(3,NUMATM)
      COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      COMMON /GRADNT/ GRAD(MAXPAR),GNORM
      COMMON /GRAVEC/ COSINE
      COMMON /PATH  / LATOM, LPARAM, REACT(200)
      COMMON /PPARAM/ CURRT
      COMMON /KLOOP / KLOOP
      COMMON /PROFIL/ PROFIL
      COMMON /KEYWRD/ KEYWRD
      DIMENSION GD(MAXPAR),XLAST(MAXPAR),MDFP(20),XDFP(20)
      DIMENSION PROFIL(200)
      CHARACTER KEYWRD*241, GETNAM*80
      STEP=READA(KEYWRD,INDEX(KEYWRD,'STEP')+5)
      NPTS=READA(KEYWRD,INDEX(KEYWRD,'POINT')+6)
C
C  THE SMALLEST VALUE IN THE PATH IS
C      REACT(1) DEGREE OR GEO(LPARAM,LATOM) RADIANS
C
      DEGREE=180.D0/3.14159265359D0
      IF(LPARAM.NE.1)STEP=STEP/DEGREE
C
C  NOW TO SWEEP THROUGH THE PATH
C
      IF(LPARAM.NE.1) THEN
         C1=DEGREE
      ELSE
         C1=1.D0
      ENDIF
C
      KLOOP=1
      CPUTOT=0.0
      CURRT=GEO(LPARAM,LATOM)
      PROFIL(1)=0.D0
      IF (INDEX(KEYWRD,'RESTART').NE.0) THEN
         MDFP(9)=0
         CALL DFPSAV(CPUTOT,XPARAM,GD,XLAST,ESCF,MDFP,XDFP)
         WRITE(6,'(//10X,'' RESTARTING AT POINT '',I3)') KLOOP
      ENDIF
C
      GEO(LPARAM,LATOM)=CURRT
      DO 10 ILOOP=KLOOP,NPTS
         CPU1=SECOND()
         CURRT=GEO(LPARAM,LATOM)
         CALL FLEPO(XPARAM, NVAR, ESCF)
         KLOOP=KLOOP+1
         CPU2=SECOND()
         CPU3=CPU2-CPU1
         CPUTOT=CPUTOT+CPU3
         PROFIL(ILOOP)=ESCF
         WRITE(6,'(/''          VARIABLE        FUNCTION'')')
         WRITE(6,'('' :'',F16.5,F16.6)')GEO(LPARAM,LATOM)*C1,ESCF
         CALL GEOUT(6)
         GEO(LPARAM,LATOM)=GEO(LPARAM,LATOM)+STEP
   10 CONTINUE
      DO 20 I=2,NPTS
   20 REACT(I)=REACT(I-1)+STEP*C1
      WRITE(6,'(/16X,''POINTS ON REACTION PATH '',
     1          /16X,''AND CORRESPONDING HEATS'',//)')
      OPEN(UNIT=12,FILE=GETNAM('FOR012'),STATUS='UNKNOWN')
      WRITE(12,30)
      CALL WRTTXT(12)
   30 FORMAT(' ARCHIVE FILE FOR PATH CALCULATION'/
     1 'A PROFIL OF COORDINATES - HEATS'/)
      WRITE(12,'(/'' TOTAL CPU TIME IN FLEPO : '',F10.3/)') CPUTOT
C
      L=NPTS/8
      M=NPTS - L*8
      IF (L.LT.1) GO TO 50
      DO 40 K=0,L-1
         WRITE(6,'(8F7.2)') (REACT(I),I=K*8+1,K*8+8)
         WRITE(6,'(8F7.2,/)') (PROFIL(I),I=K*8+1,K*8+8)
         WRITE(12,'(8F7.2)') (REACT(I),I=K*8+1,K*8+8)
   40 WRITE(12,'(8F7.2,/)') (PROFIL(I),I=K*8+1,K*8+8)
   50 IF (M.GT.0) THEN
         WRITE(6,'(8F7.2)') (REACT(I),I=L*8+1,L*8+M)
         WRITE(6,'(8F7.2,/)') (PROFIL(I),I=L*8+1,L*8+M)
         WRITE(12,'(8F7.2)') (REACT(I),I=L*8+1,L*8+M)
         WRITE(12,'(8F7.2,/)') (PROFIL(I),I=L*8+1,L*8+M)
      ENDIF
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
C.C. Cooper & Sons Insurance operated out of the Varsity Theater building for 55 years before moving to their current location. They are a family owned and operated insurance agency who takes pride in their work and their clients.

 Services
   Homeowners Insurance
   Renters Insurance
   Flood Insurance
   Auto Insurance
   Commercial Auto Insurance
   Mobile Home Insurance
   Recreational Vehicle Insurance
   Boat Insurance
   General Liability Insurance
   One Day Event Insurance

20100804 16:40:06 nbsp C. C. Cooper & Sons Insurance is the place to go!  They are always nice and willing to help you out with any questions, this is the best place to go for insurance hands down Users/dnm12

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
       subroutine eqblok ( t, n, kpm,  work1, work2,
     *		       bloks, lenblk, integs, nbloks,  b )
c  from  * a practical guide to splines *  by c. de boor
calls putit(difequ,bsplvd(bsplvb))
c  to be called in  c o l l o c
c
c******  i n p u t  ******
c  t   the knot sequence, of length n+kpm
c  n   the dimension of the approximating spline space, i.e., the order
c      of the linear system to be constructed.
c  kpm = k+m, the order of the approximating spline
c  lenblk   the maximum length of the array  bloks  as allowed by the
c	    dimension statement in  colloc .
c
c******  w o r k   a r e a s  ******
c  work1    used in  putit, of size (kpm,kpm)
c  work2    used in  putit, of size (kpm,m+1)
c
c******  o u t p u t  ******
c  bloks    the coefficient matrix of the linear system, stored in al-
c	    most block diagonal form, of size
c	       kpm*sum(integs(1,i) , i=1,...,nbloks)
c  integs   an integer array, of size (3,nbloks), describing the block
c	    structure.
c	    integs(1,i)  =  number of rows in block  i
c	    integs(2,i)  =  number of columns in block	i
c	    integs(3,i)  =  number of elimination steps which can be
c			carried out in block  i  before pivoting might
c			bring in an equation from the next block.
c  nbloks   number of blocks, equals number of polynomial pieces
c  b   the right side of the linear system, stored corresponding to the
c      almost block diagonal form, of size sum(integs(1,i) , i=1,...,
c      nbloks).
c
c******  m e t h o d  ******
c  each breakpoint interval gives rise to a block in the linear system.
c  this block is determined by the  k  colloc.equations in the interval
c  with the side conditions (if any) in the interval interspersed ap-
c  propriately, and involves the  kpm  b-splines having the interval in
c  their support. correspondingly, such a block has  nrow = k + isidel
c  rows, with  isidel = number of side conditions in this and the prev-
c  ious intervals, and	ncol = kpm  columns.
c     further, because the interior knots have multiplicity  k, we can
c  carry out (in slvblk)  k  elimination steps in a block before pivot-
c  ing might involve an equation from the next block. in the last block,
c  of course, all kpm elimination steps will be carried out (in slvblk).
c
c  see the detailed comments in the solveblok package for further in-
c  formation about the almost block diagonal form used here.
      integer integs(3,1),kpm,lenblk,n,nbloks,	 i,index,indexb,iside
     *					    ,isidel,itermx,k,left,m,nrow
      real b(1),bloks(1),t(1),work1(1),work2(1),   rho,xside
      common /side/ m, iside, xside(10)
      common /other/ itermx,k,rho(19)
      index = 1
      indexb = 1
      i = 0
      iside = 1
      do 20 left=kpm,n,k
	 i = i+1
c	 determine integs(.,i)
	 integs(2,i) = kpm
	 if (left .lt. n)		go to 14
	 integs(3,i) = kpm
	 isidel = m
					go to 16
   14	 integs(3,i) = k
c	 at this point,  iside-1  gives the number of side conditions
c	 incorporated so far. adding to this the side conditions in the
c	 current interval gives the number  isidel .
	 isidel = iside-1
   15	 if (isidel .eq. m)		go to 16
	 if (xside(isidel+1) .ge. t(left+1))
     *					go to 16
	 isidel = isidel+1
					go to 15
   16	 nrow = k + isidel
	 integs(1,i) = nrow
c	 the detailed equations for this block are generated and put
c	 together in  p u t i t .
	 if (lenblk .lt. index+nrow*kpm-1)go to 999
	 call putit(t,kpm,left,work1,work2,bloks(index),nrow,b(indexb))
	 index = index + nrow*kpm
   20	 indexb = indexb + nrow
      nbloks = i
					return
  999 print 699,lenblk
  699 format(11h **********/23h the assigned dimension,i5
     *	      ,38h for	bloks  in  colloc  is too small.)
					stop
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


module fields_arrays
  use constants, only: run_name_size
  implicit none

  private

  public :: phi, apar, bpar, phinew, aparnew, bparnew
  public :: aparold
  public :: phistar_old, phistar_new ! NDCTESTmichaelnew
  public :: gf_phi, gf_apar, gf_bpar, gf_phinew, gf_aparnew, gf_bparnew
  public :: apar_ext, time_field, response_file

  !Main fields
  complex, dimension (:,:,:), allocatable :: phi,    apar,    bpar
  complex, dimension (:,:,:), allocatable :: phinew, aparnew, bparnew
  complex, dimension (:,:,:), allocatable :: aparold
  complex, dimension (:,:,:), allocatable :: phistar_old, phistar_new ! NDCTESTmichaelnew
  complex, dimension (:,:,:), allocatable :: gf_phi, gf_apar, gf_bpar
  complex, dimension (:,:,:), allocatable :: gf_phinew, gf_aparnew, gf_bparnew
  !For antenna
  complex, dimension (:,:,:), allocatable :: apar_ext
  ! (-ntgrid:ntgrid,ntheta0,naky) replicated

  !Timing data
  real :: time_field(2)=0.

  !For response data
  character(run_name_size) :: response_file
end module fields_arrays


text:
c     Continuum calculations according to Clough et al (1994)
      subroutine avhrr51(z0, p0, t0, u0, ux, nlev,
     *co2, n2o, f11, f12, ch4,
     *od, wght, ntau)
      implicit real*8 (v)
      INCLUDE "avhrr.inc"
      real co2, n2o, ch4, f11, f12
      real z0(mm)
      real od(10,mm)
      real wght(10)
      integer ntau,nlev
      real wi(8),xi(8),p0(mq),p(mm),t0(mq),t(mm),bnu(mm),u0(mq),u(mm),
     *ucfc(mm),tau(6,mm),taucfc(6,mm),tauf(6,6,mm),delw,x(mm,8),y(mm,8),
     *sx(6,6,mm),sy(6,6,mm),f(6),fcfc(6),fac(mm),tsx(mq),tsy(mq)
      real dvabs,absrb(2030),pres(mm),war(mm),wn2(mm),wo2(mm),wco2(mm),
     *w(mm)
      data absrb /2030*0.0/
      data xi /.01986,.10167,.23723,.40828,
     *         .59172,.76277,.89833,.98014/
      data wi /.05062,.11119,.15685,.18134,
     *         .18134,.15685,.11119,.05062/
      tb=0
      nlyr = nlev-1
      if(nlyr.GT.mm)then
        write (*,*) 'error: increase mm to at least', nlyr
        write (*,*) '       in avhrr.inc'
        return
      end if
      if(nlev.GT.mq)then
        write (*,*) 'error: increase mq to at least', nlev
        write (*,*) '       in avhrr.inc'
        return
      end if
      omega=840.0
      delw=80.0
      v1abs=800.0
      v2abs=880.0
      dvabs=5.0
      nptabs=17
      mlv=nlyr
      h=3.0
      do 5 m=1,nlyr
        t(m)=(t0(m)+t0(m+1))/2.0
        p(m)=(p0(m)+p0(m+1))/2.0
        call flux_51(bnu(m),omega,t(m))
        fac(m)=z0(m+1)-z0(m)
c       u(m)=fac(m)*(u0(m)+u0(m+1))/2.0*6.023e+23/2.687e+19/18.01534
        u(m)=fac(m)*(u0(m)+u0(m+1))/2.0
        w(m)=u(m)*6.023e+23/18.01534
        ucfc(m)=p(m)*1.01325e+06/(1.3805e-16*t(m))*f11*1e-6*fac(m)*
     *          1.0e+05/6.023e+23*137.3685
        if(z0(m).ge.24.999)
     *ucfc(m)=ucfc(m)*exp((-7.5-(z0(m)-20.))/h)
        if(z0(m).lt.24.999.and.z0(m).ge.14.999)
     *ucfc(m)=ucfc(m)*exp((14.5-z0(m))/h)
        wn2(m)=p(m)*1.01325e+06/(1.3805e-16*t(m))*0.78084*fac(m)*
     *         1.0e+05
        wo2(m)=p(m)*1.01325e+06/(1.3805e-16*t(m))*0.20948*fac(m)*
     *         1.0e+05
        war(m)=p(m)*1.01325e+06/(1.3805e-16*t(m))*0.00934*fac(m)*
     *         1.0e+05
        wco2(m)=p(m)*1.01325e+06/(1.3805e-16*t(m))*co2*1e-6*fac(m)*
     *          1.0e+05
        pres(m)=p(m)*1013.25
 5    continue
      do 7 m=1,nlev
        tsx(m)=0.0
        tsy(m)=0.0
 7    continue
      call flux_51(back,omega,t0(1))
      call ck_51(ni,mlv,u,delw,f,p,t,tau)
      call ckcfc_51(ncfc,mlv,ucfc,delw,fcfc,p,t,taucfc)
      do 10 m=1,nlyr
        do k=1,nptabs
          absrb(k)=0.0
        end do
        call clough_51(pres(m),t(m),war(m),wn2(m),wo2(m),wco2(m),w(m),
     *       v1abs,v2abs,dvabs,nptabs,absrb)
        twin=0.5*(absrb(1)+absrb(17))
        do k=2,nptabs-1
          twin=twin+absrb(k)
        end do
        twin=twin/16.0
        do i=1,ni
          do ia=1,ncfc
c           tauf(i,ia,m)=twin
c           tauf(i,ia,m)=tau(i,m)+twin
c           tauf(i,ia,m)=taucfc(ia,m)
            tauf(i,ia,m)=tau(i,m)+twin+taucfc(ia,m)
          end do
        end do
 10   continue
      do 90 i=1,ni
       do 85 ia=1,ncfc
        do 30 m=1,nlyr
          sx(i,ia,m)=0.0
          sy(i,ia,m)=0.0
          do 20 j=1,8
           if(m.eq.1)then
            x(m,j)=back*exp(-tauf(i,ia,m)/xi(j))+
     *      bnu(m)*(1.0-exp(-tauf(i,ia,m)/xi(j)))
           else
            x(m,j)=x(m-1,j)*exp(-tauf(i,ia,m)/xi(j))+
     *      bnu(m)*(1.0-exp(-tauf(i,ia,m)/xi(j)))
           end if
 20       continue
 30     continue
        do 50 m=nlyr,1,-1
          do 40 j=1,8
           if(m.eq.nlyr)then
            y(m,j)=bnu(m)*(1.0-exp(-tauf(i,ia,m)/xi(j)))
           else
            y(m,j)=y(m+1,j)*exp(-tauf(i,ia,m)/xi(j))+
     *      bnu(m)*(1.0-exp(-tauf(i,ia,m)/xi(j)))
           end if
 40       continue
 50     continue
        do 70 m=1,nlyr
          do 60 j=1,8
            sx(i,ia,m)=sx(i,ia,m)+2.0*x(m,j)*wi(j)*xi(j)
            sy(i,ia,m)=sy(i,ia,m)+2.0*y(m,j)*wi(j)*xi(j)
 60       continue
 70     continue
        do 80 m=1,nlyr
           tsx(m+1)=tsx(m+1)+sx(i,ia,m)*delw*f(i)*fcfc(ia)
           tsy(m)=tsy(m)+sy(i,ia,m)*delw*f(i)*fcfc(ia)
 80     continue
        tb=tb+back*delw*f(i)*fcfc(ia)
        tsx(1)=tb
 85    continue
 90   continue
      do 95 m=1,nlev
cbm         write(6,100)m-1,tsx(m),tb-tsx(m),tsy(m)
 95   continue
      do i1=1,5
       do i2=1,2
        do m=1,nlyr
         od(1+(i2-1)+(i1-1)*2,m) =
     *tauf(i1,i2,m)
         wght(1+(i2-1)+(i1-1)*2) =
     *f(i1) * fcfc(i2)
        enddo
       enddo
      enddo
      ntau=10
 100  format(i3,3f15.7)
      end
c *********************************************************************
      subroutine flux_51(bnu,omega,t)
      INCLUDE "avhrr.inc"
      bnu=0.0
      z=3.14159
      bnu=z/1000.0*((2.*6.625e-27)*(2.997925e+10)**2
     *    *omega**3)/(exp(min(88.0,1.4388*omega/t))-1.0)
      return
      end
c *********************************************************************
      subroutine ck_51(ni,mlv,u,delw,f,p,t,tau)
      INCLUDE "avhrr.inc"
      real f(6),k(mm),tau(6,mm),delw,p(mm),t(mm),u(mm),coefk(6,3,19),
     *     stp(19),fkg(6,mm),pmb(mm)
      ni=5
      f(1)=0.771954
      f(2)=0.096366
      f(3)=0.084180
      f(4)=0.034118
      f(5)=0.013382
      k(1)=0.00183332
      do i=2,ni
        k(i)=6.0*k(i-1)
      end do
      data ( (coefk(1,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0000312,
     *0.000385,0.00179,0.00873,0.0439,0.1173,0.2507,0.4602,0.7777,
     *1.3008,2.1112,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
     *0.000e+00,0.000e+00,1.593e-05,1.684e-05,8.016e-05,3.721e-04,
     *1.368e-03,3.312e-03,6.749e-03,1.227e-02,2.046e-02,3.343e-02,
     *5.354e-02,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
     *0.000e+00,0.000e+00,2.256e-08,2.197e-07,1.096e-06,4.947e-06,
     *1.429e-05,3.257e-05,6.305e-05,1.151e-04,1.871e-04,3.031e-04,
     *4.766e-04/
      data ( (coefk(2,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.000133,0.000923,
     *0.00487,0.00248,0.0727,0.1313,0.2139,0.3371,0.5231,0.8021,1.2575,
     *1.9333,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
     *0.000e+00,6.381e-06,3.835e-05,1.914e-04,8.819e-04,2.263e-03,
     *3.927e-03,6.408e-03,9.993e-03,1.518e-02,2.252e-02,3.447e-02,
     *5.203e-02,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
     *0.000e+00,8.859e-08,4.956e-07,2.353e-06,9.859e-06,2.330e-05,
     *4.071e-05,6.659e-05,1.043e-04,1.560e-04,2.234e-04,3.312e-04,
     *4.915e-04/
      data ( (coefk(3,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.0,0.0,0.0,0.0,0.00000543,0.000120,0.000707,0.00383,0.0141,
     *0.0342,0.0627,0.1003,0.1547,0.2357,0.3587,0.5496,0.8192,1.2191,
     *1.7196,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,1.846e-07,4.780e-06,
     *2.770e-05,1.519e-04,5.093e-04,1.137e-03,2.001e-03,3.169e-03,
     *4.808e-03,7.223e-03,1.089e-02,1.659e-02,2.463e-02,3.677e-02,
     *5.074e-02,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,1.617e-09,5.863e-08,
     *3.469e-07,1.865e-06,5.869e-06,1.255e-05,2.169e-05,3.438e-05,
     *5.125e-05,7.591e-05,1.135e-04,1.711e-04,2.527e-04,3.823e-04,
     *5.181e-04/
      data ( (coefk(4,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.0000320,0.0000391,0.000143,0.000775,0.00300,0.00793,0.0164,
     *0.0282,0.0441,0.0658,0.0980,0.1469,0.2216,0.3319,0.4871,0.6852,
     *0.8953,1.1005,1.2216,
     *5.089e-06,5.999e-06,8.501e-06,4.752e-05,1.522e-04,3.379e-04,
     *6.123e-04,9.700e-04,1.434e-03,2.098e-03,3.057e-03,4.523e-03,
     *6.710e-03,1.001e-02,1.470e-02,2.046e-02,2.624e-02,3.193e-02,
     *3.556e-02,
     *1.078e-07,1.262e-07,1.343e-07,7.650e-07,2.190e-06,4.466e-06,
     *7.469e-06,1.125e-05,1.598e-05,2.304e-05,3.293e-05,4.689e-05,
     *6.966e-05,1.029e-04,1.518e-04,2.096e-04,2.609e-04,3.115e-04,
     *3.470e-04/
      data ( (coefk(5,jt,jp), jp = 1, 19), jt = 1, 3)/
     *1.6376,1.6361,1.6353,1.6354,1.6349,1.6337,1.6332,1.6308,
     *1.6297,1.6279,1.6230,1.6073,1.5710,1.5034,1.4069,1.2739,1.0877,
     *0.9120,0.7257,
     *5.102e-02,5.099e-02,5.093e-02,5.091e-02,5.086e-02,5.076e-02,
     *5.068e-02,5.047e-02,5.024e-02,4.987e-02,4.937e-02,4.858e-02,
     *4.729e-02,4.524e-02,4.225e-02,3.844e-02,3.295e-02,2.775e-02,
     *2.245e-02,
     *5.617e-04,5.614e-04,5.603e-04,5.599e-04,5.588e-04,5.570e-04,
     *5.550e-04,5.517e-04,5.475e-04,5.412e-04,5.333e-04,5.226e-04,
     *5.074e-04,4.866e-04,4.544e-04,4.171e-04,3.615e-04,3.101e-04,
     *2.585e-04/
	data stp / 0.251, 0.398, 0.631, 1.000, 1.58, 2.51,
     *	             3.98, 6.31, 10.0, 15.8, 25.1, 39.8, 63.1,
     *	             100.0, 158.0, 251.0, 398.0, 631.0, 1000.0/
      do i=1,ni
      do m=1,mlv
        ml=1
        pmb(m)=p(m)*1013.25
        if(pmb(m).lt.stp(1))then
         x1=coefk(i,1,1)+coefk(i,2,1)*(t(m)-250.0)
     *   +coefk(i,3,1)*(t(m)-250.0)**2
         fkg(i,m)=x1*pmb(m)/stp(1)
        else if (pmb(m).gt.stp(19)) then
         x1=coefk(i,1,18)+coefk(i,2,18)*(t(m)-250.0)
     *   +coefk(i,3,18)*(t(m)-250.0)**2
         x2=coefk(i,1,19)+coefk(i,2,19)*(t(m)-250.0)
     *   +coefk(i,3,19)*(t(m)-250.0)**2
         fkg(i,m)=x1+(x2-x1)/(stp(19)-stp(18))
     *  *(pmb(m)-stp(18))
        else
         do while(pmb(m).ge.stp(ml))
           ml=ml+1
         end do
         x1=coefk(i,1,ml-1)+coefk(i,2,ml-1)*(t(m)-250.0)
     *   + coefk(i,3,ml-1)*(t(m)-250.0)**2
         x2=coefk(i,1,ml)+coefk(i,2,ml)*(t(m)-250.0)
     *   + coefk(i,3,ml)*(t(m)-250.0)**2
         fkg(i,m)=x1+(x2-x1)/(stp(ml)-stp(ml-1))
     *   *(pmb(m)-stp(ml-1))
        end if
      end do
      end do
      do i=1,ni
        do m=1,mlv
          tau(i,m)=k(i)*u(m)*fkg(i,m)
c         write(6,*)m,i,u(m),k(i),fkg(i,m),tau(i,m)
        end do
      end do
      return
      end
c *********************************************************************
      subroutine ckcfc_51(ni,mlv,u,delw,f,p,t,tau)
      INCLUDE "avhrr.inc"
      real f(6),k(mm),tau(6,mm),delw,p(mm),t(mm),u(mm)
      ni=2
c     f(1)=ratio of spectral interval used vs. cfc spectral interval
c     correction for spectral interval offset
c      cor=1.000000
      cor=0.945681
      f(1)=0.625*cor
      f(2)=1.0-f(1)
      do m=1,mlv
        k(1)=5611.777*(1.0000+1.2661e-03*(t(m)-250.0)+
     *       3.5571e-06*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
        tau(2,m)=0.0
c       write(6,*)m,u(m),k(1),tau(1,m)
      end do
      return
      end
c**********************************************************************
      subroutine clough_51(pres,temp,wa,wn2,wo2,wco2,wh2o,v1abs,v2abs,
     *dvabs,nptabs,absrb)
      INCLUDE "avhrr.inc"
      implicit real*8 (v)
      dimension absrb(*)
      dimension wk(60)
      data xlosmt/2.68675e+19/
      data wk/60*0/
c
c
c   this program calculates the continuum optical depth
c         for an homogeneous layer
c
c   the following quantities must be specified:
c
c          pressure                   pave (mb)
c
c          temperature                tave ( k)
c
c          column amount
c            nitrogen                 wn2    (molec/cm**2)
c            oxygen                   wk(7)  (molec/cm**2)
c            carbon dioxide           wk(2)  (molec/cm**2)
c            water vapor              wk(1)  (molec/cm**2)
c
c          number of molecules        nmol
c
c          beginning wavenumber       v1abs (cm-1)
c
c          ending wavenumber          v2abs (cm-1)
c
c          sampling interval          dvabs (cm-1)
c
c          number of values           nptabs
c
c
c   the results are in array absorb
c
c   note that for an atmospheric layer:
c
c            wtot   = xlosmt * (pave/1013) * (273/tave) * (path length)
c
c            wbroad = the column amount for all species not explicitly provided
c
c            wk(m)  = (volume mixing ratio) * (column of dry air)
c
c
c
c      note that continua are included for
c
c             m=1      water vapor       air & self         0 - 20,000 cm-1
c             wn2      nitrogen          air             2020 -   2800 cm-1
c
c   the following is an example for a one km path (see cntnout for results)
c
      pave = pres
      tave = temp
      wk(7) = wo2
      wk(2) = wco2
      wk(1) =wh2o
      wbroad=wn2+wa
      nmol = 7
c
      jrad=1
c
      call contnm_51(jrad,pave,tave,wk,wbroad,nmol,v1abs,v2abs,dvabs,
     *nptabs,absrb)
c
      return
      end
c**********************************************************************
      subroutine xint_51(v1a,v2a,dva,a,afact,vft,dvr3,r3,n1r3,n2r3)
      INCLUDE "avhrr.inc"
c
      implicit double precision (v)
c
c     this subroutine interpolates the a array stored
c     from v1a to v2a in increments of dva using a multiplicative
c     factor afact, into the r3 array from location n1r3 to n2r3 in
c     increments of dvr3
c
      dimension a(*),r3(*)
      data onepl/1.001/, onemi/0.999/, argmin/34./
c
      recdva = 1./dva
      ilo = (v1a+dva-vft)/dvr3+1.+onemi
      ilo = max(ilo,n1r3)
      ihi = (v2a-dva-vft)/dvr3+onemi
      ihi = min(ihi,n2r3)
c
      do 10 i = ilo, ihi
         vi = vft+dvr3*float(i-1)
         j = (vi-v1a)*recdva+onepl
         vj = v1a+dva*float(j-1)
         p = recdva*(vi-vj)
         c = (3.-2.*p)*p*p
         b = 0.5*p*(1.-p)
         b1 = b*(1.-p)
         b2 = b*p
         conti = -a(j-1)*b1+a(j)*(1.-c+b2)+a(j+1)*(c+b1)-a(j+2)*b2
         r3(i) = r3(i)+conti*afact
   10 continue
c
      return
c
      end
      function radfn51 (vi,xkt)
c
      implicit double precision (v)
c
c     function radfn51 calculates the radiation term for the line shape
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c               last modification:    12 august 1991
c
c                  implementation:    r.d. worsham
c
c             algorithm revisions:    s.a. clough
c                                     r.d. worsham
c                                     j.l. moncet
c
c
c                     atmospheric and environmental research inc.
c                     840 memorial drive,  cambridge, ma   02139
c
c----------------------------------------------------------------------
c
c               work supported by:    the arm program
c                                     office of energy research
c                                     department of energy
c
c
c      source of original routine:    afgl line-by-line model
c
c                                             fascod3
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      data onepl/1.001/, onemi/0.999/, argmin/34./
c
c      in the small xviokt region 0.5 is required
c
      xvi = vi
c
      if (xkt.gt.0.0) then
c
         xviokt = xvi/xkt
c
         if (xviokt.le.0.01) then
            radfn51 = 0.5*xviokt*xvi
c
         elseif (xviokt.le.10.0) then
            expvkt = exp(-xviokt)
            radfn51 = xvi*(1.-expvkt)/(1.+expvkt)
c
         else
            radfn51 = xvi
         endif
c
      else
         radfn51 = xvi
      endif
c
      return
c
      end
c**********************************************************************
      subroutine contnm_51(jrad,pave,tave,wk,wbroad,nmol,v1abs,v2abs,
     *dvabs,nptabs,absrb)
      INCLUDE "avhrr.inc"
c
      implicit double precision (v)
c
c     subroutine contnm contains the continuum data
c     which is interpolated into the array absrb
c
c
      dimension absrb(*)
      dimension wk(60)
      dimension ch(2050),csh2o(2050),cfh2o(2050)
      dimension c(2050),c0(2050),c1(2050),c2(2050)
      dimension xfac(0:50)
c
      data pi /3.1415927/
      data planck / 6.626176e-27 /,boltz / 1.380662e-16 /,
     *     clight / 2.99792458e10 /,avog / 6.022045e23 /
      data p0 / 1013. /,t0 / 296. /
      data xlosmt / 2.68675e+19 /
c
c     these are self-continuum modification factors from 700-1200 cm-1
c
      data (xfac(i),i=0,50)/
     1    1.00000,1.01792,1.03767,1.05749,1.07730,1.09708,
     2    1.10489,1.11268,1.12047,1.12822,1.13597,1.14367,
     3    1.15135,1.15904,1.16669,1.17431,1.18786,1.20134,
     4    1.21479,1.22821,1.24158,1.26580,1.28991,1.28295,
     5    1.27600,1.26896,1.25550,1.24213,1.22879,1.21560,
     6    1.20230,1.18162,1.16112,1.14063,1.12016,1.10195,
     7    1.09207,1.08622,1.08105,1.07765,1.07398,1.06620,
     8    1.05791,1.04905,1.03976,1.02981,1.00985,1.00000,
     9    1.00000,1.00000,1.00000/
c
      radcn1=2.*planck*clight*clight*1.e-07
      radcn2=planck*clight/boltz
c
c
c     assign sccs version number to module
c
cbm      hvrcnt = '3.3'
c
      rhoave = (pave/p0)*(t0/tave)
      xkt = tave/radcn2
      wtot = wbroad
      do 10 m = 1, nmol
         wtot = wtot+wk(m)
   10 continue
      wtot = wtot*1.e-20
      patm = pave/p0
c
c     ********    water    ********
c
      call sl296_51(v1c,v2c,dvc,nptc,c0,v1abs,v2abs,dvabs,nptabs)
      call sl260_51(v1c,v2c,dvc,nptc,c1,v1abs,v2abs,dvabs,nptabs)
      call frn296_51(v1c,v2c,dvc,nptc,c2,v1abs,v2abs,dvabs,nptabs)
c
      w1 = wk(1)*1.0e-20
c
c     the factor of 1.e-20 is handled this way to avoid underflows
c
      ph2o = patm*(w1/wtot)
      rh2o = ph2o*(t0/tave)
      pfrgn = patm-ph2o
      rfrgn = pfrgn*(t0/tave)
      xkt = tave/radcn2
      tfac = (tave-t0)/(260.-t0)
c
c--------------------------------------------------------------------
c                             self

      alpha2 = 200.**2
c
      alphs2= 120.**2
      betas = 5.e-06
      v0s=1310.
      factrs= 0.15
c
c--------------------------------------------------------------------
c                             foreign
      hwsqf= 330.**2
      betaf = 8.  e-11
      v0f =1130.
      factrf = 0.97
c
      v0f2 =1900.
      hwsqf2 = 150.**2
      beta2 = 3.e-06
c
c--------------------------------------------------------------------
c
      v1h=v1c
      dvh=dvc
      npth=nptc
c
      do 20 j = 1, nptc
         vj = v1c+dvc*float(j-1)
         vs2 = (vj-v0s)**2
         sh2o = 0.
         if (c0(j).gt.0.) then
            sh2o = c0(j)*(c1(j)/c0(j))**tfac
c
         sfac = 1.
         if (vj.ge.700. .and.  vj.le.1200.) then
            jfac = (vj-700.)/10. + 0.00001
            sfac = xfac(jfac)
         endif
c
c     correction to self continuum (1 sept 85); factor of 0.78 at 1000
c                             and  .......
c
      sh2o = sfac * sh2o*(1.-0.2333*(alpha2/((vj-1050.)**2+alpha2))) *
     *                  (1.-factrs*(alphs2/(vs2+(betas*vs2**2)+alphs2)))
         endif
c
c     correction to foreign continuum
c
        vf2 = (vj-v0f)**2
        vf6 = vf2 * vf2 * vf2
        fscal  = (1.-factrf*(hwsqf/(vf2+(betaf*vf6)+hwsqf)))
        vf2 = (vj-v0f2)**2
        vf4 = vf2*vf2
        fscal = fscal* (1.- 0.6*(hwsqf2/(vf2 + beta2*vf4 + hwsqf2)))
c
        c2(j)=c2(j)*fscal
c
         c(j) = w1*(sh2o*rh2o+c2(j)*rfrgn)
c
      csh2o(j)=1.e-20 * sh2o
      cfh2o(j)=1.e-20 * c2(j)
c
c     radiation field
c
         if (jrad.eq.1) c(j) = c(j)*radfn51(vj,xkt)
   20 continue
      call xint_51(v1c,v2c,dvc,c,1.0,v1abs,dvabs,absrb,1,nptabs)
c
      return
c
      end
c**********************************************************************
      subroutine sl296_51(v1c,v2c,dvc,nptc,c,v1abs,v2abs,dvabs,nptabs)
      INCLUDE "avhrr.inc"
c
      implicit double precision (v)
c
      dimension s(305)
      dimension c(*)
c
c               06/28/82
c               units of (cm**3/mol) * 1.e-20
c
      data v1s,v2s,dvs,npts / -20.0, 3020.0, 10.0, 305/
c
      data (s(i),i=1,2)/
     *     1.1109e-01 ,1.0573e-01/
      data (s(i),i=3,52)/
     *     1.0162e-01, 1.0573e-01, 1.1109e-01, 1.2574e-01, 1.3499e-01,
     *     1.4327e-01, 1.5065e-01, 1.5164e-01, 1.5022e-01, 1.3677e-01,
     *     1.3115e-01, 1.2253e-01, 1.1271e-01, 1.0070e-01, 8.7495e-02,
     *     8.0118e-02, 6.9940e-02, 6.2034e-02, 5.6051e-02, 4.7663e-02,
     *     4.2450e-02, 3.6690e-02, 3.3441e-02, 3.0711e-02, 2.5205e-02,
     *     2.2113e-02, 1.8880e-02, 1.6653e-02, 1.4626e-02, 1.2065e-02,
     *     1.0709e-02, 9.1783e-03, 7.7274e-03, 6.7302e-03, 5.6164e-03,
     *     4.9089e-03, 4.1497e-03, 3.5823e-03, 3.1124e-03, 2.6414e-03,
     *     2.3167e-03, 2.0156e-03, 1.7829e-03, 1.5666e-03, 1.3928e-03,
     *     1.2338e-03, 1.0932e-03, 9.7939e-04, 8.8241e-04, 7.9173e-04/
      data (s(i),i=53,102)/
     *     7.1296e-04, 6.4179e-04, 5.8031e-04, 5.2647e-04, 4.7762e-04,
     *     4.3349e-04, 3.9355e-04, 3.5887e-04, 3.2723e-04, 2.9919e-04,
     *     2.7363e-04, 2.5013e-04, 2.2876e-04, 2.0924e-04, 1.9193e-04,
     *     1.7618e-04, 1.6188e-04, 1.4891e-04, 1.3717e-04, 1.2647e-04,
     *     1.1671e-04, 1.0786e-04, 9.9785e-05, 9.2350e-05, 8.5539e-05,
     *     7.9377e-05, 7.3781e-05, 6.8677e-05, 6.3993e-05, 5.9705e-05,
     *     5.5788e-05, 5.2196e-05, 4.8899e-05, 4.5865e-05, 4.3079e-05,
     *     4.0526e-05, 3.8182e-05, 3.6025e-05, 3.4038e-05, 3.2203e-05,
     *     3.0511e-05, 2.8949e-05, 2.7505e-05, 2.6170e-05, 2.4933e-05,
     *     2.3786e-05, 2.2722e-05, 2.1736e-05, 2.0819e-05, 1.9968e-05/
      data (s(i),i=103,152)/
     *     1.9178e-05, 1.8442e-05, 1.7760e-05, 1.7127e-05, 1.6541e-05,
     *     1.5997e-05, 1.5495e-05, 1.5034e-05, 1.4614e-05, 1.4230e-05,
     *     1.3883e-05, 1.3578e-05, 1.3304e-05, 1.3069e-05, 1.2876e-05,
     *     1.2732e-05, 1.2626e-05, 1.2556e-05, 1.2544e-05, 1.2604e-05,
     *     1.2719e-05, 1.2883e-05, 1.3164e-05, 1.3581e-05, 1.4187e-05,
     *     1.4866e-05, 1.5669e-05, 1.6717e-05, 1.8148e-05, 2.0268e-05,
     *     2.2456e-05, 2.5582e-05, 2.9183e-05, 3.3612e-05, 3.9996e-05,
     *     4.6829e-05, 5.5055e-05, 6.5897e-05, 7.5360e-05, 8.7213e-05,
     *     1.0046e-04, 1.1496e-04, 1.2943e-04, 1.5049e-04, 1.6973e-04,
     *     1.8711e-04, 2.0286e-04, 2.2823e-04, 2.6780e-04, 2.8766e-04/
      data (s(i),i=153,202)/
     *     3.1164e-04, 3.3640e-04, 3.6884e-04, 3.9159e-04, 3.8712e-04,
     *     3.7433e-04, 3.4503e-04, 3.1003e-04, 2.8027e-04, 2.5253e-04,
     *     2.3408e-04, 2.2836e-04, 2.4442e-04, 2.7521e-04, 2.9048e-04,
     *     3.0489e-04, 3.2646e-04, 3.3880e-04, 3.3492e-04, 3.0987e-04,
     *     2.9482e-04, 2.8711e-04, 2.6068e-04, 2.2683e-04, 1.9996e-04,
     *     1.7788e-04, 1.6101e-04, 1.3911e-04, 1.2013e-04, 1.0544e-04,
     *     9.4224e-05, 8.1256e-05, 7.3667e-05, 6.2233e-05, 5.5906e-05,
     *     5.1619e-05, 4.5140e-05, 4.0273e-05, 3.3268e-05, 3.0258e-05,
     *     2.6440e-05, 2.3103e-05, 2.0749e-05, 1.8258e-05, 1.6459e-05,
     *     1.4097e-05, 1.2052e-05, 1.0759e-05, 9.1400e-06, 8.1432e-06/
      data (s(i),i=203,252)/
     *     7.1460e-06, 6.4006e-06, 5.6995e-06, 4.9372e-06, 4.4455e-06,
     *     3.9033e-06, 3.4740e-06, 3.1269e-06, 2.8059e-06, 2.5558e-06,
     *     2.2919e-06, 2.0846e-06, 1.8983e-06, 1.7329e-06, 1.5929e-06,
     *     1.4631e-06, 1.3513e-06, 1.2461e-06, 1.1519e-06, 1.0682e-06,
     *     9.9256e-07, 9.2505e-07, 8.6367e-07, 8.0857e-07, 7.5674e-07,
     *     7.0934e-07, 6.6580e-07, 6.2580e-07, 5.8853e-07, 5.5333e-07,
     *     5.2143e-07, 4.9169e-07, 4.6431e-07, 4.3898e-07, 4.1564e-07,
     *     3.9405e-07, 3.7403e-07, 3.5544e-07, 3.3819e-07, 3.2212e-07,
     *     3.0714e-07, 2.9313e-07, 2.8003e-07, 2.6777e-07, 2.5628e-07,
     *     2.4551e-07, 2.3540e-07, 2.2591e-07, 2.1701e-07, 2.0866e-07/
      data (s(i),i=253,302)/
     *     2.0082e-07, 1.9349e-07, 1.8665e-07, 1.8027e-07, 1.7439e-07,
     *     1.6894e-07, 1.6400e-07, 1.5953e-07, 1.5557e-07, 1.5195e-07,
     *     1.4888e-07, 1.4603e-07, 1.4337e-07, 1.4093e-07, 1.3828e-07,
     *     1.3569e-07, 1.3270e-07, 1.2984e-07, 1.2714e-07, 1.2541e-07,
     *     1.2399e-07, 1.2102e-07, 1.1878e-07, 1.1728e-07, 1.1644e-07,
     *     1.1491e-07, 1.1305e-07, 1.1235e-07, 1.1228e-07, 1.1224e-07,
     *     1.1191e-07, 1.1151e-07, 1.1098e-07, 1.1068e-07, 1.1109e-07,
     *     1.1213e-07, 1.1431e-07, 1.1826e-07, 1.2322e-07, 1.3025e-07,
     *     1.4066e-07, 1.5657e-07, 1.7214e-07, 1.9449e-07, 2.2662e-07,
     *     2.6953e-07, 3.1723e-07, 3.7028e-07, 4.4482e-07, 5.3852e-07/
      data (s(i),i=303,305)/
     *     6.2639e-07, 7.2175e-07, 7.7626e-07/
c
      dvc = dvs
      v1c = v1abs-dvc
      v2c = v2abs+dvc
c
      i1 = (v1c-v1s)/dvs
      if (v1c.lt.v1s) i1 = i1-1
c
      v1c = v1s+dvs*float(i1)
      i2 = (v2c-v1s)/dvs
      nptc = i2-i1+3
      v2c = v1c+dvs*float(nptc-1)
      do 10 j = 1, nptc
         i = i1+j
         c(j) = 0.
         if ((i.lt.1).or.(i.gt.npts)) go to 10
         c(j) = s(i)
   10 continue
c
      return
c
      end
c**********************************************************************
      subroutine sl260_51(v1c,v2c,dvc,nptc,c,v1abs,v2abs,dvabs,nptabs)
      INCLUDE "avhrr.inc"
c
      implicit double precision (v)
c
      dimension s(305)
      dimension c(*)
c
c               06/28/82
c               units of (cm**3/mol) * 1.e-20
c
       data v1s,v2s,dvs,npts / -20.0, 3020.0, 10.0, 305/
c
c
      data (s(i),i=1,2)/
     *     1.7750e-01, 1.7045e-01/
      data (s(i),i=3,52)/
     *     1.6457e-01, 1.7045e-01, 1.7750e-01, 2.0036e-01, 2.1347e-01,
     *     2.2454e-01, 2.3428e-01, 2.3399e-01, 2.3022e-01, 2.0724e-01,
     *     1.9712e-01, 1.8317e-01, 1.6724e-01, 1.4780e-01, 1.2757e-01,
     *     1.1626e-01, 1.0098e-01, 8.9033e-02, 7.9770e-02, 6.7416e-02,
     *     5.9588e-02, 5.1117e-02, 4.6218e-02, 4.2179e-02, 3.4372e-02,
     *     2.9863e-02, 2.5252e-02, 2.2075e-02, 1.9209e-02, 1.5816e-02,
     *     1.3932e-02, 1.1943e-02, 1.0079e-02, 8.7667e-03, 7.4094e-03,
     *     6.4967e-03, 5.5711e-03, 4.8444e-03, 4.2552e-03, 3.6953e-03,
     *     3.2824e-03, 2.9124e-03, 2.6102e-03, 2.3370e-03, 2.1100e-03,
     *     1.9008e-03, 1.7145e-03, 1.5573e-03, 1.4206e-03, 1.2931e-03/
      data (s(i),i=53,102)/
     *     1.1803e-03, 1.0774e-03, 9.8616e-04, 9.0496e-04, 8.3071e-04,
     *     7.6319e-04, 7.0149e-04, 6.4637e-04, 5.9566e-04, 5.4987e-04,
     *     5.0768e-04, 4.6880e-04, 4.3317e-04, 4.0037e-04, 3.7064e-04,
     *     3.4325e-04, 3.1809e-04, 2.9501e-04, 2.7382e-04, 2.5430e-04,
     *     2.3630e-04, 2.1977e-04, 2.0452e-04, 1.9042e-04, 1.7740e-04,
     *     1.6544e-04, 1.5442e-04, 1.4425e-04, 1.3486e-04, 1.2618e-04,
     *     1.1817e-04, 1.1076e-04, 1.0391e-04, 9.7563e-05, 9.1696e-05,
     *     8.6272e-05, 8.1253e-05, 7.6607e-05, 7.2302e-05, 6.8311e-05,
     *     6.4613e-05, 6.1183e-05, 5.8001e-05, 5.5048e-05, 5.2307e-05,
     *     4.9761e-05, 4.7395e-05, 4.5197e-05, 4.3155e-05, 4.1256e-05/
      data (s(i),i=103,152)/
     *     3.9491e-05, 3.7849e-05, 3.6324e-05, 3.4908e-05, 3.3594e-05,
     *     3.2374e-05, 3.1244e-05, 3.0201e-05, 2.9240e-05, 2.8356e-05,
     *     2.7547e-05, 2.6814e-05, 2.6147e-05, 2.5551e-05, 2.5029e-05,
     *     2.4582e-05, 2.4203e-05, 2.3891e-05, 2.3663e-05, 2.3531e-05,
     *     2.3483e-05, 2.3516e-05, 2.3694e-05, 2.4032e-05, 2.4579e-05,
     *     2.5234e-05, 2.6032e-05, 2.7119e-05, 2.8631e-05, 3.0848e-05,
     *     3.3262e-05, 3.6635e-05, 4.0732e-05, 4.5923e-05, 5.3373e-05,
     *     6.1875e-05, 7.2031e-05, 8.5980e-05, 9.8642e-05, 1.1469e-04,
     *     1.3327e-04, 1.5390e-04, 1.7513e-04, 2.0665e-04, 2.3609e-04,
     *     2.6220e-04, 2.8677e-04, 3.2590e-04, 3.8624e-04, 4.1570e-04/
      data (s(i),i=153,202)/
     *     4.5207e-04, 4.9336e-04, 5.4500e-04, 5.8258e-04, 5.8086e-04,
     *     5.6977e-04, 5.3085e-04, 4.8020e-04, 4.3915e-04, 4.0343e-04,
     *     3.7853e-04, 3.7025e-04, 3.9637e-04, 4.4675e-04, 4.7072e-04,
     *     4.9022e-04, 5.2076e-04, 5.3676e-04, 5.2755e-04, 4.8244e-04,
     *     4.5473e-04, 4.3952e-04, 3.9614e-04, 3.4086e-04, 2.9733e-04,
     *     2.6367e-04, 2.3767e-04, 2.0427e-04, 1.7595e-04, 1.5493e-04,
     *     1.3851e-04, 1.1874e-04, 1.0735e-04, 9.0490e-05, 8.1149e-05,
     *     7.4788e-05, 6.5438e-05, 5.8248e-05, 4.8076e-05, 4.3488e-05,
     *     3.7856e-05, 3.3034e-05, 2.9592e-05, 2.6088e-05, 2.3497e-05,
     *     2.0279e-05, 1.7526e-05, 1.5714e-05, 1.3553e-05, 1.2145e-05/
      data (s(i),i=203,252)/
     *     1.0802e-05, 9.7681e-06, 8.8196e-06, 7.8291e-06, 7.1335e-06,
     *     6.4234e-06, 5.8391e-06, 5.3532e-06, 4.9079e-06, 4.5378e-06,
     *     4.1716e-06, 3.8649e-06, 3.5893e-06, 3.3406e-06, 3.1199e-06,
     *     2.9172e-06, 2.7348e-06, 2.5644e-06, 2.4086e-06, 2.2664e-06,
     *     2.1359e-06, 2.0159e-06, 1.9051e-06, 1.8031e-06, 1.7074e-06,
     *     1.6185e-06, 1.5356e-06, 1.4584e-06, 1.3861e-06, 1.3179e-06,
     *     1.2545e-06, 1.1951e-06, 1.1395e-06, 1.0873e-06, 1.0384e-06,
     *     9.9250e-07, 9.4935e-07, 9.0873e-07, 8.7050e-07, 8.3446e-07,
     *     8.0046e-07, 7.6834e-07, 7.3800e-07, 7.0931e-07, 6.8217e-07,
     *     6.5648e-07, 6.3214e-07, 6.0909e-07, 5.8725e-07, 5.6655e-07/
      data (s(i),i=253,302)/
     *     5.4693e-07, 5.2835e-07, 5.1077e-07, 4.9416e-07, 4.7853e-07,
     *     4.6381e-07, 4.5007e-07, 4.3728e-07, 4.2550e-07, 4.1450e-07,
     *     4.0459e-07, 3.9532e-07, 3.8662e-07, 3.7855e-07, 3.7041e-07,
     *     3.6254e-07, 3.5420e-07, 3.4617e-07, 3.3838e-07, 3.3212e-07,
     *     3.2655e-07, 3.1865e-07, 3.1203e-07, 3.0670e-07, 3.0252e-07,
     *     2.9749e-07, 2.9184e-07, 2.8795e-07, 2.8501e-07, 2.8202e-07,
     *     2.7856e-07, 2.7509e-07, 2.7152e-07, 2.6844e-07, 2.6642e-07,
     *     2.6548e-07, 2.6617e-07, 2.6916e-07, 2.7372e-07, 2.8094e-07,
     *     2.9236e-07, 3.1035e-07, 3.2854e-07, 3.5481e-07, 3.9377e-07,
     *     4.4692e-07, 5.0761e-07, 5.7715e-07, 6.7725e-07, 8.0668e-07/
      data (s(i),i=303,305)/
     *     9.3716e-07, 1.0797e-06, 1.1689e-06/
c
      dvc = dvs
      v1c = v1abs-dvc
      v2c = v2abs+dvc
c
      i1 = (v1c-v1s)/dvs
      if (v1c.lt.v1s) i1 = i1-1
c
      v1c = v1s+dvs*float(i1)
      i2 = (v2c-v1s)/dvs
      nptc = i2-i1+3
      v2c = v1c+dvs*float(nptc-1)
      do 10 j = 1, nptc
         i = i1+j
         c(j) = 0.
         if ((i.lt.1).or.(i.gt.npts)) go to 10
         c(j) = s(i)
   10 continue
c
      return
c
      end
c**********************************************************************
      subroutine frn296_51(v1c,v2c,dvc,nptc,c,v1abs,v2abs,dvabs,nptabs)
      INCLUDE "avhrr.inc"
c
      implicit double precision (v)
c
      dimension s(305)
      dimension c(*)
c
c               06/28/82
c               units of (cm**3/mol)*1.e-20
c
       data v1s,v2s,dvs,npts / -20.0, 3020.0, 10.0, 305/
c
      data (s(i),i=1,2)/
     *     1.2859e-02, 1.1715e-02/
      data (s(i),i=3,52)/
     *     1.1038e-02, 1.1715e-02, 1.2859e-02, 1.5326e-02, 1.6999e-02,
     *     1.8321e-02, 1.9402e-02, 1.9570e-02, 1.9432e-02, 1.7572e-02,
     *     1.6760e-02, 1.5480e-02, 1.3984e-02, 1.2266e-02, 1.0467e-02,
     *     9.4526e-03, 8.0485e-03, 6.9484e-03, 6.1416e-03, 5.0941e-03,
     *     4.4836e-03, 3.8133e-03, 3.4608e-03, 3.1487e-03, 2.4555e-03,
     *     2.0977e-03, 1.7266e-03, 1.4920e-03, 1.2709e-03, 9.8081e-04,
     *     8.5063e-04, 6.8822e-04, 5.3809e-04, 4.4679e-04, 3.3774e-04,
     *     2.7979e-04, 2.1047e-04, 1.6511e-04, 1.2993e-04, 9.3033e-05,
     *     7.4360e-05, 5.6428e-05, 4.5442e-05, 3.4575e-05, 2.7903e-05,
     *     2.1374e-05, 1.6075e-05, 1.3022e-05, 1.0962e-05, 8.5959e-06/
      data (s(i),i=53,102)/
     *     6.9125e-06, 5.3808e-06, 4.3586e-06, 3.6394e-06, 2.9552e-06,
     *     2.3547e-06, 1.8463e-06, 1.6036e-06, 1.3483e-06, 1.1968e-06,
     *     1.0333e-06, 8.4484e-07, 6.7195e-07, 5.0947e-07, 4.2343e-07,
     *     3.4453e-07, 2.7830e-07, 2.3063e-07, 1.9951e-07, 1.7087e-07,
     *     1.4393e-07, 1.2575e-07, 1.0750e-07, 8.2325e-08, 5.7524e-08,
     *     4.4482e-08, 3.8106e-08, 3.4315e-08, 2.9422e-08, 2.5069e-08,
     *     2.2402e-08, 1.9349e-08, 1.6152e-08, 1.2208e-08, 8.9660e-09,
     *     7.1322e-09, 6.1028e-09, 5.2938e-09, 4.5350e-09, 3.4977e-09,
     *     2.9511e-09, 2.4734e-09, 2.0508e-09, 1.8507e-09, 1.6373e-09,
     *     1.5171e-09, 1.3071e-09, 1.2462e-09, 1.2148e-09, 1.2590e-09/
      data (s(i),i=103,152)/
     *     1.3153e-09, 1.3301e-09, 1.4483e-09, 1.6944e-09, 2.0559e-09,
     *     2.2954e-09, 2.6221e-09, 3.2606e-09, 4.2392e-09, 5.2171e-09,
     *     6.2553e-09, 8.2548e-09, 9.5842e-09, 1.1280e-08, 1.3628e-08,
     *     1.7635e-08, 2.1576e-08, 2.4835e-08, 3.0014e-08, 3.8485e-08,
     *     4.7440e-08, 5.5202e-08, 7.0897e-08, 9.6578e-08, 1.3976e-07,
     *     1.8391e-07, 2.3207e-07, 2.9960e-07, 4.0408e-07, 5.9260e-07,
     *     7.8487e-07, 1.0947e-06, 1.4676e-06, 1.9325e-06, 2.6587e-06,
     *     3.4534e-06, 4.4376e-06, 5.8061e-06, 7.0141e-06, 8.4937e-06,
     *     1.0186e-05, 1.2034e-05, 1.3837e-05, 1.6595e-05, 1.9259e-05,
     *     2.1620e-05, 2.3681e-05, 2.7064e-05, 3.2510e-05, 3.5460e-05/
      data (s(i),i=153,202)/
     *     3.9109e-05, 4.2891e-05, 4.7757e-05, 5.0981e-05, 5.0527e-05,
     *     4.8618e-05, 4.4001e-05, 3.7982e-05, 3.2667e-05, 2.7794e-05,
     *     2.4910e-05, 2.4375e-05, 2.7316e-05, 3.2579e-05, 3.5499e-05,
     *     3.8010e-05, 4.1353e-05, 4.3323e-05, 4.3004e-05, 3.9790e-05,
     *     3.7718e-05, 3.6360e-05, 3.2386e-05, 2.7409e-05, 2.3626e-05,
     *     2.0631e-05, 1.8371e-05, 1.5445e-05, 1.2989e-05, 1.1098e-05,
     *     9.6552e-06, 8.0649e-06, 7.2365e-06, 5.9137e-06, 5.2759e-06,
     *     4.8860e-06, 4.1321e-06, 3.5918e-06, 2.7640e-06, 2.4892e-06,
     *     2.1018e-06, 1.7848e-06, 1.5855e-06, 1.3569e-06, 1.1986e-06,
     *     9.4693e-07, 7.4097e-07, 6.3443e-07, 4.8131e-07, 4.0942e-07/
      data (s(i),i=203,252)/
     *     3.3316e-07, 2.8488e-07, 2.3461e-07, 1.7397e-07, 1.4684e-07,
     *     1.0953e-07, 8.5396e-08, 6.9261e-08, 5.4001e-08, 4.5430e-08,
     *     3.2791e-08, 2.5995e-08, 2.0225e-08, 1.5710e-08, 1.3027e-08,
     *     1.0229e-08, 8.5277e-09, 6.5249e-09, 5.0117e-09, 3.9906e-09,
     *     3.2332e-09, 2.7847e-09, 2.4570e-09, 2.3359e-09, 2.0599e-09,
     *     1.8436e-09, 1.6559e-09, 1.4910e-09, 1.2794e-09, 9.8229e-10,
     *     8.0054e-10, 6.0769e-10, 4.5646e-10, 3.3111e-10, 2.4428e-10,
     *     1.8007e-10, 1.3291e-10, 9.7974e-11, 7.8271e-11, 6.3833e-11,
     *     5.4425e-11, 4.6471e-11, 4.0209e-11, 3.5227e-11, 3.1212e-11,
     *     2.8840e-11, 2.7762e-11, 2.7935e-11, 3.2012e-11, 3.9525e-11/
      data (s(i),i=253,302)/
     *     5.0303e-11, 6.8027e-11, 9.3954e-11, 1.2986e-10, 1.8478e-10,
     *     2.5331e-10, 3.4827e-10, 4.6968e-10, 6.2380e-10, 7.9106e-10,
     *     1.0026e-09, 1.2102e-09, 1.4146e-09, 1.6154e-09, 1.7510e-09,
     *     1.8575e-09, 1.8742e-09, 1.8700e-09, 1.8582e-09, 1.9657e-09,
     *     2.1204e-09, 2.0381e-09, 2.0122e-09, 2.0436e-09, 2.1213e-09,
     *     2.0742e-09, 1.9870e-09, 2.0465e-09, 2.1556e-09, 2.2222e-09,
     *     2.1977e-09, 2.1047e-09, 1.9334e-09, 1.7357e-09, 1.5754e-09,
     *     1.4398e-09, 1.4018e-09, 1.5459e-09, 1.7576e-09, 2.1645e-09,
     *     2.9480e-09, 4.4439e-09, 5.8341e-09, 8.0757e-09, 1.1658e-08,
     *     1.6793e-08, 2.2694e-08, 2.9468e-08, 3.9278e-08, 5.2145e-08/
      data (s(i),i=303,305)/
     *     6.4378e-08, 7.7947e-08, 8.5321e-08/
c
      dvc = dvs
      v1c = v1abs-dvc
      v2c = v2abs+dvc
c
      i1 = (v1c-v1s)/dvs
      if (v1c.lt.v1s) i1 = i1-1
c
      v1c = v1s+dvs*float(i1)
      i2 = (v2c-v1s)/dvs
      nptc = i2-i1+3
      v2c = v1c+dvs*float(nptc-1)
      do 10 j = 1, nptc
         i = i1+j
         c(j) = 0.
         if ((i.ge.1).and.(i.le.npts)) then
            c(j) = s(i)
         endif
   10 continue
c
      return
c
      end


text:
      subroutine doflux
c******************************************************************************
c     This routine produces flux curves
c******************************************************************************

      implicit real*8 (a-h,o-z)
      include 'Atmos.com'
      include 'Linex.com'
      include 'Dummy.com'
      include 'Pstuff.com'
      include 'Source.com'

c*****examine the parameter file
      call params

c*****open the files for standard output and summary curves-of-growth
      nf1out = 20
      lscreen = 4
      array = 'STANDARD OUTPUT'
      nchars = 15
      call infile ('output ',nf1out,'formatted  ',0,nchars,
     .             f1out,lscreen)
      nf2out = 21
      lscreen = 6
      array = 'SUMMARY FLUX OUTPUT'
      nchars = 19
      call infile ('output ',nf2out,'formatted  ',0,nchars,
     .             f2out,lscreen)
      nf5out = 26
      lscreen = lscreen + 2
      array = 'POSTSCRIPT PLOT OUTPUT'
      nchars = 22
      call infile ('output ',nf5out,'formatted  ',0,nchars,
     .             f5out,lscreen)

c*****open and read the model atmosphere
      nfmodel = 30
      lscreen = lscreen + 2
      array = 'THE MODEL ATMOSPHERE'
      nchars = 20
      call infile ('input  ',nfmodel,'formatted  ',0,nchars,
     .             fmodel,lscreen)
      call inmodel

c*****compute the flux curve
      wave = start
1     call opacit (2,wave)
      if (modprintopt .ge. 2)
     .    write(nf1out,1002) wave,(kaplam(i),i=1,ntau)
      if (scatopt .eq. 0) then
         call cdcalc (1)
         first = 0.4343*cd(1)
         flux = rinteg(xref,cd,dummy1,ntau,first)
      else
         call cdcalc_JS (1)
         flux = Flux_cont_moog
      endif
      if (flux .le. 0.1) flux = 0.
      if (iunits .eq. 1) then
         write (nf1out,1003) 1.d-4*wave,flux
      else
         write (nf1out,1004) wave,flux
      endif
      waveinv = 1.0d4/wave
      if (flux .gt. 0.) then
         fluxlog = dlog10(flux)
      else
         fluxlog = -1.0
      endif
      write (nf2out,1001) wave, flux, waveinv, fluxlog
      wave = wave + step
      if (wave .le. sstop) go to 1
      call pltflux

c****end the computations
      call finish (0)
      return

c*****format statements
1001  format (1p2d12.4,0p2f10.4)
1002  format ('  kaplam from 1 to ntau at wavelength',f11.3/
     1        (6(1pd12.4)))
1003  format ('AT WAVELENGTH/FREQUENCY =',f11.7,
     .        '  CONTINUUM FLUX/INTENSITY =', 1pd12.5)
1004  format ('AT WAVELENGTH/FREQUENCY =',f11.3,
     .        '  CONTINUUM FLUX/INTENSITY =', 1pd12.5)

      end


text:
subroutine foo()
   I = 0
   DO WHILE (I.GE.10)
      I = 7
   END DO
end subroutine

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab