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 StrmHydrograph
  USE MessageLogger          , ONLY: SetLastMessage                , &
                                     LogMessage                    , &
                                     MessageArray                  , &
                                     iFatal                        , &
                                     iInfo
  USE IOInterface            , ONLY: GenericFileType               , &
                                     iGetFileType_FromName         , &
                                     DSS                           , &
                                     UNKNOWN                       , &
                                     iDataset
  USE GeneralUtilities       , ONLY: LocateInList                  , &
                                     StripTextUntilCharacter       , &
                                     IntToText                     , &
                                     ArrangeText                   , &
                                     CleanSpecialCharacters        , &
                                     EstablishAbsolutePathFilename , &
                                     PrepareTitle                  , &
                                     UpperCase                     , &
                                     ConvertID_To_Index
  USE TimeSeriesUtilities    , ONLY: TimeStepType                  , &
                                     IncrementTimeStamp
  USE Package_Misc           , ONLY: RealTSDataInFileType          , &
                                     PrepareTSDOutputFile          , &
                                     ReadTSData                    , &
                                     f_iDataUnitType_Length        , &
                                     f_iDataUnitType_Volume
  USE Class_StrmState
  IMPLICIT NONE

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

  ! -------------------------------------------------------------
  ! --- PUBLIC ENTITIES
  ! -------------------------------------------------------------
  PRIVATE
  PUBLIC :: StrmHydrographType , &
            iHydFlow           , &
            iHydStage          , &
            iHydBoth

  ! -------------------------------------------------------------
  ! --- TYPES OF STREAM HYDROGRAPHS
  ! -------------------------------------------------------------
  INTEGER,PARAMETER :: iHydStage = 1 , &
                       iHydFlow  = 0 , &
                       iHydBoth  = 2

  ! -------------------------------------------------------------
  ! --- STREAM HYDROGRAPH PRINT DATA TYPE
  ! -------------------------------------------------------------
  TYPE StrmHydrographType
    TYPE(GenericFileType)         :: HydFile
    TYPE(RealTSDataInFileType)    :: HydFile_ForInquiry
    LOGICAL                       :: HydFile_Defined    = .FALSE.
    INTEGER                       :: iHydType           = iHydFlow
    REAL(8)                       :: FactFlow           = 1.0
    CHARACTER(LEN=10)             :: UnitFlow           = ''
    REAL(8)                       :: FactElev           = 1.0
    CHARACTER(LEN=10)             :: UnitElev           = ''
    INTEGER                       :: NHyd               = 0
    INTEGER,ALLOCATABLE           :: iHydNodes(:)
    CHARACTER(LEN=30),ALLOCATABLE :: cHydNames(:)
  CONTAINS
    PROCEDURE,PASS :: New
    PROCEDURE,PASS :: Kill
    PROCEDURE,PASS :: ReadStrmHydrograph_AtNode
    PROCEDURE,PASS :: PrintResults
    PROCEDURE,PASS :: IsOutFileDefined
    PROCEDURE,PASS :: IsHydrographAtNodeRequested
    PROCEDURE,PASS :: GetFileName
    PROCEDURE,PASS :: GetNHydrographs
    PROCEDURE,PASS :: GetHydrographNodes
    PROCEDURE,PASS :: GetCoordinates
    PROCEDURE,PASS :: GetNames
    PROCEDURE,PASS :: Transfer_To_HDF
  END TYPE StrmHydrographType

  ! -------------------------------------------------------------
  ! --- MISC. ENTITIES
  ! -------------------------------------------------------------
  INTEGER,PARAMETER                   :: ModNameLen    = 16
  CHARACTER(LEN=ModNameLen),PARAMETER :: ModName       = 'StrmHydrograph::'

CONTAINS

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

  ! -------------------------------------------------------------
  ! --- INSTANTIATE STREAM HYDROGRAPH PRINT DATA
  ! -------------------------------------------------------------
  SUBROUTINE New(StrmHyd,IsRoutedStreams,IsForInquiry,cWorkingDirectory,NStrmNodes,iStrmNodeIDs,TimeStep,InFile,iStat)
    CLASS(StrmHydrographType),INTENT(OUT) :: StrmHyd
    LOGICAL,INTENT(IN)                    :: IsRoutedStreams,IsForInquiry
    CHARACTER(LEN=*),INTENT(IN)           :: cWorkingDirectory
    INTEGER,INTENT(IN)                    :: NStrmNodes,iStrmNodeIDs(NStrmNodes)
    TYPE(TimeStepType),INTENT(IN)         :: TimeStep
    TYPE(GenericFileType)                 :: InFile
    INTEGER,INTENT(OUT)                   :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+3) :: ThisProcedure = ModName // 'New'
    INTEGER                     :: NHyd,indx,ErrorCode,iNode
    CHARACTER                   :: ALine*1000,cHydOutFile*1000,cNumber*7
    INTEGER,ALLOCATABLE         :: iHydNodeIDs(:)
    CHARACTER(:),ALLOCATABLE    :: cAbsPathFileName

    !Read data
    CALL InFile%ReadData(NHyd,iStat)  ;  IF (iStat .EQ. -1) RETURN  ;  StrmHyd%NHyd = NHyd
    IF (NHyd .EQ. 0) THEN
        DO indx=1,6
            CALL InFile%ReadData(ALine,iStat)
            IF (iStat .EQ. -1) RETURN
        END DO
        RETURN
    END IF
    CALL InFile%ReadData(StrmHyd%iHydType,iStat)  ;  IF (iStat .EQ. -1) RETURN
    CALL InFile%ReadData(StrmHyd%FactFlow,iStat)  ;  IF (iStat .EQ. -1) RETURN
    CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN  ; CALL CleanSpecialCharacters(ALine) ; ALine = StripTextUntilCharacter(ALine,'/')
    StrmHyd%UnitFlow = ADJUSTL(TRIM(ALine))
    CALL InFile%ReadData(StrmHyd%FactElev,iStat)  ;  IF (iStat .EQ. -1) RETURN
    CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN ; CALL CleanSpecialCharacters(ALine) ; ALine = StripTextUntilCharacter(ALine,'/')
    StrmHyd%UnitElev = ADJUSTL(TRIM(ALine))

    !Output file name
    CALL InFile%ReadData(cHydOutFile,iStat)  ;  IF (iStat .EQ. -1) RETURN ; cHydOutFile = StripTextUntilCharacter(cHydOutFile,'/') ; CALL CleanSpecialCharacters(cHydOutFile)
    IF (cHydOutFile .EQ. '') THEN
        IF (IsRoutedStreams) CALL LogMessage('Stream hydrograph printing is suppressed because an output file name is not specified!',iInfo,ThisProcedure)
        DO indx=1,NHyd
            CALL InFile%ReadData(iNode,iStat)
            IF (iStat .EQ. -1) RETURN
        END DO
        StrmHyd%NHyd = 0
        RETURN
    END IF

    !Allocate memory
    ALLOCATE (StrmHyd%iHydNodes(StrmHyd%NHyd)  ,  StrmHyd%cHydNames(StrmHyd%NHyd) , iHydNodeIDs(StrmHyd%NHyd))
    StrmHyd%cHydNames = ''

    !Read the node numbers and hydrograph names for hydrograph printing
    DO indx=1,NHyd
        CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN  ;  CALL CleanSpecialCharacters(ALine)  ;  ALine = StripTextUntilCharacter(ALine,'/')  ;  ALine = ADJUSTL(ALine)
        READ (ALine,*,IOSTAT=ErrorCode) iHydNodeIDs(indx)
        CALL ConvertID_To_Index(iHydNodeIDs(indx),iStrmNodeIDs,iNode)
        IF (iNode .EQ. 0) THEN
            CALL SetLastMessage('Stream node ID '//TRIM(IntToText(iHydNodeIDs(indx)))//' listed for hydrograph printing is not in the model!',iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
        cNumber = ''  ;  cNumber = IntToText(iHydNodeIDs(indx))  ;  cNumber = ADJUSTL(cNumber)
        StrmHyd%cHydNames(indx) = ADJUSTL(ALine(LEN_TRIM(cNumber)+1:))  ;  CALL CleanSpecialCharacters(StrmHyd%cHydNames(indx))
        IF (ErrorCode .NE. 0) THEN
            MessageArray(1) = 'Error in reading stream hydrograph print data!'
            MessageArray(2) = TRIM(ALine)
            CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF
        StrmHyd%iHydNodes(indx) = iNode
    END DO

    IF (IsRoutedStreams) THEN
        !Prepare hydrograph file
        CALL EstablishAbsolutePathFileName(TRIM(ADJUSTL(cHydOutFile)),cWorkingDirectory,cAbsPathFileName)
        IF (IsForInquiry) THEN
            CALL PrepStrmHydFile_ForInquiry(TimeStep,cAbsPathFileName,StrmHyd%iHydType,StrmHyd%NHyd,iHydNodeIDs,StrmHyd%HydFile_ForInquiry,iStat)
        ELSE
            CALL PrepStrmHydFile(TimeStep,cAbsPathFileName,StrmHyd%UnitFlow,StrmHyd%UnitElev,StrmHyd%iHydType,StrmHyd%NHyd,iHydNodeIDs,StrmHyd%HydFile,iStat)
        END IF
        IF (iStat .EQ. -1) RETURN

        !Set the flag that checks if the hydrograph output file is defined
        StrmHyd%HydFile_Defined = .TRUE.
    ELSE
        CALL StrmHyd%Kill()
    END IF

  END SUBROUTINE New

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

  ! -------------------------------------------------------------
  ! --- KILL STREAM HYDROGRAPH PRINT DATA OBJECT
  ! -------------------------------------------------------------
  SUBROUTINE Kill(StrmHyd)
    CLASS(StrmHydrographType) :: StrmHyd

    !Local variables
    INTEGER                  :: ErrorCode
    TYPE(StrmHydrographType) :: Dummy

    !Deallocate array attributes
    DEALLOCATE (StrmHyd%iHydNodes , StrmHyd%cHydNames , STAT=ErrorCode)

    !Close output (or input) file
    CALL StrmHyd%HydFile%Kill()
    CALL StrmHyd%HydFile_ForInquiry%Close()

    !Set attributes to their default values
    SELECT TYPE (p => StrmHyd)
        TYPE IS (StrmHydrographType)
            p = Dummy
    END SELECT

  END SUBROUTINE Kill

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

  ! -------------------------------------------------------------
  ! --- GET OUTPUT FILE NAME
  ! -------------------------------------------------------------
  SUBROUTINE GetFileName(StrmHyd,cFileName)
    CLASS(StrmHydrographType)            :: StrmHyd
    CHARACTER(:),ALLOCATABLE,INTENT(OUT) :: cFileName

    !Local variables
    INTEGER :: ErrorCode

    !Initialize
    DEALLOCATE (cFileName , STAT=ErrorCode)

    IF (StrmHyd%HydFile_Defined) THEN
        IF (StrmHyd%HydFile%iGetFileType() .NE. UNKNOWN) THEN
            CALL StrmHyd%HydFile%GetName(cFileName)
        ELSE
            CALL StrmHyd%HydFile_ForInquiry%GetFileName(cFileName)
        END IF
    END IF

  END SUBROUTINE GetFileName

  ! -------------------------------------------------------------
  ! --- GET STREAM NODE INDICES WHERE HYDROGRAPHS ARE PRINTED
  ! -------------------------------------------------------------
  PURE SUBROUTINE GetHydrographNodes(StrmHyd,iNodes)
    CLASS(StrmHydrographType),INTENT(IN) :: StrmHyd
    INTEGER,INTENT(OUT)                  :: iNodes(:)

    IF (STrmHyd%NHyd .GT. 0) iNodes = StrmHyd%iHydNodes

  END SUBROUTINE GetHydrographNodes

  ! -------------------------------------------------------------
  ! --- GET NUMBER OF HYDROGRAPHS
  ! -------------------------------------------------------------
  FUNCTION GetNHydrographs(StrmHyd) RESULT(NHydrographs)
    CLASS(StrmHydrographType) :: StrmHyd
    INTEGER                   :: NHydrographs

    NHydrographs = StrmHyd%NHyd

  END FUNCTION GetNHydrographs

  ! -------------------------------------------------------------
  ! --- GET HYDROGRAPH COORDINATES
  ! -------------------------------------------------------------
  SUBROUTINE GetCoordinates(StrmHyd,iGWNodes,GridX,GridY,XHyd,YHyd)
    CLASS(StrmHydrographType) :: StrmHyd
    INTEGER,INTENT(IN)        :: iGWNodes(:)
    REAL(8),INTENT(IN)        :: GridX(:),GridY(:)
    REAL(8),INTENT(OUT)       :: XHyd(:),YHyd(:)

    !Local variables
    INTEGER :: indx,iGWNode

    DO indx=1,StrmHyd%NHyd
        iGWNode    = iGWNodes(StrmHyd%iHydNodes(indx))
        XHyd(indx) = GridX(iGWNode)
        YHyd(indx) = GridY(iGWNode)
    END DO

  END SUBROUTINE GetCoordinates

  ! -------------------------------------------------------------
  ! --- GET HYDROGRAPH NAMES
  ! -------------------------------------------------------------
  PURE SUBROUTINE GetNames(StrmHyd,cNamesList)
    CLASS(StrmHydrographType),INTENT(IN) :: StrmHyd
    CHARACTER(LEN=*),INTENT(OUT)         :: cNamesList(:)   !Assumes array is previously dimensioned according to the number of hydrographs

    cNamesList = StrmHyd%cHydNames

  END SUBROUTINE GetNames

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

  ! -------------------------------------------------------------
  ! --- READ STREAM HYDROGRAPH FOR A STREAM NODE FOR A TIME RANGE
  ! -------------------------------------------------------------
  SUBROUTINE ReadStrmHydrograph_AtNode(StrmHyd,iNodeID,iNode,iHydType,cOutputBeginDateAndTime,cOutputEndDateAndTime,rFact_LT,rFact_VL,iDataUnitType,nActualOutput,rOutputDates,rOutputValues,iStat)
    CLASS(StrmHydrographType)   :: StrmHyd
    INTEGER,INTENT(IN)          :: iNodeID,iNode,iHydType
    CHARACTER(LEN=*),INTENT(IN) :: cOutputBeginDateAndTime,cOutputEndDateAndTime
    REAL(8),INTENT(IN)          :: rFact_LT,rFact_VL
    INTEGER,INTENT(OUT)         :: iDataUnitType,nActualOutput
    REAL(8),INTENT(OUT)         :: rOutputDates(:),rOutputValues(:)
    INTEGER,INTENT(OUT)         :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+25),PARAMETER :: ThisProcedure = ModName // 'ReadStrmHydrograph_AtNode'
    INTEGER                                :: ErrorCode,iHydIndex
    REAL(8)                                :: rEffectiveFactor

    !Convert node ID to hydrograph index
    iHydIndex = LocateInList(iNode,StrmHyd%iHydNodes)
    IF (iHydIndex .EQ. 0) THEN
        CALL SetLastMessage('Stream node ID '//TRIM(IntToText(iNodeID))//' does not have a hydrograph printed as model results!',iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF

    !Read data
    IF (StrmHyd%iHydType .EQ. iHydBoth) THEN
        IF (iHydType .EQ. iHydFlow) THEN
            CALL StrmHyd%HydFile_ForInquiry%ReadData(iHydIndex,cOutputBeginDateAndTime,cOutputEndDateAndTime,nActualOutput,rOutputValues,rOutputDates,ErrorCode,iStat)
        ELSE
            CALL StrmHyd%HydFile_ForInquiry%ReadData(2*StrmHyd%NHyd+iHydIndex,cOutputBeginDateAndTime,cOutputEndDateAndTime,nActualOutput,rOutputValues,rOutputDates,ErrorCode,iStat)
        END IF
    ELSE
        CALL StrmHyd%HydFile_ForInquiry%ReadData(iHydIndex,cOutputBeginDateAndTime,cOutputEndDateAndTime,nActualOutput,rOutputValues,rOutputDates,ErrorCode,iStat)
    END IF
    IF (iStat .EQ. -1) RETURN

    !Convert flow/elevation data to simulation units, then to output units
    IF (StrmHyd%iHydType .EQ. iHydFlow) THEN
        iDataUnitType    = f_iDataUnitType_Volume
        rEffectiveFactor = rFact_VL / StrmHyd%FactFlow
        IF (rEffectiveFactor .NE. 1d0) rOutputValues(1:nActualOutput) = rOutputValues(1:nActualOutput) * rEffectiveFactor
    ELSEIF (StrmHyd%iHydType .EQ. iHydStage) THEN
        iDataUnitType    = f_iDataUnitType_Length
        rEffectiveFactor = rFact_LT / StrmHyd%FactElev
        IF (rEffectiveFactor .NE. 1d0) rOutputValues(1:nActualOutput) = rOutputValues(1:nActualOutput) * rEffectiveFactor
    ELSE
        IF (iHydType .EQ. iHydFlow) THEN
            iDataUnitType    = f_iDataUnitType_Volume
            rEffectiveFactor = rFact_VL / StrmHyd%FactFlow
        ELSE
            iDataUnitType    = f_iDataUnitType_Length
            rEffectiveFactor = rFact_LT / StrmHyd%FactElev
        END IF
        IF (rEffectiveFactor .NE. 1d0) rOutputValues(1:nActualOutput) = rOutputValues(1:nActualOutput) * rEffectiveFactor
    END IF

    !Rewind file
    CALL StrmHyd%HydFile_ForInquiry%File%RewindFile_To_BeginningOfTSData(iStat)

  END SUBROUTINE ReadStrmHydrograph_AtNode

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** DATA WRITERS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- PRINT OUT HYDROGRAPHS
  ! -------------------------------------------------------------
  SUBROUTINE PrintResults(StrmHyd,StrmState,BottomElev,TimeStep,lEndOfSimulation)
    CLASS(StrmHydrographType)      :: StrmHyd
    TYPE(StrmStateType),INTENT(IN) :: StrmState(:)
    REAL(8),INTENT(IN)             :: BottomElev(:)
    TYPE(TimeStepType),INTENT(IN)  :: TimeStep
    LOGICAL,INTENT(IN)             :: lEndOfSimulation

    !Local variables
    INTEGER           :: NHyd
    REAL(8)           :: rValues(2*StrmHyd%NHyd)
    CHARACTER(LEN=21) :: SimulationTime

    !Initialize
    NHyd = StrmHyd%NHyd

    !Create the simulation time
    IF (TimeStep%TrackTime) THEN
        SimulationTime=ADJUSTL(TimeStep%CurrentDateAndTime)
    ELSE
        WRITE(SimulationTime,'(F10.2,1X,A10)') TimeStep%CurrentTime,ADJUSTL(TimeStep%Unit)
    END IF

    !Print-out
    SELECT CASE (StrmHyd%iHydType)
      !Flow hydrograph
      CASE (iHydFlow)
        rValues(1:NHyd) = StrmState(StrmHyd%iHydNodes)%Flow * StrmHyd%FactFlow
        CALL StrmHyd%HydFile%WriteData(SimulationTime,rValues(1:NHyd),FinalPrint=lEndOfSimulation)

      !Stage hydrograph
      CASE (iHydStage)
        rValues(1:NHyd) = (StrmState(StrmHyd%iHydNodes)%Head - BottomElev(StrmHyd%iHydNodes)) * StrmHyd%FactElev
        CALL StrmHyd%HydFile%WriteData(SimulationTime,rValues(1:NHyd),FinalPrint=lEndOfSimulation)

      !Both flow and stage hydrographs
      CASE (iHydBoth)
        rValues(1:NHyd)  = StrmState(StrmHyd%iHydNodes)%Flow * StrmHyd%FactFlow
        rValues(NHyd+1:) = (StrmState(StrmHyd%iHydNodes)%Head - BottomElev(StrmHyd%iHydNodes)) * StrmHyd%FactElev
        CALL StrmHyd%HydFile%WriteData(SimulationTime,rValues,FinalPrint=lEndOfSimulation)

    END SELECT

  END SUBROUTINE PrintResults

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

  ! -------------------------------------------------------------
  ! --- CHECK IF OUTPUT FILE IS DEFINED
  ! -------------------------------------------------------------
  PURE FUNCTION IsOutFileDefined(StrmHyd) RESULT(lDefined)
    CLASS(StrmHydrographType),INTENT(IN) :: StrmHyd
    LOGICAL                              :: lDefined

    lDefined = StrmHyd%HydFile_Defined

  END FUNCTION IsOutFileDefined

  ! -------------------------------------------------------------
  ! --- CHECK IF HYDROGRAPH AT A SPECIFIC NODE IS REQUESTED
  ! -------------------------------------------------------------
  FUNCTION IsHydrographAtNodeRequested(StrmHyd,iStrmNode) RESULT(lRequested)
    CLASS(StrmHydrographType),INTENT(IN) :: StrmHyd
    INTEGER,INTENT(IN)                   :: iStrmNode
    LOGICAL                              :: lRequested

    IF (LocateInList(iStrmNode,StrmHyd%iHydNodes) .GT. 0) THEN
        lRequested = .TRUE.
    ELSE
        lRequested = .FALSE.
    END IF

  END FUNCTION IsHydrographAtNodeRequested

  ! -------------------------------------------------------------
  ! --- NEW STREAM HYDROGRAPH OUTPUT FILE
  ! -------------------------------------------------------------
  SUBROUTINE PrepStrmHydFile(TimeStep,cOutFileName,cUnitFlow,cUnitElev,iHydType,NHyd,iHydNodes,HydFile,iStat)
    TYPE(TimeStepType),INTENT(IN) :: TimeStep
    CHARACTER(LEN=*),INTENT(IN)   :: cOutFileName,cUnitFlow,cUnitElev
    INTEGER,INTENT(IN)            :: iHydType,NHyd,iHydNodes(NHyd)
    TYPE(GenericFileType)         :: HydFile
    INTEGER,INTENT(OUT)           :: iStat

    !Local variables
    CHARACTER             :: Text*20,FormatSpec*500,FPart(1)*32,HeaderFormat(3)*500,  &
                             WorkArray(2)*3000,TitleLines(1)*325
    INTEGER               :: NColumnsOfData,NRowsOfData,indx,ErrorCode
    LOGICAL               :: OverwriteNColumnsOfData,PrintColumnNo
    CHARACTER,ALLOCATABLE :: cPart*32(:),DataUnit*10(:),DataType*10(:),Header*50(:,:)
    INTEGER,ALLOCATABLE   :: iHydNodes_Local(:)

    !Initialize
    iStat = 0

    !Open file
    CALL HydFile%New(FileName=cOutFileName,InputFile=.FALSE.,IsTSFile=.TRUE.,Descriptor='stream hydrographs output file',iStat=iStat)
    IF (iStat .EQ. -1) RETURN

    !Prepare file
    IF (iHydType .EQ. iHydBoth) THEN
        NColumnsOfData = 2 * NHyd
        ALLOCATE (CPart(NColumnsOfData) , DataUnit(NColumnsOfData) , DataType(NColumnsOfData) , Header(3,1+NColumnsOfData) , iHydNodes_Local(NColumnsOfData))
        CPart                    = ''
        CPart(1:NHyd)            = ADJUSTL('FLOW')
        CPart(NHyd+1:)           = ADJUSTL('FLOW_DEPTH')
        DataUnit                 = ''
        DataUnit(1:NHyd)         = ADJUSTL(cUnitFlow)
        DataUnit(NHyd+1:)        = ADJUSTL(cUnitElev)
        DataType                 = ''
        DataType(1:NHyd)         = 'PER-AVER'
        DataType(NHyd+1:)        = 'INST-VAL'
        iHydNodes_Local(1:NHyd)  = iHydNodes
        iHydNodes_Local(NHyd+1:) = iHydNodes
    ELSE
        NColumnsOfData = NHyd
        ALLOCATE (CPart(NColumnsOfData) , DataUnit(NColumnsOfData) , DataType(NColumnsOfData) , Header(3,1+NColumnsOfData) , iHydNodes_Local(NColumnsOfData))
        iHydNodes_Local = iHydNodes
        CPart           = ''
        DataUnit        = ''
        DataType        = ''
        IF (iHydType .EQ. iHydFlow) THEN
            CPart    = ADJUSTL('FLOW')
            DataUnit = ADJUSTL(cUnitFlow)
            DataType = 'PER-AVER'
        ELSE
            CPart    = ADJUSTL('FLOW_DEPTH')
            DataUnit = ADJUSTL(cUnitElev)
            DataType = 'INST-VAL'
        END IF
    END IF
    Text                    = IntToText(NColumnsOfData)
    NRowsOfData             = 1
    OverwriteNColumnsOfData = .TRUE.
    PrintColumnNo           = .FALSE.
    FormatSpec              = '(A21,'//TRIM(Text)//'(2X,F12.2))'
    FPart(1)      = 'STREAM_HYDROGRAPHS'

    !Prepare header lines
    WorkArray(1) = ArrangeText('STREAM HYDROGRAPH',37)
    IF (iHydType .EQ. iHydFlow) THEN
        WorkArray(2) = ArrangeText('(UNIT=',cUnitFlow,')',37)
    ELSEIF (iHydType .EQ. iHydStage) THEN
        WorkArray(2) = ArrangeText('(UNIT=',cUnitElev,')',37)
    ELSE
        WorkArray(2) = ArrangeText('(UNIT=',TRIM(ADJUSTL(cUnitFlow))//' or '//TRIM(ADJUSTL(cUnitElev)),')',37)
    END IF
    CALL PrepareTitle(TitleLines(1),WorkArray(1:2),39,42)
    Header = ''
    WRITE (Header(1,1),'(A1,7X,A13)') '*','HYDROGRAPH ID'
    WRITE (Header(2,1),'(A1,15X,A5)') '*','NODES'
    WRITE (Header(3,1),'(A1,8X,A4)') '*','TIME'
    IF (iHydType .EQ. iHydBoth) THEN
        DO indx=1,NHyd
            WRITE (Header(1,indx+1),'(I12)') indx
            WRITE (Header(2,indx+1),'(I12)') iHydNodes(indx)
        END DO
        DO indx=NHyd+1,NColumnsOfData
            WRITE (Header(1,indx+1),'(I12)') indx
            WRITE (Header(2,indx+1),'(I12)') iHydNodes(indx-NHyd)
        END DO
    ELSE
        DO indx=1,NColumnsOfData
            WRITE (Header(1,indx+1),'(I12)') indx
            WRITE (Header(2,indx+1),'(I12)') iHydNodes(indx)
        END DO
    END IF
    HeaderFormat(1)='(A21,'//TRIM(Text)//'(2X,A12))'
    HeaderFormat(2)='(A21,'//TRIM(Text)//'(2X,A12))'
    HeaderFormat(3)='(A13,8X,'//TRIM(Text)//'(A))'

    !Prepare the time series output file
    CALL PrepareTSDOutputFile(HydFile                          , &
                              NColumnsOfData                   , &
                              NRowsOfData                      , &
                              OverwriteNColumnsOfData          , &
                              FormatSpec                       , &
                              TitleLines                       , &
                              Header                           , &
                              HeaderFormat                     , &
                              PrintColumnNo                    , &
                              DataUnit                         , &
                              DataType                         , &
                              CPart                            , &
                              FPart                            , &
                              TimeStep%Unit                    , &
                              StrmNodes=iHydNodes_Local        , &
                              iStat=iStat                      )

    !Clear memory
    DEALLOCATE (cPart , DataUnit , DataType , Header , iHydNodes_Local , STAT=ErrorCode)

  END SUBROUTINE PrepStrmHydFile

  ! -------------------------------------------------------------
  ! --- NEW STREAM HYDROGRAPH FILE AS INPUT FILE FOR POST-PROCESSING
  ! -------------------------------------------------------------
  SUBROUTINE PrepStrmHydFile_ForInquiry(TimeStep,cFileName,iHydType,NHyd,iHydNodes,InFile,iStat)
    TYPE(TimeStepType),INTENT(IN) :: TimeStep
    CHARACTER(LEN=*),INTENT(IN)   :: cFileName
    INTEGER,INTENT(IN)            :: iHydType,NHyd,iHydNodes(NHyd)
    TYPE(RealTSDataInFileType)    :: InFile
    INTEGER,INTENT(OUT)           :: iStat

    !Local variables
    INTEGER               :: indx,ErrorCode
    CHARACTER             :: BPart*32,EandFParts*64
    CHARACTER,ALLOCATABLE :: cPathNames*80(:),CPart*32(:)

    !Initialize
    iStat = 0

    !Pathname parts
    IF (iHydType .EQ. iHydBoth) THEN
        ALLOCATE (cPathNames(2*NHyd) , CPart(2*NHyd))
        cPathNames     = ''
        CPart          = ''
        CPart(1:NHyd)  = 'FLOW'
        CPart(NHyd+1:) = 'FLOW_DEPTH'
    ELSE
        ALLOCATE (cPathNames(NHyd) , CPart(NHyd))
        cPathNames = ''
        CPart      = ''
        IF (iHydType .EQ. iHydFlow) THEN
            CPart = 'FLOW'
        ELSE
            CPart ='FLOW_DEPTH'
        END IF
    END IF
    EandFParts = UpperCase(TRIM(TimeStep%Unit)) // '/STREAM_HYDROGRAPHS/'

    !Instantiate file according to its type
    IF (iGetFileType_FromName(cFileName) .EQ. DSS) THEN
        !Form pathnames
        IF (iHydType .EQ. iHydBoth) THEN
            DO indx=1,NHyd
                BPart            = 'R' // TRIM(IntToText(iHydNodes(indx)))
                cPathnames(indx) = '/IWFM/' // TRIM(BPart) // '/' // TRIM(CPart(indx)) // '//' // TRIM(EandFParts)
            END DO
            DO indx=NHyd+1,2*NHyd
                BPart            = 'R' // TRIM(IntToText(iHydNodes(indx-NHyd)))
                cPathnames(indx) = '/IWFM/' // TRIM(BPart) // '/' // TRIM(CPart(indx)) // '//' // TRIM(EandFParts)
            END DO
            CALL InFile%Init(cFileName,'streamflow hydrograph at selected locations',TimeStep%TrackTime,nCol=2*NHyd,cPathNames=cPathNames,iStat=iStat)
        ELSE
            DO indx=1,NHyd
                BPart            = 'R' // TRIM(IntToText(iHydNodes(indx)))
                cPathnames(indx) = '/IWFM/' // TRIM(BPart) // '/' // TRIM(CPart(indx)) // '//' // TRIM(EandFParts)
            END DO
            CALL InFile%Init(cFileName,'streamflow hydrograph at selected locations',TimeStep%TrackTime,nCol=NHyd,cPathNames=cPathNames,iStat=iStat)
        END IF
    ELSE
        IF (iHydType .EQ. iHydBoth) THEN
            CALL InFile%Init(cFileName,'streamflow hydrograph at selected locations',BlocksToSkip=0,nCol=2*NHyd,iStat=iStat)
        ELSE
            CALL InFile%Init(cFileName,'streamflow hydrograph at selected locations',BlocksToSkip=0,nCol=NHyd,iStat=iStat)
        END IF
    END IF

    !Clear memory
    DEALLOCATE (cPathNames , CPart , STAT=ErrorCode)

  END SUBROUTINE PrepStrmHydFile_ForInquiry

  ! -------------------------------------------------------------
  ! --- TRANSFER HYDROGRAPHS FROM TEXT/DSS FILE TO HDF FILE
  ! -------------------------------------------------------------
  SUBROUTINE Transfer_To_HDF(StrmHyd,NTIME,TimeStep,iStat)
    CLASS(StrmHydrographType)      :: StrmHyd
    INTEGER,INTENT(IN)             :: NTIME
    TYPE(TimeStepType),INTENT(IN)  :: TimeStep
    INTEGER,INTENT(OUT)            :: iStat

    !Local variables
    INTEGER                  :: NColumns(1),FileReadCode,indxTime,ErrorCode
    CHARACTER                :: cDataSetName(1)*50,cHDFFileName*500
    TYPE(GenericFileType)    :: OutFile
    CHARACTER(:),ALLOCATABLE :: cFileName
    REAL(8),ALLOCATABLE      :: rData(:,:),rConvFactor(:)
    TYPE(TimeStepType)       :: TimeStep_Local

    !Initialize
    iStat = 0

    !Return if no output
    IF (.NOT. StrmHyd%HydFile_Defined) RETURN

    !Get the name of the text/DSS file
    CALL StrmHyd%HydFile_ForInquiry%GetFileName(cFileName)

    !Name for the HDF file
    cHDFFileName = TRIM(ADJUSTL(StripTextUntilCharacter(cFileName,'.',Back=.TRUE.))) // '.hdf'

    !Open output file HDF file
    CALL OutFile%New(FileName=TRIM(cHDFFileName),InputFile=.FALSE.,IsTSFile=.TRUE.,iStat=iStat)
    IF (iStat .EQ. -1) RETURN

    !Create dataset
    IF (StrmHyd%iHydType .EQ. iHydBoth) THEN
        NColumns(1) = 2 * StrmHyd%NHyd
    ELSE
        NColumns(1) = StrmHyd%NHyd
    END IF
    cDataSetName(1) = '/Stream_Hydrographs'
    CALL OutFile%CreateHDFDataSet(cPathNames=cDataSetName,NColumns=NColumns,NTime=NTIME,TimeStep=TimeStep,DataType=0d0,iStat=iStat)
    IF (iStat .EQ. -1) RETURN

    !Write the type of hydrograph
    CALL OutFile%WriteData(iDataset,'/Stream_Hydrographs','HydrographType',ScalarAttrData=StrmHyd%iHydType)

    !Conversion factor used when printing out results
    IF (StrmHyd%iHydType .EQ. iHydBoth) THEN
        ALLOCATE (rConvFactor(2*StrmHyd%NHyd))
        rConvFactor(1:StrmHyd%NHyd)  = 1d0 / StrmHyd%FactFlow
        rConvFactor(StrmHyd%NHyd+1:) = 1d0 / StrmHyd%FactElev
    ELSE
        ALLOCATE (rConvFactor(StrmHyd%NHyd))
        IF (StrmHyd%iHydType .EQ. iHydFlow) THEN
            rConvFactor = 1d0 / StrmHyd%FactFlow
        ELSE
            rConvFactor = 1d0 / StrmHyd%FactElev
        END IF
    END IF

    !Advance time to t=1
    TimeStep_Local                    = TimeStep
    TimeStep_Local%CurrentTimeStep    = TimeStep_Local%CurrentTimeStep + 1
    TimeStep_Local%CurrentDateAndTime = IncrementTimeStamp(TimeStep_Local%CurrentDateAndTime,TimeStep_Local%DELTAT_InMinutes)

    !Transfer heads to HDF file
    IF (StrmHyd%iHydType .EQ. iHydBoth) THEN
        ALLOCATE (rData(2*StrmHyd%NHyd,1))
    ELSE
        ALLOCATE (rData(StrmHyd%NHyd,1))
    END IF
    DO indxTime=1,NTIME
        !Read data
        CALL ReadTSData(TimeStep_Local,'stream hydrographs',StrmHyd%HydFile_ForInquiry,FileReadCode,iStat)
        IF (iStat .EQ. -1) RETURN

        !Transfer  values to matrix to be written to HDF file
        rData(:,1) = StrmHyd%HydFile_ForInquiry%rValues * rConvFactor

        !Store heads into HDF file
        CALL OutFile%WriteData(rData)

        !Advance time
        TimeStep_Local%CurrentTimeStep    = TimeStep_Local%CurrentTimeStep + 1
        TimeStep_Local%CurrentDateAndTime = IncrementTimeStamp(TimeStep_Local%CurrentDateAndTime,TimeStep_Local%DELTAT_InMinutes)
    END DO

    !Rewind input file
    CALL StrmHyd%HydFile_ForInquiry%File%RewindFile_To_BeginningOfTSData(iStat)

    !Close HDF file
    CALL OutFile%Kill()

    !Clear memory
    DEALLOCATE (rData,rConvFactor,STAT=ErrorCode)

  END SUBROUTINE Transfer_To_HDF

END MODULE

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE


text:
      SUBROUTINE ENDREG
*
*
*       Transformation to physical variables.
*       -------------------------------------
*
      IMPLICIT REAL*8  (A-H,O-Z)
      REAL*8  M,MIJ,CM(3)
      LOGICAL  SWITCH,GTYPE,GTYPE0
      COMMON/CREG/  M(4),X(12),XD(12),P(12),Q(12),TIME4,ENERGY,EPSR2,
     &              XR(9),W(9),R(6),TA(6),MIJ(6),CMX(10),RMAX4,TMAX,
     &              DS,TSTEP,EPS,NSTEP4,NAME4(4),KZ15,KZ27,NREG,NFN
      COMMON/TPR/   SWITCH,GTYPE,GTYPE0
      COMMON/SAVEP/  PI(12)
      COMMON/ICONF/  I1,I2,I3,I4
*
*
      DO 1 L = 1,3
          L1 = 3*(L - 1) + 1
          L2 = L1 + 1
          L3 = L2 + 1
*
          K1 = 4*(L - 1) + 1
          K2 = K1 + 1
          K3 = K2 + 1
          K4 = K3 + 1
*
          XR(L1) = Q(K1)**2 - Q(K2)**2 - Q(K3)**2 + Q(K4)**2
          XR(L2) = 2.D0*(Q(K1)*Q(K2) - Q(K3)*Q(K4))
          XR(L3) = 2.D0*(Q(K1)*Q(K3) + Q(K2)*Q(K4))
          R(L) = Q(K1)**2 + Q(K2)**2 + Q(K3)**2 + Q(K4)**2
          A = 0.5D0/R(L)
          W(L1) = (Q(K1)*P(K1) - Q(K2)*P(K2) - Q(K3)*P(K3) +
     &                                         Q(K4)*P(K4))*A
          W(L2) = (Q(K2)*P(K1) + Q(K1)*P(K2) - Q(K4)*P(K3) -
     &                                         Q(K3)*P(K4))*A
          W(L3) = (Q(K3)*P(K1) + Q(K4)*P(K2) + Q(K1)*P(K3) +
     &                                         Q(K2)*P(K4))*A
    1 CONTINUE
*
      IP1 = I1*3 - 3
      IP2 = I2*3 - 3
      IP3 = I3*3 - 3
      IP4 = I4*3 - 3
*
      DO 2 K = 1,3
          PI(IP1+K) = -W(K  )
          PI(IP2+K) =  W(K  ) - W(K+3)
          PI(IP3+K) =  W(K+3) - W(K+6)
          PI(IP4+K) =  W(K+6)
          X(IP1+K) = -XR(K)
          X(IP2+K) = 0.0D0
*       Note terms X(IP1+K) = 0 & X(IP2+K) = XR(K) in early formulation.
          X(IP3+K) = X(IP2+K) + XR(K+3)
          X(IP4+K) = X(IP3+K) + XR(K+6)
    2 CONTINUE
*
*       Skip reduction to centre of mass after switching.
      IF (SWITCH) RETURN
*
*       Obtain velocities from momenta and form c.m. coordinates.
      CM(1) = 0.0D0
      CM(2) = 0.0D0
      CM(3) = 0.0D0
      SMASS = M(1) + M(2) + M(3) + M(4)
      DO 4 I = 1,4
          LI = 3*(I - 1)
          DO 3 K = 1,3
              XD(LI+K) = PI(LI+K)/M(I)
              CM(K) = CM(K) + M(I)*X(LI+K)/SMASS
    3     CONTINUE
    4 CONTINUE
*
*       Express coordinates in the c.m. frame.
      L = 0
      DO 6 I = 1,4
          DO 5 K = 1,3
              L = L + 1
              X(L) = X(L) - CM(K)
    5     CONTINUE
    6 CONTINUE
*
      RETURN
*
      END

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
! This file is part of flox.
! SPDX-Identifier: Apache-2.0

module flox_ast_printer
  use flox_ast, only : lox_visitor, lox_expr, lox_ast, &
    & lox_block, lox_expr_stmt, lox_print, lox_var, lox_if, lox_while, lox_fun, lox_return, &
    & lox_assign, lox_logical, lox_binary, lox_grouping, lox_literal, lox_unary, lox_call
  implicit none
  private

  public :: lox_ast_printer

  type, extends(lox_visitor) :: lox_ast_printer
    character(len=:), allocatable :: string
  contains
    procedure :: visit_ast
    procedure :: visit_assign
    procedure :: visit_logical
    procedure :: visit_binary
    procedure :: visit_grouping
    procedure :: visit_literal
    procedure :: visit_unary
    procedure :: visit_call
    procedure :: visit_block
    procedure :: visit_expr_stmt
    procedure :: visit_print
    procedure :: visit_var
    procedure :: visit_if
    procedure :: visit_while
    procedure :: visit_fun
    procedure :: visit_return
  end type lox_ast_printer

contains

  recursive subroutine visit_assign(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_assign), intent(in) :: expr

    call parenthesize(self, "= "//expr%name%val, expr%value)
  end subroutine visit_assign

  recursive subroutine visit_logical(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_logical), intent(in) :: expr

    call parenthesize(self, expr%operator%val, expr%left, expr%right)
  end subroutine visit_logical

  recursive subroutine visit_binary(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_binary), intent(in) :: expr

    call parenthesize(self, expr%operator%val, expr%left, expr%right)
  end subroutine visit_binary

  recursive subroutine visit_grouping(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_grouping), intent(in) :: expr

    call parenthesize(self, "group", expr%expression)
  end subroutine visit_grouping

  recursive subroutine visit_literal(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_literal), intent(in) :: expr

    self%string = self%string // expr%object%val
  end subroutine visit_literal

  recursive subroutine visit_unary(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_unary), intent(in) :: expr

    call parenthesize(self, expr%operator%val, expr%right)
  end subroutine visit_unary

  recursive subroutine visit_call(self, expr)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_call), intent(in) :: expr

    integer :: iarg

    self%string = self%string // "(call "
    call expr%callee%accept(self)
    do iarg = 1, expr%narg
      self%string = self%string // " "
      associate(arg => expr%args(iarg)%expr)
        call arg%accept(self)
      end associate
    end do
    self%string = self%string // ")"
  end subroutine visit_call

  recursive subroutine visit_ast(self, ast)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_ast), intent(in) :: ast

    integer :: istmt

    self%string = "(root"
    do istmt = 1, ast%nstmt
      self%string = self%string // " "
      associate(entry => ast%stmts(istmt)%stmt)
        call entry%accept(self)
      end associate
    end do
    self%string = self%string // ")"
  end subroutine visit_ast

  recursive subroutine visit_block(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_block), intent(in) :: stmt

    integer :: istmt

    self%string = self%string // "(block"
    do istmt = 1, stmt%nstmt
      self%string = self%string // " "
      associate(entry => stmt%stmts(istmt)%stmt)
        call entry%accept(self)
      end associate
    end do
    self%string = self%string // ")"
  end subroutine visit_block

  recursive subroutine visit_expr_stmt(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_expr_stmt), intent(in) :: stmt

    call parenthesize(self, "stmt", stmt%expression)
  end subroutine visit_expr_stmt

  recursive subroutine visit_print(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_print), intent(in) :: stmt

    call parenthesize(self, "print", stmt%expression)
  end subroutine visit_print

  recursive subroutine visit_var(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_var), intent(in) :: stmt

    if (allocated(stmt%expression)) then
      call parenthesize(self, "var "//stmt%name%val, stmt%expression)
    else
      self%string = self%string // "(var " // stmt%name%val // ")"
    end if
  end subroutine visit_var

  recursive subroutine visit_if(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_if), intent(in) :: stmt

    self%string = self%string // "(if "
    call parenthesize(self, "condition", stmt%condition)
    self%string = self%string // " "
    call stmt%then_branch%accept(self)
    if (allocated(stmt%else_branch)) then
      self%string = self%string // " "
      call stmt%else_branch%accept(self)
    end if
    self%string = self%string // ")"
  end subroutine visit_if

  recursive subroutine visit_while(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_while), intent(in) :: stmt

    self%string = self%string // "(while "
    call parenthesize(self, "condition", stmt%condition)
    self%string = self%string // " "
    call stmt%body%accept(self)
    self%string = self%string // ")"
  end subroutine visit_while

  recursive subroutine visit_fun(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_fun), intent(in) :: stmt

    integer :: ipar

    self% string = self%string // "(fun " // stmt%name%val // " "
    do ipar = 1, size(stmt%params)
      self%string = self%string // stmt%params(ipar)%val // " "
    end do
    call stmt%body%accept(self)
    self%string = self%string // ")"
  end subroutine visit_fun

  recursive subroutine visit_return(self, stmt)
    class(lox_ast_printer), intent(inout) :: self
    class(lox_return), intent(in) :: stmt

    if (allocated(stmt%expression)) then
      call parenthesize(self, "return", stmt%expression)
    else
      self%string = self%string // "(return)"
    end if
  end subroutine visit_return

  recursive subroutine parenthesize(self, name, expr1, expr2)
    class(lox_ast_printer), intent(inout) :: self
    character(len=*), intent(in) :: name
    class(lox_expr), intent(in) :: expr1
    class(lox_expr), intent(in), optional :: expr2

    self%string = self%string // "(" // name // " "
    call expr1%accept(self)
    if (present(expr2)) then
      self%string = self%string // " "
      call expr2%accept(self)
    end if
    self%string = self%string // ")"
  end subroutine parenthesize

end module flox_ast_printer

max_stars_repo_name: awvwgk/flox

max_issues_repo_name: awvwgk/flox

max_forks_repo_name: awvwgk/flox


text:
! { dg-do compile }
! { dg-options "-c -Wall" }
!
! PR fortran/40041
! cf. also PR fortran/20373

subroutine valid_one
  REAL :: a
  INTEGER :: n
  INTRINSIC   ABS, MAX
  a(n) = MAX(ABS(2),ABS(3),n)
end subroutine

subroutine valid_two
  IMPLICIT NONE
  REAL :: a
  INTEGER :: n
  INTRINSIC   ABS, MAX
  a(n) = MAX(ABS(2),ABS(3),n)
end subroutine

subroutine warnings_one
  REAL :: a
  INTEGER :: n
  REAL :: ABS ! { dg-warning "Type specified for intrinsic function" }
  REAL :: MAX ! { dg-warning "Type specified for intrinsic function" }
  INTRINSIC   ABS, MAX
  a(n) = MAX(ABS(2),ABS(3),n)
end subroutine

subroutine warnings_two
  IMPLICIT NONE
  REAL :: a
  INTEGER :: n
  INTRINSIC ABS ! { dg-warning "Type specified for intrinsic function" }
  INTRINSIC MAX ! { dg-warning "Type specified for intrinsic function" }
  REAL :: ABS
  REAL :: MAX
  a(n) = MAX(ABS(2),ABS(3),n)
end subroutine

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!RUN: %flang -gdwarf-4 -S -emit-llvm %s -o - | FileCheck %s

!CHECK-LABEL: distinct !DIGlobalVariable(name: "alc_arr"
!CHECK-SAME: type: [[TYPE:![0-9]+]]
!CHECK: [[TYPE]] = !DICompositeType(tag: DW_TAG_array_type,
!CHECK-SAME: elements: [[ELEMENTS:![0-9]+]], dataLocation: !DIExpression(DW_OP_push_object_address, DW_OP_deref), allocated: !DIExpression(DW_OP_push_object_address, DW_OP_deref)
!CHECK: [[ELEMENTS]] = !{[[ELEMENT:![0-9]+]]}
!CHECK: [[ELEMENT]] = !DISubrange(lowerBound: !DIExpression(DW_OP_push_object_address, DW_OP_plus_uconst, 96, DW_OP_deref), upperBound: !DIExpression(DW_OP_push_object_address, DW_OP_plus_uconst, 136, DW_OP_deref), stride: !DIExpression(DW_OP_push_object_address, DW_OP_plus_uconst, 128, DW_OP_deref, DW_OP_push_object_address, DW_OP_plus_uconst, 40, DW_OP_deref, DW_OP_mul))

module mod_allocatable_arr
  integer, allocatable :: alc_arr(:)
end module

program main
  use mod_allocatable_arr
  allocate (alc_arr(10))
  alc_arr = 99
  print *, alc_arr
end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
module fcl_kinetic_energy
#include <messenger.h>

    use mod_kinds,                  only: ik, rk
    use type_evaluator,             only: evaluator_t
    use type_chidg_worker,          only: chidg_worker_t
    use mod_constants,              only: HALF, ONE, TWO
    use type_integral_cache,        only: integral_cache_t
    use mod_functional_operators,   only: integrate_volume

    use DNAD_D

    implicit none

    !>  This functional computes the integral of the kientic energy over a volume
    !!
    !!  Ek = 1/2*u**2
    !!
    !!  @author Matteo Ugolotti
    !!  @date   12/17/2017
    !!
    !!
    !------------------------------------------------------------------------------------------
    type,  extends(evaluator_t), public   :: kinetic_energy_t

        ! Quantities needed for the computation
        ! Used to check completeness of information provided
        integer(ik)     :: min_ref_geom = 1
        integer(ik)     :: min_aux_geom = 0
        integer(ik)     :: max_ref_geom = 100
        integer(ik)     :: max_aux_geom = 0

    contains

        procedure, public   :: init
        procedure, public   :: check
        procedure, public   :: compute_functional
        procedure, public   :: finalize_functional

    end type kinetic_energy_t
    !******************************************************************************************

contains

    !>  Initialize the functional
    !!
    !!
    !!  @author Matteo Ugolotti
    !!  @date   05/11/2017
    !!
    !!
    !------------------------------------------------------------------------------------------
    subroutine init(self)
        class(kinetic_energy_t), intent(inout)   :: self

        ! Start defining the evaluator information
        call self%set_name("Kinetic Energy")
        call self%set_eval_type("Functional")
        call self%set_int_type("VOLUME INTEGRAL")

        call self%add_integral("kinetic energy")

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

    !>  This procedure checks that all the information have been provided by the user to fully
    !!  compute the functional
    !!
    !!  @author Matteo Ugolotti
    !!  @date   12/17/2017
    !!
    !!
    !------------------------------------------------------------------------------------------
    subroutine check(self)
        class(kinetic_energy_t), intent(inout)   :: self

        integer(ik)                     :: aux_geoms, ref_geoms
        character(len=:), allocatable   :: usr_msg_r, usr_msg_a
        logical                         :: ref_geoms_exceed, aux_geoms_exceed

        !
        ! Check that the functional has all the information needed
        !
        ref_geoms = self%n_ref_geom()
        aux_geoms = self%n_aux_geom()

        ref_geoms_exceed = (ref_geoms < self%min_ref_geom .or. ref_geoms > self%max_ref_geom)
        aux_geoms_exceed = (aux_geoms < self%min_aux_geom .or. aux_geoms > self%max_aux_geom)

        usr_msg_r = "fcl_kinetic_energy: wrong number of reference geometries. The minimum number of reference geoemtries required is 1"
        usr_msg_a = "fcl_kinetic_energy: This functional does not require any auxiliary geometry. Please, use `chidg edit` to remove the input auxiliary geometry"

        if (ref_geoms_exceed) call chidg_signal(FATAL, usr_msg_r)
        if (aux_geoms_exceed) call chidg_signal(FATAL, usr_msg_a)

    end subroutine check
    !******************************************************************************************

    !>  Compute Kinetic Energy Integral of a single element.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   12/17/2017
    !!
    !!  param[inout]       worker
    !!  param[inout]       cache        ! Storage for integrals computed in each element/face
    !!
    !------------------------------------------------------------------------------------------
    subroutine compute_functional(self,worker,cache)
        class(kinetic_energy_t),    intent(inout)   :: self
        type(chidg_worker_t),       intent(inout)   :: worker
        type(integral_cache_t),     intent(inout)   :: cache

        integer(ik)             :: idomain_g, ielement_g, iface, ierr

        type(AD_D), allocatable, dimension(:) :: density,                   &
                                                 mom_1, mom_2, mom_3,       &
                                                 u_1, u_2, u_3, norm_vel,   &
                                                 umag, Ek,                  &
                                                 norm_1, norm_2, norm_3,    &
                                                 unorm_1, unorm_2, unorm_3, &
                                                 da, r
        real(rk),   allocatable, dimension(:) :: weights
        type(AD_D)                            :: energy_integral

        ! Get primary fields
        !
        density     = worker%get_field('Density'   ,'value','element')
        mom_1       = worker%get_field('Momentum-1','value','element')
        mom_2       = worker%get_field('Momentum-2','value','element')
        mom_3       = worker%get_field('Momentum-3','value','element')

        !
        ! Account for cylindrical coordinates
        !
        r = worker%coordinate('1','face interior')
        if (worker%coordinate_system() == 'Cylindrical') then
            mom_2 = mom_2/r
        else if (worker%coordinate_system() == 'Cartesian') then

        else
            call chidg_signal(FATAL,'bad coordinate system')
        end if

        !
        ! Compute velocities
        !
        u_1 = mom_1/density
        u_2 = mom_2/density
        u_3 = mom_3/density
        umag = sqrt(u_1*u_1 + u_2*u_2 + u_3*u_3)

        !
        ! Compute Kinetic Energy
        !
        Ek = HALF * density * umag**TWO

        !
        ! Compute volume integral
        !
        energy_integral = integrate_volume(worker,Ek)

        !
        ! Store in cache
        !
        call cache%push_back(energy_integral,'kinetic energy')

    end subroutine compute_functional
    !******************************************************************************************

    !>  Finalize subroutine carries out the final computation before storing the functional
    !!  It comes after the communcation among processore and input argument cache has
    !!  the over all integrals (in this case only 'kinetic energy') over the reference
    !!  geometry.
    !!
    !!  @auhtor Matteo Ugolotti
    !!  @date   12/17/2017
    !!
    !!  param[inout]       cache     Storage for integrals, this contains the overall
    !!                               functional (ie after parallel communication).
    !!
    !---------------------------------------------------------------------------------------------
    function finalize_functional(self,worker,cache) result (res)
        class(kinetic_energy_t),    intent(inout)   :: self
        type(chidg_worker_t),       intent(in)      :: worker
        type(integral_cache_t),     intent(inout)   :: cache

        type(AD_D)      :: res

        res = cache%get_value(worker%mesh,'kinetic energy')

    end function finalize_functional
    !*********************************************************************************************

end module fcl_kinetic_energy

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


      program CO2system_test
! **********************************************************************
! *                                                                    *
! *   Test program of CO2system.f                                      *
! *                                                                    *
! **********************************************************************
!
      USE mod_geochem

      implicit none

      real(8) T,S
      real(8) pH, d11BT, d11B_BOH4
      real(8) :: DOsatu, ssO2flux, U10, DO
!
      integer i,j
!      integer istep
!
      open(52,file='./output/debug2.txt')!!!!!!!!!!!!!!!!!!!for debug
      open(51,file='./output/result2.txt')

      T=28.0d0+273.15
      S=34.5
      U10=1.0d0

      DO=300.0d0

!      do i=1,100
!        pH = 7.0d0+0.03*dble(i)
!        d11BT= 39.5d0
!        d11B_BOH4= d11B_BOH4_frompHd11BT(pH,d11BT, T, S)
!        write(51,*) pH, d11B_BOH4
!      enddo

!      write(*,*) densSW(25.0d0,35.0d0),densW(25.0d0)

      DOsatu=O2satu(T,S)
      ssO2flux = Flux_O2(DO, DOsatu, U10, T, S )  ! sea to air is positive

      write(*,*) T, DOsatu, ssO2flux

      end program CO2system_test

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


text:
../aes/rtl/AesControl.v
../aes/rtl/AesDecipher.v
../aes/rtl/AesEncipher.v
../aes/rtl/AesInvMixColumns.v
../aes/rtl/AesInvShiftRows.v
../aes/rtl/AesInvSubBytes.v
../aes/rtl/AesKeyExpansion.v
../aes/rtl/AesMixColumns.v
../aes/rtl/AesShiftRows.v
../aes/rtl/AesSubBytes.v
../aes/rtl/Aes.v
../aes/rtl/SPSR.v

./rtl/axis_interface.sv
./rtl/krnl_aes_axi_ctrl_slave.v
./rtl/aes_wrapper.sv
./rtl/krnl_aes.sv

./ip_generation/clk_gen/clk_gen_clk_wiz.v
./ip_generation/clk_gen/clk_gen.v

./ip_generation/axi_clock_converter/sim/axi_clock_converter.v
./ip_generation/axi_clock_converter/hdl/axi_clock_converter_v2_1_vl_rfs.v

./ip_generation/axis_clock_converter/sim/axis_clock_converter.v
./ip_generation/axis_clock_converter/hdl/axis_clock_converter_v1_1_vl_rfs.v
./ip_generation/axis_clock_converter/hdl/axis_infrastructure_v1_1_vl_rfs.v

./ip_generation/axi_vip_mst/sim/axi_vip_mst_pkg.sv
./ip_generation/axi_vip_mst/sim/axi_vip_mst.sv
./ip_generation/axi_vip_mst/hdl/axi_vip_v1_1_vl_rfs.sv
./ip_generation/axi_vip_mst/hdl/axi_infrastructure_v1_1_vl_rfs.v

./ip_generation/axis_vip_mst/sim/axis_vip_mst_pkg.sv
./ip_generation/axis_vip_mst/sim/axis_vip_mst.sv
./ip_generation/axis_vip_slv/sim/axis_vip_slv_pkg.sv
./ip_generation/axis_vip_slv/sim/axis_vip_slv.sv
./ip_generation/axis_vip_mst/hdl/axi4stream_vip_v1_1_vl_rfs.sv

./tbench/tb_krnl_aes.sv

$XILINX_VIVADO/data/verilog/src/glbl.v

max_stars_repo_name: jlamperez/Vitis-Tutorials

max_issues_repo_name: jlamperez/Vitis-Tutorials

max_forks_repo_name: jlamperez/Vitis-Tutorials


text:
! This file is part of mctc-lib.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

module test_read_xyz
   use mctc_env_testing, only : new_unittest, unittest_type, error_type, check
   use mctc_io_read_xyz
   use mctc_io_structure
   implicit none
   private

   public :: collect_read_xyz

contains

!> Collect all exported unit tests
subroutine collect_read_xyz(testsuite)

   !> Collection of tests
   type(unittest_type), allocatable, intent(out) :: testsuite(:)

   testsuite = [ &
      & new_unittest("valid1-xyz", test_valid1_xyz), &
      & new_unittest("valid2-xyz", test_valid2_xyz), &
      & new_unittest("valid3-xyz", test_valid3_xyz), &
      & new_unittest("valid4-xyz", test_valid4_xyz), &
      & new_unittest("valid5-xyz", test_valid5_xyz), &
      & new_unittest("invalid1-xyz", test_invalid1_xyz, should_fail=.true.), &
      & new_unittest("invalid2-xyz", test_invalid2_xyz, should_fail=.true.), &
      & new_unittest("invalid3-xyz", test_invalid3_xyz, should_fail=.true.), &
      & new_unittest("invalid4-xyz", test_invalid4_xyz, should_fail=.true.), &
      & new_unittest("invalid5-xyz", test_invalid5_xyz, should_fail=.true.), &
      & new_unittest("invalid6-xyz", test_invalid6_xyz, should_fail=.true.) &
      & ]

end subroutine collect_read_xyz

subroutine test_valid1_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "9", &
      "WATER27, (H2O)3", &
      "O     1.1847029    1.1150792   -0.0344641 ", &
      "H     0.4939088    0.9563767    0.6340089 ", &
      "H     2.0242676    1.0811246    0.4301417 ", &
      "O    -1.1469443    0.0697649    1.1470196 ", &
      "H    -1.2798308   -0.5232169    1.8902833 ", &
      "H    -1.0641398   -0.4956693    0.3569250 ", &
      "O    -0.1633508   -1.0289346   -1.2401808 ", &
      "H     0.4914771   -0.3248733   -1.0784838 ", &
      "H    -0.5400907   -0.8496512   -2.1052499 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

   call check(error, allocated(struc%comment), "Comment line should be preserved")
   if (allocated(error)) return
   call check(error, struc%comment, "WATER27, (H2O)3")
   if (allocated(error)) return
   call check(error, struc%nat, 9, "Number of atoms does not match")
   if (allocated(error)) return
   call check(error, struc%nid, 2, "Number of species does not match")
   if (allocated(error)) return

end subroutine test_valid1_xyz

subroutine test_valid2_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "24", &
      "", &
      "C          1.07317        0.04885       -0.07573", &
      "N          2.51365        0.01256       -0.07580", &
      "C*         3.35199        1.09592       -0.07533", &
      "N          4.61898        0.73028       -0.07549", &
      "C*         4.57907       -0.63144       -0.07531", &
      "C          3.30131       -1.10256       -0.07524", &
      "C          2.98068       -2.48687       -0.07377", &
      "18O        1.82530       -2.90038       -0.07577", &
      "N          4.11440       -3.30433       -0.06936", &
      "C*         5.45174       -2.85618       -0.07235", &
      "O          6.38934       -3.65965       -0.07232", &
      "N          5.66240       -1.47682       -0.07487", &
      "C          7.00947       -0.93648       -0.07524", &
      "C          3.92063       -4.74093       -0.06158", &
      "D          0.73398        1.08786       -0.07503", &
      "D          0.71239       -0.45698        0.82335", &
      "D          0.71240       -0.45580       -0.97549", &
      "H          2.99301        2.11762       -0.07478", &
      "H          7.76531       -1.72634       -0.07591", &
      "H          7.14864       -0.32182        0.81969", &
      "H          7.14802       -0.32076       -0.96953", &
      "H          2.86501       -5.02316       -0.05833", &
      "H          4.40233       -5.15920        0.82837", &
      "H          4.40017       -5.16929       -0.94780"
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

   call check(error, .not.allocated(struc%comment), "Empty comment line should not be saved")
   if (allocated(error)) return
   call check(error, struc%nat, 24, "Number of atoms does not match")
   if (allocated(error)) return
   call check(error, struc%nid, 7, "Number of species does not match")
   if (allocated(error)) return

end subroutine test_valid2_xyz

subroutine test_valid3_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "24", &
      "#", &
      "c  1.07317  0.04885 -0.07573 -0.05445590", &
      "n  2.51365  0.01256 -0.07580 -0.00457526", &
      "c  3.35199  1.09592 -0.07533  0.08391889", &
      "n  4.61898  0.73028 -0.07549 -0.27870751", &
      "c  4.57907 -0.63144 -0.07531  0.11914924", &
      "c  3.30131 -1.10256 -0.07524 -0.02621044", &
      "c  2.98068 -2.48687 -0.07377  0.26115960", &
      "o  1.82530 -2.90038 -0.07577 -0.44071824", &
      "n  4.11440 -3.30433 -0.06936 -0.10804747", &
      "c  5.45174 -2.85618 -0.07235  0.30411699", &
      "o  6.38934 -3.65965 -0.07232 -0.44083760", &
      "n  5.66240 -1.47682 -0.07487 -0.07457706", &
      "c  7.00947 -0.93648 -0.07524 -0.04790859", &
      "c  3.92063 -4.74093 -0.06158 -0.03738239", &
      "h  0.73398  1.08786 -0.07503  0.06457802", &
      "h  0.71239 -0.45698  0.82335  0.08293905", &
      "h  0.71240 -0.45580 -0.97549  0.08296802", &
      "h  2.99301  2.11762 -0.07478  0.05698136", &
      "h  7.76531 -1.72634 -0.07591  0.09025556", &
      "h  7.14864 -0.32182  0.81969  0.07152988", &
      "h  7.14802 -0.32076 -0.96953  0.07159003", &
      "h  2.86501 -5.02316 -0.05833  0.08590674", &
      "h  4.40233 -5.15920  0.82837  0.06906357", &
      "h  4.40017 -5.16929 -0.94780  0.06926350"
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

   call check(error, struc%nat, 24, "Number of atoms does not match")
   if (allocated(error)) return
   call check(error, struc%nid, 4, "Number of species does not match")
   if (allocated(error)) return

end subroutine test_valid3_xyz

subroutine test_valid4_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "3", &
      "WATER27, H2O", &
      "O     1.1847029    1.1150792   -0.0344641 ", &
      "H     0.4939088    0.9563767    0.6340089 ", &
      "H     2.0242676    1.0811246    0.4301417 ", &
      "3", &
      "WATER27, H2O", &
      "O    -1.1469443    0.0697649    1.1470196 ", &
      "H    -1.2798308   -0.5232169    1.8902833 ", &
      "H    -1.0641398   -0.4956693    0.3569250 ", &
      "3", &
      "WATER27, H2O", &
      "O    -0.1633508   -1.0289346   -1.2401808 ", &
      "H     0.4914771   -0.3248733   -1.0784838 ", &
      "H    -0.5400907   -0.8496512   -2.1052499 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   if (.not.allocated(error)) then
      call read_xyz(struc, unit, error)
   end if
   close(unit)
   if (allocated(error)) return

   call check(error, struc%nat, 3, "Number of atoms does not match")
   if (allocated(error)) return
   call check(error, struc%nid, 2, "Number of species does not match")
   if (allocated(error)) return

end subroutine test_valid4_xyz

subroutine test_valid5_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "3", &
      "WATER27, H2O", &
      "8     1.1847029    1.1150792   -0.0344641 ", &
      "1     0.4939088    0.9563767    0.6340089 ", &
      "1     2.0242676    1.0811246    0.4301417 ", &
      "3", &
      "WATER27, H2O", &
      "8    -1.1469443    0.0697649    1.1470196 ", &
      "1    -1.2798308   -0.5232169    1.8902833 ", &
      "1    -1.0641398   -0.4956693    0.3569250 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   if (.not.allocated(error)) then
      call read_xyz(struc, unit, error)
   end if
   close(unit)
   if (allocated(error)) return

   call check(error, struc%nat, 3, "Number of atoms does not match")
   if (allocated(error)) return
   call check(error, struc%nid, 2, "Number of species does not match")
   if (allocated(error)) return

end subroutine test_valid5_xyz

subroutine test_invalid1_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "9", &
      "WATER27, (H2O)3", &
      "O     1.1847029    1.1150792   -0.0344641 ", &
      "H     0.4939088    0.9563767    0.6340089 ", &
      "H     2.0242676    1.0811246    0.4301417 ", &
      "O    -1.1469443    0.0697649    1.1470196 ", &
      "H    -1.2798308   -0.5232169    1.8902833 ", &
      "H    -1.0641398   -0.4956693    0.3569250 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

end subroutine test_invalid1_xyz

subroutine test_invalid2_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "", &
      "WATER27, (H2O)3", &
      "O     1.1847029    1.1150792   -0.0344641 ", &
      "H     0.4939088    0.9563767    0.6340089 ", &
      "H     2.0242676    1.0811246    0.4301417 ", &
      "O    -1.1469443    0.0697649    1.1470196 ", &
      "H    -1.2798308   -0.5232169    1.8902833 ", &
      "H    -1.0641398   -0.4956693    0.3569250 ", &
      "O    -0.1633508   -1.0289346   -1.2401808 ", &
      "H     0.4914771   -0.3248733   -1.0784838 ", &
      "H    -0.5400907   -0.8496512   -2.1052499 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

end subroutine test_invalid2_xyz

subroutine test_invalid3_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      & "120"
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

end subroutine test_invalid3_xyz

subroutine test_invalid4_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "-3", &
      "H2O", &
      "O     1.1847029    1.1150792   -0.0344641 ", &
      "H     0.4939088    0.9563767    0.6340089 ", &
      "H     2.0242676    1.0811246    0.4301417 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

end subroutine test_invalid4_xyz

subroutine test_invalid5_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "3", &
      "H2O", &
      "***   1.1847029    1.1150792   -0.0344641 ", &
      "***   0.4939088    0.9563767    0.6340089 ", &
      "***   2.0242676    1.0811246    0.4301417 "
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

end subroutine test_invalid5_xyz

subroutine test_invalid6_xyz(error)

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   type(structure_type) :: struc
   integer :: unit

   open(status='scratch', newunit=unit)
   write(unit, '(a)') &
      "3", &
      "H2O", &
      "1.1847029    1.1150792   -0.0344641      O", &
      "0.4939088    0.9563767    0.6340089      H", &
      "2.0242676    1.0811246    0.4301417      H"
   rewind(unit)

   call read_xyz(struc, unit, error)
   close(unit)
   if (allocated(error)) return

end subroutine test_invalid6_xyz

end module test_read_xyz

max_stars_repo_name: loriab/mctc-lib

max_issues_repo_name: loriab/mctc-lib

max_forks_repo_name: loriab/mctc-lib


text:
C
C
*+ cln_setdef

       subroutine cln_setdef(full,status)
C      ----------------------------------
C
C reset the CLEAN defaults
C Given:
C   logical flag indicating whether a full initialisation is required
       logical       full
C Returned:
C    error status
       integer       status
C
C Initialise the control information for CLEAN.  If full is set then
C the results information is initialised as well.
*-
       include '../include/clean_record.inc'
       include '/mrao/include/maplib_redtape.inc'

       integer  n

       if (status.ne.0) return

C reset values for parameters
       lnotbx = .false.
       do n=1,4
         box_window(n) = 0
       end do
       if (full) then
         ltrbm = .false.
         lprhat = .false.
         prhat = 0.0
         lclean = .true.
         lrest = .true.
         fllim = 0.001
         fract = 0.2
         itlim = 100
         lsort = .false.
         sort_gate  = 0.0
         sort_depth = 100
         sort_inner = 100
         number_cc  = 0
         times_cleaned = 0
         number_it_done = 0
         max_on_resid = 0.0
         how_ended = 0
       end if

       do n=1,4
         box_window(n) = 0
       end do

       if (map_defined) then

         call ldredt(redt_map,status)
         call iuv_load(clean_window,status)
         number_search = 1
         do n=1,9
           call iuv_load(search_window(1,n),status)
         end do

         call ldredt(redt_beam,status)
         beam_pa = 0.0
         if (hpfbwu.gt.0.0) then
             beam_size_u = hpfbwu
             beam_size_v = hpfbwv
         else
           if (freq.gt.100.0 .and. freq.lt.40000.0) then
             beam_size_u = 2.0*5000.0/freq
           else
             beam_size_u = 1.0
           end if
           beam_size_v = beam_size_u
           beam_size_v = beam_size_v/sin(decmap)
         end if
         call iuv_load(beam_window,status)

       end if

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
program main
    use iso_c_binding
    use slam_pg_connector

    implicit none

!    character, parameter :: conninfo(*) = "host=localhost port=5432 dbname=tracklets sslmode=disable"
    TYPE(C_PTR)         :: conn
    TYPE(C_PTR)         :: res
    INTEGER(C_INT)      :: fields
    INTEGER(C_INT)      :: tuples
    INTEGER(C_INT)      :: i,j
    character(kind=c_char, len=:), allocatable :: c_tmp
    !character(len=100) :: c_tmp

    ! Make a connection to the database
    conn = connect_db("host=localhost port=5432 dbname=tracklets sslmode=disable"//C_NULL_CHAR)
    !conn = connect_db("host=localhost port=5432 dbname=tracklets sslmode=disable")

    ! Check to see that the backend connection was successfully made
    if (connStatus(conn) .eq. CONN_OK) then
        write(*,*) "Connection to database established."
    else
        write (*,*) "Connection to database failed: ", get_db_error_message(conn)
        call finish(conn)
    end if

    res = exec(conn, "select * from objects;"//C_NULL_CHAR)
    !res = exec(conn, "select * from objects;")
    if (resultStatus(res) .eq. RES_TUPLES_OK) then
        ! first, print out the table collumn attribute names
        fields = nfields(res);
        do i = 0, fields -1
            c_tmp = fname(res, i)
            write (*,'(A15)', advance='no') c_tmp
            deallocate(c_tmp)
        end do
        write (*,'(A)',advance='yes') ""

        ! next, print out the rows of data
        tuples = ntuples(res)
        do i = 0, tuples -1
            do j = 0, fields -1
                c_tmp = get_value(res, i, j)
                write (*,'(A15)',advance='no') get_value(res, i, j)
                deallocate(c_tmp)
            end do
            write (*,'(A)',advance='yes') ""
        end do
    else
        write (*,*) "select * failed: ", get_db_error_message(conn)
        call clear(res)
        call finish(conn)
    end if

    call clear(res)
    ! close the connection to the database and cleanup
    call finish(conn)

end program main

max_stars_repo_name: Space-Systems/libslam

max_issues_repo_name: Space-Systems/libslam

max_forks_repo_name: Space-Systems/libslam


text:
           program demo_any
           implicit none
           logical l
              l = any([.true., .true., .true.])
              print *, l
              call section
              contains
                subroutine section
                integer a(2,3), b(2,3)
                  a = 1
                  b = 1
                  b(2,2) = 2
                  print *, any(a .eq. b, 1)
                  print *, any(a .eq. b, 2)
                end subroutine section
           end program demo_any

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages


text:
*DECK CV
      REAL FUNCTION CV (XVAL, NDATA, NCONST, NORD, NBKPT, BKPT, W)
C***BEGIN PROLOGUE  CV
C***PURPOSE  Evaluate the variance function of the curve obtained
C            by the constrained B-spline fitting subprogram FC.
C***LIBRARY   SLATEC
C***CATEGORY  L7A3
C***TYPE      SINGLE PRECISION (CV-S, DCV-D)
C***KEYWORDS  ANALYSIS OF COVARIANCE, B-SPLINE,
C             CONSTRAINED LEAST SQUARES, CURVE FITTING
C***AUTHOR  Hanson, R. J., (SNLA)
C***DESCRIPTION
C
C     CV( ) is a companion function subprogram for FC( ).  The
C     documentation for FC( ) has complete usage instructions.
C
C     CV( ) is used to evaluate the variance function of the curve
C     obtained by the constrained B-spline fitting subprogram, FC( ).
C     The variance function defines the square of the probable error
C     of the fitted curve at any point, XVAL.  One can use the square
C     root of this variance function to determine a probable error band
C     around the fitted curve.
C
C     CV( ) is used after a call to FC( ).  MODE, an input variable to
C     FC( ), is used to indicate if the variance function is desired.
C     In order to use CV( ), MODE must equal 2 or 4 on input to FC( ).
C     MODE is also used as an output flag from FC( ).  Check to make
C     sure that MODE = 0 after calling FC( ), indicating a successful
C     constrained curve fit.  The array SDDATA, as input to FC( ), must
C     also be defined with the standard deviation or uncertainty of the
C     Y values to use CV( ).
C
C     To evaluate the variance function after calling FC( ) as stated
C     above, use CV( ) as shown here
C
C          VAR=CV(XVAL,NDATA,NCONST,NORD,NBKPT,BKPT,W)
C
C     The variance function is given by
C
C          VAR=(transpose of B(XVAL))*C*B(XVAL)/MAX(NDATA-N,1)
C
C     where N = NBKPT - NORD.
C
C     The vector B(XVAL) is the B-spline basis function values at
C     X=XVAL.  The covariance matrix, C, of the solution coefficients
C     accounts only for the least squares equations and the explicitly
C     stated equality constraints.  This fact must be considered when
C     interpreting the variance function from a data fitting problem
C     that has inequality constraints on the fitted curve.
C
C     All the variables in the calling sequence for CV( ) are used in
C     FC( ) except the variable XVAL.  Do not change the values of these
C     variables between the call to FC( ) and the use of CV( ).
C
C     The following is a brief description of the variables
C
C     XVAL    The point where the variance is desired.
C
C     NDATA   The number of discrete (X,Y) pairs for which FC( )
C             calculated a piece-wise polynomial curve.
C
C     NCONST  The number of conditions that constrained the B-spline in
C             FC( ).
C
C     NORD    The order of the B-spline used in FC( ).
C             The value of NORD must satisfy 1 < NORD < 20 .
C
C             (The order of the spline is one more than the degree of
C             the piece-wise polynomial defined on each interval.  This
C             is consistent with the B-spline package convention.  For
C             example, NORD=4 when we are using piece-wise cubics.)
C
C     NBKPT   The number of knots in the array BKPT(*).
C             The value of NBKPT must satisfy NBKPT .GE. 2*NORD.
C
C     BKPT(*) The real array of knots.  Normally the problem data
C             interval will be included between the limits BKPT(NORD)
C             and BKPT(NBKPT-NORD+1).  The additional end knots
C             BKPT(I),I=1,...,NORD-1 and I=NBKPT-NORD+2,...,NBKPT, are
C             required by FC( ) to compute the functions used to fit
C             the data.
C
C     W(*)    Real work array as used in FC( ).  See FC( ) for the
C             required length of W(*).  The contents of W(*) must not
C             be modified by the user if the variance function is
C             desired.
C
C***REFERENCES  R. J. Hanson, Constrained least squares curve fitting
C                 to discrete data using B-splines, a users guide,
C                 Report SAND78-1291, Sandia Laboratories, December
C                 1978.
C***ROUTINES CALLED  BSPLVN, SDOT
C***REVISION HISTORY  (YYMMDD)
C   780801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CV
      DIMENSION BKPT(NBKPT), W(*), V(40)
C***FIRST EXECUTABLE STATEMENT  CV
      ZERO = 0.
      MDG = NBKPT - NORD + 3
      MDW = NBKPT - NORD + 1 + NCONST
      IS = MDG*(NORD+1) + 2*MAX(NDATA,NBKPT) + NBKPT + NORD**2
      LAST = NBKPT - NORD + 1
      ILEFT = NORD
   10 IF (.NOT.(XVAL.GE.BKPT(ILEFT+1) .AND. ILEFT.LT.LAST-1)) GO TO 20
      ILEFT = ILEFT + 1
      GO TO 10
   20 CALL BSPLVN(BKPT, NORD, 1, XVAL, ILEFT, V(NORD+1))
      ILEFT = ILEFT - NORD + 1
      IP = MDW*(ILEFT-1) + ILEFT + IS
      N = NBKPT - NORD
      DO 30 I=1,NORD
         V(I) = SDOT(NORD,W(IP),1,V(NORD+1),1)
         IP = IP + MDW
   30 CONTINUE
      CV = MAX(SDOT(NORD,V,1,V(NORD+1),1),ZERO)
C
C     SCALE THE VARIANCE SO IT IS AN UNBIASED ESTIMATE.
      CV = CV/MAX(NDATA-N,1)
      RETURN
      END


text:
c program bestscale

      include 'VICMAIN_FOR'
      subroutine main44

      parameter (nlmax=1100,nsmax=1100)
      real*4 inbuf(nsmax,nlmax),obuf(nsmax,nlmax)
      real*4 scale(2),refbuf(nsmax,nlmax)
      real*4 cor(100),sc(100)
      integer*4 count,def,status,unit,ounit
      real*8 r,sumx,sumy,sumx2,sumy2,sumxy

      call xvpcnt('INP',nin)
      call xvpcnt('OUT',nout)

c parameters
      call xveaction('SA',' ')
      call xvpcnt('INP',nin)
      call xvpcnt('OUT',nout)
      if(nin.ne.nout)then
        write(*,*)'inputs must match outputs in number'
        call abend
      endif
      call xvparm('SCALE',scale,count,def,2)

c open input1
      call xvunit(unit,'INP',1,status,' ')
      call xvopen(unit,status,'U_FORMAT','REAL',' ')
      call xvget(unit,status,'NL',nl,'NS',ns,' ')
      if(ns.gt.nsmax)then
          call xvmessage('Line length too long',' ')
          call abend
      endif
      if(nl.gt.nlmax)then
          call xvmessage('Line column too long',' ')
          call abend
      endif

c open output1
      call xvunit(ounit,'OUT',1,status,' ')
      call xvopen(ounit,status,'U_FORMAT','REAL',
     +              'OP','WRITE',' ')

c copy first input to first output
      do j=1,nl
        call xvread(unit,refbuf(1,j),status,'LINE',j,' ')
        call xvwrit(ounit,refbuf(1,j),status,' ')
      enddo
      call xvclose(unit,status,'CLOS_ACT','FREE',' ')
      call xvclose(ounit,status,'CLOS_ACT','FREE',' ')

      sinc=(scale(2)-scale(1))/50. ! can go to 100
      cx=ns/2.0
      cy=nl/2.0
      inc=5
      mode=0 ! bilinear
      write(*,*)'Image   scale   correlation'
      write(*,*)1,1.0,1.0
      inp=1
10    inp=inp+1

c read next image
      call xvselpi(inp) ! use this input label for the output.
      call xvunit(unit,'INP',inp,status,' ')
      call xvopen(unit,status,'U_FORMAT','REAL',' ')
      do j=1,nl
        call xvread(unit,inbuf(1,j),status,'LINE',j,' ')
      enddo

c compute correlations between each image and refimage
c as a function of scale.
      k=0
      do s=scale(1),scale(2),sinc
        k=k+1
        n=0
        sumx=0.d0
        sumx2=0.d0
        sumy=0.d0
        sumy2=0.d0
        sumxy=0.d0
        do j=1,nl,inc
          y=(j-cy)*s+cy
          do i=1,ns,inc
            x=(i-cx)*s+cx
            call interpolate(inbuf,nlmax,nsmax,nl,ns,mode,
     +                         x,y,dn,ind)
            if(ind.eq.0)then
              n=n+1
              sumx=sumx+refbuf(i,j)
              sumx2=sumx2+refbuf(i,j)*refbuf(i,j)
              sumy=sumy+dn
              sumy2=sumy2+dn*dn
              sumxy=sumxy+refbuf(i,j)*dn
            endif
          enddo
        enddo
        r=(sumxy-sumx*sumy/n)**2/
     +    ((sumx2-sumx**2/n)*(sumy2-sumy**2/n))
        cor(k)=dsqrt(r)
        sc(k)=s
c        write(*,*)k,sc(k),cor(k)
      enddo

c obtain the best estimate for the scale
      best=-1.0
      do i=1,k
        if(cor(i).gt.best)then
          best=cor(i)
          n=i
        endif
      enddo
      if(n.eq.1)then
        s=sc(1)
      else if(n.eq.k)then
        s=sc(k)
      else
        s=sc(n)-0.5*
     +    ((sc(n)-sc(n-1))**2*(cor(n)-cor(n+1))-(sc(n)-sc(n+1))**2*
     +    (cor(n)-cor(n-1)))  /
     +    ((sc(n)-sc(n-1))*(cor(n)-cor(n+1))-(sc(n)-sc(n+1))*
     +    (cor(n)-cor(n-1)))
      endif
      write(*,*)inp,s,cor(n)

c generate image at best scale
      do j=1,nl
        y=(j-cy)*s+cy
        do i=1,ns
          x=(i-cx)*s+cx
          call interpolate(inbuf,nlmax,nsmax,nl,ns,mode,
     +                         x,y,dn,ind)
          obuf(i,j)=dn
        enddo
      enddo

c write next output
      call xvunit(ounit,'OUT',inp,status,' ')
      call xvopen(ounit,status,'U_FORMAT','REAL',
     +              'OP','WRITE',' ')
      do j=1,nl
        call xvwrit(ounit,obuf(1,j),status,' ')
      enddo
      call xvclose(unit,status,'CLOS_ACT','FREE',' ')
      call xvclose(ounit,status,'CLOS_ACT','FREE',' ')

      if(inp.lt.nin)goto 10

      return
      end

      subroutine interpolate(inbuf,nlbuf,nsbuf,nl,ns,mode,
     +  x,y,dn,ind)
c inbuf of dimension nsbuf by nlbuf, input , real, ns by nl data.
c mode is the interpolation mode.
c mode=0 bilinear interpolation, 4 neighbors.
c mode=1 nearest neighbor.
c mode=2 2nd order, 6 neighbors.
c mode=3 sampling theorem, 15 by 15 convolution.
c mode=4 bicubic spline
c x,y are real coordinates of input point, input, real.
c dn is interpolated intensity, returned, real.
c ind=0 if dn computed ok. ind=1 if x,y out of bounds.
c    returned, integer.
c If ind=0 on return, dn is set to zero.

      parameter (kerw=9)
      real*4 inbuf(nsbuf,nlbuf),tw,bw,lw,rw
      real*4 kerx(-kerw:kerw),kery(-kerw:kerw)
      real*4 ya(4),y1(4),y2(4),y12(4),c(4,4),cl(16),xa(16)
      real*4 wt(16,16)
      DATA WT/1.,0.,-3.,2.,4*0.,-3.,0.,9.,-6.,2.,0.,-6.,
     *  4.,8*0.,3.,0.,-9.,6.,-2.,0.,6.,-4.,10*0.,9.,-6.,
     *  2*0.,-6.,4.,2*0.,3.,-2.,6*0.,-9.,6.,2*0.,6.,-4.,
     *  4*0.,1.,0.,-3.,2.,-2.,0.,6.,-4.,1.,0.,-3.,2.,8*0.,
     *  -1.,0.,3.,-2.,1.,0.,-3.,2.,10*0.,-3.,2.,2*0.,3.,
     *  -2.,6*0.,3.,-2.,2*0.,-6.,4.,2*0.,3.,-2.,0.,1.,-2.,
     *  1.,5*0.,-3.,6.,-3.,0.,2.,-4.,2.,9*0.,3.,-6.,3.,0.,
     *  -2.,4.,-2.,10*0.,-3.,3.,2*0.,2.,-2.,2*0.,-1.,1.,
     *  6*0.,3.,-3.,2*0.,-2.,2.,5*0.,1.,-2.,1.,0.,-2.,4.,
     *  -2.,0.,1.,-2.,1.,9*0.,-1.,2.,-1.,0.,1.,-2.,1.,10*0.,
     *  1.,-1.,2*0.,-1.,1.,6*0.,-1.,1.,2*0.,2.,-2.,2*0.,-1.,1./

      ind=1
      dn=0.0

      if(mode.eq.0)then ! bilinear
        i=x
        j=y
        if(i.lt.1)return
        if(j.lt.1)return
        if(i.ge.ns)return
        if(j.ge.nl)return
        ind=0
        lw=x-i
        rw=1.0-lw
        tw=y-j
        bw=1.0-tw
        top=lw*inbuf(i+1,j)+rw*inbuf(i,j)
        bot=lw*inbuf(i+1,j+1)+rw*inbuf(i,j+1)
        dn=bw*top+tw*bot

      else if(mode.eq.1)then ! nearest neighbor
        i=nint(x)
        j=nint(y)
        if(i.lt.1)return
        if(j.lt.1)return
        if(i.gt.ns)return
        if(j.gt.nl)return
        ind=0
        dn=inbuf(i,j)

      else if(mode.eq.2)then ! second order
        i=nint(x)
        j=nint(y)
        if(i.lt.2)return
        if(j.lt.2)return
        if(i.ge.ns)return
        if(j.ge.nl)return
        ind=0
        p=x-i
        q=y-j
        dn=q*(q-1.0)*0.5*inbuf(i,j-1)+
     +     p*(p-1.0)*0.5*inbuf(i-1,j)+
     +    (1.0+p*q-p*p-q*q)*inbuf(i,j)+
     +    p*(p-2.0*q+1.0)*0.5*inbuf(i+1,j)+
     +    q*(q-2.0*p+1.0)*0.5*inbuf(i,j+1)+
     +    p*q*inbuf(i+1,j+1)

      else if(mode.eq.3)then ! Sampling Theorem
        ii=nint(x)
        jj=nint(y)
        if(ii-kerw.lt.1)return
        if(jj-kerw.lt.1)return
        if(ii+kerw.gt.ns)return
        if(jj+kerw.gt.nl)return
        pi=3.14159
        ind=0
        xx=x-ii
        do i=-kerw,kerw
          d=pi*(i-xx)
          if(d.ne.0.0)then
            kerx(i)=sin(d)/d
          else
            kerx(i)=1.0
          endif
        enddo
        yy=y-jj
        do j=-kerw,kerw
          d=pi*(j-yy)
          if(d.ne.0.0)then
            kery(j)=sin(d)/d
          else
            kery(j)=1.0
          endif
        enddo
        dn=0.0
        do j=-kerw,kerw
          do i=-kerw,kerw
            dn=dn+kery(j)*kerx(i)*inbuf(ii+i,jj+j)
          enddo
        enddo

      else if(mode.eq.4)then ! bicubic spline
        i=x
        j=y
        if(i.lt.2)return
        if(j.lt.2)return
        if(i.gt.ns-2)return
        if(j.gt.nl-2)return
        ind=0
        ya(1)=inbuf(i,j)
        ya(2)=inbuf(i+1,j)
        ya(3)=inbuf(i+1,j+1)
        ya(4)=inbuf(i,j+1)
        y1(1)=(inbuf(i+1,j)-inbuf(i-1,j))/2.0
        y1(2)=(inbuf(i+2,j)-inbuf(i,j))/2.0
        y1(3)=(inbuf(i+2,j+1)-inbuf(i,j+1))/2.0
        y1(4)=(inbuf(i+1,j+1)-inbuf(i-1,j+1))/2.0
        y2(1)=(inbuf(i,j+1)-inbuf(i,j-1))/2.0
        y2(2)=(inbuf(i+1,j+1)-inbuf(i+1,j-1))/2.0
        y2(3)=(inbuf(i+1,j+2)-inbuf(i+1,j))/2.0
        y2(4)=(inbuf(i,j+2)-inbuf(i,j))/2.0
        y12(1)=(inbuf(i+1,j+1)-inbuf(i+1,j-1)-
     +          inbuf(i-1,j+1)+inbuf(i-1,j-1))/4.0
        y12(2)=(inbuf(i+2,j+1)-inbuf(i+2,j-1)-
     +          inbuf(i,j+1)+inbuf(i,j-1))/4.0
        y12(3)=(inbuf(i+2,j+2)-inbuf(i+2,j)-
     +          inbuf(i,j+2)+inbuf(i,j))/4.0
        y12(4)=(inbuf(i+1,j+2)-inbuf(i+1,j)-
     +          inbuf(i-1,j+2)+inbuf(i-1,j))/4.0
        ii=i
        jj=j
        DO 11 I=1,4
          Xa(I)=Ya(I)
          Xa(I+4)=Y1(I)
          Xa(I+8)=Y2(I)
          Xa(I+12)=Y12(I)
11      CONTINUE
        DO 13 I=1,16
          XX=0.
          DO 12 K=1,16
            XX=XX+WT(I,K)*Xa(K)
12        CONTINUE
          CL(I)=XX
13      CONTINUE
        L=0
        DO 15 I=1,4
          DO 14 J=1,4
            L=L+1
            C(I,J)=CL(L)
14        CONTINUE
15      CONTINUE
        T=x-ii
        U=y-jj
        dn=0.
        DO 16 I=4,1,-1
          dn=T*dn+((C(I,4)*U+C(I,3))*U+C(I,2))*U+C(I,1)
16      CONTINUE

      else
        write(*,*)'INTERPOLATE: invalid mode'
        stop
      endif
      return
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


module FarmClass
    use, intrinsic :: iso_fortran_env
    use SoilClass
    use SoybeanClass
    implicit none

    type :: Farm_
        type(Soybean_),allocatable :: Soybean(:,:)
        type(Soil_) ::  Soil
        integer(int32) :: num_of_ridge
        integer(int32) :: num_of_plant_per_ridge
        real(real64)   :: width_of_ridge
        real(real64)   :: width_of_plant_per_ridge
        real(real64)   :: length_of_farm, width_of_farm
        real(real64)   :: soil_depth
        real(real64)   :: seed_depth
        integer(int32) :: total_num_of_plant
        real(real64) :: plant_density
        real(real64) :: plant_density_m
        real(real64) :: total_weight_of_seed
        real(real64) :: total_area
        real(real64) :: g_per_100seed
        real(real64) :: locale(2)

        real(real64) :: Water_kg
    contains
        procedure :: init => initFarm
        procedure :: sowing => initFarm
        procedure :: fertilize => fertilizeFarm
        procedure :: diagnosis => diagnosisFarm
        procedure :: grow   => growFarm
        procedure :: export => exportFarm
    end type

contains

! ############################################
subroutine initFarm(obj,crop_name,num_of_ridge, num_of_plant_per_ridge,width_of_ridge,&
    width_of_plant_per_ridge,length_of_farm, width_of_farm,soil_depth,seed_depth,g_per_100seed,&
    meter,single,Variety)

    class(Farm_),intent(inout)  :: obj
    integer(int32),optional,intent(in) :: num_of_ridge, num_of_plant_per_ridge
    real(real64),optional,intent(in) :: width_of_ridge, width_of_plant_per_ridge,soil_depth
    real(real64),optional,intent(in) :: length_of_farm, width_of_farm,seed_depth,g_per_100seed
    character(*),intent(in) :: crop_name,Variety
    logical,optional,intent(in)      :: meter,single
    integer(int32) :: i,j,k,l,n,m

    ! input default value

    obj%num_of_ridge                = 10      ! 10 ridge
    obj%num_of_plant_per_ridge      = 10      ! 10 plant per a ridge
    obj%width_of_ridge              = 10.0d0  ! 10.0 cm
    obj%width_of_plant_per_ridge    = 10.0d0  ! 10.0 cm
    obj%length_of_farm              = 100.0d0 ! 1.0 m
    obj%width_of_farm               = 100.0d0 ! 1.0 m

    obj%soil_depth               = input(default=-100.0d0, option=soil_depth)
    obj%seed_depth               = input(default=-3.0d0, option=seed_depth)
    obj%g_per_100seed           = input(default=30.0d0, option = g_per_100seed )

    if(present(single) )then
        if(single .eqv. .true.)then
            obj%num_of_ridge                = 1      ! 1 ridge
            obj%num_of_plant_per_ridge      = 1      ! 1 plant per a ridge
            obj%width_of_ridge              = 1.0d0  ! 1.0 cm
            obj%width_of_plant_per_ridge    = 1.0d0  ! 1.0 cm
            obj%length_of_farm              = 1.0d0 ! 1.0 m
            obj%width_of_farm               = 1.0d0 ! 1.0 m
        endif
    endif

    ! determine area default = 1 m^2 = (100 cm)^2
    if(present(length_of_farm) .and. present(width_of_farm)  )then
        obj%length_of_farm = length_of_farm
        obj%width_of_farm  = width_of_farm

        ! if num_of_ridge and num_of_plant_per_ridge are imported
        if(present(num_of_ridge) .and. present(num_of_plant_per_ridge) )then
            obj%num_of_ridge             = num_of_ridge
            obj%num_of_plant_per_ridge   = num_of_plant_per_ridge

            obj%width_of_ridge           = dble(  int( obj%width_of_farm    /dble(obj%num_of_ridge             ) ) )
            obj%width_of_plant_per_ridge = dble(  int( obj%length_of_farm   /dble(obj%num_of_plant_per_ridge   ) ) )
        elseif(present(width_of_ridge) .and. present(width_of_plant_per_ridge) )then
            obj%width_of_ridge           = width_of_ridge
            obj%width_of_plant_per_ridge = width_of_plant_per_ridge

            obj%num_of_ridge             = int(   dble(  obj%width_of_farm    /dble(obj%width_of_ridge          )   ) )
            obj%num_of_plant_per_ridge   = int(   dble(  obj%length_of_farm   /dble(obj%width_of_plant_per_ridge)   ) )

        else
            print *, "initFarm #1 >> please give enough information about sowing condition >> default value is set."
        endif

    else
        print *, "initFarm #2 >> please give enough information about sowing condition >> default value is set."
    endif

    ! compute total number of plant
    obj%total_area          = obj%length_of_farm * obj%width_of_farm
    obj%total_num_of_plant  = obj%num_of_plant_per_ridge*obj%num_of_ridge
    obj%plant_density       = dble(obj%total_num_of_plant)/dble(obj%total_area) ! plant per cm^2
    obj%plant_density_m     = obj%plant_density*100.0d0*100.0d0 ! plant per m^2
    obj%total_weight_of_seed= obj%g_per_100seed* dble(obj%total_num_of_plant)/100.0d0

    ! compute soil profile
    call obj%soil%init()

    if(allocated(obj%soybean) )then
        deallocate(obj%soybean)
    endif
    allocate(obj%soybean(obj%num_of_ridge, obj%num_of_plant_per_ridge) )

    if( trim(crop_name) == "soybean" .or. trim(crop_name) == "Soybean" )then
        do i=1,obj%num_of_ridge
            do j=1, obj%num_of_plant_per_ridge
                call obj%soybean(i,j)%sowing(x=obj%width_of_ridge*dble(i-1)  ,&
                    y=obj%width_of_plant_per_ridge*dble(j-1) ,z=obj%seed_depth,Variety=Variety )
            enddo
        enddo
    else
        print *, "Sorry, crop_name :: ",crop_name ,"is not implemented yet."
    endif

end subroutine
! ############################################

! ############################################
subroutine fertilizeFarm(obj,N_kg,P_kg,K_kg,Ca_kg,Mg_kg,S_kg,Fe_kg,&
    Mn_kg,B_kg,Zn_kg,Mo_kg,Cu_kg,Cl_kg)
    class(Farm_),intent(inout) :: obj
    ! ================
    real(real64),optional,intent(in) :: N_kg
    real(real64),optional,intent(in) :: P_kg
    real(real64),optional,intent(in) :: K_kg
    real(real64),optional,intent(in) :: Ca_kg
    real(real64),optional,intent(in) :: Mg_kg
    real(real64),optional,intent(in) :: S_kg
    ! ================
    real(real64),optional,intent(in) :: Fe_kg
    real(real64),optional,intent(in) :: Mn_kg
    real(real64),optional,intent(in) :: B_kg
    real(real64),optional,intent(in) :: Zn_kg
    real(real64),optional,intent(in) :: Mo_kg
    real(real64),optional,intent(in) :: Cu_kg
    real(real64),optional,intent(in) :: Cl_kg
    ! ================

    call obj%Soil%fertilize(N_kg=N_kg,P_kg=P_kg,K_kg=K_kg,Ca_kg=Ca_kg,Mg_kg=Mg_kg,S_kg=S_kg,Fe_kg=Fe_kg,&
    Mn_kg=Mn_kg,B_kg=B_kg,Zn_kg=Zn_kg,Mo_kg=Mo_kg,Cu_kg=Cu_kg,Cl_kg=Cl_kg)

end subroutine
! ############################################

! ############################################
subroutine exportFarm(obj,FileName,withSTL,withMesh,TimeStep)
    class(Farm_),intent(inout)::obj
    integer(int32) :: obj_id,plant_id,i,j,tstep
    integer(int32),optional,intent(in) :: TimeStep
    character(*),intent(in) :: FileName
    logical,optional,intent(in) :: withSTL,withMesh
    character(200) :: id

    ! Visualize soybean-field
    tstep=input(default=0,option=TimeStep)
    ! initialize
    obj_id = 1 + tstep
    plant_id  = 0

    ! export soybeans
    do i=1,obj%num_of_ridge
        do j=1,obj%num_of_plant_per_ridge
            plant_id = plant_id + 1
            id=trim(  adjustl(fstring( plant_id ) ))
            if( present(withSTL) )then
                if(withSTL .eqv. .true. )then
                    call obj%soybean(i,j)%export(FileName=FileName//trim(id)//".geo",SeedID=obj_id,withSTL=withSTL)
                endif
            endif
            if( present(withMesh) )then
                if(withMesh .eqv. .true. )then
                    call obj%soybean(i,j)%export(FileName=FileName//trim(id)//".geo",SeedID=obj_id,withMesh=withMesh)
                endif
            endif
            call obj%soybean(i,j)%export(FileName=FileName//trim(id)//".geo",SeedID=obj_id )
        enddo
    enddo
    print *, "Total "//trim(id)//" plants are exported."

    ! export soil
    call obj%soil%export(FileName=FileName//trim(id)//"soil",format=".geo",objID=obj_id)

end subroutine
! ############################################

! ########################################
subroutine diagnosisFarm(obj,FileName)
    class(Farm_),intent(inout) :: obj
    character(*),optional,intent(in)::FileName

    call obj%Soil%diagnosis(FileName=FileName)

end subroutine
! ########################################

! ########################################
subroutine growFarm(obj,dt,temp,crop_name)
    class(Farm_),intent(inout) :: obj
    character(*),intent(in) :: crop_name
    real(real64),intent(in) :: dt,temp
    integer(int32) :: i,j

    if(crop_name=="Soybean" .or. crop_name=="soybean")then
        do i=1,size(obj%Soybean,1)
            do j=1,size(obj%Soybean,2)
                call obj%Soybean(i,j)%grow(dt=dt,Temp=temp)
            enddo
        enddo
    endif

end subroutine
! ########################################

end module


text:
      subroutine qlI4sub2mo(xpi,musq,Ival)
      implicit none
C---  Calculates the box with two opposite internal masses

C     Uses the ordering for the routine xpi wanted by FF
C     psq(1) lies between msq(1) and msq(2) and so on
C     xpi(1-4) = msq(1),msq(2),msq(3),msq(4)
C     xpi(5-8) = psq(1),psq(2),psq(3),psq(4)
C     xpi(9-10) = s12,s23
C     xpi(11) = +xpi(5)+xpi(6)+xpi(7)+xpi(8)-xpi(9)-xpi(10)
C     xpi(12) = -xpi(5)+xpi(6)-xpi(7)+xpi(8)+xpi(9)+xpi(10)
C     xpi(13) = +xpi(5)-xpi(6)+xpi(7)-xpi(8)+xpi(9)+xpi(10)

      include 'qlconstants.f'
      double precision xpi(13),musq,Y(4,4),Yalt(4,4)
      double complex Ival(-2:0)
      logical qlnonzero,qlzero
      integer ier

      If (qlnonzero(xpi(1)) .or. qlnonzero(xpi(3))) then
      write(6,*)'Error in qlI4sub2mo,qlnonzero m1sq/m3sq',xpi(1),xpi(3)
      stop
      endif

      call qlYcalc(xpi,Y,Yalt)
C----case 14
      if    ((qlzero(Y(1,1)))
     . .and. (qlzero(Y(3,3)))
     . .and. (qlzero(Y(1,2)))
     . .and. (qlzero(Y(1,4)))
     . .and. (qlzero(Y(2,3)))
     . .and. (qlzero(Y(3,4)))) then
      call qlbox14(Y,musq,Ival)

C----case 15
      elseif((qlzero(Y(1,1)))
     . .and. (qlzero(Y(3,3)))
     . .and. (qlzero(Y(1,2)))
     . .and. (qlzero(Y(1,4))))then
      call qlbox15(Y,musq,Ival)

      elseif((qlzero(Y(1,1)))
     . .and. (qlzero(Y(3,3)))
     . .and. (qlzero(Y(2,3)))
     . .and. (qlzero(Y(3,4))))then
      call qlbox15(Yalt,musq,Ival)

      else
      Ival(-2)=czip
      Ival(-1)=czip
      xpi(11) = +xpi(5)+xpi(6)+xpi(7)+xpi(8)-xpi(9)-xpi(10)
      xpi(12) = -xpi(5)+xpi(6)-xpi(7)+xpi(8)+xpi(9)+xpi(10)
      xpi(13) = +xpi(5)-xpi(6)+xpi(7)-xpi(8)+xpi(9)+xpi(10)
      call qlI4fin(Ival(0),xpi,ier)
      endif

      return
      end