text:
#define type(x) TYPE(x), target

Module matProp_Mod
  ! Module for loading and validating material properties

  Private

  Type matProps
    ! Stores the set of properties

    Character(len=80) :: name ! Material name, required

    ! Material properties, grouped by feature
    ! (Required inputs)
    Double Precision :: E1, E2, G12, v12, v23                            ! Req. Transversely isotropic.
    ! (Optional inputs)
    Double Precision :: YT, SL, GYT, GSL, eta_BK, YC, alpha0             ! Opt. Matrix failure
    Double Precision :: E3, G13, G23, v13                                ! Opt. Orthotropic. Defaults are calculated based on transverse isotropy
    Double Precision :: cte(3)                                           ! Opt. CTEs (Required if a temperature change is applied in the analysis). Defaults to zero
    Double Precision :: aPL, nPL                                         ! Opt. Shear nonlinearity
    Double Precision :: XT, fXT, GXT, fGXT                               ! Opt. Longitudinal tensile damage (max strain failure criterion, bilinear softening)
    Double Precision :: XC, fXC, GXC, fGXC                               ! Opt. Longitudinal compressive damage (max strain failure criterion, bilinear softening)
    Double Precision :: w_kb                                             ! Opt. LaRC15 kink band model (Also requires XC, shear nonlinearity).
    Double Precision :: cl                                           ! Opt. fiber nonlinearity (defaults to zero, which is no fiber nonlinearity)
    Double Precision :: mu                                               ! Opt. Friction. Defaults to zero
    Double Precision :: es(4), gs(4)

    ! min and max values for acceptable range
    Double Precision, private :: modulus_min, modulus_max
    Double Precision, private :: poissons_min, poissons_max
    Double Precision, private :: strength_min, strength_max
    Double Precision, private :: toughness_min, toughness_max
    Double Precision, private :: eta_BK_min, eta_BK_max
    Double Precision, private :: cte_min, cte_max
    Double Precision, private :: aPL_min, aPL_max, nPL_min, nPL_max
    Double Precision, private :: w_kb_min, w_kb_max
    Double Precision, private :: cl_min, cl_max
    Double Precision, private :: alpha0_min, alpha0_max
    Double Precision, private :: mu_min, mu_max
    Double Precision, private :: schapery_min, schapery_max

    ! Flags that indicate if values have been set
    Logical, private :: E1_def, E2_def, G12_def, v12_def, v23_def
    Logical, private :: YT_def, SL_def, GYT_def, GSL_def, eta_BK_def
    Logical, private :: E3_def, G13_def, G23_def, v13_def
    Logical, private :: cte_def(3)
    Logical, private :: aPL_def, nPL_def
    Logical, private :: XT_def, fXT_def, GXT_def, fGXT_def
    Logical, private :: XC_def, fXC_def, GXC_def, fGXC_def
    Logical, private :: YC_def, w_kb_def, alpha0_def
    Logical, private :: cl_def
    Logical, private :: mu_def
    Logical, private :: es_def(4), gs_def(4)

    ! Calculated properties
    Double Precision :: v21, v31, v32
    Double Precision :: etaL, etaT, ST, phic, alpha0_deg

    Double Precision :: thickness

    ! Flags for features
    Logical :: matrixDam
    Logical :: shearNonlinearity12
    Logical :: shearNonlinearity13
    Logical :: schapery
    Logical :: fiberTenDam
    Logical :: fiberCompDamBL
    Logical :: fiberCompDamFKT
    Logical :: friction

  End Type matProps

  Type materialList
    ! An array of matProps

    Type(matProps), Allocatable :: materials(:)

  End Type materialList

  ! Public interface
  Public :: matProps
  Public :: materialList
  Public :: loadMatProps
  Public :: checkForSnapBack
  Public :: initializePhi0
  Public :: consistencyChecks

  ! Reference to object for singleton
  type(matProps) :: m
  type(materialList), Save :: user

Contains

  Type(matProps) Function loadMatProps(materialName, nprops, props)
    ! Loads material properties into module variable m

    Use forlog_Mod

    ! Arguments
    Character(len=*), intent(IN) :: materialName
    ! Logical, intent(IN) :: issueWarnings
    Integer, intent(IN) :: nprops
    Double Precision :: props(nprops)

    ! Locals
    Character(len=256) :: outputDir, jobName, fileName
    Integer :: lenOutputDir, lenJobName
    Logical :: fileExists

    Double Precision :: valueDbl

    Double Precision, Parameter :: zero=0.d0, one=1.d0, two=2.d0

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

    ! Initializations
    Call initializeMinMaxValues()

    ! Record material name
    m%name = trim(materialName)

#ifndef PYEXT
    ! Get the output directory (location to search for .props file)
    Call VGETOUTDIR(outputDir, lenOutputDir)

    ! Get the jobname
    Call VGETJOBNAME(jobName, lenJobName)

    ! Look to see if a material properties file exists
    ! First look for: jobName_materialName.props
    fileName = trim(outputDir) // '/' // trim(jobName) // '_' // trim(materialName) // '.props'
    Inquire(FILE=fileName, EXIST=fileExists)
    ! Try looking for: materialName.props
    If (.NOT. fileExists) Then
      fileName = trim(outputDir) // '/' // trim(materialName) // '.props'
      Inquire(FILE=fileName, EXIST=fileExists)
    End If
#else
    fileName = trim(materialName) // '.props'
    Inquire(FILE=fileName, EXIST=fileExists)
#endif

    ! Handle cases
    If (fileExists .AND. nprops > 8) Then
      Call log%warn("loadMatProps: Found material property file " // trim(fileName) // " and properties in input deck; using properties in the input deck.")
      Call loadPropertiesFromInp(nprops, props)
    Else If (fileExists) Then
      Call log%info("loadMatProps: Using material property file " // trim(fileName))
      Call loadPropertiesFromFile(trim(fileName), nprops, props)
    Else
      Call loadPropertiesFromInp(nprops, props)
      Call log%info("loadMatProps: Material property file not found. Looking for: " // trim(fileName))
      Call log%info("loadMatProps: Reading material properties from the input deck instead.")
    End If

    ! Feature flags
    If (.NOT. m%friction) m%mu = zero

    ! ! Checks that a consistent set of properties has been defined
    ! Call consistencyChecks(issueWarnings)

    ! Calculated properties
    m%v21 = m%E2*m%v12/m%E1
    m%v31 = m%E3*m%v13/m%E1
    m%v32 = m%E3*m%v23/m%E2
    m%cte(3) = m%cte(2)

    ! TODO - additional admissibility checks

    ! Return a reference to the matProps object
    loadMatProps = m

    Return
  End Function loadMatProps

  Subroutine loadPropertiesFromFile(fileName, nprops, props)
    ! Populates m with the material properties in the specified .props file

    Use forlog_Mod

    !Arguments
    Character(len=*), intent(IN) :: fileName
    Integer, intent(IN) :: nprops
    Double Precision :: props(nprops)

    ! Locals
    Integer, Parameter :: unit=108
    Integer :: iostat
    Character(len=256) :: line, key, value, tmp
    Character(len=30) :: featureFlags
    Integer :: commentTokenPos, equalTokenPos
    Double Precision, Parameter :: zero=0.d0, one=1.d0, two=2.d0
    ! -------------------------------------------------------------------- !

    Call log%debug("loadMatProps: props file name: " // fileName)

    ! Try to load material properties from a '.props' file
    Open(UNIT=unit, FILE=fileName, STATUS='old', ACTION='read', position='rewind', IOSTAT=iostat)

    If (iostat /= 0) Call log%error("loadMatProps: Unable to access the .props file")
    ReadLines: Do

      ! Read the next line in the file
      Read(unit,'(A255)',IOSTAT=iostat) line

      If (iostat > 0) Then
        Call log%error("loadMatProps: Unknown error reading file")
        Return

      Else If (iostat < 0) Then
        Call log%debug("loadMatProps: Reached end of props file")
        Exit ReadLines

      Else  ! Parse the line in the file
        ! Skip blank lines
        If (Len(Trim(line)) == 0) Cycle ReadLines

        ! Ignore comment lines (token: //)
        commentTokenPos = Index(line, '//')
        If (commentTokenPos == 1) Then
          Cycle ReadLines

        Else If (commentTokenPos > 1) Then
          ! This line has a trailing comment, retain everything before the comment
          line = line(1:commentTokenPos-1)
        End If

        ! Split the line into a key and value
        equalTokenPos = Index(line, '=')

        If (equalTokenPos == 0) Call log%error("loadMatProps: Expected [name] = [value] format not found. Line image: " // trim(line))

        ! Parse key and value from the line
        key = line(1:equalTokenPos-1)
        value = line(equalTokenPos+1:)

        ! Check if the key is a known property. If it's a known property, check validity, and store it
        Select Case (trim(key))

          Case ('E1')
            Call verifyAndSaveProperty_str(trim(key), value, m%modulus_min, m%modulus_max, m%E1, m%E1_def)

          Case ('E2')
            Call verifyAndSaveProperty_str(trim(key), value, m%modulus_min, m%modulus_max, m%E2, m%E2_def)

          Case ('G12')
            Call verifyAndSaveProperty_str(trim(key), value, m%modulus_min, m%modulus_max, m%G12, m%G12_def)

          Case ('v12')
            Call verifyAndSaveProperty_str(trim(key), value, m%poissons_min, m%poissons_max, m%v12, m%v12_def)

          Case ('v23')
            Call verifyAndSaveProperty_str(trim(key), value, m%poissons_min, m%poissons_max, m%v23, m%v23_def)

          Case ('YT')
            Call verifyAndSaveProperty_str(trim(key), value, m%strength_min, m%strength_max, m%YT, m%YT_def)

          Case ('SL')
            Call verifyAndSaveProperty_str(trim(key), value, m%strength_min, m%strength_max, m%SL, m%SL_def)

          Case ('GYT')
            Call verifyAndSaveProperty_str(trim(key), value, m%toughness_min, m%toughness_max, m%GYT, m%GYT_def)

          Case ('GSL')
            Call verifyAndSaveProperty_str(trim(key), value, m%toughness_min, m%toughness_max, m%GSL, m%GSL_def)

          Case ('eta_BK')
            Call verifyAndSaveProperty_str(trim(key), value, m%eta_BK_min, m%eta_BK_max, m%eta_BK, m%eta_BK_def)

          Case ('E3')
            Call verifyAndSaveProperty_str(trim(key), value, m%modulus_min, m%modulus_max, m%E3, m%E3_def)

          Case ('G13')
            Call verifyAndSaveProperty_str(trim(key), value, m%modulus_min, m%modulus_max, m%G13, m%G13_def)

          Case ('G23')
            Call verifyAndSaveProperty_str(trim(key), value, m%modulus_min, m%modulus_max, m%G23, m%G23_def)

          Case ('v13')
            Call verifyAndSaveProperty_str(trim(key), value, m%poissons_min, m%poissons_max, m%v13, m%v13_def)

          Case ('alpha11')
            Call verifyAndSaveProperty_str(trim(key), value, m%cte_min, m%cte_max, m%cte(1), m%cte_def(1))

          Case ('alpha22')
            Call verifyAndSaveProperty_str(trim(key), value, m%cte_min, m%cte_max, m%cte(2), m%cte_def(2))

          Case ('alpha_PL')
            Call verifyAndSaveProperty_str(trim(key), value, m%aPL_min, m%aPL_max, m%aPL, m%aPL_def)

          Case ('n_PL')
            Call verifyAndSaveProperty_str(trim(key), value, m%nPL_min, m%nPL_max, m%nPL, m%nPL_def)

          Case ('XT')
            Call verifyAndSaveProperty_str(trim(key), value, m%strength_min, m%strength_max, m%XT, m%XT_def)

          Case ('fXT')
            Call verifyAndSaveProperty_str(trim(key), value, zero, one, m%fXT, m%fXT_def)

          Case ('GXT')
            Call verifyAndSaveProperty_str(trim(key), value, m%toughness_min, m%toughness_max, m%GXT, m%GXT_def)

          Case ('fGXT')
            Call verifyAndSaveProperty_str(trim(key), value, zero, one, m%fGXT, m%fGXT_def)

          Case ('XC')
            Call verifyAndSaveProperty_str(trim(key), value, m%strength_min, m%strength_max, m%XC, m%XC_def)

          Case ('fXC')
            Call verifyAndSaveProperty_str(trim(key), value, Tiny(zero), one, m%fXC, m%fXC_def)

          Case ('GXC')
            Call verifyAndSaveProperty_str(trim(key), value, m%toughness_min, m%toughness_max, m%GXC, m%GXC_def)

          Case ('fGXC')
            Call verifyAndSaveProperty_str(trim(key), value, zero, one, m%fGXC, m%fGXC_def)

          Case ('YC')
            Call verifyAndSaveProperty_str(trim(key), value, m%strength_min, m%strength_max, m%YC, m%YC_def)

          Case ('w_kb')
            Call verifyAndSaveProperty_str(trim(key), value, m%w_kb_min, m%w_kb_max, m%w_kb, m%w_kb_def)

          Case ('cl')
            Call verifyAndSaveProperty_str(trim(key), value, m%cl_min, m%cl_max, m%cl, m%cl_def)

          Case ('alpha0')
            Call verifyAndSaveProperty_str(trim(key), value, m%alpha0_min, m%alpha0_max, m%alpha0, m%alpha0_def)

          Case ('mu')
            Call verifyAndSaveProperty_str(trim(key), value, m%mu_min, m%mu_max, m%mu, m%mu_def)

          Case ('es0')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%es(1), m%es_def(1))
          Case ('es1')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%es(2), m%es_def(2))
          Case ('es2')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%es(3), m%es_def(3))
          Case ('es3')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%es(4), m%es_def(4))

          Case ('gs0')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%gs(1), m%gs_def(1))
          Case ('gs1')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%gs(2), m%gs_def(2))
          Case ('gs2')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%gs(3), m%gs_def(3))
          Case ('gs3')
            Call verifyAndSaveProperty_str(trim(key), value, m%schapery_min, m%schapery_max, m%gs(4), m%gs_def(4))

          Case Default
            Call log%error("loadMatProps: Property not recognized: " // trim(key))
        End Select
      End If
    End Do ReadLines

    ! Close the .props file
    Close(unit)

    ! Require at least the first three material properties to be specified in the inp deck
    If (nprops < 3) Call log%error("loadMatProps: Must define the first three properties in the inp deck")

    ! Load material properties from input deck (at most, first 8)
    If (nprops > 8) Then
      Call loadPropertiesFromInp(8, props)
    Else
      Call loadPropertiesFromInp(nprops, props)
    End If

    Return
  End Subroutine loadPropertiesFromFile

  Subroutine loadPropertiesFromInp(nprops, props)
    Use forlog_Mod

    !Arguments
    Integer, intent(IN) :: nprops
    Double Precision :: props(nprops)

    ! Locals
    Character(len=30) :: tmp, featureFlags
    Integer :: numMissingLeadingZeros
    Double Precision, Parameter :: zero=0.d0, one=1.d0, two=2.d0
    ! -------------------------------------------------------------------- !

    Do i=1, nprops
      Select Case (i)
        Case (1)
          ! Convert to string (remove decimal and trailing zeros)
          write (tmp, *) props(i)
          featureFlags = Adjustl(tmp(1:Index(tmp, '.')-1))
          numMissingLeadingZeros = 6 - Len_trim(featureFlags)
          If (numMissingLeadingZeros > 0) Then
            Do j=1, numMissingLeadingZeros
              featureFlags = '0' // featureFlags
            End Do
          End If

          ! Parse each position into a feature flag
          ! Position 1: Matrix damage
          ! Position 2: Peak-peak nonlinearity (shear plasticity or schapery micro-damage)
          ! Position 3: Fiber tensile damage
          ! Position 4: Fiber compression damage, CDM [Mutually exclusive with Position 5]
          ! Position 5: Fiber compression damage, DGD [Mutually exclusive with Position 4]
          ! Position 6: Friction

          ! Default friction to enabled
          If (Len_trim(featureFlags) < 6) Call log%error("loadMatProps: All 6 feature flags must be specified, only found " // trim(str(Len_trim(featureFlags))))

          ! Set flags
          Do j=1,Len_trim(featureFlags)
            Select Case (j)
              Case (1)
                If (featureFlags(j:j) == '1') Then
                  m%matrixDam = .TRUE.
                  Call log%info("loadMatProps: Matrix damage ENABLED")
                Else
                  m%matrixDam = .FALSE.
                  Call log%info("loadMatProps: Matrix damage DISABLED")
                End If

              Case (2)
                If (featureFlags(j:j) == '1') Then
                  m%shearNonlinearity12 = .TRUE.
                  m%shearNonlinearity13 = .FALSE.
                  m%schapery = .FALSE.
                  Call log%info("loadMatProps: Shear nonlinearity (1-2 plane) ENABLED")
                Else If (featureFlags(j:j) == '2') Then
                  m%shearNonlinearity12 = .FALSE.
                  m%shearNonlinearity13 = .FALSE.
                  m%schapery = .TRUE.
                  Call log%info("loadMatProps: Schapery micro-damage ENABLED")
                Else If (featureFlags(j:j) == '3') Then
                  m%shearNonlinearity12 = .TRUE.
                  m%shearNonlinearity13 = .TRUE.
                  m%schapery = .FALSE.
                  Call log%info("loadMatProps: Shear nonlinearity (3-D) ENABLED")
                Else If (featureFlags(j:j) == '4') Then
                  m%shearNonlinearity12 = .FALSE.
                  m%shearNonlinearity13 = .TRUE.
                  m%schapery = .FALSE.
                  Call log%info("loadMatProps: Shear nonlinearity (1-3 plane) ENABLED")
                Else
                  m%shearNonlinearity12 = .FALSE.
                  m%shearNonlinearity13 = .FALSE.
                  m%schapery = .FALSE.
                  Call log%info("loadMatProps: pre-peak nonlinearity DISABLED")
                End If

              Case (3)
                If (featureFlags(j:j) == '1') Then
                  m%fiberTenDam = .TRUE.
                  Call log%info("loadMatProps: fiber tensile damage ENABLED")
                Else
                  m%fiberTenDam = .FALSE.
                  Call log%info("loadMatProps: fiber tensile damage DISABLED")
                End If

              Case (4)
                If (featureFlags(j:j) == '1') Then
                  m%fiberCompDamBL = .TRUE.
                  m%fiberCompDamFKT = .FALSE.
                  Call log%info("loadMatProps: fiber comp. damage (CDM) ENABLED")
                Else If (featureFlags(j:j) == '2') Then
                  Call log%error("loadMatProps: fiber comp. model 2 not implemented. TODO")
                Else If (featureFlags(j:j) == '3') Then
                  m%fiberCompDamBL = .FALSE.
                  m%fiberCompDamFKT = .TRUE.
                  Call log%info("loadMatProps: fiber comp. damage (DGD) ENABLED")
                Else
                  m%fiberCompDamBL = .FALSE.
                  m%fiberCompDamFKT = .FALSE.
                  Call log%info("loadMatProps: fiber comp. damage DISABLED")
                End If

              Case (5)
                If (featureFlags(j:j) == '1') Then
                  Call log%info("loadMatProps: feature flag position 5 is not used.")
                End If

              Case (6)
                If (featureFlags(j:j) == '1') Then
                  m%friction = .TRUE.
                  Call log%info("loadMatProps: friction ENABLED")
                Else
                  m%friction = .FALSE.
                  Call log%info("loadMatProps: friction DISABLED")
                End If

              Case Default
                Call log%error("loadMatProps: Unknown position found in feature flags")
            End Select
          End Do
        Case (2)
        Case (3)
          m%thickness = props(i)
        ! props(4:8) are reserved for use in the future as needed
        Case (4)
        Case (5)
        Case (6)
        Case (7)
        Case (8)

        Case (9)
          Call verifyAndSaveProperty_double('E1', props(i), m%modulus_min, m%modulus_max, m%E1, m%E1_def)

        Case (10)
          Call verifyAndSaveProperty_double('E2', props(i), m%modulus_min, m%modulus_max, m%E2, m%E2_def)

        Case (11)
          Call verifyAndSaveProperty_double('G12', props(i), m%modulus_min, m%modulus_max, m%G12, m%G12_def)

        Case (12)
          Call verifyAndSaveProperty_double('v12', props(i), m%poissons_min, m%poissons_max, m%v12, m%v12_def)

        Case (13)
          Call verifyAndSaveProperty_double('v23', props(i), m%poissons_min, m%poissons_max, m%v23, m%v23_def)

        Case (14)
          Call verifyAndSaveProperty_double('YT', props(i), m%strength_min, m%strength_max, m%YT, m%YT_def)

        Case (15)
          Call verifyAndSaveProperty_double('SL', props(i), m%strength_min, m%strength_max, m%SL, m%SL_def)

        Case (16)
          Call verifyAndSaveProperty_double('GYT', props(i), m%toughness_min, m%toughness_max, m%GYT, m%GYT_def)

        Case (17)
          Call verifyAndSaveProperty_double('GSL', props(i), m%toughness_min, m%toughness_max, m%GSL, m%GSL_def)

        Case (18)
          Call verifyAndSaveProperty_double('eta_BK', props(i), m%eta_BK_min, m%eta_BK_max, m%eta_BK, m%eta_BK_def)

        Case (19)
          Call verifyAndSaveProperty_double('YC', props(i), m%strength_min, m%strength_max, m%YC, m%YC_def)

        Case (20)
          Call verifyAndSaveProperty_double('alpha0', props(i), m%alpha0_min, m%alpha0_max, m%alpha0, m%alpha0_def)

        Case (21)
          Call verifyAndSaveProperty_double('E3', props(i), m%modulus_min, m%modulus_max, m%E3, m%E3_def)

        Case (22)
          Call verifyAndSaveProperty_double('G13', props(i), m%modulus_min, m%modulus_max, m%G13, m%G13_def)

        Case (23)
          Call verifyAndSaveProperty_double('G23', props(i), m%modulus_min, m%modulus_max, m%G23, m%G23_def)

        Case (24)
          Call verifyAndSaveProperty_double('v13', props(i), m%poissons_min, m%poissons_max, m%v13, m%v13_def)

        Case (25)
          Call verifyAndSaveProperty_double('alpha11', props(i), m%cte_min, m%cte_max, m%cte(1), m%cte_def(1))

        Case (26)
          Call verifyAndSaveProperty_double('alpha22', props(i), m%cte_min, m%cte_max, m%cte(2), m%cte_def(2))

        Case (27)
          Call verifyAndSaveProperty_double('alpha_PL', props(i), m%aPL_min, m%aPL_max, m%aPL, m%aPL_def)

        Case (28)
          Call verifyAndSaveProperty_double('n_PL', props(i), zero, m%nPL_max, m%nPL, m%nPL_def)

        Case (29)
          Call verifyAndSaveProperty_double('XT', props(i), m%strength_min, m%strength_max, m%XT, m%XT_def)

        Case (30)
          Call verifyAndSaveProperty_double('fXT', props(i), zero, one, m%fXT, m%fXT_def)

        Case (31)
          Call verifyAndSaveProperty_double('GXT', props(i), m%toughness_min, m%toughness_max, m%GXT, m%GXT_def)

        Case (32)
          Call verifyAndSaveProperty_double('fGXT', props(i), zero, one, m%fGXT, m%fGXT_def)

        Case (33)
          Call verifyAndSaveProperty_double('XC', props(i), m%strength_min, m%strength_max, m%XC, m%XC_def)

        Case (34)
          Call verifyAndSaveProperty_double('fXC', props(i), Tiny(zero), one, m%fXC, m%fXC_def)

        Case (35)
          Call verifyAndSaveProperty_double('GXC', props(i), m%toughness_min, m%toughness_max, m%GXC, m%GXC_def)

        Case (36)
          Call verifyAndSaveProperty_double('fGXC', props(i), zero, one, m%fGXC, m%fGXC_def)

        Case (37)
          Call verifyAndSaveProperty_double('cl', props(i), m%cl_min, m%cl_max, m%cl, m%cl_def)

        Case (38)
          Call verifyAndSaveProperty_double('w_kb', props(i), m%w_kb_min, m%w_kb_max, m%w_kb, m%w_kb_def)

        Case (39)

        Case (40)
          Call verifyAndSaveProperty_double('mu', props(i), m%mu_min, m%mu_max, m%mu, m%mu_def)

        ! TODO: add schapery micro-damage input properties to the input file parser -- currently only works with
        ! external properties file.

        Case Default
          Call log%error("loadMatProps: Unknown property #"//trim(str(i))//" found. Expecting nprops <= 40")
      End Select

    End Do

    Return
  End Subroutine loadPropertiesFromInp

  Subroutine initializeMinMaxValues()

    ! Locals
    Double Precision, Parameter :: zero=0.d0, one=1.d0, two=2.d0
    Double Precision, parameter :: pi = 3.141592653589793239
    ! -------------------------------------------------------------------- !

    m%modulus_min = Tiny(zero)
    m%modulus_max = Huge(zero)

    m%poissons_min = zero
    m%poissons_max = one

    m%strength_min = Tiny(zero)
    m%strength_max = Huge(zero)

    m%toughness_min = Tiny(zero)
    m%toughness_max = Huge(zero)

    m%eta_BK_min = one
    m%eta_BK_max = Huge(zero)

    m%cte_min = -one
    m%cte_max = one

    m%aPL_min = zero
    m%aPL_max = Huge(zero)

    m%nPL_min = zero
    m%nPL_max = Huge(zero)

    m%w_kb_min = Tiny(zero)
    m%w_kb_max = Huge(zero)

    m%cl_min = zero
    m%cl_max = 1000.d0

    m%alpha0_min = zero
    m%alpha0_max = pi/two

    m%mu_min = zero
    m%mu_max = one

    m%schapery_min = -Huge(zero)
    m%schapery_max = Huge(zero)

    Return
  End Subroutine initializeMinMaxValues

  Subroutine verifyAndSaveProperty_str(key, value, min, max, saveTo, flag)
    ! Checks if the value is within the specified bounds. Prints an error message
    ! which kills the analysis if a value is out of bounds. Otherwise, save the property
    ! and sets the flag to indicate that the property has been read in.

    Use forlog_Mod

    !Arguments
    Character(len=*), intent(IN) :: key, value
    Double Precision, intent(IN) :: min, max
    Double Precision, intent(OUT) :: saveTo
    Logical, intent(OUT) :: flag

    ! Locals
    Double Precision :: valueDbl
    ! -------------------------------------------------------------------- !

    ! Convert to double
    Read(value,*) valueDbl

    ! Verify that the value is within the specified bounds
    If (valueDbl < min) Then
      Call log%error(" PROPERTY ERROR " // trim(key) // " cannot be less than " // trim(str(min)) // ". Found value: " // trim(str(valueDbl)))

    Else If (valueDbl > max) Then
      Call log%error(" PROPERTY ERROR " // trim(key) // " cannot be greater than " // trim(str(max)) // ". Found value: " // trim(str(valueDbl)))

    Else
      Call log%debug(" loadMatProps: Loaded " // trim(key) // " = " // trim(str(valueDbl)))
    End If

    ! Save the value and set the flag
    saveTo = valueDbl
    flag = .TRUE.

    Return
  End Subroutine verifyAndSaveProperty_str

  Subroutine verifyAndSaveProperty_double(key, value, min, max, saveTo, flag)
    ! Checks if the value is within the specified bounds. Prints an error message
    ! which kills the analysis if a value is out of bounds. Otherwise, save the property
    ! and sets the flag to indicate that the property has been read in.

    Use forlog_Mod

    !Arguments
    Character(len=*), intent(IN) :: key
    Double Precision, intent(IN) :: value, min, max
    Double Precision, intent(OUT) :: saveTo
    Logical, intent(OUT) :: flag

    ! Locals
    Double Precision, Parameter :: zero=0.d0
    ! -------------------------------------------------------------------- !

    ! Verify that the value is within the specified bounds
    If (value < min) Then
      If (value == zero) Then
        flag = .FALSE.
        Call log%info(" loadMatProps: assuming " // trim(key) // " is not defined")
        Return
      Else
        Call log%error(" PROPERTY ERROR " // trim(key) // " cannot be less than " // trim(str(min)) // ". Found value: " // trim(str(value)))
      End If

    Else If (value > max) Then
      Call log%error(" PROPERTY ERROR " // trim(key) // " cannot be greater than " // trim(str(max)) // ". Found value: " // trim(str(value)))

    Else
      Call log%debug(" loadMatProps: Loaded " // trim(key) // " = " // trim(str(value)))
    End If

    ! Save the value and set the flag
    saveTo = value
    flag = .TRUE.

    Return
  End Subroutine verifyAndSaveProperty_double

  Subroutine consistencyChecks(m, issueWarnings)
    ! Checks that a consistent set of properties has been defined

    Use forlog_Mod

    ! Arguments
    Type(matProps), intent(INOUT) :: m
    Logical, intent(IN) :: issueWarnings

    ! Locals
    Double Precision, Parameter :: zero=0.d0, one=1.d0, two=2.d0, four=4.d0
    ! -------------------------------------------------------------------- !

    ! Check that all elastic properties have been defined for transverse isotropy
    If (.NOT. m%E1_def) Call log%error('PROPERTY ERROR: Must define a value for E1')
    If (.NOT. m%E2_def) Call log%error('PROPERTY ERROR: Must define a value for E2')
    If (.NOT. m%G12_def) Call log%error('PROPERTY ERROR: Must define a value for G12')
    If (.NOT. m%v12_def) Call log%error('PROPERTY ERROR: Must define a value for v12')
    If (.NOT. m%v23_def) Call log%error('PROPERTY ERROR: Must define a value for v23')

    ! Check if orthotropic elastic properties have been defined
    If (m%E3_def .OR. m%G13_def .OR. m%G23_def) Then
      If (.NOT. m%E3_def) Call log%error('PROPERTY ERROR: Some orthotropic elastic properties are missing. Must define a value for E3.')
      If (.NOT. m%G13_def) Call log%error('PROPERTY ERROR: Some orthotropic elastic properties are missing. Must define a value for G13.')
      If (.NOT. m%G23_def) Call log%error('PROPERTY ERROR: Some orthotropic elastic properties are missing. Must define a value for G23.')
      If (.NOT. m%v13_def) Call log%error('PROPERTY ERROR: Some orthotropic elastic properties are missing. Must define a value for v13.')
      Call log%info('PROPERTY: Orthotropic constants are defined directly')
    Else
      ! Compute these properties assuming transverse isotropy
      Call log%info('PROPERTY: Assuming transverse isotropy')
      m%E3  = m%E2
      m%G13 = m%G12
      m%v13 = m%v12
      m%G23 = m%E2/two/(one + m%v23)
    End If

    ! Check if matrix damage properties have been specified
    If (m%matrixDam) Then
      If (.NOT. m%YT_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for YT.')
      If (.NOT. m%SL_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for SL_def.')
      If (.NOT. m%GYT_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for GYT.')
      If (.NOT. m%GSL_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for GSL.')
      If (.NOT. m%eta_BK_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for eta_BK.')
      If (.NOT. m%YC_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for YC.')
      If (.NOT. m%alpha0_def) Call log%error('PROPERTY ERROR: Some matrix damage properties are missing. Must define a value for alpha0.')
      m%etaL = -m%SL*COS(two*m%alpha0)/(m%YC*COS(m%alpha0)*COS(m%alpha0))
      m%etaT = -one/TAN(two*m%alpha0)
      m%ST   = m%YC*COS(m%alpha0)*(SIN(m%alpha0) + COS(m%alpha0)/TAN(two*m%alpha0))

      Call log%info('PROPERTY: Matrix damage is enabled')
    Else
      Call log%info('PROPERTY: Matrix damage is disabled')
    End If

    ! Check if CTEs have been defined
    If (m%cte_def(1) .OR. m%cte_def(2)) Then
      If (.NOT. m%cte_def(1)) Call log%error('PROPERTY ERROR: Some CTE properties are missing. Must define a value for alpha11.')
      If (.NOT. m%cte_def(2)) Call log%error('PROPERTY ERROR: Some CTE properties are missing. Must define a value for alpha22.')
      Call log%info('PROPERTY: CTEs have been defined')
    Else
      If (issueWarnings) Call log%warn('PROPERTY: CTEs are being set to zero')
      m%cte(1) = zero
      m%cte(2) = zero
    End If

    ! Check if shear nonlinearity properties have been defined
    If (m%shearNonlinearity12 .OR. m%shearNonlinearity13) Then
      If (.NOT. m%aPL_def) Call log%error('PROPERTY ERROR: Some shear-nonlinearity properties are missing. Must define a value for alpha_PL.')
      If (.NOT. m%nPL_def) Call log%error('PROPERTY ERROR: Some shear-nonlinearity properties are missing. Must define a value for n_PL.')
      Call log%info('PROPERTY: Shear-nonlinearity properties have been defined')
    Else
      Call log%info('PROPERTY: Shear-nonlinearity is disabled')
    End If

    ! Check if all Schapery micro-damage material properties have been defined
    If (m%schapery) Then
      If (.NOT. (all(m%es_def) .AND. all(m%gs_def))) Call log%error('PROPERTY ERROR: Some Schapery micro-damage are missing.')
      Call log%info('PROPERTY: Schapery micro-damage properties have been defined')
    Else
      Call log%info('PROPERTY: Schapery micro-damage is disabled')
      ! Default values that will cause no pre-peak nonlinearity
      m%es = zero
      m%es(1) = one
      m%gs = zero
      m%gs(1) = one
    End If

    ! Check if fiber tensile damage properties have been defined
    If (m%fiberTenDam) Then
      If (.NOT. m%XT_def) Call log%error('PROPERTY ERROR: Some fiber tensile damage properties are missing. Must define a value for XT.')
      If (.NOT. m%fXT_def) Call log%error('PROPERTY ERROR: Some fiber tensile damage properties are missing. Must define a value for fXT.')
      If (.NOT. m%GXT_def) Call log%error('PROPERTY ERROR: Some fiber tensile damage properties are missing. Must define a value for GXT.')
      If (.NOT. m%fGXT_def) Call log%error('PROPERTY ERROR: Some fiber tensile damage properties are missing. Must define a value for fGXT.')

      Call log%info('PROPERTY: fiber tensile damage properties have been defined')
    Else
      Call log%info('PROPERTY: fiber tensile damage is disabled')
    End If

    ! Check if the CDM fiber compression damage properties have been defined
    If (m%fiberCompDamBL) Then
      If (.NOT. m%XC_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (BL model). Must define a value for XC.')
      If (.NOT. m%fXC_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (BL model). Must define a value for fXC.')
      If (.NOT. m%GXC_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (BL model). Must define a value for GXC.')
      If (.NOT. m%fGXC_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (BL model). Must define a value for fGXC.')

      Call log%info('PROPERTY: fiber compression damage BL (model 1) properties have been defined')
    Else
      Call log%info('PROPERTY: fiber compression damage BL (model 1) is disabled')
    End If

    ! Check if the DGD fiber compression damage properties have been defined
    If (m%fiberCompDamFKT) Then
      If (.NOT. m%XC_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (FKT model). Must define a value for XC.')
      If (.NOT. m%YC_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (FKT model). Must define a value for YC.')
      If (.NOT. m%w_kb_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (FKT model). Must define a value for w_kb.')
      If (.NOT. m%alpha0_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (FKT model). Must define a value for alpha0.')
      If (.NOT. m%SL_def) Call log%error('PROPERTY ERROR: Some fiber compression damage properties are missing (FKT model). Must define a value for SL.')

      ! Make sure shear nonlinearity is enabled
      If (.NOT. m%shearNonlinearity12) Call log%error('PROPERTY ERROR: In-plane shear-nonlinearity must be enabled with fiber compression damage FKT.')

      Call log%info('PROPERTY: fiber compression damage FKT (model 3) properties have been defined')
    Else
      Call log%info('PROPERTY: fiber compression damage FKT (model 3) is disabled')
    End If

    ! check if fiber nonlinearity has been defined
    If (m%cl_def) Then
      Call log%info('PROPERTY: fiber nonlinearity has been defined')
    Else
      m%cl = zero
    End If

    Return
  End Subroutine consistencyChecks

  Subroutine checkForSnapBack(m, Lc, elementNumber)
    ! Issues a warning if an element greater than the snap back size

    Use forlog_Mod

    ! Arguments
    Type(matProps), intent(IN) :: m
    Double Precision, intent(IN) :: Lc(3)                    ! Element length
    Integer, intent(IN) :: elementNumber                     ! Element number

    ! Locals
    Double Precision :: Lc_fT, Lc_fC, Lc_mT, Lc_SL, Lc_YC    ! Maximum element size for each damage mode
    Double Precision, Parameter :: zero=0.d0, one=1.d0, two=2.d0
    ! -------------------------------------------------------------------- !

    ! Fiber Tension
    If (m%fiberTenDam) Then
      If (two*m%GXT*m%E1 < Lc(1)*m%XT**2) Then
        Lc_fT = two*m%GXT*m%E1/m%XT**2
        Call log%warn("Element " // trim(str(elementNumber)) // " has size " // trim(str(Lc(1))) // " which is > fiber tension snap-back threshold. Set element size < " // trim(str(Lc_fT)))
      End If

      If (two*m%GXT*m%fGXT*m%E1 < Lc(1)*(m%XT*m%fXT)**2) Then
        Lc_fT = two*m%GXT*m%fGXT*m%E1/(m%XT*m%fXT)**2
        Call log%warn("Snap-back in 1st part of fiber tensile softening law. Adjust bilinear ratios or decrease the element size.")
      End If

      If (two*m%GXT*(one - m%fGXT)*m%E1 < Lc(1)*(m%XT*(one - m%fXT))**2) Then
        Lc_fT = two*m%GXT*(one - m%fGXT)*m%E1/(m%XT*(one - m%fXT))**2
        Call log%warn("Snap-back in 2nd part of fiber tensile softening law. Adjust bilinear ratios or decrease the element size.")
      End If
    End If

    ! Fiber Compression
    If (m%fiberCompDamBL) Then
      If (two*m%GXC*m%E1 < Lc(1)*m%XC**2) Then
        Lc_fC = two*m%GXC*m%E1/m%XC**2
        Call log%warn("Element " // trim(str(elementNumber)) // " has size " // trim(str(Lc(1))) // " which is > fiber compression snap-back threshold 1. Set element size < " // trim(str(Lc_fC)))
      End If

      If (two*m%GXC*m%fGXC*m%E1 < Lc(1)*(m%XC*m%fXC)**2) Then
        Lc_fT = two*m%GXC*m%fGXC*m%E1/(m%XC*m%fXC)**2
        Call log%warn("Snap-back in 1st part of fiber compression softening law. Adjust bilinear ratios or decrease the element size.")
      End If

      If (two*m%GXC*(one - m%fGXC)*m%E1 < Lc(1)*(m%XC*(one - m%fXC))**2) Then
        Lc_fT = two*m%GXC*(one - m%fGXC)*m%E1/(m%XC*(one - m%fXC))**2
        Call log%warn("Snap-back in 2nd part of fiber compression softening law. Adjust bilinear ratios or decrease the element size.")
      End If
    End If

    If (m%matrixDam) Then
      ! Matrix Tension
      If (two*m%GYT*m%E2 < Lc(2)*m%YT**2) Then
        Lc_mT = two*m%GYT*m%E2/m%YT/m%YT
        Call log%warn("Element " // trim(str(elementNumber)) // " has size " // trim(str(Lc(2))) // " which is > matrix mode I snap-back threshold. Set element size < " // trim(str(Lc_mT)))
      End If

      ! Matrix Shear
      If (two*m%GSL*m%G12 < Lc(2)*m%SL**2) Then
        Lc_SL = two*m%GSL*m%G12/m%SL**2
        Call log%warn("Element " // trim(str(elementNumber)) // " has size " // trim(str(Lc(2))) // " which is > matrix mode II snap-back threshold. Set element size < " // trim(str(Lc_SL)))
      End If

      ! Matrix Compression
      If (two*m%GSL*m%E2 < Lc(2)*m%YC**2*COS(m%alpha0)) Then
        Lc_YC = two*m%GSL*m%E2/(m%YC**2*COS(m%alpha0))
        Call log%warn("Element " // trim(str(elementNumber)) // " has size " // trim(str(Lc(2))) // " which is > matrix compression snap-back threshold. Set element size < " // trim(str(Lc_YC)))
      End If
    End If

    Return
  End Subroutine checkForSnapBack

  Function initializePhi0(phi0, G12, XC, aPL, nPL, Lc, position)
    ! Determines the phi0 for shear instability

    Use forlog_Mod

    ! Arguments
    Double Precision, intent(IN) :: phi0                         ! Guess for phi0
    Double Precision, intent(IN) :: G12, XC                      ! Material properties
    Double Precision, intent(IN) :: aPL, nPL                     ! Shear nonlinearity
    Double Precision, intent(IN) :: Lc(3)                        ! Element characteristic lengths
    Double Precision, intent(IN) :: position(3)                  ! Position
    Double Precision :: initializePhi0

    ! Locals
    Double Precision :: shearInstabilityPhi0
    Double Precision :: randomNumbers(1000)
    Double Precision :: rn
    Integer :: index, rnsize
    Double Precision, parameter :: zero=0.d0, half=0.5d0, one=1.d0, two=2.d0

    ! Common
    Common randomNumbers
    ! -------------------------------------------------------------------- !

    shearInstabilityPhi0 = (nPL-one)/G12*((G12-XC)/(XC*nPL*aPL**(one/nPL)))**(nPL/(nPL-one))

    If (phi0 == zero) Then  ! Special case to use phi0 corresponding to instability
      ! Calculate phi0
      initializePhi0 = shearInstabilityPhi0

    Else If (phi0 <= half) Then   ! Use initial condition value of phi0
      initializePhi0 = phi0

    Else If (phi0 == one) Then   ! 1-D variation
      ! Set index to be integer count of row
      index = NINT(position(1)/Lc(1))

      ! Wrap in a loop if there are more rows than randomNumbers is long
      rnsize = SIZE(randomNumbers, 1)
      If (index > rnsize) Then
        index = MOD(index, rnsize)
        Call log%warn('Random fiber misalignments are being wrapped in a loop. Increase randomNumberCount in vexternaldb.')
      End If

      ! Get the random number (0 to 1)
      rn = randomNumbers(index)

      ! Scale -phi to +phi
      initializePhi0 = rn*two*shearInstabilityPhi0 - shearInstabilityPhi0

    Else
      Call log%error('Invalid phi0. Found phi0 = '  // trim(str(phi0)) // ' See README.')

    End If

    Return
  End Function initializePhi0

End Module matProp_Mod

ext: for


text:
      SUBROUTINE FD01AD( JP, L, LAMBDA, XIN, YIN, EFOR, XF, EPSBCK,
     $                   CTETA, STETA, YQ, EPOS, EOUT, SALPH, IWARN,
     $                   INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To solve the least-squares filtering problem recursively in time.
C     Each subroutine call implements one time update of the solution.
C     The algorithm uses a fast QR-decomposition based approach.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JP      CHARACTER*1
C             Indicates whether the user wishes to apply both prediction
C             and filtering parts, as follows:
C             = 'B':  Both prediction and filtering parts are to be
C                     applied;
C             = 'P':  Only the prediction section is to be applied.
C
C     Input/Output Parameters
C
C     L       (input) INTEGER
C             The length of the impulse response of the equivalent
C             transversal filter model.  L >= 1.
C
C     LAMBDA  (input) DOUBLE PRECISION
C             Square root of the forgetting factor.
C             For tracking capabilities and exponentially stable error
C             propagation, LAMBDA < 1.0 (strict inequality) should
C             be used.  0.0 < LAMBDA <= 1.0.
C
C     XIN     (input) DOUBLE PRECISION
C             The input sample at instant n.
C             (The situation just before and just after the call of
C             the routine are denoted by instant (n-1) and instant n,
C             respectively.)
C
C     YIN     (input) DOUBLE PRECISION
C             If JP = 'B', then YIN must contain the reference sample
C             at instant n.
C             Otherwise, YIN is not referenced.
C
C     EFOR    (input/output) DOUBLE PRECISION
C             On entry, this parameter must contain the square root of
C             exponentially weighted forward prediction error energy
C             at instant (n-1).  EFOR >= 0.0.
C             On exit, this parameter contains the square root of the
C             exponentially weighted forward prediction error energy
C             at instant n.
C
C     XF      (input/output) DOUBLE PRECISION array, dimension (L)
C             On entry, this array must contain the transformed forward
C             prediction variables at instant (n-1).
C             On exit, this array contains the transformed forward
C             prediction variables at instant n.
C
C     EPSBCK  (input/output) DOUBLE PRECISION array, dimension (L+1)
C             On entry, the leading L elements of this array must
C             contain the normalized a posteriori backward prediction
C             error residuals of orders zero through L-1, respectively,
C             at instant (n-1), and EPSBCK(L+1) must contain the
C             square-root of the so-called "conversion factor" at
C             instant (n-1).
C             On exit, this array contains the normalized a posteriori
C             backward prediction error residuals, plus the square root
C             of the conversion factor at instant n.
C
C     CTETA   (input/output) DOUBLE PRECISION array, dimension (L)
C             On entry, this array must contain the cosines of the
C             rotation angles used in time updates, at instant (n-1).
C             On exit, this array contains the cosines of the rotation
C             angles at instant n.
C
C     STETA   (input/output) DOUBLE PRECISION array, dimension (L)
C             On entry, this array must contain the sines of the
C             rotation angles used in time updates, at instant (n-1).
C             On exit, this array contains the sines of the rotation
C             angles at instant n.
C
C     YQ      (input/output) DOUBLE PRECISION array, dimension (L)
C             On entry, if JP = 'B', then this array must contain the
C             orthogonally transformed reference vector at instant
C             (n-1). These elements are also the tap multipliers of an
C             equivalent normalized lattice least-squares filter.
C             Otherwise, YQ is not referenced and can be supplied as
C             a dummy array (i.e., declare this array to be YQ(1) in
C             the calling program).
C             On exit, if JP = 'B', then this array contains the
C             orthogonally transformed reference vector at instant n.
C
C     EPOS    (output) DOUBLE PRECISION
C             The a posteriori forward prediction error residual.
C
C     EOUT    (output) DOUBLE PRECISION
C             If JP = 'B', then EOUT contains the a posteriori output
C             error residual from the least-squares filter at instant n.
C
C     SALPH   (output) DOUBLE PRECISION array, dimension (L)
C             The element SALPH(i), i=1,...,L, contains the opposite of
C             the i-(th) reflection coefficient for the least-squares
C             normalized lattice predictor (whose value is -SALPH(i)).
C
C     Warning Indicator
C
C     IWARN   INTEGER
C             = 0:  no warning;
C             = 1:  an element to be annihilated by a rotation is less
C                   than the machine precision (see LAPACK Library
C                   routine DLAMCH).
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The output error EOUT at instant n, denoted by EOUT(n), is the
C     reference sample minus a linear combination of L successive input
C     samples:
C
C                           L-1
C        EOUT(n) = YIN(n) - SUM h_i * XIN(n-i),
C                           i=0
C
C     where YIN(n) and XIN(n) are the scalar samples at instant n.
C     A least-squares filter uses those h_0,...,h_{L-1} which minimize
C     an exponentially weighted sum of successive output errors squared:
C
C         n
C        SUM [LAMBDA**(2(n-k)) * EOUT(k)**2].
C        k=1
C
C     Each subroutine call performs a time update of the least-squares
C     filter using a fast least-squares algorithm derived from a
C     QR decomposition, as described in references [1] and [2] (the
C     notation from [2] is followed in the naming of the arrays).
C     The algorithm does not compute the parameters h_0,...,h_{L-1} from
C     the above formula, but instead furnishes the parameters of an
C     equivalent normalized least-squares lattice filter, which are
C     available from the arrays SALPH (reflection coefficients) and YQ
C     (tap multipliers), as well as the exponentially weighted input
C     signal energy
C
C         n                                              L
C        SUM [LAMBDA**(2(n-k)) * XIN(k)**2] = EFOR**2 + SUM XF(i)**2.
C        k=1                                            i=1
C
C     For more details on reflection coefficients and tap multipliers,
C     references [2] and [4] are recommended.
C
C     REFERENCES
C
C     [1]  Proudler, I. K., McWhirter, J. G., and Shepherd, T. J.
C          Fast QRD based algorithms for least-squares linear
C          prediction.
C          Proceedings IMA Conf. Mathematics in Signal Processing
C          Warwick, UK, December 1988.
C
C     [2]  Regalia, P. A., and Bellanger, M. G.
C          On the duality between QR methods and lattice methods in
C          least-squares adaptive filtering.
C          IEEE Trans. Signal Processing, SP-39, pp. 879-891,
C          April 1991.
C
C     [3]  Regalia, P. A.
C          Numerical stability properties of a QR-based fast
C          least-squares algorithm.
C          IEEE Trans. Signal Processing, SP-41, June 1993.
C
C     [4]  Lev-Ari, H., Kailath, T., and Cioffi, J.
C          Least-squares adaptive lattice and transversal filters:
C          A unified geometric theory.
C          IEEE Trans. Information Theory, IT-30, pp. 222-236,
C          March 1984.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires O(L) operations for each subroutine call.
C     It is backward consistent for all input sequences XIN, and
C     backward stable for persistently exciting input sequences,
C     assuming LAMBDA < 1.0 (see [3]).
C     If the condition of the signal is very poor (IWARN = 1), then the
C     results are not guaranteed to be reliable.
C
C     FURTHER COMMENTS
C
C     1.  For tracking capabilities and exponentially stable error
C         propagation, LAMBDA < 1.0 should be used.  LAMBDA is typically
C         chosen slightly less than 1.0 so that "past" data are
C         exponentially forgotten.
C     2.  Prior to the first subroutine call, the variables must be
C         initialized. The following initial values are recommended:
C
C         XF(i) = 0.0,        i=1,...,L
C         EPSBCK(i) = 0.0     i=1,...,L
C         EPSBCK(L+1) = 1.0
C         CTETA(i) = 1.0      i=1,...,L
C         STETA(i) = 0.0      i=1,...,L
C         YQ(i) = 0.0         i=1,...,L
C
C         EFOR = 0.0          (exact start)
C         EFOR = "small positive constant" (soft start).
C
C         Soft starts are numerically more reliable, but result in a
C         biased least-squares solution during the first few iterations.
C         This bias decays exponentially fast provided LAMBDA < 1.0.
C         If sigma is the standard deviation of the input sequence
C         XIN, then initializing EFOR = sigma*1.0E-02 usually works
C         well.
C
C     CONTRIBUTOR
C
C     P. A. Regalia (October 1994).
C     Release 4.0: V. Sima, Katholieke Univ. Leuven, Belgium, Dec. 1999.
C
C     REVISIONS
C
C     -
C
C     KEYWORDS
C
C     Kalman filtering, least-squares estimator, optimal filtering,
C     orthogonal transformation, recursive estimation, QR decomposition.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ONE, ZERO
      PARAMETER         ( ONE = 1.0D0, ZERO = 0.0D0 )
C     .. Scalar Arguments ..
      CHARACTER         JP
      INTEGER           INFO, IWARN, L
      DOUBLE PRECISION  EFOR, EOUT, EPOS, LAMBDA, XIN, YIN
C     .. Array Arguments ..
      DOUBLE PRECISION  CTETA(*), EPSBCK(*), SALPH(*), STETA(*), XF(*),
     $                  YQ(*)
C     .. Local Scalars ..
      LOGICAL           BOTH
      INTEGER           I
      DOUBLE PRECISION  CTEMP, EPS, FNODE, NORM, TEMP, XFI, YQI
C     .. External Functions ..
      LOGICAL           LSAME
      DOUBLE PRECISION  DLAMCH, DLAPY2, DNRM2
      EXTERNAL          DLAMCH, DLAPY2, DNRM2, LSAME
C     .. External Subroutines ..
      EXTERNAL          DLARTG, XERBLA
C     .. Intrinsic Functions
      INTRINSIC         ABS, SQRT
C     .. Executable statements ..
C
C     Test the input scalar arguments.
C
      BOTH  = LSAME( JP, 'B' )
      IWARN = 0
      INFO  = 0
C
      IF( .NOT.BOTH .AND. .NOT.LSAME( JP, 'P' ) ) THEN
         INFO = -1
      ELSE IF( L.LT.1 ) THEN
         INFO = -2
      ELSE IF( ( LAMBDA.LE.ZERO ) .OR. ( LAMBDA.GT.ONE  ) ) THEN
         INFO = -3
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'FD01AD', -INFO )
         RETURN
      END IF
C
C     Computation of the machine precision EPS.
C
      EPS = DLAMCH( 'Epsilon' )
C
C     Forward prediction rotations.
C
      FNODE = XIN
C
      DO 10  I = 1, L
         XFI   = XF(I) * LAMBDA
         XF(I) = STETA(I) * FNODE + CTETA(I) * XFI
         FNODE = CTETA(I) * FNODE - STETA(I) * XFI
   10 CONTINUE
C
      EPOS = FNODE * EPSBCK(L+1)
C
C     Update the square root of the prediction energy.
C
      EFOR = EFOR * LAMBDA
      TEMP = DLAPY2( FNODE, EFOR )
      IF ( TEMP.LT.EPS ) THEN
         FNODE = ZERO
         IWARN = 1
      ELSE
         FNODE = FNODE * EPSBCK(L+1)/TEMP
      END IF
      EFOR = TEMP
C
C     Calculate the reflection coefficients and the backward prediction
C     errors.
C
      DO 20 I = L, 1, -1
         IF ( ABS( XF(I) ).LT.EPS )
     $      IWARN = 1
         CALL DLARTG( TEMP, XF(I), CTEMP, SALPH(I), NORM )
         EPSBCK(I+1) = CTEMP * EPSBCK(I) - SALPH(I) * FNODE
         FNODE = CTEMP * FNODE + SALPH(I) * EPSBCK(I)
         TEMP  = NORM
   20 CONTINUE
C
      EPSBCK(1) = FNODE
C
C     Update to new rotation angles.
C
      NORM = DNRM2( L, EPSBCK, 1 )
      TEMP = SQRT( ( ONE + NORM )*( ONE - NORM ) )
      EPSBCK(L+1) = TEMP
C
      DO 30 I = L, 1, -1
         IF ( ABS( EPSBCK(I) ).LT.EPS )
     $      IWARN = 1
         CALL DLARTG( TEMP, EPSBCK(I), CTETA(I), STETA(I), NORM )
         TEMP = NORM
   30 CONTINUE
C
C     Joint process section.
C
      IF ( BOTH) THEN
         FNODE = YIN
C
         DO 40  I = 1, L
            YQI   = YQ(I) * LAMBDA
            YQ(I) = STETA(I) * FNODE + CTETA(I) * YQI
            FNODE = CTETA(I) * FNODE - STETA(I) * YQI
   40    CONTINUE
C
         EOUT = FNODE * EPSBCK(L+1)
      END IF
C
      RETURN
C *** Last line of FD01AD ***
      END

ext: f

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


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 GSLIBV (STRING)
C=======================================================================

C***********************************************************************

C     *** GSLIBV *** Returns current version number of the suplib library
C                    Get SupLIB Version number
C     Parameters:
C     STRING - OUT - string containing version number of suplib

C     Version Number Format:
C           n1.n2.n3

C     where n1 is the major version number
C           n2 is the minor version or change capability number
C           n3 is the bug-fix number
C***********************************************************************

      CHARACTER*32 STRING

C     Version Format: 'major.minor.bug_fix'
      STRING = '1.2.0'

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2015 NCAR/RAL
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! 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, either version 3 of the License, or
! (at your option) any later version.
!
! 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, see <http://www.gnu.org/licenses/>.

module read_icond_module
USE nrtype
USE netcdf
implicit none
private
public::read_icond
public::read_icond_nlayers
contains

 ! ************************************************************************************************
 ! public subroutine read_icond_nlayers: read model initial conditions file for number of snow/soil layers
 ! ************************************************************************************************
 subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message)
 ! --------------------------------------------------------------------------------------------------------
 ! modules
 USE nrtype
 USE var_lookup,only:iLookIndex                        ! variable lookup structure
 USE globalData,only:gru_struc                         ! gru-hru mapping structures
 USE netcdf_util_module,only:nc_file_close             ! close netcdf file
 USE netcdf_util_module,only:nc_file_open              ! close netcdf file
 USE netcdf_util_module,only:netcdf_err                ! netcdf error handling
 USE data_types,only:gru_hru_intVec                    ! actual data
 USE data_types,only:var_info                          ! metadata
 implicit none

 ! --------------------------------------------------------------------------------------------------------
 ! variable declarations
 ! dummies
 character(*)        ,intent(in)     :: iconFile       ! name of input (restart) file
 integer(i4b)        ,intent(in)     :: nGRU           ! total # of GRUs in run domain
 type(var_info)      ,intent(in)     :: indx_meta(:)   ! metadata
 integer(i4b)        ,intent(out)    :: err            ! error code
 character(*)        ,intent(out)    :: message        ! returned error message

 ! locals
 integer(i4b)             :: ncID                       ! netcdf file id
 integer(i4b)             :: dimID                      ! netcdf file dimension id
 integer(i4b)             :: fileHRU                    ! number of HRUs in netcdf file
 integer(i4b)             :: snowID, soilID             ! netcdf variable ids
 integer(i4b)             :: iGRU, iHRU                 ! loop indexes
 integer(i4b),allocatable :: snowData(:)                ! number of snow layers in all HRUs
 integer(i4b),allocatable :: soilData(:)                ! number of soil layers in all HRUs
 character(len=256)       :: cmessage                   ! downstream error message

 ! --------------------------------------------------------------------------------------------------------
 ! initialize error message
 err=0
 message = 'read_icond_nlayers/'

 ! open netcdf file
 call nc_file_open(iconFile,nf90_nowrite,ncid,err,cmessage);
 if (err/=0) then; message=trim(message)//trim(cmessage); return; end if

 ! get number of HRUs in file
 err = nf90_inq_dimid(ncID,"hru",dimId);               if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if
 err = nf90_inquire_dimension(ncID,dimId,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if

 ! allocate sotrage for reading from file
 allocate(snowData(fileHRU))
 allocate(soilData(fileHRU))
 snowData = 0
 soilData = 0

 ! get variable ids
 err = nf90_inq_varid(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),snowid); call netcdf_err(err,message)
 err = nf90_inq_varid(ncid,trim(indx_meta(iLookIndex%nSoil)%varName),soilid); call netcdf_err(err,message)

 ! get data
 err = nf90_get_var(ncid,snowid,snowData); call netcdf_err(err,message)
 err = nf90_get_var(ncid,soilid,soilData); call netcdf_err(err,message)

 ! assign to index structure - gru by hru
 do iGRU = 1,nGRU
  do iHRU = 1,gru_struc(iGRU)%hruCount
   gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc)
   gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc)
  end do
 end do

 ! close file
 call nc_file_close(ncid,err,cmessage)
 if(err/=0)then;message=trim(message)//trim(cmessage);return;end if

 ! cleanup
 deallocate(snowData,soilData)

 end subroutine read_icond_nlayers

 ! ************************************************************************************************
 ! public subroutine read_icond: read model initial conditions
 ! ************************************************************************************************
 subroutine read_icond(iconFile,                      & ! name of initial conditions file
                       nGRU,                          & ! number of GRUs
                       prog_meta,                     & ! metadata
                       progData,                      & ! model prognostic (state) variables
                       indxData,                      & ! layer index data
                       err,message)                     ! error control
 ! --------------------------------------------------------------------------------------------------------
 ! modules
 USE nrtype
 USE var_lookup,only:iLookVarType                       ! variable lookup structure
 USE var_lookup,only:iLookProg                          ! variable lookup structure
 USE var_lookup,only:iLookIndex                         ! variable lookup structure
 USE globalData,only:gru_struc                          ! gru-hru mapping structures
 USE globaldata,only:iname_soil,iname_snow              ! named variables to describe the type of layer
 USE netcdf_util_module,only:nc_file_close              ! close netcdf file
 USE netcdf_util_module,only:nc_file_open               ! close netcdf file
 USE netcdf_util_module,only:netcdf_err                 ! netcdf error handling
 USE data_types,only:gru_hru_doubleVec                  ! actual data
 USE data_types,only:gru_hru_intVec                     ! actual data
 USE data_types,only:var_info                           ! metadata
 USE get_ixName_module,only:get_varTypeName             ! to access type strings for error messages
 implicit none

 ! --------------------------------------------------------------------------------------------------------
 ! variable declarations
 ! dummies
 character(*)           ,intent(in)     :: iconFile     ! name of netcdf file containing the initial conditions
 integer(i4b)           ,intent(in)     :: nGRU         ! number of grouped response units in simulation domain
 type(var_info)         ,intent(in)     :: prog_meta(:) ! prognostic metadata
 type(gru_hru_doubleVec),intent(inout)  :: progData     ! prognostic vars
 type(gru_hru_intVec)   ,intent(inout)  :: indxData     ! layer indexes
 integer(i4b)           ,intent(out)    :: err          ! error code
 character(*)           ,intent(out)    :: message      ! returned error message

 ! locals
 character(len=256)                     :: cmessage     ! downstream error message
 integer(i4b)                           :: fileHRU      ! number of HRUs in file
 integer(i4b)                           :: iVar         ! loop index
 integer(i4b)                           :: iGRU         ! loop index
 integer(i4b)                           :: iHRU         ! loop index
 integer(i4b)                           :: dimID        ! varible dimension ids
 integer(i4b)                           :: ncVarID      ! variable ID in netcdf file
 character(256)                         :: dimName      ! not used except as a placeholder in call to inq_dim function
 integer(i4b)                           :: dimLen       ! data dimensions
 integer(i4b)                           :: ncID         ! netcdf file ID
 real(dp),allocatable                   :: varData(:,:) ! variable data storage
 integer(i4b)                           :: nSoil, nSnow, nToto ! # layers
 integer(i4b),parameter                 :: nBand=2      ! number of spectral bands

 character(len=32),parameter            :: scalDimName   ='scalarv'  ! dimension name for scalar data
 character(len=32),parameter            :: midSoilDimName='midSoil'  ! dimension name for soil-only layers
 character(len=32),parameter            :: midTotoDimName='midToto'  ! dimension name for layered varaiables
 character(len=32),parameter            :: ifcTotoDimName='ifcToto'  ! dimension name for layered varaiables

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

 ! Start procedure here
 err=0; message="read_icond/"

 ! --------------------------------------------------------------------------------------------------------
 ! (1) read the file
 ! --------------------------------------------------------------------------------------------------------
 ! open netcdf file
 call nc_file_open(iconFile,nf90_nowrite,ncID,err,cmessage)
 if (err/=0) then; message=trim(message)//trim(cmessage); return; end if

 ! get number of HRUs in file
 err = nf90_inq_dimid(ncID,"hru",dimID);               if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if
 err = nf90_inquire_dimension(ncID,dimID,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if

 ! loop through prognostic variables
 do iVar = 1,size(prog_meta)

  ! skip variables that are computed later
  if(prog_meta(iVar)%varName=='scalarCanopyWat'           .or. &
     prog_meta(iVar)%varName=='spectralSnowAlbedoDiffuse' .or. &
     prog_meta(iVar)%varName=='scalarSurfaceTemp'         .or. &
     prog_meta(iVar)%varName=='mLayerVolFracWat'          .or. &
     prog_meta(iVar)%varName=='mLayerHeight'                   ) cycle

  ! get variable id
  err = nf90_inq_varid(ncID,trim(prog_meta(iVar)%varName),ncVarID); call netcdf_err(err,message)
  if(err/=0)then
   message=trim(message)//': problem with getting variable id, var='//trim(prog_meta(iVar)%varName)
   return
  endif

  ! get variable dimension IDs
  select case (prog_meta(iVar)%varType)
   case (iLookVarType%scalarv); err = nf90_inq_dimid(ncID,trim(scalDimName)   ,dimID); call netcdf_err(err,message)
   case (iLookVarType%midSoil); err = nf90_inq_dimid(ncID,trim(midSoilDimName),dimID); call netcdf_err(err,message)
   case (iLookVarType%midToto); err = nf90_inq_dimid(ncID,trim(midTotoDimName),dimID); call netcdf_err(err,message)
   case (iLookVarType%ifcToto); err = nf90_inq_dimid(ncID,trim(ifcTotoDimName),dimID); call netcdf_err(err,message)
   case default
    message=trim(message)//"unexpectedVariableType[name='"//trim(prog_meta(iVar)%varName)//"';type='"//trim(get_varTypeName(prog_meta(iVar)%varType))//"']"
    err=20; return
  end select

  ! check errors
  if(err/=0)then
   message=trim(message)//': problem with dimension ids, var='//trim(prog_meta(iVar)%varName)
   return
  endif

  ! get the dimension length
  err = nf90_inquire_dimension(ncID,dimID,dimName,dimLen); call netcdf_err(err,message)
  if(err/=0)then; message=trim(message)//': problem getting the dimension length'; return; endif

  ! iniitialize the variable data
  allocate(varData(fileHRU,dimLen),stat=err)
  if(err/=0)then; message=trim(message)//'problem allocating variable data'; return; endif

  ! get data
  err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message)
  if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif

  ! check data are not set to the fill value
  if( any( abs(varData - nf90_fill_double) < epsilon(varData) ) )then
   message=trim(message)//"data set to the fill value (name='"//trim(prog_meta(iVar)%varName)//"')"
   err=20; return
  endif

  ! store data in prognostics structure
  ! loop through GRUs
  do iGRU = 1,nGRU
   do iHRU = 1,gru_struc(iGRU)%hruCount

    ! get ther number of layers
    nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow
    nSoil = gru_struc(iGRU)%hruInfo(iHRU)%nSoil
    nToto = nSnow + nSoil

    ! put the data into data structures
    select case (prog_meta(iVar)%varType)
     case (iLookVarType%scalarv); progData%gru(iGRU)%hru(iHRU)%var(iVar)%dat(1        ) = varData(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc,1        )
     case (iLookVarType%midSoil); progData%gru(iGRU)%hru(iHRU)%var(iVar)%dat(1:nSoil  ) = varData(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc,1:nSoil  )
     case (iLookVarType%midToto); progData%gru(iGRU)%hru(iHRU)%var(iVar)%dat(1:nToto  ) = varData(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc,1:nToto  )
     case (iLookVarType%ifcToto); progData%gru(iGRU)%hru(iHRU)%var(iVar)%dat(0:nToto  ) = varData(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc,1:nToto+1)
     case default
      message=trim(message)//"unexpectedVariableType[name='"//trim(prog_meta(iVar)%varName)//"';type='"//trim(get_varTypeName(prog_meta(iVar)%varType))//"']"
      err=20; return
    end select

    ! initialize the spectral albedo
    progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1:nBand) = progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1)

   end do ! iHRU
  end do ! iGRU

  ! deallocate storage vector for next variable
  deallocate(varData, stat=err)
  if(err/=0)then; message=trim(message)//'problem deallocating variable data'; return; endif

 end do ! iVar

 ! --------------------------------------------------------------------------------------------------------
 ! (2) set number of layers
 ! --------------------------------------------------------------------------------------------------------
 do iGRU = 1,nGRU
  do iHRU = 1,gru_struc(iGRU)%hruCount

   ! save the number of layers
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%nSnow)%dat(1)   = gru_struc(iGRU)%hruInfo(iHRU)%nSnow
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%nSoil)%dat(1)   = gru_struc(iGRU)%hruInfo(iHRU)%nSoil
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%nLayers)%dat(1) = gru_struc(iGRU)%hruInfo(iHRU)%nSnow + gru_struc(iGRU)%hruInfo(iHRU)%nSoil

   ! initalize the indices for midSnow, midSoil, midToto, and ifcToto
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%midSnowStartIndex)%dat(1) = 1
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%midSoilStartIndex)%dat(1) = 1
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%midTotoStartIndex)%dat(1) = 1
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%ifcSnowStartIndex)%dat(1) = 1
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%ifcSoilStartIndex)%dat(1) = 1
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%ifcTotoStartIndex)%dat(1) = 1

   ! set layer type
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%layerType)%dat(1:gru_struc(iGRU)%hruInfo(iHRU)%nSnow) = iname_snow
   indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%layerType)%dat((gru_struc(iGRU)%hruInfo(iHRU)%nSnow+1):(gru_struc(iGRU)%hruInfo(iHRU)%nSnow+gru_struc(iGRU)%hruInfo(iHRU)%nSoil)) = iname_soil

  end do
 end do

 end subroutine read_icond

end module read_icond_module

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


! { dg-do run }
! { dg-add-options ieee }
! { dg-skip-if "NaN not supported" { spu-*-* } }
!
! List-directed part of PR fortran/43298
! and follow up to PR fortran/34319.
!
! Check handling of "NAN(alphanum)"
!
character(len=200) :: str
real :: r
complex :: z

! read_real:

r = 1.0
str = 'INfinity' ; read(str,*) r
if (r < 0 .or. r /= r*1.1) STOP 1

r = 1.0
str = '-INF' ; read(str,*) r
if (r > 0 .or. r /= r*1.1) STOP 2

r = 1.0
str = '+INF' ; read(str,*) r
if (r < 0 .or. r /= r*1.1) STOP 3

r = 1.0
str = '-inFiniTY' ; read(str,*) r
if (r > 0 .or. r /= r*1.1) STOP 4

r = 1.0
str = 'NAN' ; read(str,*) r
if (.not. isnan(r)) STOP 5

r = 1.0
str = '-NAN' ; read(str,*) r
if (.not. isnan(r)) STOP 6

r = 1.0
str = '+NAN' ; read(str,*) r
if (.not. isnan(r)) STOP 7

r = 1.0
str = 'NAN(0x111)' ; read(str,*) r
if (.not. isnan(r)) STOP 8

r = 1.0
str = '-NAN(123)' ; read(str,*) r
if (.not. isnan(r)) STOP 9

r = 1.0
str = '+NAN(0xFFE)' ; read(str,*) r
if (.not. isnan(r)) STOP 10

! parse_real

z = cmplx(-2.0,-4.0)
str = '(0.0,INfinity)' ; read(str,*) z
if (aimag(z) < 0 .or. aimag(z) /= aimag(z)*1.1) STOP 11

z = cmplx(-2.0,-4.0)
str = '(-INF,0.0)' ; read(str,*) z
if (real(z) > 0 .or. real(z) /= real(z)*1.1) STOP 12

z = cmplx(-2.0,-4.0)
str = '(0.0,+INF)' ; read(str,*) z
if (aimag(z) < 0 .or. aimag(z) /= aimag(z)*1.1) STOP 13

z = cmplx(-2.0,-4.0)
str = '(-inFiniTY,0.0)' ; read(str,*) z
if (real(z) > 0 .or. real(z) /= real(z)*1.1) STOP 14

z = cmplx(-2.0,-4.0)
str = '(NAN,0.0)' ; read(str,*) z
if (.not. isnan(real(z))) STOP 15

z = cmplx(-2.0,-4.0)
str = '(0.0,-NAN)' ; read(str,*) z
if (.not. isnan(aimag(z))) STOP 16

z = cmplx(-2.0,-4.0)
str = '(+NAN,0.0)' ; read(str,*) z
if (.not. isnan(real(z))) STOP 17

z = cmplx(-2.0,-4.0)
str = '(NAN(0x111),0.0)' ; read(str,*) z
if (.not. isnan(real(z))) STOP 18

z = cmplx(-2.0,-4.0)
str = '(0.0,-NaN(123))' ; read(str,*) z
if (.not. isnan(aimag(z))) STOP 19

z = cmplx(-2.0,-4.0)
str = '(+nan(0xFFE),0.0)' ; read(str,*) z
if (.not. isnan(real(z))) STOP 20

end


text:
C$Procedure      BODMAT ( Return transformation matrix for a body )

      SUBROUTINE BODMAT ( BODY, ET, TIPM )

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

      IMPLICIT NONE

      INTEGER               BODY
      DOUBLE PRECISION      ET
      DOUBLE PRECISION      TIPM   ( 3,3 )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     BODY       I   ID code of body.
C     ET         I   Epoch of transformation.
C     TIPM       O   Transformation from Inertial to PM for BODY at ET.
C
C$ Detailed_Input
C
C     BODY        is the integer ID code of the body for which the
C                 transformation is requested. Bodies are numbered
C                 according to the standard NAIF numbering scheme.
C
C     ET          is the epoch at which the transformation is
C                 requested. (This is typically the epoch of
C                 observation minus the one-way light time from
C                 the observer to the body at the epoch of
C                 observation.)
C
C$ Detailed_Output
C
C     TIPM        is the transformation matrix from Inertial to body
C                 Equator and Prime Meridian.  The X axis of the PM
C                 system is directed to the intersection of the
C                 equator and prime meridian. The Z axis points north.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If data required to define the body-fixed frame associated
C        with BODY are not found in the binary PCK system or the kernel
C        pool, the error SPICE(FRAMEDATANOTFOUND) is signaled. In
C        the case of IAU style body-fixed frames, the absence of
C        prime meridian polynomial data (which are required) is used
C        as an indicator of missing data.
C
C     2) If the test for exception (1) passes, but in fact requested
C        data are not available in the kernel pool, the error will be
C        signaled by routines in the call tree of this routine.
C
C     3) If the kernel pool does not contain all of the data required
C        to define the number of nutation precession angles
C        corresponding to the available nutation precession
C        coefficients, the error SPICE(INSUFFICIENTANGLES) is
C        signaled.
C
C     4) If the reference frame REF is not recognized, a routine
C        called by BODMAT will diagnose the condition and invoke the
C        SPICE error handling system.
C
C     5) If the specified body code BODY is not recognized, the
C        error is diagnosed by a routine called by BODMAT.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine is related to the more general routine TIPBOD
C     which returns a matrix that transforms vectors from a
C     specified inertial reference frame to body equator and
C     prime meridian coordinates.  TIPBOD accepts an input argument
C     REF that allows the caller to specify an inertial reference
C     frame.
C
C     The transformation represented by BODMAT's output argument TIPM
C     is defined as follows:
C
C        TIPM = [W] [DELTA] [PHI]
C                 3        1     3
C
C     If there exists high-precision binary PCK kernel information
C     for the body at the requested time, these angles, W, DELTA
C     and PHI are computed directly from that file.  The most
C     recently loaded binary PCK file has first priority followed
C     by previously loaded binary PCK files in backward time order.
C     If no binary PCK file has been loaded, the text P_constants
C     kernel file is used.
C
C     If there is only text PCK kernel information, it is
C     expressed in terms of RA, DEC and W (same W as above), where
C
C        RA    = PHI - HALFPI()
C        DEC   = HALFPI() - DELTA
C
C     RA, DEC, and W are defined as follows in the text PCK file:
C
C           RA  = RA0  + RA1*T  + RA2*T*T   + a  sin theta
C                                              i          i
C
C           DEC = DEC0 + DEC1*T + DEC2*T*T  + d  cos theta
C                                              i          i
C
C           W   = W0   + W1*d   + W2*d*d    + w  sin theta
C                                              i          i
C
C     where:
C
C           d = days past J2000.
C
C           T = Julian centuries past J2000.
C
C           a , d , and w  arrays apply to satellites only.
C            i   i       i
C
C           theta  = THETA0 * THETA1*T are specific to each planet.
C                i
C
C     These angles -- typically nodal rates -- vary in number and
C     definition from one planetary system to the next.
C
C$ Examples
C
C     In the following code fragment, BODMAT is used to rotate
C     the position vector (POS) from a target body (BODY) to a
C     spacecraft from inertial coordinates to body-fixed coordinates
C     at a specific epoch (ET), in order to compute the planetocentric
C     longitude (PCLONG) of the spacecraft.
C
C        CALL BODMAT ( BODY, ET, TIPM )
C        CALL MXV    ( TIPM, POS, POS )
C        CALL RECLAT ( POS, RADIUS, PCLONG, LAT )
C
C     To compute the equivalent planetographic longitude (PGLONG),
C     it is necessary to know the direction of rotation of the target
C     body, as shown below.
C
C        CALL BODVCD ( BODY, 'PM', 3, DIM, VALUES )
C
C        IF ( VALUES(2) .GT. 0.D0 ) THEN
C           PGLONG = PCLONG
C        ELSE
C           PGLONG = TWOPI() - PCLONG
C        END IF
C
C     Note that the items necessary to compute the transformation
C     TIPM must have been loaded into the kernel pool (by one or more
C     previous calls to FURNSH).
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     1)  Refer to the NAIF_IDS required reading file for a complete
C         list of the NAIF integer ID codes for bodies.
C
C$ Author_and_Institution
C
C     N.J. Bachman    (JPL)
C     B.V. Semenov    (JPL)
C     W.L. Taber      (JPL)
C     I.M. Underwood  (JPL)
C     K.S. Zukor      (JPL)
C
C$ Version
C
C-    SPICELIB Version 4.2.0, 27-JUL-2016 (BVS)
C
C        Updated to use the 3x3 top-left corner of the 6x6 matrix
C        returned by TISBOD instead of fetching kernel data and doing
C        computations in-line.
C
C-    SPICELIB Version 4.1.1, 01-FEB-2008 (NJB)
C
C        The routine was updated to improve the error messages created
C        when required PCK data are not found. Now in most cases the
C        messages are created locally rather than by the kernel pool
C        access routines. In particular missing binary PCK data will
C        be indicated with a reasonable error message.
C
C-    SPICELIB Version 4.1.0, 25-AUG-2005 (NJB)
C
C        Updated to remove non-standard use of duplicate arguments
C        in MXM call.
C
C         Calls to ZZBODVCD have been replaced with calls to
C         BODVCD.
C
C-     SPICELIB Version 4.0.0, 12-FEB-2004 (NJB)
C
C         Code has been updated to support satellite ID codes in the
C         range 10000 to 99999 and to allow nutation precession angles
C         to be associated with any object.
C
C         Implementation changes were made to improve robustness
C         of the code.
C
C-     SPICELIB Version 3.2.0, 22-MAR-1995 (KSZ)
C
C        Gets TSIPM matrix from PCKMAT (instead of Euler angles
C        from PCKEUL.)
C
C-     SPICELIB Version 3.0.0, 10-MAR-1994 (KSZ)
C
C        Ability to get Euler angles from binary PCK file added.
C        This uses the new routine PCKEUL.
C
C-     SPICELIB Version 2.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 2.0.0, 04-SEP-1991 (NJB)
C
C         Updated to handle P_constants referenced to different epochs
C         and inertial reference frames.
C
C         The header was updated to specify that the inertial reference
C         frame used by BODMAT is restricted to be J2000.
C
C-    SPICELIB Version 1.0.0, 31-JAN-1990 (WLT) (IMU)
C
C-&

C$ Index_Entries
C
C     fetch transformation matrix for a body
C     transformation from j2000 position to bodyfixed
C     transformation from j2000 to bodyfixed coordinates
C
C-&

C$ Revisions
C
C-    SPICELIB Version 4.2.0, 02-MAR-2016 (BVS)
C
C        Updated to use the 3x3 top-left corner of the 6x6 matrix
C        returned by TISBOD instead of fetching kernel data and doing
C        computations in-line.
C
C-    SPICELIB Version 4.1.0, 25-AUG-2005 (NJB)
C
C        Updated to remove non-standard use of duplicate arguments
C        in MXM call.
C
C         Calls to ZZBODVCD have been replaced with calls to
C         BODVCD.
C
C-     SPICELIB Version 4.0.0, 12-FEB-2004 (NJB)
C
C         Code has been updated to support satellite ID codes in the
C         range 10000 to 99999 and to allow nutation precession angles
C         to be associated with any object.
C
C         Calls to deprecated kernel pool access routine RTPOOL
C         were replaced by calls to GDPOOL.
C
C         Calls to BODVAR have been replaced with calls to
C         ZZBODVCD.
C
C-     SPICELIB Version 3.2.0, 22-MAR-1995 (KSZ)
C
C        BODMAT now get the TSIPM matrix from PCKMAT, and
C        unpacks TIPM from it.  Also the calculated but unused
C        variable LAMBDA was removed.
C
C-     SPICELIB Version 3.0.0, 10-MAR-1994 (KSZ)
C
C        BODMAT now uses new software to check for the
C        existence of binary PCK files, search the for
C        data corresponding to the requested body and time,
C        and return the appropriate Euler angles, using the
C        new routine PCKEUL.  Otherwise the code calculates
C        the Euler angles from the P_constants kernel file.
C
C-     SPICELIB Version 2.0.0, 04-SEP-1991 (NJB)
C
C         Updated to handle P_constants referenced to different epochs
C         and inertial reference frames.
C
C         The header was updated to specify that the inertial reference
C         frame used by BODMAT is restricted to be J2000.
C
C         BODMAT now checks the kernel pool for presence of the
C         variables
C
C            BODY#_CONSTANTS_REF_FRAME
C
C         and
C
C            BODY#_CONSTANTS_JED_EPOCH
C
C         where # is the NAIF integer code of the barycenter of a
C         planetary system or of a body other than a planet or
C         satellite.  If either or both of these variables are
C         present, the P_constants for BODY are presumed to be
C         referenced to the specified inertial frame or epoch.
C         If the epoch of the constants is not J2000, the input
C         time ET is converted to seconds past the reference epoch.
C         If the frame of the constants is not J2000, the rotation from
C         the P_constants' frame to body-fixed coordinates is
C         transformed to the rotation from J2000 coordinates to
C         body-fixed coordinates.
C
C         For efficiency reasons, this routine now duplicates much
C         of the code of BODEUL so that it doesn't have to call BODEUL.
C         In some cases, BODEUL must covert Euler angles to a matrix,
C         rotate the matrix, and convert the result back to Euler
C         angles.  If this routine called BODEUL, then in such cases
C         this routine would convert the transformed angles back to
C         a matrix.  That would be a bit much....
C
C
C-    Beta Version 1.1.0, 16-FEB-1989 (IMU) (NJB)
C
C        Examples section completed.  Declaration of unused variable
C        FOUND removed.
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local variables
C
      DOUBLE PRECISION      TSIPM   ( 6, 6 )

      INTEGER               I
      INTEGER               J

C
C     Standard SPICE Error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'BODMAT' )
      END IF

C
C     Get 6x6 state transformation from TISBOD. If succeeded, pull out
C     left-top 3x3 matrix.
C
      CALL TISBOD ( 'J2000', BODY, ET, TSIPM )

      IF ( FAILED() ) THEN
         CALL CHKOUT ( 'BODMAT' )
         RETURN
      END IF

      DO  I = 1, 3
         DO   J  = 1, 3
            TIPM  ( I, J ) = TSIPM  ( I, J )
         END DO
      END DO

      CALL CHKOUT ( 'BODMAT' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
module covlocal
!$$$  module documentation block
!
! module: covlocal                     functions for Gaspari-Cohn localization,
!                                      blending northern, tropical and southern
!                                      hemisphere parameter values.
!
! prgmmr: whitaker         org: esrl/psd               date: 2009-02-23
!
! abstract: two functions, one defining the Gaspari-Cohn polynomial for
!  compact localization (taper) and one to smoothly blend parameter
!  values defined separately for the northern hemisphere, southern
!  hemisphere and tropical regions as a function of latitude.
!
! Public Functions:
!  taper: Compute the localization function defined by Gaspari and Cohn, 1999, QJRMS,
!   p. 723-757, eqn 4.10.
!  latval: blend parameters defined as separate NH, SH and tropical constants.
!   The blending is done linearly in latitude band defined by the
!   parameters latbound and delat.
!
! Public Variables: None
!
! program history log:
!   2009-02-23  Initial version.
!
! attributes:
!   language: f95
!
!$$$

use kinds, only : r_single

private
public :: taper, latval

real(r_single) a1,a2,a3,a4,a5,a6,a7,a8,a9
parameter(a1 = -8.0_r_single)
parameter(a2 = 8.0_r_single)
parameter(a3 = 5.0_r_single)
parameter(a4 = 20.0_r_single/3.0_r_single)
parameter(a5 = 1.0_r_single)
parameter(a6 = 8.0_r_single/3.0_r_single)
parameter(a7 = -10.0_r_single)
parameter(a8 =  4.0_r_single)
parameter(a9 = -1.0_r_single/3.0_r_single)

contains

real(r_single) function taper(r)
 ! Gaspari-Cohn taper function.
 ! r should be positive, and normalized so taper = 0 at r = 1
 ! very close to exp(-(r/c)**2), where c = 0.388
 implicit none
 real(r_single), intent(in) :: r
 if(r < a5)then
   if(r > 0.5_r_single)then
      taper = ( ( ( ( a6*r -a2 )*r +a3 )*r +a4 )*r +a7)*r + a8 + a9/r
   else
      taper = ( ( ( a1*r +a2)*r +a3 )*r -a4)*r*r + a5
   end if
 else
    taper = 0._r_single
 end if
end function taper

real(r_single) function latval(deglat,valnh,valtr,valsh)
 ! blend valnh, valtr and valsh (nh, tropical and sh values)
 ! so they vary linearly from one to another in the latitude
 ! bands defined by latbound and delat.
 ! Inputs:
 !  delat: latitude (in degrees)
 !  valnh,valtr,valsh: NH, tropical and SH parameter values.
 ! Returns:  the parameter value corresponding the specified
 !           latitude.
 use params, only: latboundpp,latboundpm,latboundmp,latboundmm, delatinv
 !  latbound:  the latitude that defines the poleward extent
 !             of the "tropics".
 !  delat:  the transition zone (in degree latitude) over
 !          which the parameters change linearly from one
 !          value to another.
 implicit none
 real(r_single), intent(in) :: deglat,valnh,valtr,valsh
 if (deglat > latboundpp) then
    latval = valnh
 else if (deglat >= latboundpm) then
    latval = ((latboundpp-deglat)*valtr + (deglat-latboundpm)*valnh)*delatinv
 else if (deglat > latboundmp) then
    latval = valtr
 else if (deglat >= latboundmm) then
    latval = ((latboundpp+deglat)*valtr + (-deglat+latboundmp)*valsh)*delatinv
 else
    latval = valsh
 end if
end function latval

end module covlocal

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text:
program ch0907
    ! Simple e Edit Descriptor Example
    implicit none

    integer :: i
    real :: r1 = 1.23456
    real :: r2 = 1.23456

    do i = 1, 10
        print 100, i, r1, r2
        r1 = r1/10.
        r2 = r2*10.
    end do

100 format (' ', i3, ' ', e12.4, ' ', e12.4)
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
      real*8 function tmr()

      real*8 second

c      tmr = second
      tmr = 0

      return
      end


text:
Name: Steve Pop

Phone Number: (916) 4477878

Office: 2801 J St., Sacramento

Email: spop (at) golyon.com

Personality: Steve is a realtor with Lyon Real Estate. He operates out of Sacramento, but often works in Davis and throughout Yolo County. Hes been an agent for over 20 years, and is familiar with every aspect of the real estate business. Once you get in touch with him through his office number listed here, hes always quickly reachable by cell phone or email.

20090206 14:52:06 nbsp My fiance and I have been working with Steve Pop in the downtown office. He is the most phenomenal real estate agent Ive ever met! We are purchasing our first home, but actually went through two agents before Steve. The first agent tried to unload his garbage onto us, and the second agent was never available to show us houses we found. In a fast moving market (such as we have now in Sacramento and Yolo County), you just cant afford to have an agent who isnt at the top of his game.  Please note that the other agent was not with Lyon.

When our second agent refused for the umpteenth time to show us a property we liked, we called Lyon and asked for someone who would always be available to show us places. They connected us to Steve and we saw the property that day. He was incredibly informative and actually went beyond the usual agent. he had the insight of an inspector and the financial knowledge of a broker: he not only knew where to look for dry rot or cracked foundations, he could also tell us how much we could raise the value of the house by, say, adding wood floors, or what kind of steps we should take to lock in our interest rate.

And when we were outbid (by an investor paying in cash) on the first house, he was there to show us the second. And when we lost the second house because people bid tens of thousands over the askingSteve was there to show us the third. Honestly, rain or shine, this guy would meet us anywhere.

In the end, we even switched from Sacramento County to Yolo county and he didnt even blink (even though this lengthened his drive to the houses we wanted to see). Both his flexibility and availability got us our first home...the exact home we had imagined.

And at the final stages of closing, he is still here...calling occasionally to ask if we need anything or check on our stress levels. He has said time and time again that we could call him in the middle of the night if we start to freak out. Hes an angel.

If youre really ready to buy a house and are ready to go at it full onthis is your man. I would recommend him a thousand times over. Actually, I made this yelp account so people wouldnt have to fumble around with bad agents like we had in the beginning.

I have no doubt that if you need someone to sell your house, you couldnt ask for a better agent. He will go to amazing lengths to make sure youre a happy client.

Amazing! Fast, helpful, caring, and fun too! This couldnt have been a better experience. I mean...I was ignored by some agents for only having about a 200k budget, but Steve respected me and treated me like I had a 2 million dollar budget. Amazing! Id definitely recommend him to anyone looking to buy, sell, or flip a house anywhere in the area. Users/natalies
 This is not yelp.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search