text:
!***********************************************************************
!  Integrated Water Flow Model (IWFM)
!  Copyright (C) 2005-2021
!  State of California, Department of Water Resources
!
!  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 2
!  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.
!  (http://www.gnu.org/copyleft/gpl.html)
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
!  For tecnical support, e-mail: IWFMtechsupport@water.ca.gov
!***********************************************************************
MODULE Class_AppLake_v40
  USE Class_Version                , ONLY: ReadVersion
  USE MessageLogger                , ONLY: SetLastMessage                , &
                                           EchoProgress                  , &
                                           MessageArray                  , &
                                           iFatal
  USE GeneralUtilities             , ONLY: StripTextUntilCharacter       , &
                                           IntToText                     , &
                                           FirstLocation                 , &
                                           ShellSort                     , &
                                           LocateInList                  , &
                                           EstablishAbsolutePathFilename , &
                                           CleanSpecialCharacters        , &
                                           ConvertID_To_Index
  USE TimeSeriesUtilities          , ONLY: TimeStepType
  USE IOInterface                  , ONLY: GenericFileType
  USE Package_Misc                 , ONLY: PairedDataType                , &
                                           f_iFlowDest_Outside           , &
                                           f_iFlowDest_StrmNode          , &
                                           f_iFlowDest_Lake              , &
                                           f_iLakeComp
  USE Package_Discretization       , ONLY: AppGridType                   , &
                                           StratigraphyType
  USE Package_Budget               , ONLY: BudgetHeaderType
  USE Package_ComponentConnectors  , ONLY: LakeGWConnectorType           , &
                                           StrmLakeConnectorType         , &
                                           f_iStrmToLakeFlow             , &
                                           f_iBypassToLakeFlow           , &
                                           f_iLakeToStrmFlow
  USE Class_Lake                   , ONLY: LakeType                      , &
                                           ReadInitialLakeElevs
  USE Class_MaxLakeElevFile        , ONLY: MaxLakeElevFileType
  USE Package_PrecipitationET      , ONLY: ETType                        , &
                                           PrecipitationType
  USE Package_Matrix               , ONLY: MatrixType
  USE Class_BaseAppLake            , ONLY: BaseAppLakeType               , &
                                           PrepareLakeBudgetHeader       , &
                                           GenerateRatingTable
  IMPLICIT NONE

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** VARIABLE DEFINITIONS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- PUBLIC ENTITIES
  ! -------------------------------------------------------------
  PRIVATE
  PUBLIC :: AppLake_v40_Type

  ! -------------------------------------------------------------
  ! --- APPLICATION LAKES DATA TYPE
  ! -------------------------------------------------------------
  TYPE,EXTENDS(BaseAppLakeType) :: AppLake_v40_Type
      PRIVATE
      INTEGER,ALLOCATABLE                   :: iColMaxElev(:)        !Pointer to data column in the maximum lake elevation data file
      TYPE(MaxLakeElevFileType),ALLOCATABLE :: MaxLakeElevFile
  CONTAINS
      PROCEDURE,PASS :: SetStaticComponent              => AppLake_v40_SetStaticComponent
      PROCEDURE,PASS :: SetStaticComponentFromBinFile   => AppLake_v40_SetStaticComponentFromBinFile
      PROCEDURE,PASS :: SetDynamicComponent             => AppLake_v40_SetDynamicComponent
      PROCEDURE,PASS :: SetAllComponents                => AppLake_v40_SetAllComponents
      PROCEDURE,PASS :: SetAllComponentsWithoutBinFile  => AppLake_v40_SetAllComponentsWithoutBinFile
      PROCEDURE,PASS :: KillImplementation              => AppLake_v40_Kill
      PROCEDURE,PASS :: GetVersion                      => AppLake_v40_GetVersion
      PROCEDURE,PASS :: Simulate                        => AppLake_v40_Simulate
      PROCEDURE,PASS :: ReadTSData                      => AppLake_v40_ReadTSData
      PROCEDURE,PASS :: CheckExternalTSDataPointers     => AppLake_v40_CheckExternalTSDataPointers
      PROCEDURE,PASS :: ConvertTimeUnit                 => AppLake_v40_ConvertTimeUnit
  END TYPE AppLake_v40_Type

  ! -------------------------------------------------------------
  ! --- LAKE PACKAGE VERSION RELATED DATA
  ! -------------------------------------------------------------
  INTEGER,PARAMETER                    :: iVersion    = 40
  INTEGER,PARAMETER                    :: iLenVersion = 8
  CHARACTER(LEN=iLenVersion),PARAMETER :: cVersion    ='4.0.0000'
  INCLUDE 'AppLake_v40_Revision.fi'

  ! -------------------------------------------------------------
  ! --- MISC. DATA
  ! -------------------------------------------------------------
  INTEGER,PARAMETER                   :: ModNameLen = 19
  CHARACTER(LEN=ModNameLen),PARAMETER :: ModName    = 'Class_AppLake_v40::'

CONTAINS

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** CONSTRUCTORS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- NEW LAKE FROM RAW DATA (GENERALLY CALLED IN PRE-PROCESSOR)
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_SetStaticComponent(AppLake,cFileName,Stratigraphy,AppGrid,StrmLakeConnector,LakeGWConnector,iStat)
    CLASS(AppLake_v40_Type),INTENT(OUT)   :: AppLake
    CHARACTER(LEN=*),INTENT(IN)           :: cFileName
    TYPE(StratigraphyType),INTENT(IN)     :: Stratigraphy
    TYPE(AppGridType),INTENT(IN)          :: AppGrid
    TYPE(StrmLakeConnectorType)           :: StrmLakeConnector
    TYPE(LakeGWConnectorType),INTENT(OUT) :: LakeGWConnector
    INTEGER,INTENT(OUT)                   :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+30) :: ThisProcedure = ModName // 'AppLake_v40_SetStaticComponent'
    INTEGER                      :: NLakes,ErrorCode,indxLake,DummyArray(5),ID,iDestType,NElements, &
                                    indxElem,iDestID,iElem,indxLake1,iElemIDs(AppGrid%NElements),   &
                                    iDest
    INTEGER,ALLOCATABLE          :: iElems_Work(:),iLakeIDs(:)
    INTEGER,PARAMETER            :: f_iDestTypes(3) = [f_iFlowDest_Outside , f_iFlowDest_StrmNode , f_iFlowDest_Lake]
    CHARACTER(:),ALLOCATABLE     :: cVersion
    TYPE(GenericFileType)        :: InFile

    !Initialize
    iStat    = 0
    iElemIDs = AppGrid%AppElement%ID

    !Return if filename is empty
    IF (cFileName .EQ. '') RETURN

    !Echo progress
    CALL EchoProgress('Instantiating static component of application lakes')

    !Open file
    CALL InFile%New(FileName=cFileName,InputFile=.TRUE.,IsTSFile=.FALSE.,Descriptor='pre-processor lake data file',iStat=iStat)
    IF (iStat .EQ. -1) RETURN

    !Read the first line that holds version number
    CALL ReadVersion(InFile,'LAKE',cVersion,iStat)
    IF (iStat .EQ. -1) RETURN

    !Number of lakes
    CALL InFile%ReadData(NLakes,iStat)  ;  IF(iStat .EQ. -1) RETURN
    AppLake%NLakes = NLakes

    !Allocate memory
    ALLOCATE (AppLake%Lakes(NLakes) , STAT=ErrorCode)
    IF (ErrorCode .NE. 0) THEN
        CALL SetLastMessage('Error in allocating memory for lakes!',iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF
    IF (NLakes .EQ. 0) RETURN

    !Read lake elements and outflow destination
    DO indxLake=1,NLakes
        CALL InFile%ReadData(DummyArray,iStat)  ;  IF(iStat .EQ. -1) RETURN

        !Lake ID
        ID                         = DummyArray(1)
        AppLake%Lakes(indxLake)%ID = ID

        !Make sure lake ID is not used more than once
        DO indxLake1=1,indxLake-1
            IF (ID .EQ. AppLake%Lakes(indxLake1)%ID) THEN
                CALL SetLastMessage('Lake ID '//TRIM(IntToText(ID))//' is used more than once!',iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
        END DO

        !Outflow destination type
        iDestType = DummyArray(2)
        IF (.NOT. ANY(iDestType.EQ.f_iDestTypes)) THEN
            CALL SetLastMessage('Outflow destination type for lake '//TRIM(IntToText(ID))//' is not recognized!',iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
        AppLake%Lakes(indxLake)%OutflowDestType = iDestType

        !Outflow destination
        AppLake%Lakes(indxLake)%OutflowDest = DummyArray(3)

        !Lake elements
        NElements                         = DummyArray(4)
        AppLake%Lakes(indxLake)%NElements = NElements
        DEALLOCATE (iElems_Work , STAT=ErrorCode)
        ALLOCATE (AppLake%Lakes(indxLake)%Elements(NElements) , iElems_Work(NElements))
        iElems_Work(1) = DummyArray(5)
        DO indxElem=2,NElements
            CALL InFile%ReadData(iElems_Work(indxElem),iStat)
            IF(iStat .EQ. -1) RETURN
        END DO
        CALL ConvertID_To_Index(iElems_Work,iElemIDs,AppLake%Lakes(indxLake)%Elements)
        IF (ANY(AppLake%Lakes(indxLake)%Elements.EQ.0)) THEN
            CALL SetLastMessage('One or more elements listed for lake '//TRIM(IntToText(ID))//' are not in the model!',iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
        CALL ShellSort(AppLake%Lakes(indxLake)%Elements)

        !Make sure lake elements are not listed for more than one lake
        DO indxElem=1,NElements
            iElem = AppLake%Lakes(indxLake)%Elements(indxElem)
            DO indxLake1=1,indxLake-1
                IF (ANY(iElem .EQ. AppLake%Lakes(indxLake1)%Elements)) THEN
                    CALL SetLastMessage('Element '//TRIM(IntToText(iElemIDs(iElem)))//' listed for lake '//TRIM(IntToText(ID))//' is also listed for lake '//TRIM(IntToText(AppLake%Lakes(indxLake1)%ID))//'!',iFatal,ThisProcedure)
                    iStat = -1
                    RETURN
                END IF
            END DO
        END DO

        !Lake nodes
        CALL AppLake%Lakes(indxLake)%CompileLakeNodes(AppGrid,iStat)
        IF (iStat .EQ. -1) RETURN

        !Rating table
        CALL GenerateRatingTable(AppGrid,AppLake%Lakes(indxLake)%Elements,AppLake%Lakes(indxLake)%Nodes,Stratigraphy%GSElev,AppLake%Lakes(indxLake)%RatingTable,iStat)
        IF (iStat .EQ. -1) RETURN

        !Lake node areas
        CALL AppLake%Lakes(indxLake)%ComputeLakeNodeAreas(AppGrid)

        !Add GW nodes to LakeGWConnector
        CALL LakeGWConnector%AddGWNodes(indxLake,AppLake%Lakes(indxLake)%Elements,AppGrid,Stratigraphy)

    END DO

    !Process outflow destinations
    ALLOCATE (iLakeIDs(NLakes))
    iLakeIDs = AppLake%Lakes%ID
    DO indxLake=1,NLakes
        ID      = iLakeIDs(indxLake)
        iDestID = AppLake%Lakes(indxLake)%OutflowDest
        SELECT CASE (AppLake%Lakes(indxLake)%OutflowDestType)
            !If flow to stream, make sure that stream node is modeled
            CASE (f_iFlowDest_StrmNode)
                AppLake%Lakes(indxLake)%OutflowDest = iDestID
                CALL StrmLakeConnector%AddData(f_iLakeToStrmFlow,indxLake,iDestID)

            !If flow to lake, make sure downstream lake is modeled
            CASE (f_iFlowDest_Lake)
                IF (iDestID .EQ. ID) THEN
                    CALL SetLastMessage('Outflow from lake '//TRIM(IntToText(ID))//' cannot flow into itself!',iFatal,ThisProcedure)
                    iStat = -1
                    RETURN
                END IF
                iDest = LocateInList(iDestID,iLakeIDs)
                IF (iDest .EQ. 0) THEN
                    CALL SetLastMessage('Lake '//TRIM(IntToText(iDestID))//' as outflow destination for lake '//TRIM(IntToText(ID))//' is not in the model!',iFatal,ThisProcedure)
                    iStat = -1
                    RETURN
                END IF
                AppLake%Lakes(indxLake)%OutflowDest = iDest
        END SELECT
    END DO

    !Close lake data file
    CALL InFile%Kill()

  END SUBROUTINE AppLake_v40_SetStaticComponent

  ! -------------------------------------------------------------
  ! --- NEW STATIC LAKE DATA FROM PRE-PROCESSOR BINARY FILE
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_SetStaticComponentFromBinFile(AppLake,BinFile,iStat)
    CLASS(AppLake_v40_Type),INTENT(OUT) :: AppLake
    TYPE(GenericFileType)               :: BinFile
    INTEGER,INTENT(OUT)                 :: iStat

    !Initialize
    iStat = 0

    !Read the preprocessed data for lakes
    CALL AppLake%ReadPreprocessedData(BinFile,iStat)

  END SUBROUTINE AppLake_v40_SetStaticComponentFromBinFile

  ! -------------------------------------------------------------
  ! --- INSTANTIATE THE DYNAMIC COMPONENT OF LAKE DATA (GENERALLY CALLED IN SIMULATION)
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_SetDynamicComponent(AppLake,IsForInquiry,cFileName,cWorkingDirectory,TimeStep,NTIME,AppGrid,LakeGWConnector,iStat)
    CLASS(AppLake_v40_Type)       :: AppLake
    LOGICAL,INTENT(IN)            :: IsForInquiry
    CHARACTER(LEN=*),INTENT(IN)   :: cFileName,cWorkingDirectory
    TYPE(TimeStepType),INTENT(IN) :: TimeStep
    INTEGER,INTENT(IN)            :: NTIME
    TYPE(AppGridType),INTENT(IN)  :: AppGrid
    TYPE(LakeGWConnectorType)     :: LakeGWConnector
    INTEGER,INTENT(OUT)           :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+31) :: ThisProcedure = ModName // 'AppLake_v40_SetDynamicComponent'
    TYPE(GenericFileType)        :: LakeDataFile
    CHARACTER                    :: ALine*2000,TimeUnitConductance*6,cLakeBudgetFileName*2000
    TYPE(BudgetHeaderType)       :: BudHeader
    REAL(8)                      :: FactK,DummyArray(6),FactL,CLAKE,DLAKE,FactC
    INTEGER                      :: indxLake,ID,indx,iLoc,iLakeIDs(AppLake%NLakes),iLake
    CHARACTER(:),ALLOCATABLE     :: cVersion,cAbsPathFileName
    LOGICAL                      :: lProcessed(AppLake%NLakes)

    !Initialize
    iStat = 0

    !Return if no file name is specified
    IF (cFileName .EQ. '') RETURN

    !Echo progress
    CALL EchoProgress('Instantiating dynamic component of application lakes')

    !Open lake data file
    CALL LakeDataFile%New(cFileName,InputFile=.TRUE.,iStat=iStat)
    IF (iStat .EQ. -1) RETURN

    !Read first line that stores the version number
    CALL ReadVersion(LakeDataFile,'LAKE',cVersion,iStat)
    IF (iStat .EQ. -1) RETURN

    !Maximum lake elevation data file
    CALL LakeDataFile%ReadData(ALine,iStat)  ;  IF(iStat .EQ. -1) RETURN
    ALine = StripTextUntilCharacter(ALine,'/')
    CALL CleanSpecialCharacters(ALine)
    IF (ALine .EQ. '') THEN
        CALL SetLastMessage('Maximum lake elevations file must be specified when lakes are simulated!',iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF
    ALLOCATE (AppLake%MaxLakeElevFile , AppLake%iColMaxElev(AppLake%NLakes))
    CALL EstablishAbsolutePathFileName(TRIM(ADJUSTL(ALine)),cWorkingDirectory,cAbsPathFileName)
    CALL AppLake%MaxLakeElevFile%New(cAbsPathFileName,cWorkingDirectory,TimeStep,iStat)
    IF (iStat .EQ. -1) RETURN

    !Lake budget raw file name
    CALL LakeDataFile%ReadData(cLakeBudgetFileName,iStat)  ;  IF(iStat .EQ. -1) RETURN
    cLakeBudgetFileName = StripTextUntilCharacter(cLakeBudgetFileName,'/')
    CALL CleanSpecialCharacters(cLakeBudgetFileName)
    CALL EstablishAbsolutePathFileName(TRIM(ADJUSTL(cLakeBudgetFileName)),cWorkingDirectory,cAbsPathFileName)
    cLakeBudgetFileName = cAbsPathFileName

    !Final lake elevations output file
    CALL LakeDataFile%ReadData(ALine,iStat)  ;  IF(iStat .EQ. -1) RETURN  ;  ALine = ADJUSTL(StripTextUntilCharacter(ALine,'/'))  ;  CALL CleanSpecialCharacters(ALine)
    IF (ALine .NE. '') THEN
        ALLOCATE (AppLake%FinalElevFile)
        CALL EstablishAbsolutePathFileName(TRIM(ADJUSTL(ALine)),cWorkingDirectory,cAbsPathFileName)
        IF (IsForInquiry) THEN
            CALL AppLake%FinalElevFile%New(FileName=cAbsPathFileName,InputFile=.TRUE.,IsTSFile=.FALSE.,Descriptor='end-of-simulation lake elevations data',iStat=iStat)
        ELSE
            CALL AppLake%FinalElevFile%New(FileName=cAbsPathFileName,InputFile=.FALSE.,IsTSFile=.FALSE.,Descriptor='end-of-simulation lake elevations data',iStat=iStat)
        END IF
        IF (iStat .EQ. -1) RETURN
        AppLake%lFinalElevFile_Defined = .TRUE.
    END IF

    !Read factors
    CALL LakeDataFile%ReadData(FactK,iStat)  ;  IF(iStat .EQ. -1) RETURN
    CALL LakeDataFile%ReadData(ALine,iStat)  ;  IF(iStat .EQ. -1) RETURN
    CALL CleanSpecialCharacters(ALine)
    TimeUnitConductance = ADJUSTL(StripTextUntilCharacter(ALine,'/'))
    CALL LakeDataFile%ReadData(FactL,iStat)  ;  IF(iStat .EQ. -1) RETURN

    !Read lake parameters
    lProcessed = .FALSE.
    iLakeIDs   = AppLake%Lakes%ID
    DO indxLake=1,AppLake%NLakes
        CALL LakeDataFile%ReadData(ALine,iStat)  ;  IF(iStat .EQ. -1) RETURN
        READ (ALine,*) DummyArray

        !Make sure that lake ID number is recognized
        ID = INT(DummyArray(1))
        CALL ConvertID_To_Index(ID,iLakeIDs,iLake)
        IF (iLake .EQ. 0) THEN
            CALL SetLastMessage('Lake ID '//TRIM(IntToText(ID))//' listed for lake parameters is not recognized!',iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF

        !Make sure lake data was not entered previously
        IF (lProcessed(iLake)) THEN
            CALL SetLastMessage('Parameters for lake '//TRIM(IntToText(ID))//' are entered more than once!',iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
        lProcessed(iLake) = .TRUE.

        !Process parameters
        CLAKE                               = DummyArray(2)*FactK
        DLAKE                               = DummyArray(3)*FactL
        AppLake%iColMaxElev(iLake)       = INT(DummyArray(4))
        AppLake%Lakes(iLake)%iColET      = INT(DummyArray(5))
        AppLake%Lakes(iLake)%iColPrecip  = INT(DummyArray(6))
        FactC                               = CLAKE / DLAKE

        !Compile lake-gw connector
        CALL LakeGWConnector%SetConductance(AppGrid,iLake,TimeUnitConductance,FactC,iStat)
        IF (iStat .EQ. -1) RETURN

        !Make sure that iColMaxElev is consistent with the data columns in the file
        IF (AppLake%iColMaxElev(iLake) .GT. AppLake%MaxLakeElevFile%iSize) THEN
            MessageArray(1) = 'Maximum lake elevation data column for lake '//TRIM(IntToText(ID))//' is greater than the'
            MessageArray(2) = 'available data columns in the Maximum Lake Elevations Data File!'
            CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF

        !Extract the name of the lake from ALine
        CALL CleanSpecialCharacters(ALine)
        ALine = ADJUSTL(ALine)
        DO indx=1,SIZE(DummyArray)
          iLoc  = FirstLocation(' ',ALine)
          ALine = ADJUSTL(ALine(iLoc+1:))
        END DO
        AppLake%Lakes(iLake)%cName = StripTextUntilCharacter(ADJUSTL(ALine),'/')

    END DO

    !Instantiate the lake budget file
    IF (cLakeBudgetFileName .NE. '') THEN
        ALLOCATE (AppLake%LakeBudRawFile)
        IF (IsForInquiry) THEN
            CALL AppLake%LakeBudRawFile%New(TRIM(cLakeBudgetFileName),iStat)
            IF (iStat .EQ. -1) RETURN
        ELSE
            BudHeader = PrepareLakeBudgetHeader(AppLake%Lakes,NTIME,TimeStep,AppLake%GetVersion())
            CALL AppLake%LakeBudRawFile%New(TRIM(cLakeBudgetFileName),BudHeader,iStat)
            IF (iStat .EQ. -1) RETURN
            CALL BudHeader%Kill()
        END IF
        AppLake%LakeBudRawFile_Defined = .TRUE.
    END IF

    !Initial lake elevations
    CALL ReadInitialLakeElevs(LakeDataFile,AppLake%Lakes,iLakeIDs,iStat)

    !Close file
    CALL LakeDataFile%Kill()

  END SUBROUTINE AppLake_v40_SetDynamicComponent

  ! -------------------------------------------------------------
  ! --- INSTANTIATE COMPLETE LAKE DATA
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_SetAllComponents(AppLake,IsForInquiry,cFileName,cSimWorkingDirectory,TimeStep,NTIME,AppGrid,BinFile,LakeGWConnector,iStat)
    CLASS(AppLake_v40_Type),INTENT(OUT) :: AppLake
    LOGICAL,INTENT(IN)                  :: IsForInquiry
    CHARACTER(LEN=*),INTENT(IN)         :: cFileName,cSimWorkingDirectory
    TYPE(TimeStepType),INTENT(IN)       :: TimeStep
    INTEGER,INTENT(IN)                  :: NTIME
    TYPE(AppGridType),INTENT(IN)        :: AppGrid
    TYPE(GenericFileType)               :: BinFile
    TYPE(LakeGWConnectorType)           :: LakeGWConnector
    INTEGER,INTENT(OUT)                 :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+28) :: ThisProcedure = ModName // 'AppLake_v40_SetAllComponents'

    !Initialize
    iStat = 0

    !Read the preprocessed data for lakes
    CALL AppLake%ReadPreprocessedData(BinFile,iStat)
    IF (iStat .EQ. -1) RETURN

    !Return if there are no lakes defined in Pre-Processor
    IF (AppLake%NLakes .EQ. 0) RETURN

    !Set the dynamic part of AppLake
    CALL AppLake%SetDynamicComponent(IsForInquiry,cFileName,cSimWorkingDirectory,TimeStep,NTIME,AppGrid,LakeGWConnector,iStat)
    IF (iStat .EQ. -1) RETURN

    !Make sure that if static part is defined, so is the dynamic part
    IF (AppLake%NLakes .GT. 0) THEN
        IF (cFileName .EQ. '') THEN
            MessageArray(1) = 'For proper simulation of lakes, relevant lake data files must'
            MessageArray(2) = 'be specified when lakes are defined in Pre-Processor.'
            CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
    END IF

  END SUBROUTINE AppLake_v40_SetAllComponents

  ! -------------------------------------------------------------
  ! --- INSTANTIATE COMPLETE LAKE DATA WITHOUT INTERMEDIATE BINARY FILE
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_SetAllComponentsWithoutBinFile(AppLake,IsForInquiry,cPPFileName,cSimFileName,cSimWorkingDirectory,AppGrid,Stratigraphy,TimeStep,NTIME,StrmLakeConnector,LakeGWConnector,iStat)
    CLASS(AppLake_v40_Type),INTENT(OUT)   :: AppLake
    LOGICAL,INTENT(IN)                    :: IsForInquiry
    CHARACTER(LEN=*),INTENT(IN)           :: cPPFileName,cSimFileName,cSimWorkingDirectory
    TYPE(AppGridType),INTENT(IN)          :: AppGrid
    TYPE(StratigraphyType),INTENT(IN)     :: Stratigraphy
    TYPE(TimeStepType),INTENT(IN)         :: TimeStep
    INTEGER,INTENT(IN)                    :: NTIME
    TYPE(StrmLakeConnectorType)           :: StrmLakeConnector
    TYPE(LakeGWConnectorType),INTENT(OUT) :: LakeGWConnector
    INTEGER,INTENT(OUT)                   :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+43) :: ThisProcedure = ModName // 'AppLake_v40_SetAllComponentsWithoutBinFile'

    !Initialize
    iStat = 0

    !Instantiate the static components of the AppLake data
    CALL AppLake%SetStaticComponent(cPPFileName,Stratigraphy,AppGrid,StrmLakeConnector,LakeGWConnector,iStat)
    IF (iStat .EQ. -1) RETURN

    !Instantiate the dynamic component of the AppLake data
    CALL AppLake%SetDynamicComponent(IsForInquiry,cSimFileName,cSimWorkingDirectory,TimeStep,NTIME,AppGrid,LakeGWConnector,iStat)
    IF (iStat .EQ. -1) RETURN

    !Make sure that if static part is defined, so is the dynamic part
    IF (AppLake%NLakes .GT. 0) THEN
        IF (AppLake%MaxLakeElevFile%iSize .EQ. 0) THEN
            MessageArray(1) = 'For proper simulation of lakes, relevant lake data files must'
            MessageArray(2) = 'be specified when lakes are defined in Pre-Processor.'
            CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
    END IF

  END SUBROUTINE AppLake_v40_SetAllComponentsWithoutBinFile

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** DESTRUCTOR
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- KILL LAKE DATA OBJECT
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_Kill(AppLake)
    CLASS(AppLake_v40_Type) :: AppLake

    !Local variables
    INTEGER :: ErrorCode

    DEALLOCATE (AppLake%iColMaxElev , STAT=ErrorCode)

    !Kill maximum lake elevation data file object
    IF (ALLOCATED(AppLake%MaxLakeElevFile)) THEN
        CALL AppLake%MaxLakeElevFile%Kill()
        DEALLOCATE (AppLake%MaxLakeElevFile , STAT=ErrorCode)
    END IF

  END SUBROUTINE AppLake_v40_Kill

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** GETTERS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- GET VERSION NUMBER
  ! -------------------------------------------------------------
  FUNCTION AppLake_v40_GetVersion(AppLake) RESULT(cVrs)
    CLASS(AppLake_v40_Type)  :: AppLake
    CHARACTER(:),ALLOCATABLE :: cVrs

    IF (.NOT. AppLake%Version%IsDefined())   &
        AppLake%Version = AppLAke%Version%New(iLenVersion,cVersion,cRevision)

    cVrs = AppLake%Version%GetVersion()

  END FUNCTION AppLake_v40_GetVersion

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** DATA READERS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- READ TIME SERIES  DATA
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_ReadTSData(AppLake,TimeStep,ET,Precip,iStat)
    CLASS(AppLake_v40_Type)            :: AppLake
    TYPE(TimeStepType),INTENT(IN)      :: TimeStep
    TYPE(ETType),INTENT(IN)            :: ET
    TYPE(PrecipitationType),INTENT(IN) :: Precip
    INTEGER,INTENT(OUT)                :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+22) :: ThisProcedure = ModName // 'AppLake_v40_ReadTSData'
    INTEGER                      :: indxLake,ID
    REAL(8)                      :: rLowGSElev,rMaxElev

    !Initialzie
    iStat = 0

    !Lake precipitation
    IF (Precip%IsUpdated()) AppLake%Lakes%PrecipRate = Precip%GetValues(AppLake%Lakes%iColPrecip)

    !Lake potential evapotranspiration
    IF (ET%IsUpdated()) AppLake%Lakes%ETp_Rate = ET%GetValues(AppLake%Lakes%iColET)

    !Maximum lake elevation
    CALL AppLake%MaxLakeElevFile%ReadTSData(TimeStep,iStat)
    IF (iStat .EQ. -1) RETURN

    !Make sure that maximum lake elevation is greater than the minimum ground surface elevation
    IF (AppLake%MaxLakeElevFile%lUpdated) THEN
        AppLake%Lakes%MaxElev = AppLake%MaxLakeElevFile%rValues(AppLake%iColMaxElev)
        DO indxLake=1,AppLake%NLakes
            rLowGSElev = AppLake%Lakes(indxLake)%RatingTable%XPoint(1)
            rMaxElev   = AppLake%Lakes(indxLake)%MaxElev
            IF (rLowGSElev .GT. rMaxElev) THEN
                ID              = AppLake%Lakes(indxLake)%ID
                MessageArray(1) = 'Maximum lake elevation at lake '//TRIM(IntToText(ID))//' is lower than lake bottom!'
                WRITE (MessageArray(2),'(A,F6.2)') 'Maximum lake elevation = ',rMaxElev
                WRITE (MessageArray(3),'(A,F6.2)') 'Lake bottom elevation  = ',rLowGSElev
                CALL SetLastMessage(MessageArray(1:3),iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
        END DO
    END IF

  END SUBROUTINE AppLake_v40_ReadTSData

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** MISC. METHODS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- SIMULATE LAKES
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_Simulate(AppLake,GSElevs,GWHeads,Runoff,ReturnFlow,LakeGWConnector,StrmLakeConnector,Matrix)
    CLASS(AppLake_v40_Type)              :: AppLake
    REAL(8),INTENT(IN)                   :: GSElevs(:),GWHeads(:,:),Runoff(:),ReturnFlow(:)
    TYPE(LakeGWConnectorType),INTENT(IN) :: LakeGWConnector
    TYPE(StrmLakeConnectorType)          :: StrmLakeConnector
    TYPE(MatrixType)                     :: Matrix

    !Local variables
    INTEGER                           :: indxLake,NLakes,iDest,iNodeIDs(1)
    REAL(8),DIMENSION(AppLake%NLakes) :: StrmInflows,LakeGWFlow_AtMaxLakeElev
    REAL(8)                           :: Elev,rOutflow,rUpdateValues(1),rUpdateRHS(AppLake%NLakes),OtherInflow
    INTEGER,PARAMETER                 :: iCompIDs(1) = [f_iLakeComp]

    !Initialize
    NLakes                     =  AppLake%NLakes
    AppLake%Lakes%InflowUplake =  0.0
    AppLake%Lakes%Outflow      =  0.0
    CALL StrmLakeConnector%ResetLakeToStrmFlows()

    !Compute inflows from streams
    DO indxLake=1,NLakes
      StrmInflows(indxLake) = StrmLakeConnector%GetFlow(f_iStrmToLakeFlow,indxLake)  &
                            + StrmLakeConnector%GetFlow(f_iBypassToLakeFlow,indxLake)
    END DO

    !Lake-gw interaction at maximum lake elevations
    LakeGWFlow_AtMaxLakeElev= LakeGWConnector%ComputeLakeGWFlow(NLakes,AppLake%Lakes%MaxElev,GWHeads,GSElevs)

    !Actual ET
    CALL AppLake%ComputeLakeETa(GSElevs,StrmInflows)

    !Compile lake equations
    DO indxLake=1,NLakes
      !Initialize
      Elev        = AppLake%Lakes(indxLake)%Elev                                    !Lake elevation
      iNodeIDs(1) = indxLake                                                        !Lake to lake connectivity node
      OtherInflow = StrmInflows(indxLake) + Runoff(indxLake) + ReturnFlow(indxLake) !Inflows into lake

      !Lake storage
      AppLake%Lakes(indxLake)%Storage = MAX(AppLake%Lakes(indxLake)%RatingTable%Evaluate(Elev) , 0.0)

      !RHS vector (compute lake outflow if necessary)
      rOutflow                        = RHSLake(AppLake%Lakes(indxLake),AppLake%Lakes(indxLake)%MaxElev,OtherInflow) + LakeGWFlow_AtMaxLakeElev(indxLake)
      rOutflow                        = -MIN(rOutflow , 0.0)
      rUpdateRHS(indxLake)            = RHSLake(AppLake%Lakes(indxLake),Elev,OtherInflow) + rOutflow
      AppLake%Lakes(indxLake)%Outflow = rOutflow

      !Send outflow to destination
      IF (rOutflow .GT. 0.0) THEN
          iDest = AppLake%Lakes(indxLake)%OutflowDest
          SELECT CASE (AppLake%Lakes(indxLake)%OutflowDestType)
              CASE (f_iFlowDest_Lake)
                  AppLake%Lakes(iDest)%InflowUpLake = AppLake%Lakes(iDest)%InflowUpLake + rOutflow
              CASE (f_iFlowDest_StrmNode)
                  CALL StrmLakeConnector%SetFlow(f_iLakeToStrmFlow,indxLake,iDest,rOutflow)
          END SELECT
      END IF

      !Update COEFF matrix
      rUpdateValues(1) = AppLake%Lakes(indxLake)%RatingTable%Derivative(Elev)
      CALL Matrix%UpdateCOEFF(f_iLakeComp,indxLake,1,iCompIDs,iNodeIDs,rUpdateValues)

    END DO

    !Update RHS vector
    CALL Matrix%UpdateRHS(f_iLakeComp,1,rUpdateRHS)

  END SUBROUTINE AppLake_v40_Simulate

  ! -------------------------------------------------------------
  ! --- COMPUTE RIGHT-HAND-SIDE OF A LAKE EQUATION
  ! -------------------------------------------------------------
  FUNCTION RHSLake(Lake,HLake,OtherInflow) RESULT(RHS)
    TYPE(LakeType),INTENT(IN) :: Lake
    REAL(8),INTENT(IN)        :: HLake,OtherInflow
    REAL(8)                   :: RHS

    !Local variables
    REAL(8) :: Storage

    !Compute storage at HLake
    Storage = MAX(Lake%RatingTable%Evaluate(HLake) , 0.0)

    RHS =  Storage                      & !Storage
         - Lake%Storage_P               & !Storage at the previous time step
         - Lake%PrecipRate * Lake%Area  & !Precipitation
         + Lake%ETa                     & !Evaporation
         - OtherInflow                  & !All other inflows to lake
         - Lake%InflowUpLake              !Inflow from upstream lake

  END FUNCTION RHSLake

  ! -------------------------------------------------------------
  ! --- MAKE SURE THAT POINTED TIME-SERIES DATA HAVE ENOUGH COLUMNS
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_CheckExternalTSDataPointers(AppLake,Precip,ET,iStat)
    CLASS(AppLake_v40_Type),INTENT(IN) :: AppLake
    TYPE(PrecipitationType),INTENT(IN) :: Precip
    TYPE(ETType),INTENT(IN)            :: ET
    INTEGER,INTENT(OUT)                :: iStat

    !Local variables
    CHARACTER(Len=ModNameLen+39) :: ThisProcedure = ModName // 'AppLake_v40_CheckExternalTSDataPointers'
    INTEGER                      :: iLake(1),ID

    !Initialize
    iStat = 0

    !Check precip columns
    IF (Precip%GetNDataColumns() .LT. MAXVAL(AppLake%Lakes%iColPrecip)) THEN
        iLake = MAXLOC(AppLake%Lakes%iColPrecip)
        ID    = AppLake%Lakes(iLake(1))%ID
        MessageArray(1) = 'Precipitation data column for lake '//TRIM(IntToText(ID))//' is greater than the'
        MessageArray(2) = 'available data columns in the Precipitation Data file!'
        CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF

    !Check ET columns
    IF (ET%GetNDataColumns() .LT. MAXVAL(AppLake%Lakes%iColET)) THEN
        iLake = MAXLOC(AppLake%Lakes%iColET)
        ID    = AppLake%Lakes(iLake(1))%ID
        MessageArray(1) = 'Evapotranspiration data column for lake '//TRIM(IntToText(ID))//' is greater than the'
        MessageArray(2) = 'available data columns in the Evapotranspiration Data file!'
        CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF

  END SUBROUTINE AppLake_v40_CheckExternalTSDataPointers

  ! -------------------------------------------------------------
  ! --- CONVERT TIME UNIT OF LAKE RELATED ENTITIES
  ! -------------------------------------------------------------
  SUBROUTINE AppLake_v40_ConvertTimeUnit(AppLake,NewUnit)
    CLASS(AppLake_v40_Type)     :: AppLake
    CHARACTER(LEN=*),INTENT(IN) :: NewUnit

    !No time unit to be converted in this version

  END SUBROUTINE AppLake_v40_ConvertTimeUnit

END MODULE

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE


text:
C=======================================================================
C  PATH, Subroutine
C
C  Program to read the PATH from DSSATPRO file
C-----------------------------------------------------------------------
C  Revision history
C
C  01/01/1991 GH  Written
C  05/28/1993 PWW Header revision and minor changes
C  09/17/2007 JIL Added codes for IXIM maize model
C  08/09/2012 GH  Add codes for CSCAS cassava model
C  04/16/2013 CHP/KAD Added codes for SALUS model
!  05/09/2013 CHP/FR/JZW Added N-wheat module
C  06/03/2015 LPM Added codes for CSYCA CIAT cassava model
C  06/18/2015 GH  Add error code for configuration file issues
C  09/26/2017 WP  Add DSSAT_HOME environment variable configuration
C-----------------------------------------------------------------------
C  INPUT  : PROCOD,PFLAG
C
C  LOCAL  : BLANK ERRKEY FILEPR LINE LUNPR ERRNUM PATHL I K FEXIST
C
C  OUTPUT : PATHC,NAMEF
C-----------------------------------------------------------------------
C  Called : SECLI SEWTH IPEXP
C
C  Calls  : ERROR
C=======================================================================

      SUBROUTINE PATH (PROCOD,DSSATP,PATHC,PFLAG,NAMEF)

      USE ModuleDefs
      IMPLICIT     NONE

      CHARACTER*1  BLANK
      CHARACTER*3  PROCOD
      CHARACTER*6  ERRKEY
      CHARACTER*12 NAMEF
      CHARACTER*102 DSSATP
      CHARACTER*80 LINE,PATHC

      INTEGER      LUNPR,ERRNUM,PATHL,I,K,PFLAG
      LOGICAL      FEXIST

!     PARAMETER (LUNPR  = 15)
      PARAMETER (ERRKEY = 'PATH  ')
      PARAMETER (BLANK  = ' ')

      PATHC = BLANK
      CALL GETLUN('DSPRO',LUNPR)
      INQUIRE (FILE = DSSATP,EXIST = FEXIST)
      IF (FEXIST) THEN
         OPEN (LUNPR,FILE=DSSATP,STATUS = 'OLD',IOSTAT=ERRNUM)
      ENDIF
      IF (.NOT. FEXIST .OR. ERRNUM .NE. 0) RETURN

      DO I = 1, 500
         READ (LUNPR,'(A80)',IOSTAT=ERRNUM) LINE
         IF (LINE(1:3) .EQ. PROCOD) THEN
!            PATHC  = LINE(5:6) // LINE(8:80)
            call get_next_string(line,4,pathc)
            if(index(line,trim(pathc))<7)then
               call get_next_string(line,7,pathc)
               pathc = trim(adjustl(line(5:6)))//trim(pathc)
            end if

C-SUN       PATHC  = LINE(8:80)
            PATHL  = INDEX (PATHC,BLANK)
            IF (PATHL .EQ. 1) THEN
                 CALL ERROR (ERRKEY,3,DSSATP,I)
            ENDIF
            IF (PATHC(PATHL-1:PATHL-1) /= SLASH) THEN
              WRITE (PATHC(PATHL:PATHL),'(A1)') SLASH
            ELSE
              PATHL = PATHL - 1
            ENDIF
            NAMEF  = PATHC (PATHL+1:PATHL+13)
            IF (PATHL .LT. 80) THEN
               DO K = (PATHL+1),80
                  IF (PATHC(K:K) .NE. BLANK) THEN
                     PATHC(K:K) = BLANK
                  ENDIF
               END DO
            ENDIF
            CLOSE (LUNPR)
            RETURN
         ENDIF
      END DO

      IF (PFLAG .EQ. 1) THEN
         WRITE(*,600) PROCOD, DSSATP
 600     FORMAT(' CODE ',A3,' not found in file ',A)
         CALL ERROR (ERRKEY,1,DSSATP,I)
      ENDIF

      CLOSE (LUNPR)
      END SUBROUTINE PATH

C=======================================================================
C  PATHD, Subroutine
C
C  Program to read the PATH of DSSATPRO file
C    Assume in current directory first.
C    If not in current directory, DSSATPRO file is stored with the EXEs
C-----------------------------------------------------------------------
C  Revision history
C
C  1. Written                                       G.H         8-10-93
C  2  Modified by
C  3. Header revision and minor changes             P.W.W.      5-28-93
! 07/17/2008 CHP DSSATPRO name specified in ModuleDefs
C-----------------------------------------------------------------------
C  INPUT  : INPUTX
C
C  LOCAL  : BLANK ERRKEY FILEPR LINE LUNPR ERRNUM PATHL I K FEXIST
C
C  OUTPUT : DSSATP
C-----------------------------------------------------------------------
C  Called : INPUT
C
C-----------------------------------------------------------------------
C                         DEFINITIONS
C
C  FEXIST :
C=======================================================================

      SUBROUTINE PATHD (DSSATP,INPUTX,IP)

      USE ModuleDefs
      IMPLICIT     NONE

!     CHARACTER*12 DSSATF
      CHARACTER*102 DSSATP
      CHARACTER*120 INPUTX
      CHARACTER(len=255) :: DSSAT_HOME

      INTEGER      I
      INTEGER    IP

      LOGICAL      FEXIST

!     DSSATP(1:12) = DSSATF
      DSSATP(1:12) = DSSATPRO
      INQUIRE (FILE = DSSATP,EXIST = FEXIST)
      IF (.NOT. FEXIST .AND. IP .GT. 12) THEN
         DO I = IP, 0, -1
           IF (INPUTX(I:I) .EQ. SLASH .OR. INPUTX(I:I) .EQ. "/")GO TO 10
         END DO
   10    CONTINUE
!        DSSATP(1:I+12) = INPUTX(1:I) // DSSATF
         DSSATP(1:I+12) = INPUTX(1:I) // DSSATPRO
      ENDIF

      INQUIRE (FILE = DSSATP,EXIST = FEXIST)
      IF (.NOT. FEXIST) THEN
        CALL get_environment_variable("DSSAT_HOME", DSSAT_HOME)
        IF(TRIM(DSSAT_HOME) .NE. '') THEN
            STDPATH = TRIM(DSSAT_HOME)
        ENDIF
        DSSATP = trim(STDPATH) // TRIM(DSSATPRO)
      ENDIF

      INQUIRE (FILE = DSSATP,EXIST = FEXIST)
      IF (.NOT. FEXIST) THEN
        CALL ERROR ('PATH  ',2,DSSATP,0)
      ENDIF

      RETURN
      END SUBROUTINE PATHD

C=======================================================================
C  MODEL_NAME, Subroutine
C
C  Program to read the Module name from DSSATPRO file
C-----------------------------------------------------------------------
C  Revision history
C
C  02/21/2006 GH  Written
!  10/25/2006 CHP Added CRMODEL, the model name from FILEX which will
!                   override MODEL from DSSATPRO, if valid.
!  09/17/2007 JIL Added codes for IXIM maize model
!  04/17/2013 CHP Added exception for crop-model matching for
!                   SALUS generic crop model
C-----------------------------------------------------------------------
C  INPUT  : PROCOD,PFLAG
C
C  LOCAL  : BLANK ERRKEY FILEPR LINE LUNPR ERRNUM PATHL I K FEXIST
C
C  OUTPUT : PATHC,NAMEF
C-----------------------------------------------------------------------
C  Called : SECLI SEWTH IPEXP
C
C  Calls  : ERROR
C=======================================================================

      SUBROUTINE MODEL_NAME (CROP, DSSATP, CRMODEL, MODEL)

      USE ModuleDefs
      IMPLICIT     NONE

      CHARACTER*1   UPCASE
      CHARACTER*2   CROP
      CHARACTER*3   PROCOD
      CHARACTER*6   ERRKEY
      CHARACTER*8   MODEL, CRMODEL
      CHARACTER*78  MSG(4)
      CHARACTER*80  LINE
      CHARACTER*102 DSSATP
!     CHARACTER*120 PATHX

      INTEGER      EXE_POS,LUNPR,LINPR,ERRNUM,ISECT,I,J   !, IPX
      INTEGER      ValidModel, LENSTRING

      LOGICAL FEXIST

      PARAMETER (LUNPR = 25)
      PARAMETER (ERRKEY = 'MODELN')
!      PARAMETER (EXE_STRING = 'EXE')

      IF (LENSTRING(CRMODEL) > 0) THEN
        DO I = 1,8
          CRMODEL(I:I)= UPCASE(CRMODEL(I:I))
        ENDDO
!       CHP 5/22/2008 Check for model version being integer value
!       Pick off name of executable, if necessary.
        DO I = 6,8
          J = ICHAR(CRMODEL(I:I))
          IF (J < 48 .OR. J > 57) THEN
!            Fill in model version from executable name
             CRMODEL(6:8) = ModelVerTxt
            EXIT
          ENDIF
        ENDDO
        ERRNUM = ValidModel(CROP, CRMODEL)
        IF (ERRNUM == 0) THEN
          MODEL = CRMODEL
          RETURN
        ELSE
          MSG(1) = "Invalid model specified."
          WRITE(MSG(2),'("Crop: ",A2,"; Model: ",A8)') CROP, CRMODEL
          MSG(3) = "Use default model from DSSATPRO."
          CALL WARNING(3,ERRKEY,MSG)
        ENDIF
      ELSE
        ERRNUM = 8
      ENDIF
!     If no model in FILEX, or invalid, use DSSATPRO default model

      INQUIRE (FILE = DSSATP,EXIST = FEXIST)
      IF (FEXIST) THEN
         OPEN (LUNPR,FILE=TRIM(DSSATP),STATUS = 'OLD',IOSTAT=ERRNUM)
         IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,1,DSSATP,0)
      ELSE
         CALL ERROR (ERRKEY,2,DSSATP,0)
      ENDIF

      PROCOD = 'M' // CROP

      REWIND(LUNPR)
      LINPR = 0

 100  CONTINUE
      CALL IGNORE (LUNPR, LINPR, ISECT,LINE)
      IF (ISECT .EQ. 0) THEN
         MSG(1) = "Code not found in DSSATPRO file."
         MSG(2) = "File: " // DSSATP(1:72)
         WRITE(MSG(3),'("Code = ",A3)') PROCOD
         MSG(4) = "Program will end."
         CALL WARNING(4,ERRKEY,MSG)
         CALL ERROR (ERRKEY,3,DSSATP,LINPR)
      ELSEIF (LINE(1:3) .NE. PROCOD) THEN
         GO TO 100
      ELSE
         EXE_POS = INDEX(LINE,'DSCSM')
         if(exe_pos==0) exe_pos = index(line,'dscsm')
!         MODEL = LINE((EXE_POS+4):(EXE_POS+11))
         call get_next_string(line,8,model)
         exe_pos = index(line,model)
         call get_next_string(line,exe_pos,model)
         DO I = 1,5
           MODEL(I:I)= UPCASE(MODEL(I:I))
         ENDDO
      ENDIF

      CLOSE (LUNPR)

!     Check MODEL for validity.
      ERRNUM = ValidModel(CROP, MODEL)
      IF (ERRNUM /= 0) THEN
        MSG(1) = "No valid model found in DSSATPRO for this crop."
        WRITE(MSG(2),'("Crop: ",A2)') CROP
        MSG(3) = "Simulation will stop."
        CALL WARNING(3,ERRKEY,MSG)
        CALL ERROR (ERRKEY,ERRNUM,DSSATP,0)
      ENDIF

      RETURN
      END SUBROUTINE MODEL_NAME

!======================================================================
      FUNCTION ValidModel(CROP, MODEL)

      IMPLICIT NONE
      SAVE

      INTEGER ValidModel
      CHARACTER*2, INTENT(IN) :: CROP
      CHARACTER*8, INTENT(IN) :: MODEL

      INTEGER, PARAMETER :: MaxNum = 100
      INTEGER I
      CHARACTER*2 CropName(MaxNum)
      CHARACTER*5 ModelName(MaxNum)
!      CHARACTER*78 MSG(2)

      LOGICAL FIRST
      DATA FIRST /.TRUE./

      ValidModel = 0      !0 indicates no error

!     Blank model = error 6
      IF (MODEL(1:3) .EQ. '   ') THEN
        ValidModel = 6
        RETURN
      ENDIF

!!    Check for valid crop modeling approaches
!      IF ((INDEX(MODEL(3:5),'CER') .EQ. 0) .AND.   !CERES
!!    &    (INDEX(MODEL(3:5),'ORZ') .EQ. 0) .AND.   !ORYZA
!     &    (INDEX(MODEL(3:5),'IXM') .EQ. 0) .AND.   !IXIM
!     &    (INDEX(MODEL(3:5),'GRO') .EQ. 0) .AND.   !CROPGRO
!     &    (INDEX(MODEL(3:5),'FRM') .EQ. 0) .AND.   !FORAGE
!     &    (INDEX(MODEL(3:5),'CSM') .EQ. 0) .AND.   !CROPSIM (Cereal)
!     &    (INDEX(MODEL(3:5),'CAS') .EQ. 0) .AND.   !CSCAS (Cassava)
!     &    (INDEX(MODEL(3:5),'YCA') .EQ. 0) .AND.   !CSYCA (CIAT-Cassava)
!     &    (INDEX(MODEL(3:5),'SIM') .EQ. 0) .AND.   !CROPSIM (Cassava)
!     &    (INDEX(MODEL(3:5),'SUB') .EQ. 0) .AND.   !SUBSTOR
!     &    (INDEX(MODEL(3:5),'CAN') .EQ. 0) .AND.   !CANEGRO
!     &    (INDEX(MODEL(3:5),'CSP') .EQ. 0) .AND.   !CASUPRO
!     &    (INDEX(MODEL(3:5),'SAM') .EQ. 0) .AND.   !SAMUCA
!     &    (INDEX(MODEL(3:5),'ALO') .EQ. 0) .AND.   !ALOHA
!     &    (INDEX(MODEL(3:5),'ARO') .EQ. 0) .AND.   !AROIDS
!     &    (INDEX(MODEL(3:5),'CRP') .EQ. 0) .AND.   !CropSim cassava
!     &    (INDEX(MODEL(3:5),'APS') .EQ. 0) .AND.   !APSIM N-wheat
!     &    (INDEX(MODEL(3:5),'OIL') .EQ. 0) .AND.   !OILCROP
!     &    (INDEX(MODEL(3:5),'LUS') .EQ. 0)         !SALUS
!     &    ) THEN
!!       Invalid model = error 4
!        ValidModel = 4
!        RETURN
!      ENDIF
!
!!     Check for valid crop modules
!      IF ((INDEX(MODEL(1:5),'CSCER') .EQ. 0) .AND. !Wheat and Barley
!     &    (INDEX(MODEL(1:5),'CSCRP') .EQ. 0) .AND. !Wheat and barley
!     &    (INDEX(MODEL(1:5),'CSCAS') .EQ. 0) .AND. !Cassava
!     &    (INDEX(MODEL(1:5),'CSYCA') .EQ. 0) .AND. !Cassava CIAT
!     &    (INDEX(MODEL(1:5),'WHAPS') .EQ. 0) .AND. !APSIM N-wheat
!     &    (INDEX(MODEL(1:5),'TFAPS') .EQ. 0) .AND. !Tef based on N-wheat
!     &    (INDEX(MODEL(1:5),'CRGRO') .EQ. 0) .AND. !CROPGRO (All
!!                         grain legumes, grasses, vegetables and cotton
!     &    (INDEX(MODEL(1:5),'PRFRM') .EQ. 0) .AND. !FORAGE
!     &    (INDEX(MODEL(1:5),'MZCER') .EQ. 0) .AND. !Maize CERES
!     &    (INDEX(MODEL(1:5),'MZIXM') .EQ. 0) .AND. !Maize IXIM
!     &    (INDEX(MODEL(1:5),'MLCER') .EQ. 0) .AND. !Millet
!     &    (INDEX(MODEL(1:5),'PIALO') .EQ. 0) .AND. !Aloha Pineapple
!     &    (INDEX(MODEL(1:5),'PTSUB') .EQ. 0) .AND. !Potato
!     &    (INDEX(MODEL(1:5),'RICER') .EQ. 0) .AND. !CERES-Rice
!!    &    (INDEX(MODEL(1:5),'RIORZ') .EQ. 0) .AND. !ORYZA-Rice
!     &    (INDEX(MODEL(1:5),'SGCER') .EQ. 0) .AND. !Sorghum
!     &    (INDEX(MODEL(1:5),'SCCSP') .EQ. 0) .AND. !Sugarcane CASUPRO
!     &    (INDEX(MODEL(1:5),'SCCAN') .EQ. 0) .AND. !Sugarcane CaneGro
!     &    (INDEX(MODEL(1:5),'SCSAM') .EQ. 0) .AND. !Sugarcane SAMUCA
!     &    (INDEX(MODEL(1:5),'BSCER') .EQ. 0) .AND. !Sugarbeet VSH
!     &    (INDEX(MODEL(1:5),'SWCER') .EQ. 0) .AND. !Sweet corn
!     &    (INDEX(MODEL(1:5),'TNARO') .EQ. 0) .AND. !Tanier
!     &    (INDEX(MODEL(1:5),'TRARO') .EQ. 0)       !Taro
!!    &    (INDEX(MODEL(1:5),'TRARO') .EQ. 0) .AND. !Taro
!!    &    (INDEX(MODEL(1:5),'SUOIL') .EQ. 0)       !Oilcrop
!     &    ) THEN
!!       Invalid module name = error 5
!        ValidModel = 5
!        RETURN
!      ENDIF

!     Check that crop and model are compatible based on codes in
!     Simulation.CDE
      IF (FIRST) THEN
        CALL ReadCropModels(MaxNum, ModelName, CropName)
        FIRST = .FALSE.
      ENDIF

!!     Generic SALUS model can be used for any crop,
!!       identified in cultivar data
!      IF (MODEL(1:5) == 'SALUS') THEN
!        ValidModel = 0
!        RETURN
!      ENDIF
!
      ValidModel = 7  !Assume crop does NOT match model
      DO I = 1, MaxNum
        IF (CropName(I) == CROP .AND. ModelName(I) == MODEL(1:5)) THEN
!         Model matches crop -- return with no error
          ValidModel = 0
          RETURN
        ENDIF
      ENDDO

!      IF (ValidModel /= 0) THEN
!        WRITE(MSG(1),100) CROP, MODEL
!  100   FORMAT("Crop ",A2," can not be simulated with model ",A5)
!        MSG(2)="Check SIMULATION.CDE file for valid crop models."
!        CALL WARNING(2,"MODCHECK",MSG)
!      ENDIF

      RETURN
      END FUNCTION ValidModel
!======================================================================
!  MODEL and CROP should be modified when model versions change
!     or when a crop specific model is created.

!  MODEL   CROP  Description
!  BSCER    BS   CERES-Sugarbeet
!  CRGRO    BM   CROPGRO-Bermuda
!  CRGRO    BH   CROPGRO-Bahia
!  CRGRO    BN   CROPGRO-Dry bean
!  CRGRO    BR   CROPGRO-Brachiaria
!  CRGRO    CB   CROPGRO-Cabbage
!  CRGRO    CH   CROPGRO-Chickpea
!  CRGRO    CI   CROPGRO-Chia
!  CRGRO    CN   CROPGRO-Canola
!  CRGRO    CO   CROPGRO-Cotton
!  CRGRO    CP   CROPGRO-Cowpea
!  CRGRO    FA   CROPGRO-Fallow
!  CRGRO    FB   CROPGRO-Faba bean
!  CRGRO    GB   CROPGRO-Grean bean
!  CRGRO    LT   CROPGRO-Lentil
!  CRGRO    NP   CROPGRO-Napier grass
!  CRGRO    PN   CROPGRO-Peanut
!  CRGRO    PP   CROPGRO-Pigeonpea
!  CRGRO    PR   CROPGRO-Bellpepper
!  CRGRO    QU   CROPGRO-Quinoa
!  CRGRO    SF   CROPGRO-Safflower
!  CRGRO    SB   CROPGRO-Soybean
!  CRGRO    SU   CROPGRO-Sunflower
!  CRGRO    TM   CROPGRO-Tomato
!  CRGRO    VB   CROPGRO-Velvetbean
!  CSCER    BA   CROPSIM-CERES-Barley
!  CSCER    WH   CROPSIM-CERES-Wheat
!  CSCRP    BA   CSCRP-Barley
!  CSCAS    CS   CSCAS-Cassava
!  CSYCA    CS   MANIHOT-Cassava
!  CSCRP    WH   CSCRP-Wheat
!  MLCER    ML   CERES-Millet
!  MZCER    MZ   CERES-Maize
!  MZIXM    MZ   IXIM-Maize
!  PTSUB    PT   SUBSTOR-Potato
!  RICER    RI   CERES-Rice
!  SCCAN    SC   CANEGRO-Sugarcane
!  SCCSP    SC   CASUPRO-Sugarcane
!  SCSAM    SC   SAMUCA-Sugarcane
!  SGCER    SG   CERES-Sorghum
!  SWCER    SW   CERES-Sweetcorn
!  PIALO    PI   ALOHA-Pineapple
!  TRARO    TR   AROIDS-Taro
!  TNARO    TN   AROIDS-Tanier
!  TFAPS    TF   NWHEAT-Teff
!  WHAPS    WH   NWHEAT-Wheat
!  PRFRM    BM   FORAGE-Bermudagrass
!  PRFRM    BR   FORAGE-Brachiaria
!  PRFRM    BH   FORAGE-Bahiagrass
!  PRFRM    AL   FORAGE-Alfalfa
!  SUOIL    SU   OILCROP-Sunflower
!  CRGRO    SR   CRGRO-Strawberry
!  CRGRO    SQ   CRGRO-WinterSquash
!  CRGRO    CI   CRGRO-Chia
!-----------------------------------------------------------------------

ext: for

max_stars_repo_name: xiangyi-wang/dssat-csm-os

max_issues_repo_name: xiangyi-wang/dssat-csm-os

max_forks_repo_name: xiangyi-wang/dssat-csm-os


text:
use DictionaryClass
use RandomClass
implicit none

type(Dictionary_) :: df,df_copy
type(Random_)     :: random

df = dict()
!
do i_i=1,100
    call df%update("2022/01/01/"+zfill(i_i,4), random%gauss(mu=0.0d0,sigma=1.0d0) )
enddo

! to show
call df%show()

call df%to_csv("filename")
call df%to_json("filename")

df_copy = to_dict("filename.csv")

do i_i=20,30
    call df_copy%update("2022/01/01/"+zfill(i_i,4),0.0d0)
enddo

call df_copy%to_csv("filename2")

end

max_stars_repo_name: kazulagi/faba

max_issues_repo_name: kazulagi/faba

max_forks_repo_name: kazulagi/faba


text:
program mean_mask
  use stdlib_stats, only: mean
  implicit none
  real :: y(1:2, 1:3) = reshape([1., 2., 3., 4., 5., 6.], [2, 3])

  print*, mean(y, y > 2.)
  print*, mean(y, 1, y > 2.)
  print*, mean(y, 2, y > 2.)
end program mean_mask

max_stars_repo_name: awvwgk/stdlib-docs

max_issues_repo_name: awvwgk/stdlib-docs

max_forks_repo_name: awvwgk/stdlib-docs


text:
      subroutine solsbg (n, t, cc, f0, wk, hl0, bd, ipbd, v, lr, ier)
c-----------------------------------------------------------------------
c This routine applies one or two inverse preconditioner matrices
c to the array v, using the interaction-only block-diagonal Jacobian
c with block-grouping, and Gauss-Seidel applied to the diffusion terms.
c When lr = 1 or 3, it calls gs for a Gauss-Seidel approximation
c to ((I-hl0*Jd)-inverse)*v, and stores the result in v.
c When lr = 2 or 3, it computes ((I-hl0*dg/dc)-inverse)*v, using LU
c factors of the blocks in bd, and pivot information in ipbd.
c In both cases, the array v is overwritten with the solution.
c-----------------------------------------------------------------------
      integer n, ipbd, lr, ier
      double precision t, cc, f0, wk, hl0, bd, v
      dimension cc(*), f0(*), wk(*), bd(*), ipbd(*), v(*)
      dimension n(1)
      integer mp, mq, mpsq, itmax,
     2        meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr
      integer ibd, ig0, igm1, igx, igy, iip, iv, jx, jy
      double precision uround, srur
c
      common /pcom2/ uround, srur, mp, mq, mpsq, itmax
      common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp,
     1   jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20)
c
      ier = 0
c
      if (lr.eq.0 .or. lr.eq.1 .or. lr.eq.3) call gs (n, hl0, v, wk)
      if (lr.eq.0 .or. lr.eq.2 .or. lr.eq.3) then
        iv = 1
        do 20 jy = 1,meshy
          igy = jigy(jy)
          ig0 = (igy - 1)*ngx
          do 10 jx = 1,meshx
            igx = jigx(jx)
            igm1 = igx - 1 + ig0
            ibd = 1 + igm1*mpsq
            iip = 1 + igm1*mp
            call dgesl (bd(ibd), mp, mp, ipbd(iip), v(iv), 0)
            iv = iv + mp
 10         continue
 20       continue
        endif
c
      return
c------------  end of subroutine solsbg  -------------------------------
      end

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Right smack dab in the middle of the intersection of 3rd Street and University Avenue lie these crazy poles called bollards that prevent you from driving in a sane fashion (especially when trying to get to XCafe RomaX 3rd & U Cafe by campus or XNavins Copy ShopX Davis Copy Shop).  The poles are filled and chained to the street, but can be removed when needed.  During Picnic Day the City of Davis city takes them out of their holes so police/parade traffic can get through. The poles were removed March 1214 for the Downtown Davis Criterium bicycle race organized by the Cal Aggie Cycling team.

The poles form a diagonal from the northeast corner to the southwest corner.

In 2004, the owners of several businesses sitting on the start of 3rd Street investigated getting the poles removed.  The poles presence, they maintained, hampered business on their area of the street.

The poles are surely there to keep traffic out of the area. 3rd Street would likely have much more traffic and fast driving on it if this Driving in Davis traffic slowing measure wasnt there. As is, the street is conductive to walkers and bikers.

Despite the reduced car traffic, it is very unpleasant for cyclists. Navigating through the bollards with a dozen other bikes in close proximity can be hair raising. The odd dip in the road and poor pavement generates great discomfort in the groin. It also accumulates standing water in the wet seasons, leading to serious cases of the Freshman Stripe.

Photos

Sometimes the poles get abused or vandalized.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      function ignbin( n, pp )

c*********************************************************************72
c
cc IGNBIN generates a binomial random deviate.
c
c  Discussion:
c
c    This procedure generates a single random deviate from a binomial
c    distribution whose number of trials is N and whose
c    probability of an event in each trial is P.
c
c    The previous version of this program relied on the assumption that
c    local memory would be preserved between calls.  It set up data
c    one time to be preserved for use over multiple calls.  In the
c    interests of portability, this assumption has been removed, and
c    the "setup" data is recomputed on every call.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    29 March 2013
c
c  Author:
c
c    Original FORTRAN77 version by Barry Brown, James Lovato.
c    Modifications by John Burkardt.
c
c  Reference:
c
c    Voratas Kachitvichyanukul, Bruce Schmeiser,
c    Binomial Random Variate Generation,
c    Communications of the ACM,
c    Volume 31, Number 2, February 1988, pages 216-222.
c
c  Parameters:
c
c    Input, integer N, the number of binomial trials, from which a
c    random deviate will be generated.
c    0 < N.
cf
c    Input, real*8 PP, the probability of an event in each trial of the
c    binomial distribution from which a random deviate is to be generated.
c    0.0 < PP < 1.0.
c
c    Output, integer IGNBIN, a random deviate from the distribution.
c
      implicit none

      real*8 al
      real*8 alv
      real*8 amaxp
      real*8 c
      real*8 f
      real*8 f1
      real*8 f2
      real*8 ffm
      real*8 fm
      real*8 g
      integer i
      integer ignbin
      integer ix
      integer ix1
      integer k
      integer m
      integer mp
      real*8 pp
      integer n
      real*8 p
      real*8 p1
      real*8 p2
      real*8 p3
      real*8 p4
      real*8 q
      real*8 qn
      real*8 r
      real*8 t
      real*8 u
      real*8 v
      real*8 w
      real*8 w2
      real*8 x
      real*8 x1
      real*8 x2
      real*8 xl
      real*8 xll
      real*8 xlr
      real*8 xm
      real*8 xnp
      real*8 xnpq
      real*8 xr
      real*8 ynorm
      real*8 z
      real*8 z2

      if (pp == 0.0D+00 .or. n == 0) then
         ignbin=0
         return
      endif
c$$$      if ( pp < 0.0D+00 .or. 1.0D+00 <= pp ) then
c$$$        write ( *, '(a)' ) ' '
c$$$        write ( *, '(a)' ) 'IGNBIN - Fatal error!'
c$$$        write ( *, '(a)' ) '  PP is out of range.'
c$$$        print*, pp
c$$$        print*,n
c$$$        return
c$$$        !stop 1
c$$$      end if

      p = min ( pp, 1.0D+00 - pp )
      q = 1.0D+00 - p
      xnp = dble ( n ) * p

      if ( xnp .lt. 30.0D+00 ) then
        qn = q ** n
        r = p / q
        g = r * real ( n + 1 )
        go to 20
      end if
c
c  The calculation of this data was originally intended to be
c  done once, then saved for later calls.
c
      ffm = xnp + p
      m = ffm
      fm = m
      xnpq = xnp * q
      p1 = nint ( 2.195D+00 * sqrt ( xnpq ) - 4.6D+00 * q ) + 0.5D+00
      xm = fm + 0.5D+00
      xl = xm - p1
      xr = xm + p1
      c = 0.134D+00 + 20.5D+00 / ( 15.3D+00 + fm )
      al = ( ffm - xl ) / ( ffm - xl * p )
      xll = al * ( 1.0D+00 + 0.5D+00 * al )
      al = ( xr - ffm ) / ( xr * q )
      xlr = al * ( 1.0D+00 + 0.5D+00 * al )
      p2 = p1 * ( 1.0D+00 + c + c )
      p3 = p2 + c / xll
      p4 = p3 + c / xlr
c
c  Generate a variate.
c
10    continue

c      u = ran1(iseed) * p4
c      v = ran1(iseed)
      u = rand() * p4
      v = rand()
c
c  Triangle
c
      if ( u .lt. p1 ) then
        ix = xm - p1 * v + u
        if ( 0.5D+00 .lt. pp ) then
          ix = n - ix
        end if
        ignbin = ix
        return
      end if
c
c  Parallelogram
c
      if ( u .le. p2 ) then

        x = xl + ( u - p1 ) / c
        v = v * c + 1.0D+00 - dabs ( xm - x ) / p1

        if ( v .le. 0.0D+00 .or. 1.0D+00 .lt. v ) then
          go to 10
        end if

        ix = x

      else if ( u .le. p3 ) then

        ix = xl + log ( v ) / xll
        if ( ix .lt. 0 ) then
          go to 10
        end if
        v = v * ( u - p2 ) * xll

      else

        ix = xr - log ( v ) / xlr
        if ( n .lt. ix ) then
          go to 10
        end if
        v = v * ( u - p3 ) * xlr

      end if

      k = abs ( ix - m )

      if ( k .le. 20 .or. xnpq / 2.0 - 1.0 .le. k ) then

        f = 1.0D+00
        r = p / q
        g = ( n + 1 ) * r

        if ( m .lt. ix ) then
          mp = m + 1
          do i = mp, ix
            f = f * ( g / i - r )
          end do
        else if ( ix .lt. m ) then
          ix1 = ix + 1
          do i = ix1, m
            f = f / ( g / i - r )
          end do
        end if

        if ( v .le. f ) then
          if ( 0.5D+00 .lt. pp ) then
            ix = n - ix
          end if
          ignbin = ix
          return
        end if

      else

        amaxp = ( k / xnpq ) * ( ( k * ( k / 3.0D+00
     &    + 0.625D+00 ) + 0.1666666666666D+00 ) / xnpq + 0.5D+00 )
        ynorm = - real ( k * k ) / ( 2.0D+00 * xnpq )
        alv = log ( v )

        if ( alv .lt. ynorm - amaxp ) then
          if ( 0.5D+00 .lt. pp ) then
            ix = n - ix
          end if
          ignbin = ix
        return
        end if

        if ( ynorm + amaxp .lt. alv ) then
          go to 10
        end if

        x1 = real ( ix + 1 )
        f1 = fm + 1.0D+00
        z = real ( n + 1 ) - fm
        w = real ( n - ix + 1 )
        z2 = z * z
        x2 = x1 * x1
        f2 = f1 * f1
        w2 = w * w

        t = xm * log ( f1 / x1 ) + ( n - m + 0.5D+00 ) * log ( z / w )
     &    + real ( ix - m ) * log ( w * p / ( x1 * q ))
     &    + ( 13860.0D+00 - ( 462.0D+00 - ( 132.0D+00
     &    - ( 99.0D+00 - 140.0D+00
     &    / f2 ) / f2 ) / f2 ) / f2 ) / f1 / 166320.0D+00
     &    + ( 13860.0D+00 - ( 462.0D+00 - ( 132.0D+00
     &    - ( 99.0D+00 - 140.0D+00
     &    / z2 ) / z2 ) / z2 ) / z2 ) / z / 166320.0D+00
     &    + ( 13860.0D+00 - ( 462.0D+00 - ( 132.0D+00
     &    - ( 99.0D+00 - 140.0D+00
     &    / x2 ) / x2 ) / x2 ) / x2 ) / x1 / 166320.0D+00
     &    + ( 13860.0D+00 - ( 462.0D+00 - ( 132.0D+00
     &    - ( 99.0D+00 - 140.0D+00
     &    / w2 ) / w2 ) / w2 ) / w2 ) / w / 166320.0D+00

        if ( alv .le. t ) then
          if ( 0.5D+00 .lt. pp ) then
            ix = n - ix
          end if
          ignbin = ix
          return
        end if

      end if

      go to 10
c
c  Mean less than 30.
c
20    continue

      ix = 0
      f = qn

c     u = ran1(iseed)
      u=rand()

30    continue

      if ( u .lt. f ) then
        if ( 0.5D+00 .lt. pp ) then
          ix = n - ix
        end if
        ignbin = ix
        return
      end if

      if ( ix .le. 110 ) then
        u = u - f
        ix = ix + 1
        f = f * ( g / dble ( ix ) - r )
        go to 30
      end if

      go to 20

      return
      end

max_stars_repo_name: predsci/PSI-CARM

max_issues_repo_name: predsci/PSI-CARM

max_forks_repo_name: predsci/PSI-CARM


text:
*     FB01SD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX, PMAX
      PARAMETER        ( NMAX = 20, MMAX = 20, PMAX = 20 )
      INTEGER          LDAINV, LDB, LDC, LDQINV, LDRINV, LDSINV
      PARAMETER        ( LDAINV = NMAX, LDB = NMAX, LDC = PMAX,
     $                   LDQINV = MMAX, LDRINV = PMAX, LDSINV = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( NMAX*(NMAX + 2*MMAX) + 3*MMAX,
     $                                (NMAX + PMAX)*(NMAX + 1) + 2*NMAX,
     $                                 3*NMAX ) )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, INFO, ISTEP, J, M, N, P
      CHARACTER*1      JOBX, MULTAB, MULTRC
*     .. Local Arrays ..
      DOUBLE PRECISION AINV(LDAINV,NMAX), B(LDB,MMAX), C(LDC,NMAX),
     $                 DIAG(MMAX), DWORK(LDWORK), E(PMAX),
     $                 QINV(LDQINV,MMAX), RINV(LDRINV,PMAX),
     $                 RINVY(PMAX), SINV(LDSINV,NMAX), X(NMAX), Z(MMAX)
      INTEGER          IWORK(NMAX)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         DCOPY, FB01SD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, M, P, JOBX, TOL, MULTAB, MULTRC
      IF ( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( AINV(I,J), J = 1,N ), I = 1,N )
         IF ( P.LE.0 .OR. P.GT.PMAX ) THEN
            WRITE ( NOUT, FMT = 99991 ) P
         ELSE
            READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
            IF ( LSAME( MULTRC, 'N' ) ) READ ( NIN, FMT = * )
     $                               ( ( RINV(I,J), J = 1,P ), I = 1,P )
            IF ( M.LE.0 .OR. M.GT.MMAX ) THEN
               WRITE ( NOUT, FMT = 99992 ) M
            ELSE
               READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
               READ ( NIN, FMT = * ) ( ( QINV(I,J), J = 1,M ), I = 1,M )
               READ ( NIN, FMT = * ) ( ( SINV(I,J), J = 1,N ), I = 1,N )
               READ ( NIN, FMT = * ) ( Z(J), J = 1,M )
               READ ( NIN, FMT = * ) ( X(J), J = 1,N )
               READ ( NIN, FMT = * ) ( RINVY(J), J = 1,P )
*              Save the strict upper triangle of QINV in its strict
*              lower triangle and the diagonal in the array DIAG.
               DO 10 I = 2, M
                  CALL DCOPY( I, QINV(1,I), 1, QINV(I,1), LDQINV )
   10          CONTINUE
               CALL DCOPY( M, QINV, LDQINV+1, DIAG, 1 )
*              Perform three iterations of the (Kalman) filter recursion
*              (in square root information form).
               ISTEP = 1
   20          CONTINUE
                  CALL FB01SD( JOBX, MULTAB, MULTRC, N, M, P, SINV,
     $                         LDSINV, AINV, LDAINV, B, LDB, RINV,
     $                         LDRINV, C, LDC, QINV, LDQINV, X, RINVY,
     $                         Z, E, TOL, IWORK, DWORK, LDWORK, INFO )
                  ISTEP = ISTEP + 1
                  IF ( INFO.EQ.0 .AND. ISTEP.LE.3 ) THEN
*                    Restore the upper triangle of QINV.
                     DO 30 I = 2, M
                        CALL DCOPY( I, QINV(I,1), LDQINV, QINV(1,I), 1 )
   30                CONTINUE
                     CALL DCOPY( M, DIAG, 1, QINV, LDQINV+1 )
                     GO TO 20
                  END IF
*
               IF ( INFO.NE.0 ) THEN
                  WRITE ( NOUT, FMT = 99998 ) INFO
               ELSE
                  WRITE ( NOUT, FMT = 99997 )
                  DO 40 I = 1, N
                     WRITE ( NOUT, FMT = 99996 ) ( SINV(I,J), J = 1,N )
   40             CONTINUE
                  IF ( LSAME( JOBX, 'X' ) ) THEN
                     WRITE ( NOUT, FMT = 99995 )
                     DO 50 I = 1, N
                        WRITE ( NOUT, FMT = 99994 ) I, X(I)
   50                CONTINUE
                  END IF
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' FB01SD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from FB01SD = ',I2)
99997 FORMAT (' The inverse of the square root of the state covariance',
     $       ' matrix is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The components of the estimated filtered state are ',
     $       //'   k       X(k)',/)
99994 FORMAT (I4,3X,F8.4)
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' M is out of range.',/' M = ',I5)
99991 FORMAT (/' P is out of range.',/' P = ',I5)
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text:
The Davisville Advertiser was the first newspapers newspaper in the town history history of Davis.  It lived up to its name by being mostly advertisements. Davis Enterprise Modern papers arent so obvious about their intent. The weekly publication represented the collaborative efforts of its two editors/publishers, Bonham and Orr.  It was initially launched with the financial support of LJ Plant. Defunct today, the paper had a very short run starting from December 4, 1869 and ending on May 7, 1870.

During its first few months of existence, the paper portrayed the slowly burgeoning Davisville in an optimistic light, one in which the town seemed to radiate with a nearly limitless potential for growth.  However, as the year wore on, the publishers quickly became jaded with dwindling community support and the presence of domineering land barons.  Unable to pay their bills, the two publishers closed up shop in frustration, lashing out against the townspeople in a highly bitter final issue.

The original incarnation of this newspaper is available on microfilm for viewing in Shields Library. Users/WilliamLewis has scans of all the microfilm of this publication and will probably upload them all at some point.

An attempt was made in 1892 to revive the Davisville Advertiser name under a new publisher.  This venture failed miserably, and no copies are known to have survived.


module test_module
  use constants
  use szelliptic
  use szlocal
  use szfuncs
  implicit none

contains

    subroutine test_get_k
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Compute the age of the universe two ways: !
        ! (i) carlson symmetric forms               !
        ! (ii) splines                              !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use cosmo_params, only : ct2
        implicit none
        real(dp) :: r,k,f,fint,p,q,ul,kint,err1
        complex(dp), dimension(3) :: rts
        type(szshell) :: shell
        integer i
        open(unit=3,file='k.txt',form='formatted',status='replace')
        write(3,fmt=*) 'reference (correct) t_age [Mpc]: ', ct2
        write(3,fmt='(9X,A,6X,A,7X,A,3X,A)') &
                'r', 'k(carlson)', 'k(spline)', 'rel. error'
        do i=0,40
            r=1._dp + (200._dp-1._dp)/40 * i
            ul=r
            !new method (carlson)
            call init_shell(r,shell)
            f=sqrt(C0)*szellip_1111(shell%szc,ul)
            !NB. last saved roots does not necessarily
            !correspond to the output of zbrent
            k=get_k(r)

            !spline method
            if (allocated(kspl%ddf)) then
                call splint(kspl%x,kspl%f,kspl%ddf,kspl%n,r,kint)
                p=-kint/OH2
                q=2._dp*Mfunc(r)/OH2
                call fcubic_roots(p,q,rts)
                fint=sqrt(C0)*szellip_1111(rts,ul)
                err1=abs(1._dp-fint/f)
            end if
            write(3,fmt='(f10.4,2f16.9,E13.3)') r, f, fint, err1
        end do
        close(3)
    end subroutine test_get_k

    subroutine test_get_R
        implicit none
        type(szlocal_class) :: szloc
        real(dp) :: t,r
        integer i
        t=-40._dp
        write(*,*) 'Fix t =', t, ', vary r'
        do i=0,40
            r=1._dp + (200._dp-1._dp)/40 * i
            szloc=init_szlocal_class(t,r)
            call get_R(szloc)
            write(*,'(f16.9,f16.9)') r,szloc%f%R
        end do
        r=20._dp
        write(*,*) '*********************************'
        write(*,*) 'Fix r =', r, ', vary t'
        do i=0,40
            t=0._dp + (-400._dp)/40 * i
            szloc=init_szlocal_class(t,r)
            call get_R(szloc)
            write(*,'(f16.9,f16.9)') t,szloc%f%R
        end do
    end subroutine test_get_R

    subroutine test_get_kprime
        use utils, only : dfridr
        implicit none
        type(szlocal_class) :: szloc
        real(dp) :: h,err,r,kp,kp1,err1
        integer i
        real(dp), parameter :: t=0._dp
        h = 0.1_dp
        open(unit=13,file='kp.txt',form='formatted',status='replace')
        write(13,'(9X,A,5X,A,7X,A,3X,A)') &
                'r','k''(numerical)','k''(carlson)','rel. error'
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            kp=dfridr(k_of_r,r,h,err)
            szloc=init_szlocal_class(t,r)
            call get_kprime(szloc)
            err1=abs(1._dp-kp/szloc%f%kp)
            write(13,'(F10.4,E18.8,E18.8,E13.3)') r,kp,szloc%f%kp,err1
        end do
        close(13)
    contains
        real(dp) function k_of_r(r)
            implicit none
            real(dp), intent(in) :: r
            k_of_r=get_k(r)
        end function k_of_r
    end subroutine test_get_kprime

    subroutine test_get_Rprime
        use utils, only : dfridr
        implicit none
        type(szlocal_class) :: szloc
        real(dp) :: t,h,r,err,err1
        real(dp) :: Rprime,Rprime1
        integer i
        h=0.1_dp
        t=-100._dp
        open(unit=11,file='Rp.txt',form='formatted',status='replace')
        write(11,'(9X,A,3X,A,3X,A,3X,A)') &
                'r', 'R''(numerical)', 'R''(carlson)', 'rel. error'
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            Rprime=dfridr(R_of_r,r,h,err)
            szloc=init_szlocal_class(t,r)
            call get_Rprime(szloc)
            err1=abs(1._dp-Rprime/szloc%f%Rp)
            write(11,fmt='(F10.4,F16.10,F14.10,E13.3)') r, Rprime, Rprime1, err1
        end do
        close(11)
    contains
        real(dp) function R_of_r(r)
            implicit none
            real(dp), intent(in) :: r
            szloc=init_szlocal_class(t,r)
            call get_R(szloc)
            R_of_r=szloc%f%R
        end function R_of_r
    end subroutine test_get_Rprime

    subroutine test_Rdotprime
        use utils, only : dfridr
        implicit none
        type(szlocal_class) :: szloc
        real(dp) :: t,h,r,err,err1
        real(dp) :: Rdp,Rdp1
        integer i
        h=0.1_dp
        t=-100._dp
        open(unit=14,file='Rdp.txt',form='formatted',status='replace')
        write(14,'(9X,A,3X,A,3X,A,3X,A)') &
                'r', 'Rdot''(numerical)', 'Rdot''(carlson)', 'rel. error'
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            Rdp=dfridr(Rdot_of_r,r,h,err)
            szloc=init_szlocal_class(t,r)
            call get_Rdotprime(szloc)
            err1=abs(1._dp-Rdp/szloc%f%Rdp)
            write(14,'(F10.4,E19.8,E17.8,E13.3)') r, Rdp, Rdp1, err1
        end do
        close(14)
    contains
        real(dp) function Rdot_of_r(r)
            real(dp), intent(in) :: r
            szloc=init_szlocal_class(t,r)
            call get_Rdot(szloc)
            Rdot_of_r=szloc%f%Rd
        end function Rdot_of_r
    end subroutine test_Rdotprime

    subroutine test_Hubble
        use cosmo_params, only : H0
        implicit none
        type(szlocal_class) :: szloc
        real(dp) :: t,r,H
        integer :: i,j
        open(unit=33,file='H.txt',form='formatted',status='replace')
        write(33,'(A)') 'Hubble function normalised to present H0'
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            write(33,'(f10.4)',advance='no') r
            do j=0,5
                t=-100._dp/5*j
                szloc=init_szlocal_class(t,r)
                call Hubble(szloc,H)
                write(33,'(f11.6)',advance='no') H/H0
            end do
            write(33,*)
        end do
    end subroutine test_Hubble

    subroutine test_discriminant
        implicit none
        type(szshell) :: shell
        real(dp) :: r,p,q,Dt
        integer i
        open(unit=9,file='discrmt.txt',form='formatted',status='replace')
        write(unit=9,fmt='(9X,A,16X,A,18X,A,7X,A)') 'r','p','q','Discriminant'
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            call init_shell(r,shell)
            p=shell%szc%p
            q=shell%szc%q
            Dt=-4._dp*p**3-27._dp*q**2
            write(unit=9,fmt='(f10.4,f17.9,f19.9,E19.9)') r,p,q,Dt
        end do
        close(9)
    end subroutine test_discriminant

    subroutine test_p3
        implicit none
        type(szshell) :: shell
        real(dp) :: r
        integer i
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            call init_shell(r,shell)
            write(*,*) r,shell%szc%p3
        end do
    end subroutine test_p3

    subroutine test_formula_S9
        use utils, only : rd,gaussint
        implicit none
        real(dp), dimension(4) :: a,b !quartic
        real(dp), dimension(:,:), allocatable :: am,bm
        real(dp), dimension(size(a)) :: X,Y
        real(dp), dimension(size(a),size(a)) :: U
        real(dp), dimension(size(a),size(a)) :: d
        real(dp), dimension(size(a),size(a)) :: d1
        real(dp) :: r,p,r2,p2 !upper and lower limits
        real(dp) :: yout,yout_integral
        complex(dp) :: a1,a2,a3
        integer :: na,nb
        a=(/0._dp,1._dp,1._dp,1._dp/)
        b=(/1._dp,1._dp,-1._dp,0._dp/)
        r=1._dp/SQRT2 !in u
        p=1._dp/SQRT3
        na=size(a)
        nb=size(b)
        if (na /= nb) STOP 'ensure a,b are same shape'
        r2=r**2 !in t
        p2=p**2
        X=sqrt(a+b*r2)
        Y=sqrt(a+b*p2)
        U=0._dp
        U(1,4)=(X(1)*X(4)*Y(2)*Y(3)+Y(1)*Y(4)*X(2)*X(3))/(r2-p2)
        allocate(am(na,1))
        allocate(bm(1,na))
        am=reshape(a,(/na,1/))
        bm=reshape(b,(/1,na/))
        d1=matmul(am,bm)
        d=d1-transpose(d1)
        a3=dcmplx(U(1,4)**2)
        a1=a3+1._dp
        a2=a3-1._dp
        yout=THIRD*rd(a1,a2,a3)+r*p/U(1,4)
        yout_integral=gaussint(fn_S9,p,r)
        write(*,*) yout,yout_integral
    contains
        real(dp) function fn_S9(ua)
            implicit none
            real(dp), intent(in) :: ua
            fn_S9=ua**2*1._dp/sqrt(1._dp-ua**4)
        end function fn_S9
    end subroutine test_formula_S9

    subroutine test_formula_213
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Evaluate integral (2.13) in Carlson 1987.                            !
        ! This routine differs from test_formula_213 as the full, general      !
        ! formula is used and quantities have not been expanded and simplified !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use szlocal
        implicit none
        type(szshell) :: shell
        complex(dp), dimension(4) :: a,b
        complex(dp) :: yout,yout2
        real(dp) :: r,ul,ll
        integer, dimension(3) :: p
        r=15._dp
        ul=r
        ll=0._dp
        call init_shell(r,shell)
        a(1) = dcmplx(0._dp)
        a(2:4) = -shell%szc%roots
        b = (/1._dp,1._dp,1._dp,1._dp /)
        p = (/1,2,3/)
        yout = f213(a,b,ul,ll)
        yout2 = f213_2(shell,ul,p)
        write(*,*) yout
        write(*,*) yout2
        ! ul=1._dp/SQRT2
        ! ll=0._dp
        ! a = (/0._dp,1._dp,1._dp,2._dp/)
        ! b = (/1._dp,-1._dp,1._dp,1._dp/)
        ! yout = f213(a,b,ul,ll)
        ! write(*,*) yout
    end subroutine test_formula_213

    complex(dp) function f213(a,b,ulim,llim)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Use the general formula to evaluate eq (2.13) !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use utils, only : rd,rf,rf_2
        implicit none
        integer, parameter :: n=4
        complex(dp), dimension(n), intent(in) :: a,b
        real(dp), intent(in) :: ulim,llim
        complex(dp), dimension(n) :: X,Y
        complex(dp), dimension(n,1) :: am
        complex(dp), dimension(1,n) :: bm
        complex(dp), dimension(n,n) :: d1,d
        complex(dp) :: U12,U13,U14,U122,U132,U142
        complex(dp) :: s1,s2,s3,t1,t2,t3
        X = sqrt(a+b*ulim)
        Y = sqrt(a+b*llim)
        U12 = (X(1)*X(2)*Y(3)*Y(4) + Y(1)*Y(2)*X(3)*X(4)) / (ulim-llim)
        U13 = (X(1)*X(3)*Y(2)*Y(4) + Y(1)*Y(3)*X(2)*X(4)) / (ulim-llim)
        U14 = (X(1)*X(4)*Y(2)*Y(3) + Y(1)*Y(4)*X(2)*X(3)) / (ulim-llim)
        U122 = U12**2
        U132 = U13**2
        U142 = U14**2
        am = reshape(a, (/n,1/) )
        bm = reshape(b, (/1,n/) )
        d1 = matmul(am,bm)
        d = d1 - transpose(d1)
        s1 = b(1) / d(1,4) - 2._dp * b(2) / d(2,4) - 2._dp * b(3) / d(3,4)
        t1 = 2._dp / 9._dp * s1 * d(1,2) * d(1,3) * rd(U122,U132,U142)
        s2 = b(4) * d(1,2) * d(1,3) / (d(1,4) * d(2,4) * d(3,4))
        t2 = -2._dp / 3._dp * s2 * rf(U122,U132,U142)
        s3 = -2._dp / 3._dp * b(4) / (d(2,4) * d(3,4))
        t3 = s3 * X(1) * X(2) * X(3) / X(4)**3
        f213 = t1 + t2 + t3
    end function f213

    complex(dp) function f213_2(shell,ulim,perm)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Only cyclic permutations of (1,2,3) allowed         !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use szelliptic
        use utils, only : rd,rf,rd_2
        implicit none
        type(szshell), intent(in) :: shell
        real(dp), intent(in) :: ulim
        integer, dimension(3), intent(in) :: perm
        complex(dp) :: x2,x3,x4,a2,a3,a4
        complex(dp) :: s,t1,t2,t3
        real(dp) :: xi
        x2 = shell%szc%roots(perm(1))
        x3 = shell%szc%roots(perm(2))
        x4 = shell%szc%roots(perm(3))
        a2 = 1._dp/ulim - 1._dp/x2
        a3 = 1._dp/ulim - 1._dp/x3
        a4 = 1._dp/ulim - 1._dp/x4
        xi = shell%szc%xi !sqrt(-x1.x2.x3)
        s = TWOTHIRD * xi * shell%szc%f(perm(3)) / x4**2
        t1 = THIRD * (1._dp / x2 + 1._dp / x3 + 1._dp / x4 &
                      - 3._dp * x4 / (x2 * x3)) * rd(a2,a3,a4)
        t2 = rf(a2,a3,a4)
        ! t3 = -sqrt(a2 * a3 / a4**3)
        t3 = -sqrt(a2) * sqrt(a3) / sqrt(a4)**3
        f213_2 = s * (t1 + t2 + t3)
    end function f213_2

    complex(dp) function f207(shell,ulim,perm)
        use szelliptic
        use utils, only : rd
        implicit none
        real(dp), intent(in) :: ulim
        type(szshell), intent(in) :: shell
        integer, dimension(3), intent(in) :: perm
        complex(dp) :: x2,x3,x4,xi,a2,a3,a4
        x2 = shell%szc%roots(perm(1))
        x3 = shell%szc%roots(perm(2))
        x4 = shell%szc%roots(perm(3))
        xi = shell%szc%xi
        a2 = 1._dp/ulim - 1._dp/x2
        a3 = 1._dp/ulim - 1._dp/x3
        a4 = 1._dp/ulim - 1._dp/x4
        f207 = -TWOTHIRD / (x4 * xi) * rd(a2,a3,a4)
    end function f207

    subroutine test_pX555
        use utils, only : gaussint
        implicit none
        type(szshell) :: shell
        real(dp) :: yout,yout1
        real(dp) :: r,ll,ul
        reaL(dp) :: a1,a2,b2,z2
        integer i
        ll=0._dp
        open(unit=17,file='p5555.txt',form='formatted',status='replace')
        open(unit=18,file='p3555.txt',form='formatted',status='replace')
        open(unit=19,file='p1555.txt',form='formatted',status='replace')
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            ul=r
            call init_shell(r,shell)
            a1=real(shell%szc%roots(1))
            a2=real(shell%szc%roots(2))
            b2=aimag(shell%szc%roots(2))
            z2=a2**2+b2**2
! All integrals ought to be real valued
            yout=szellip_5555(shell%szc,ul)
            yout1=gaussint(f5555,ll,ul)
            write(17,fmt='(F10.4,E17.9,E17.9)') r,yout,yout1
            yout=szellip_3555(shell%szc,ul)
            yout1=gaussint(f3555,ll,ul)
            write(18,fmt='(F10.4,E17.9,E17.9)') r,yout,yout1
            yout=szellip_1555(shell%szc,ul)
            yout1=gaussint(f1555,ll,ul)
            write(19,fmt='(F10.4,E17.9,E17.9)') r,yout,yout1
        end do
        close(19)
        close(18)
        close(17)
! The first root is real, the others complex conjugate pair
    contains
        real(dp) function f5555(x)
            implicit none
            real(dp), intent(in) :: x
            real(dp) q
! Check branch ambiguity if complex valued
            q=(x-a1)*(x**2-2._dp*a2*x+z2)
            if (q < 0._dp) STOP 'f5555: integrand not real'
            f5555=sqrt(x/q)**5
        end function f5555
        real(dp) function f3555(x)
            implicit none
            real(dp), intent(in) :: x
            real(dp) q
! Check branch ambiguity if complex valued
            q=(x-a1)*(x**2-2._dp*a2*x+z2)
            if (q < 0._dp) STOP 'f3555: integrand not real'
            f3555=sqrt(x)**3/sqrt(q)**5
        end function f3555
        real(dp) function f1555(x)
            implicit none
            real(dp), intent(in) :: x
            real(dp) q
! Check branch ambiguity if complex valued
            q=(x-a1)*(x**2-2._dp*a2*x+z2)
            if (q < 0._dp) STOP 'f1555: integrand not real'
            f1555=sqrt(x)/sqrt(q)**5
        end function f1555
    end subroutine test_pX555

    subroutine test_get_kpprime
        use utils, only : dfridr
        implicit none
        type(szlocal_class) :: szloc
        real(dp) :: h,err,r,kpp,kpp1,err1
        integer i
        real(dp), parameter :: t=0._dp
        h = 0.1_dp
        open(unit=13,file='kpp.txt',form='formatted',status='replace')
        write(13,'(9X,A,4X,A,6X,A,3X,A)') &
                'r','k''''(numerical)','k''''(carlson)','rel. error'
        do i=0,40
            r=1._dp+(200._dp-1._dp)/40*i
            kpp=dfridr(kp_of_r,r,h,err)
            szloc=init_szlocal_class(t,r)
            call get_kpprime(szloc)
            err1=abs(1._dp-kpp/szloc%f%kpp)
            write(13,'(F10.4,E18.8,E18.8,E13.3)') r,kpp,szloc%f%kpp,err1
        end do
        close(13)
    contains
        real(dp) function kp_of_r(r)
            implicit none
            real(dp), intent(in) :: r
            type(szlocal_class) :: szloc
            call get_kprime(szloc)
            kp_of_r=szloc%f%kp
        end function kp_of_r
    end subroutine test_get_kpprime

    subroutine test_rayshoot_dr
        use cosmo_params, only : r_obs,theta_obs
        use sznullgeo
        implicit none
        real(dp) :: RA,DEC,red,obs_pos(3)
        RA=0.78539816339744828_dp
        DEC=0.10210642238260403_dp
        obs_pos=(/r_obs,theta_obs,0.543*PI/)
        call rayshoot_dr(obs_pos,RA,DEC,red,iexit=0,val_exit=-400._dp)
        print *,'redshift', red
    end subroutine test_rayshoot_dr

end module test_module


text:
SUBROUTINE locpt (x0, y0, x, y, n, l, m)
!-----------------------------------------------------------------------
! GIVEN A POLYGONAL LINE CONNECTING THE VERTICES (X(I),Y(I)) (I = 1,...,N)
! TAKEN IN THIS ORDER.  IT IS ASSUMED THAT THE POLYGONAL PATH IS A LOOP,
! WHERE (X(N),Y(N)) = (X(1),Y(1)) OR THERE IS AN ARC FROM (X(N),Y(N)) TO
! (X(1),Y(1)).  N.B. The polygon may cross itself any number of times.

! (X0,Y0) IS AN ARBITRARY POINT AND L AND M ARE VARIABLES.
! On output, L AND M ARE ASSIGNED THE FOLLOWING VALUES ...

!    L = -1   IF (X0,Y0) IS OUTSIDE THE POLYGONAL PATH
!    L =  0   IF (X0,Y0) LIES ON THE POLYGONAL PATH
!    L =  1   IF (X0,Y0) IS INSIDE THE POLYGONAL PATH

! M = 0 IF (X0,Y0) IS ON OR OUTSIDE THE PATH.  IF (X0,Y0) IS INSIDE THE
! PATH THEN M IS THE WINDING NUMBER OF THE PATH AROUND THE POINT (X0,Y0).

! Fortran 66 version by A.H. Morris
! Converted to ELF90 compatibility by Alan Miller, 15 February 1997

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

IMPLICIT NONE
REAL, INTENT(IN)     :: x0, y0, x(:), y(:)
INTEGER, INTENT(IN)  :: n
INTEGER, INTENT(OUT) :: l, m

!     Local variables
INTEGER :: i, n0
REAL    :: angle, eps, pi, pi2, sum, theta, theta1, thetai, tol, u, v

!     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
!            SMALLEST NUMBER SUCH THAT 1.0 + EPS > 1.0

eps = EPSILON(1.0)

!-----------------------------------------------------------------------
n0 = n
IF (x(1) == x(n) .AND. y(1) == y(n)) n0 = n - 1
pi = ATAN2(0.0, -1.0)
pi2 = 2.0*pi
tol = 4.0*eps*pi
l = -1
m = 0

u = x(1) - x0
v = y(1) - y0
IF (u == 0.0 .AND. v == 0.0) GO TO 20
IF (n0 < 2) RETURN
theta1 = ATAN2(v, u)

sum = 0.0
theta = theta1
DO i = 2, n0
  u = x(i) - x0
  v = y(i) - y0
  IF (u == 0.0 .AND. v == 0.0) GO TO 20
  thetai = ATAN2(v, u)

  angle = ABS(thetai - theta)
  IF (ABS(angle - pi) < tol) GO TO 20
  IF (angle > pi) angle = angle - pi2
  IF (theta > thetai) angle = -angle
  sum = sum + angle
  theta = thetai
END DO

angle = ABS(theta1 - theta)
IF (ABS(angle - pi) < tol) GO TO 20
IF (angle > pi) angle = angle - pi2
IF (theta > theta1) angle = -angle
sum = sum + angle

!     SUM = 2*PI*M WHERE M IS THE WINDING NUMBER

m = ABS(sum)/pi2 + 0.2
IF (m == 0) RETURN
l = 1
IF (sum < 0.0) m = -m
RETURN

!     (X0, Y0) IS ON THE BOUNDARY OF THE PATH

20 l = 0
RETURN
END SUBROUTINE locpt

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
module rstm_ssglrrw_artificial_viscosity_operator
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ONE,TWO,HALF

    use type_operator,          only: operator_t
    use type_properties,        only: properties_t
    use type_chidg_worker,      only: chidg_worker_t
    use DNAD_D
    use mod_rstm_ssglrrw, only: rstm_ssglrrw_avc
    use ieee_arithmetic
    implicit none

    private

    !> Volume flux for Fluid laplacian_av Terms.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   9/23/2016
    !!
    !!
    !------------------------------------------------------------------------------
    type, extends(operator_t), public :: rstm_ssglrrw_artificial_viscosity_operator_t

    contains

        procedure   :: init
        procedure   :: compute

    end type rstm_ssglrrw_artificial_viscosity_operator_t
    !******************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   9/23/2016
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(rstm_ssglrrw_artificial_viscosity_operator_t),   intent(inout)      :: self

        ! Set operator name
        call self%set_name('RSTMSSGLRRW Artificial Viscosity Volume Operator')

        ! Set operator type
        call self%set_operator_type('Volume Diffusive Operator')

        ! Set operator equations
        call self%add_primary_field('Density * Reynolds-11'   )
        call self%add_primary_field('Density * Reynolds-22'   )
        call self%add_primary_field('Density * Reynolds-33'   )
        call self%add_primary_field('Density * Reynolds-12'   )
        call self%add_primary_field('Density * Reynolds-13'   )
        call self%add_primary_field('Density * Reynolds-23'   )

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

    !> Volume flux routine for Fluid laplacian_av Terms.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   9/23/2016
    !!
    !!
    !!------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(rstm_ssglrrw_artificial_viscosity_operator_t), intent(inout)   :: self
        type(chidg_worker_t),                   intent(inout)   :: worker
        class(properties_t),                    intent(inout)   :: prop

        type(AD_D), allocatable, dimension(:) ::    &
            mu_neg, grad1, grad2, grad3, &
            flux_1, flux_2, flux_3

        real(rk),   allocatable, dimension(:)   :: r

        integer(ik) :: p
        real(rk)    :: h(3), sigma
        real(rk),   allocatable :: h_smooth(:, :)

        !
        ! Get Model fields:
        !   Second Coefficient of Viscosity
        !   Thermal Conductivity
        !

        h_smooth = worker%h_smooth('element')

        h = worker%element_size('interior')
        p = worker%solution_order('interior')
        h = h/real(p+1,rk)
        h_smooth = h_smooth/real(p+1,rk)
        sigma = 0.5_rk

        mu_neg = worker%get_field('RSTM AV-11'   , 'value', 'element')
        grad1 = worker%get_field('Density * Reynolds-11'   , 'grad1', 'element')
        grad2 = worker%get_field('Density * Reynolds-11'   , 'grad2', 'element')
        grad3 = worker%get_field('Density * Reynolds-11'   , 'grad3', 'element')
        !flux_1 = -(h(1)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad1
        !flux_2 = -(h(2)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad2
        !flux_3 = -(h(3)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad3

        flux_1 = -(h_smooth(:,1)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad1
        flux_2 = -(h_smooth(:,2)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad2
        flux_3 = -(h_smooth(:,3)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad3

        call worker%integrate_volume_flux('Density * Reynolds-11','Diffusion',flux_1,flux_2,flux_3)

        mu_neg = worker%get_field('RSTM AV-22'   , 'value', 'element')
        grad1 = worker%get_field('Density * Reynolds-22'   , 'grad1', 'element')
        grad2 = worker%get_field('Density * Reynolds-22'   , 'grad2', 'element')
        grad3 = worker%get_field('Density * Reynolds-22'   , 'grad3', 'element')
        !flux_1 = -(h(1)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad1
        !flux_2 = -(h(2)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad2
        !flux_3 = -(h(3)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad3

        flux_1 = -(h_smooth(:,1)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad1
        flux_2 = -(h_smooth(:,2)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad2
        flux_3 = -(h_smooth(:,3)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad3

        call worker%integrate_volume_flux('Density * Reynolds-22','Diffusion',flux_1,flux_2,flux_3)

        mu_neg = worker%get_field('RSTM AV-33'   , 'value', 'element')
        grad1 = worker%get_field('Density * Reynolds-33'   , 'grad1', 'element')
        grad2 = worker%get_field('Density * Reynolds-33'   , 'grad2', 'element')
        grad3 = worker%get_field('Density * Reynolds-33'   , 'grad3', 'element')
        !flux_1 = -(h(1)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad1
        !flux_2 = -(h(2)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad2
        !flux_3 = -(h(3)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad3

        flux_1 = -(h_smooth(:,1)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad1
        flux_2 = -(h_smooth(:,2)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad2
        flux_3 = -(h_smooth(:,3)**TWO*rstm_ssglrrw_avc+mu_neg*sigma)*grad3

        call worker%integrate_volume_flux('Density * Reynolds-33','Diffusion',flux_1,flux_2,flux_3)

        grad1 = worker%get_field('Density * Reynolds-12'   , 'grad1', 'element')
        grad2 = worker%get_field('Density * Reynolds-12'   , 'grad2', 'element')
        grad3 = worker%get_field('Density * Reynolds-12'   , 'grad3', 'element')
        !flux_1 = -(h(1)**TWO*rstm_ssglrrw_avc)*grad1
        !flux_2 = -(h(2)**TWO*rstm_ssglrrw_avc)*grad2
        !flux_3 = -(h(3)**TWO*rstm_ssglrrw_avc)*grad3

        flux_1 = -(h_smooth(:,1)**TWO*rstm_ssglrrw_avc)*grad1
        flux_2 = -(h_smooth(:,2)**TWO*rstm_ssglrrw_avc)*grad2
        flux_3 = -(h_smooth(:,3)**TWO*rstm_ssglrrw_avc)*grad3
        call worker%integrate_volume_flux('Density * Reynolds-12','Diffusion',flux_1,flux_2,flux_3)

        grad1 = worker%get_field('Density * Reynolds-13'   , 'grad1', 'element')
        grad2 = worker%get_field('Density * Reynolds-13'   , 'grad2', 'element')
        grad3 = worker%get_field('Density * Reynolds-13'   , 'grad3', 'element')
        !flux_1 = -(h(1)**TWO*rstm_ssglrrw_avc)*grad1
        !flux_2 = -(h(2)**TWO*rstm_ssglrrw_avc)*grad2
        !flux_3 = -(h(3)**TWO*rstm_ssglrrw_avc)*grad3

        flux_1 = -(h_smooth(:,1)**TWO*rstm_ssglrrw_avc)*grad1
        flux_2 = -(h_smooth(:,2)**TWO*rstm_ssglrrw_avc)*grad2
        flux_3 = -(h_smooth(:,3)**TWO*rstm_ssglrrw_avc)*grad3
        call worker%integrate_volume_flux('Density * Reynolds-13','Diffusion',flux_1,flux_2,flux_3)

        grad1 = worker%get_field('Density * Reynolds-23'   , 'grad1', 'element')
        grad2 = worker%get_field('Density * Reynolds-23'   , 'grad2', 'element')
        grad3 = worker%get_field('Density * Reynolds-23'   , 'grad3', 'element')
        !flux_1 = -(h(1)**TWO*rstm_ssglrrw_avc)*grad1
        !flux_2 = -(h(2)**TWO*rstm_ssglrrw_avc)*grad2
        !flux_3 = -(h(3)**TWO*rstm_ssglrrw_avc)*grad3

        flux_1 = -(h_smooth(:,1)**TWO*rstm_ssglrrw_avc)*grad1
        flux_2 = -(h_smooth(:,2)**TWO*rstm_ssglrrw_avc)*grad2
        flux_3 = -(h_smooth(:,3)**TWO*rstm_ssglrrw_avc)*grad3
        call worker%integrate_volume_flux('Density * Reynolds-23','Diffusion',flux_1,flux_2,flux_3)

    end subroutine compute
    !*********************************************************************************************************

end module rstm_ssglrrw_artificial_viscosity_operator

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG