text:
!**********************************************************************************************************************************
! LICENSING
! Copyright (C) 2015-2016  National Renewable Energy Laboratory
!
!    This file is part of Ambient Wind and Array Effects model for FAST.Farm.
!
! 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.
!
!**********************************************************************************************************************************
! File last committed: $Date$
! (File) Revision #: $Rev$
! URL: $HeadURL$
!**********************************************************************************************************************************
MODULE AWAE_IO

   use NWTC_Library
   use AWAE_Types

   implicit none

   type(ProgDesc), parameter  :: AWAE_Ver = ProgDesc( 'AWAE', '', '' )
   character(*),   parameter  :: AWAE_Nickname = 'AWAE'

   public :: AWAE_IO_InitGridInfo
   public :: ReadLowResWindFile

   contains

subroutine HiResWindCheck(n, nt, nX, nY, nZ, dX, dY, dZ, X0, Y0, Z0, dims, gridSpacing, origin, callingRoutine, errMsg, errStat)
   integer(IntKi),             intent(in   ) :: n               !< high-resolution time step number (0-based)
   integer(IntKi),             intent(in   ) :: nt              !< turbine number
   integer(IntKi),             intent(in   ) :: nX              !< number of grid points in the X-direction for turbine 1 at high-res time step 0
   integer(IntKi),             intent(in   ) :: nY              !< number of grid points in the Y-direction for turbine 1 at high-res time step 0
   integer(IntKi),             intent(in   ) :: nZ              !< number of grid points in the Z-direction for turbine 1 at high-res time step 0
   real(ReKi),                 intent(in   ) :: dX              !< space between grid points in the X-direction for turbine 1 at high-res time step 0
   real(ReKi),                 intent(in   ) :: dY              !< space between grid points in the Y-direction for turbine 1 at high-res time step 0
   real(ReKi),                 intent(in   ) :: dZ              !< space between grid points in the Z-direction for turbine 1 at high-res time step 0
   real(ReKi),                 intent(in   ) :: X0              !< starting X-location of the grid for turbine 1 at high-res time step 0 (m)
   real(ReKi),                 intent(in   ) :: Y0              !< starting Y-location of the grid for turbine 1 at high-res time step 0 (m)
   real(ReKi),                 intent(in   ) :: Z0              !< starting Z-location of the grid for turbine 1 at high-res time step 0 (m)
   integer(IntKi),             intent(in   ) :: dims(3)         !< dimensions of the grid for turbine nt at high-res time step n (m)
   real(ReKi),                 intent(in   ) :: gridSpacing(3)  !< spacing between grid points for turbine nt at high-res time step n (m)
   real(ReKi),                 intent(in   ) :: origin(3)       !< starting coordinates of the grid for turbine nt at high-res time step n (m)
   character(*),               intent(in   ) :: callingRoutine  !< string containing the name of the calling routine.
   integer(IntKi),             intent(  out) :: errStat         !< Error status of the operation
   character(*),               intent(  out) :: errMsg          !< Error message if errStat /= ErrID_None

      ! grid must have two points in each direction
   if ( (dims(1) < 2) .or. (dims(2) < 2) .or. (dims(3) < 2) ) then
      call SetErrStat ( ErrID_Fatal, 'The high resolution grid dimensions must contain a minimum of 2 nodes in each spatial direction. Turbine #'//trim(num2lstr(nt))//', time step '//trim(num2lstr(n)), errStat, errMsg, callingRoutine )
      return
   end if

      ! All turbines and all time steps must have the same grid dimensions due to array allocation assumptions
   if ( ( dims(1) .ne. nX ) .or. ( dims(2) .ne. nY ) .or. ( dims(3) .ne. nZ ) ) then
      call SetErrStat ( ErrID_Fatal, 'The high resolution grid dimensions for turbine #'//trim(num2lstr(nt))//' and high-res time step '//trim(num2lstr(n))//' do not match turbine #1 and time step 0.', errStat, errMsg, callingRoutine )
      return
   end if

      ! spacing must be consistent for a given turbine across all time steps
   if ( ( gridSpacing(1) .ne. dX ) .or. ( gridSpacing(2) .ne. dY ) .or. ( gridSpacing(3) .ne. dZ ) ) then
      call SetErrStat ( ErrID_Fatal, 'The high resolution grid spacing for turbine #'//trim(num2lstr(nt))//' and high-res time step '//trim(num2lstr(n))//' do not match time step 0.', errStat, errMsg, callingRoutine )
      return
   end if

      ! verify origin of any given turbine is not changing with time step.
   if ( ( origin(1) .ne. X0 ) .or. ( origin(2) .ne. Y0 ) .or. ( origin(3) .ne. Z0 ) ) then
      call SetErrStat ( ErrID_Fatal, 'The high resolution grid origin for turbine #'//trim(num2lstr(nt))//' and high-res time step '//trim(num2lstr(n))//' do not match time step 0.', errStat, errMsg, callingRoutine )
      return
   end if

end subroutine HiResWindCheck

subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg )
   integer(IntKi),             intent(in   ) :: n            !<  Low-resolution time step increment
   integer(IntKi),             intent(in   ) :: WrDisSkp1    !<  Number of low resolution time step increments per one output increment
   type(AWAE_ParameterType),   intent(in   ) :: p            !< Parameters
   type(AWAE_OutputType),   intent(in   ) :: y            !< Outputs
   type(AWAE_MiscVarType),     intent(inout) :: m            !< Misc/optimization variables
   integer(IntKi),             intent(  out) :: errStat      !< Error status of the operation
   character(*),               intent(  out) :: errMsg       !< Error message if errStat /= ErrID_None

   CHARACTER(1024)                         :: FileName
   INTEGER(IntKi)                          :: Un                   ! unit number of opened file
   INTEGER(IntKi)                          :: ErrStat2                        ! Temporary Error status
   CHARACTER(ErrMsgLen)                    :: ErrMsg2                         ! Temporary Error message
   CHARACTER(*),   PARAMETER               :: RoutineName = 'WriteDisWindFiles'
   INTEGER(IntKi)                          :: nt, n_out
   REAL(ReKi)                              :: t_out
   character(p%VTK_tWidth)                 :: Tstr  ! string for current VTK write-out step (padded with zeros)

   n_out = n/WrDisSkp1
   t_out = n*p%DT_low

   ! TimeStamp
   write(Tstr, '(i' // trim(Num2LStr(p%VTK_tWidth)) //'.'// trim(Num2LStr(p%VTK_tWidth)) // ')') n_out ! TODO use n instead..

   FileName = trim(p%OutFileVTKRoot)//".Low.Dis."//trim(Tstr)//".vtk"
   call WrVTK_SP_header( FileName, "Low resolution disturbed wind for time = "//trim(num2lstr(t_out))//" seconds.", Un, errStat2, errMsg2 )
      call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName)
      if (ErrStat >= AbortErrLev) return
   call WrVTK_SP_vectors3D( Un, "Velocity", (/p%nX_low,p%nY_low,p%nZ_low/), (/p%X0_low,p%Y0_low,p%Z0_low/), (/p%dX_low,p%dY_low,p%dZ_low/), m%Vdist_low, errStat2, errMsg2 )
      call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName)
      if (ErrStat >= AbortErrLev) return

   do nt= 1,p%NumTurbines
       ! We are only writing out the first of the high res data for a given low res time step

      FileName = trim(p%OutFileVTKRoot)//".HighT"//trim(num2lstr(nt))//".Dis."//trim(Tstr)//".vtk"
      call WrVTK_SP_header( FileName, "High resolution disturbed wind for time = "//trim(num2lstr(t_out))//" seconds.", Un, errStat2, errMsg2 )
         call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName)
         if (ErrStat >= AbortErrLev) return
      call WrVTK_SP_vectors3D( Un, "Velocity", (/p%nX_high,p%nY_high,p%nZ_high/), (/p%X0_high(nt),p%Y0_high(nt),p%Z0_high(nt)/), (/p%dX_high(nt),p%dY_high(nt),p%dZ_high(nt)/), y%Vdist_high(nt)%data(:,:,:,:,0), errStat2, errMsg2 )
         call SetErrStat(ErrStat2, errMsg2, ErrStat, ErrMsg, RoutineName)
         if (ErrStat >= AbortErrLev) return

   end do

end subroutine WriteDisWindFiles

!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine
!!
subroutine ReadLowResWindFile(n, p, Vamb_Low, errStat, errMsg)
   integer(IntKi),                 intent(in   )  :: n            !< Current simulation timestep increment (zero-based)
   type(AWAE_ParameterType),       intent(in   )  :: p            !< Parameters
   real(SiKi),                     intent(inout)  :: Vamb_Low(:,0:,0:,0:)         !< Array which will contain the low resolution grid of ambient wind velocities
   integer(IntKi),                 intent(  out)  :: errStat      !< Error status of the operation
   character(*),                   intent(  out)  :: errMsg       !< Error message if errStat /= ErrID_None

   integer(IntKi)           :: dims(3)              !  dimension of the 3D grid (nX,nY,nZ)
   real(ReKi)               :: origin(3)            !  the lower-left corner of the 3D grid (X0,Y0,Z0)
   real(ReKi)               :: gridSpacing(3)       !  spacing between grid points in each of the 3 directions (dX,dY,dZ)
   integer(IntKi)           :: Un                   !  unit number of opened file
   character(1024)          :: FileName             ! Name of output file
   character(1024)          :: descr                ! Line describing the contents of the file
   character(1024)          :: vecLabel             ! descriptor of the vector data

   errStat = ErrID_None
   errMsg  = ""

   FileName = trim(p%WindFilePath)//trim(PathSep)//"Low"//trim(PathSep)//"Amb.t"//trim(num2lstr(n))//".vtk"
   Un = 0; ! Initialization different from -1, important to prevent file closing
   call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, ErrStat, ErrMsg )
      if (ErrStat >= AbortErrLev) return
   call ReadVTK_SP_vectors( FileName, Un, dims, Vamb_Low, ErrStat, ErrMsg )

!==============================================================================

end subroutine ReadLowResWindFile

!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine
!!
subroutine ReadHighResWindFile(nt, n, p, Vamb_high, errStat, errMsg)

   integer(IntKi),                 intent(in   )  :: nt
   integer(IntKi),                 intent(in   )  :: n                       !< high-res time increment
   type(AWAE_ParameterType),       intent(in   )  :: p            !< Parameters
   real(SiKi),                     intent(inout)  :: Vamb_high(:,0:,0:,0:)         !< Array which will contain the low resolution grid of ambient wind velocities
   integer(IntKi),                 intent(  out)  :: errStat      !< Error status of the operation
   character(*),                   intent(  out)  :: errMsg       !< Error message if errStat /= ErrID_None

   integer(IntKi)           :: dims(3)              !  dimension of the 3D grid (nX,nY,nZ)
   real(ReKi)               :: origin(3)            !  the lower-left corner of the 3D grid (X0,Y0,Z0)
   real(ReKi)               :: gridSpacing(3)       !  spacing between grid points in each of the 3 directions (dX,dY,dZ)
   integer(IntKi)           :: Un                   !  unit number of opened file
   character(1024)          :: FileName             ! Name of output file
   character(1024)          :: descr                ! Line describing the contents of the file
   character(1024)          :: vecLabel             ! descriptor of the vector data

   errStat = ErrID_None
   errMsg  = ""

   FileName = trim(p%WindFilePath)//trim(PathSep)//"HighT"//trim(num2lstr(nt))//trim(PathSep)//"Amb.t"//trim(num2lstr(n))//".vtk"
   Un = 0; ! Initialization different from -1, important to prevent file closing
   call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, ErrStat, ErrMsg )
      if (ErrStat >= AbortErrLev) return
   call ReadVTK_SP_vectors( FileName, Un, dims, Vamb_high, ErrStat, ErrMsg )

!==============================================================================

end subroutine ReadHighResWindFile

subroutine AWAE_IO_InitGridInfo(InitInp, p, InitOut, errStat, errMsg)

   type(AWAE_InitInputType),    intent(in   ) :: InitInp        !< Input data for initialization routine
   type(AWAE_ParameterType),    intent(inout) :: p              !< Parameters
   type(AWAE_InitOutputType),   intent(  out) :: InitOut        !< Output for initialization routine
   integer(IntKi),              intent(  out) :: errStat
   character(*),                intent(  out) :: errMsg

   integer(IntKi)                             :: NumGrid_high
   integer(IntKi)                             :: errStat2      ! temporary error status of the operation
   character(ErrMsgLen)                       :: errMsg2       ! temporary error message
   character(*), parameter                    :: RoutineName = 'AWAE_IO_InitGridInfo'
   integer(IntKi)                             :: nXYZ_low, nx_low, ny_low, nz_low, nXYZ_high, nx_high, ny_high, nz_high
   integer(IntKi)                             :: dims(3)              ! dimension of the 3D grid (nX,nY,nZ)
   real(ReKi)                                 :: origin(3)            ! the lower-left corner of the 3D grid (X0,Y0,Z0)
   real(ReKi)                                 :: gridSpacing(3)       ! spacing between grid points in each of the 3 directions (dX,dY,dZ)
   character(1024)                            :: FileName             ! Name of output file
   character(1024)                            :: descr                ! Line describing the contents of the file
   character(1024)                            :: vecLabel             ! descriptor of the vector data
   integer(IntKi)                             :: Un                   ! file unit
   integer(IntKi)                             :: n, nt, nh, n_high_low, nhigh

   errStat = ErrID_None
   errMsg  = ""

   !---------------------------------------------------------------------------
   ! Parse time 0.0, low res wind input file to gather the grid
   !    information and set data associated with the low res grid
   !---------------------------------------------------------------------------

   if ( p%Mod_AmbWind == 1 ) then
      FileName = trim(p%WindFilePath)//trim(PathSep)//"Low"//trim(PathSep)//"Amb.t0.vtk"
      Un = -1 ! Set to force closing of file on return
      call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, ErrStat, ErrMsg )
         if (ErrStat >= AbortErrLev) return

      if ( (dims(1) < 2) .or. (dims(2) < 2) .or. (dims(3) < 2) ) then
         call SetErrStat ( ErrID_Fatal, 'The low resolution grid dimensions most contain a minimum of 2 nodes in each spatial direction.', errStat, errMsg, RoutineName )
         return
      end if
   else
      ! Using InflowWind, so data has been passed in via the InitInp data structure
      origin(1) = InitInp%InputFileData%X0_low
      origin(2) = InitInp%InputFileData%Y0_low
      origin(3) = InitInp%InputFileData%Z0_low
      dims(1)   = InitInp%InputFileData%nX_low
      dims(2)   = InitInp%InputFileData%nY_low
      dims(3)   = InitInp%InputFileData%nZ_low
      gridSpacing(1) = InitInp%InputFileData%dX_low
      gridSpacing(2) = InitInp%InputFileData%dY_low
      gridSpacing(3) = InitInp%InputFileData%dZ_low

   end if

   if ( (gridSpacing(1) <= 0.0_ReKi) .or. (gridSpacing(2) <= 0.0_ReKi) .or. (gridSpacing(3) <= 0.0_ReKi) ) &
      call SetErrStat ( ErrID_Fatal, 'The low resolution spatial resolution for Turbine 1 must be greater than zero in each spatial direction. ', errStat, errMsg, RoutineName )
      if (ErrStat >= AbortErrLev) return

   p%X0_low           = origin(1)
   p%Y0_low           = origin(2)
   p%Z0_low           = origin(3)
   p%nX_low           = dims(1)
   p%nY_low           = dims(2)
   p%nZ_low           = dims(3)
   p%dX_low           = gridSpacing(1)
   p%dY_low           = gridSpacing(2)
   p%dZ_low           = gridSpacing(3)

   InitOut%X0_Low     = origin(1)
   InitOut%Y0_low     = origin(2)
   InitOut%Z0_low     = origin(3)
   InitOut%nX_Low     = dims(1)
   InitOut%nY_low     = dims(2)
   InitOut%nZ_low     = dims(3)
   InitOut%dX_low     = gridSpacing(1)
   InitOut%dY_low     = gridSpacing(2)
   InitOut%dZ_low     = gridSpacing(3)

   p%NumGrid_low      = p%nX_Low*p%nY_Low*p%nZ_Low

   p%dXYZ_Low = gridSpacing
   p%dpol = (gridSpacing(1)+gridSpacing(2)+gridSpacing(3))/3.0_ReKi
   p%n_rp_max = ceiling(pi*((p%C_Meander*(p%r(p%NumRadii-1)+p%dpol))/p%dpol)**2.0_ReKi)
      ! Grid runs from (X0_low, Y0_low, Z0_low) to (X0_low + (p%nX_Low-1)*dX_low, Y0_low+ (p%nY_Low-1)*dY_low, Z0_low+ (p%nZ_Low-1)*dZ_low)
      ! (0,0,0) to (180,180,180)

   allocate( p%Grid_low(3,p%NumGrid_low),stat=errStat2)
      if (errStat2 /= 0) then
         call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for Grid_low.', errStat, errMsg, RoutineName )
         return
      end if

   nXYZ_low = 0
   do nz_low=0, p%nZ_low-1
      do ny_low=0, p%nY_low-1
         do nx_low=0, p%nX_low-1
            nXYZ_low = nXYZ_low + 1
            p%Grid_low(1,nXYZ_low) = origin(1) + nx_low*gridSpacing(1)
            p%Grid_low(2,nXYZ_low) = origin(2) + ny_low*gridSpacing(2)
            p%Grid_low(3,nXYZ_low) = origin(3) + nz_low*gridSpacing(3)
         end do
      end do
   end do

   if ( (InitInp%InputFileData%ChkWndFiles) .and. (p%Mod_AmbWind == 1) ) then
      do n=1,p%NumDT-1  ! We have already checked the first low res time step

         FileName = trim(p%WindFilePath)//trim(PathSep)//"Low"//trim(PathSep)//"Amb.t"//trim(num2lstr(n))//".vtk"
         Un = -1 ! Set to force closing of file on return
         call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, ErrStat, ErrMsg )
            if (ErrStat >= AbortErrLev) return

            ! verify dims, origin, gridSpacing match the first input file
         if ( ( dims(1) .ne. p%nX_low ) .or. ( dims(2) .ne. p%nY_low ) .or. ( dims(3) .ne. p%nZ_low ) ) then
            call SetErrStat ( ErrID_Fatal, 'The low resolution grid dimensions for time step '//trim(num2lstr(n))//' do not match time step 0.', errStat, errMsg, RoutineName )
            return
         end if
         if ( ( origin(1) .ne. InitOut%X0_Low ) .or. ( origin(2) .ne. InitOut%Y0_Low ) .or. ( origin(3) .ne. InitOut%Z0_Low ) ) then
            call SetErrStat ( ErrID_Fatal, 'The low resolution grid origins for time step '//trim(num2lstr(n))//' do not match time step 0.', errStat, errMsg, RoutineName )
            return
         end if
         if ( ( gridSpacing(1) .ne. p%dX_low ) .or. ( gridSpacing(2) .ne. p%dY_low ) .or. ( gridSpacing(3) .ne. p%dZ_low ) ) then
            call SetErrStat ( ErrID_Fatal, 'The low resolution grid spacing for time step '//trim(num2lstr(n))//' do not match time step 0.', errStat, errMsg, RoutineName )
            return
         end if

      end do
   end if

   allocate( InitOut%X0_high(p%NumTurbines), InitOut%Y0_high(p%NumTurbines), InitOut%Z0_high(p%NumTurbines), stat=errStat2)
      if (errStat2 /= 0) call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for InitOut origin arrays.', errStat, errMsg, RoutineName )
   allocate( InitOut%dX_high(p%NumTurbines), InitOut%dY_high(p%NumTurbines), InitOut%dZ_high(p%NumTurbines), stat=errStat2)
      if (errStat2 /= 0) call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for InitOut spatial increment arrays.', errStat, errMsg, RoutineName )
   allocate( p%X0_high(p%NumTurbines), p%Y0_high(p%NumTurbines), p%Z0_high(p%NumTurbines), stat=errStat2)
      if (errStat2 /= 0) call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for p origin arrays.', errStat, errMsg, RoutineName )
   allocate( p%dX_high(p%NumTurbines), p%dY_high(p%NumTurbines), p%dZ_high(p%NumTurbines), stat=errStat2)
      if (errStat2 /= 0) call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for p spatial increment arrays.', errStat, errMsg, RoutineName )
   if (ErrStat >= AbortErrLev) return

   if ( p%Mod_AmbWind == 1 ) then

      !---------------------------------------------------------------------------
      ! Parse turbine 1, 1st timestep, high res wind input file to gather the grid
      !    information and set data associated with turbine 1
      !---------------------------------------------------------------------------

      FileName = trim(p%WindFilePath)//trim(PathSep)//"HighT1"//trim(PathSep)//"Amb.t0.vtk"
      Un = -1 ! Set to force closing of file on return
      call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, errStat2, errMsg2 )
         call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName )
         if (errStat >= AbortErrLev) return
   else

      ! Using InflowWind, so data has been passed in via the InitInp data structure
      origin(1) = InitInp%InputFileData%X0_high(1)
      origin(2) = InitInp%InputFileData%Y0_high(1)
      origin(3) = InitInp%InputFileData%Z0_high(1)
      dims(1)   = InitInp%InputFileData%nX_high
      dims(2)   = InitInp%InputFileData%nY_high
      dims(3)   = InitInp%InputFileData%nZ_high
      gridSpacing(1) = InitInp%InputFileData%dX_high(1)
      gridSpacing(2) = InitInp%InputFileData%dY_high(1)
      gridSpacing(3) = InitInp%InputFileData%dZ_high(1)
      p%dt_high = InitInp%InputFileData%dt_high

   end if

   if ( (gridSpacing(1) <= 0.0_ReKi) .or. (gridSpacing(2) <= 0.0_ReKi) .or. (gridSpacing(3) <= 0.0_ReKi) ) &
      call SetErrStat ( ErrID_Fatal, 'The high resolution spatial resolution for Turbine 1 must be greater than zero in each spatial direction. ', errStat, errMsg, RoutineName )
      if (errStat >= AbortErrLev ) return

   p%nX_high          = dims(1)
   p%nY_high          = dims(2)
   p%nZ_high          = dims(3)
   p%X0_high(1)       = origin(1)
   p%Y0_high(1)       = origin(2)
   p%Z0_high(1)       = origin(3)
   p%dX_high(1)       = gridSpacing(1)
   p%dY_high(1)       = gridSpacing(2)
   p%dZ_high(1)       = gridSpacing(3)
   NumGrid_high       = p%nX_high*p%nY_high*p%nZ_high

   InitOut%X0_high(1) = origin(1)
   InitOut%Y0_high(1) = origin(2)
   InitOut%Z0_high(1) = origin(3)
   InitOut%dX_high(1) = gridSpacing(1)
   InitOut%dY_high(1) = gridSpacing(2)
   InitOut%dZ_high(1) = gridSpacing(3)

   if ( p%Mod_AmbWind == 1 ) then
         ! Just using this to make sure dims are >=2 points in each direction
      call HiResWindCheck(0, 1, p%nX_high, p%nY_high, p%nZ_high, p%dX_high(1), p%dY_high(1), p%dZ_high(1), p%X0_high(1), p%Y0_high(1), p%Z0_high(1), dims, gridSpacing, origin, RoutineName, errMsg2, errStat2)
         call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName )
         if (errStat >= AbortErrLev) return
   end if

   allocate( p%Grid_high(3,NumGrid_high,p%NumTurbines ),stat=errStat2)
      if (errStat2 /= 0) then
         call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for Grid_high.', errStat, errMsg, RoutineName )
         return
      end if

   nXYZ_high = 0
   do nz_high=0, p%nZ_high-1
      do ny_high=0, p%nY_high-1
         do nx_high=0, p%nX_high-1
            nXYZ_high = nXYZ_high + 1
            p%Grid_high(1,nXYZ_high,1) = InitOut%X0_high(1) + nx_high*InitOut%dX_high(1)
            p%Grid_high(2,nXYZ_high,1) = InitOut%Y0_high(1) + ny_high*InitOut%dY_high(1)
            p%Grid_high(3,nXYZ_high,1) = InitOut%Z0_high(1) + nz_high*InitOut%dZ_high(1)
         end do
      end do
   end do

   !---------------------------------------------------------------------------
   ! Parse the remaining turbine's 1st timestep, high res wind input files to
   !    gather the grid information and set data associated with those turbines
   !---------------------------------------------------------------------------

   do nt = 2, p%NumTurbines

      if ( p%Mod_AmbWind == 1 ) then
         FileName = trim(p%WindFilePath)//trim(PathSep)//"HighT"//trim(num2lstr(nt))//trim(PathSep)//"Amb.t0.vtk"
         Un = -1 ! Set to force closing of file on return
         call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, ErrStat2, ErrMsg2 )
            call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName )
            if (ErrStat >= AbortErrLev) return
      else
         ! Using InflowWind, so data has been passed in via the InitInp data structure
         origin(1) = InitInp%InputFileData%X0_high(nt)
         origin(2) = InitInp%InputFileData%Y0_high(nt)
         origin(3) = InitInp%InputFileData%Z0_high(nt)
         gridSpacing(1) = InitInp%InputFileData%dX_high(nt)
         gridSpacing(2) = InitInp%InputFileData%dY_high(nt)
         gridSpacing(3) = InitInp%InputFileData%dZ_high(nt)

      end if

      if ( (gridSpacing(1) <= 0.0_ReKi) .or. (gridSpacing(2) <= 0.0_ReKi) .or. (gridSpacing(3) <= 0.0_ReKi) ) &
         call SetErrStat ( ErrID_Fatal, 'The high resolution spatial resolution for Turbine '//trim(num2lstr(nt))//' must be greater than zero in each spatial direction. ', errStat, errMsg, RoutineName )

      InitOut%X0_high(nt) = origin(1)
      InitOut%Y0_high(nt) = origin(2)
      InitOut%Z0_high(nt) = origin(3)

      InitOut%dX_high(nt) = gridSpacing(1)
      InitOut%dY_high(nt) = gridSpacing(2)
      InitOut%dZ_high(nt) = gridSpacing(3)
      p%X0_high(nt) = origin(1)
      p%Y0_high(nt) = origin(2)
      p%Z0_high(nt) = origin(3)
      p%dX_high(nt) = gridSpacing(1)
      p%dY_high(nt) = gridSpacing(2)
      p%dZ_high(nt) = gridSpacing(3)

      if ( p%Mod_AmbWind == 1 ) then
            ! Using this to make sure dims are >=2 points in each direction, and number of grid points in each direction matches turbine 1
         call HiResWindCheck(0, nt, p%nX_high, p%nY_high, p%nZ_high, p%dX_high(nt), p%dY_high(nt), p%dZ_high(nt), p%X0_high(nt), p%Y0_high(nt), p%Z0_high(nt), dims, gridSpacing, origin, RoutineName, errMsg2, errStat2)
            call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName )
            if (ErrStat >= AbortErrLev) return
      end if

      nXYZ_high = 0
      do nz_high=0, p%nZ_high-1
         do ny_high=0, p%nY_high-1
            do nx_high=0, p%nX_high-1
               nXYZ_high = nXYZ_high + 1
               p%Grid_high(1,nXYZ_high,nt) = InitOut%X0_high(nt) + nx_high*InitOut%dX_high(nt)
               p%Grid_high(2,nXYZ_high,nt) = InitOut%Y0_high(nt) + ny_high*InitOut%dY_high(nt)
               p%Grid_high(3,nXYZ_high,nt) = InitOut%Z0_high(nt) + nz_high*InitOut%dZ_high(nt)
            end do
         end do
      end do

   end do

   InitOut%nx_high = p%nx_high
   InitOut%ny_high = p%ny_high
   InitOut%nz_high = p%nz_high

   if ( (InitInp%InputFileData%ChkWndFiles) .and. (p%Mod_AmbWind == 1) ) then
      do nt=1,p%NumTurbines
         do n=0,p%NumDT-1  ! We have already checked the first high-res files associated with n=0, but need to check the remaining, so for simplicity of code we will repeat the check on the first file

               ! We only have one high res input for for the very last low res timestep
            if ( n == (p%NumDT-1) ) then
               n_high_low = 1
            else
               n_high_low = p%n_high_low
            end if

            do nh=1,n_high_low
               nhigh = nh+n*p%n_high_low-1
               FileName = trim(p%WindFilePath)//trim(PathSep)//"HighT"//trim(num2lstr(nt))//trim(PathSep)//"Amb.t"//trim(num2lstr(nhigh))//".vtk"  !TODO: Should the turbine numbers be padding with leading zero(es)?
               Un = -1 ! Set to force closing of file on return
               call ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel, Un, ErrStat2, ErrMsg2 )
                  call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName )
                  if (ErrStat >= AbortErrLev) return

               call HiResWindCheck(nhigh, nt, p%nX_high, p%nY_high, p%nZ_high, p%dX_high(nt), p%dY_high(nt), p%dZ_high(nt), p%X0_high(nt), p%Y0_high(nt), p%Z0_high(nt), dims, gridSpacing, origin, RoutineName, errMsg2, errStat2)
                  call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName )
                  if (errStat >= AbortErrLev ) return

            end do
         end do
      end do
   end if

! End simulated read of low and high res ambient wind files
!==============================================================================

end subroutine AWAE_IO_InitGridInfo

!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE AWAE_PrintSum(  p, u, y, ErrStat, ErrMsg )
! This routine generates the summary file, which contains a summary of input file options.

      ! passed variables
   !TYPE(AWAE_InitInput),        INTENT(IN)  :: InputFileData                        ! Input-file data
   type(AWAE_ParameterType),    intent(in)  :: p                                    ! Parameters
   type(AWAE_InputType),        intent(in)  :: u                                    ! inputs
   type(AWAE_OutputType),       intent(in)  :: y                                    ! outputs
   integer(IntKi),            intent(out) :: ErrStat
   character(*),              intent(out) :: ErrMsg

      ! Local variables.

   INTEGER(IntKi)               :: I                                               ! Index for the nodes.
   INTEGER(IntKi)               :: UnSu                                            ! I/O unit number for the summary output file

   CHARACTER(*), PARAMETER      :: FmtDat    = '(A,T35,1(:,F13.3))'                ! Format for outputting mass and modal data.
   CHARACTER(*), PARAMETER      :: FmtDatT   = '(A,T35,1(:,F13.8))'                ! Format for outputting time steps.

   CHARACTER(30)                :: OutPFmt                                         ! Format to print list of selected output channels to summary file
   CHARACTER(100)               :: Msg                                             ! temporary string for writing appropriate text to summary file

   errStat = ErrID_None
   errMsg  = ""

   ! Open the summary file and give it a heading.

   CALL GetNewUnit( UnSu, ErrStat, ErrMsg )
   CALL OpenFOutFile ( UnSu, TRIM( p%OutFileRoot )//'.sum', ErrStat, ErrMsg )
   IF ( ErrStat >= AbortErrLev ) RETURN

   CLOSE(UnSu)

RETURN
END SUBROUTINE AWAE_PrintSum
!----------------------------------------------------------------------------------------------------------------------------------

END MODULE AWAE_IO

max_stars_repo_name: ndevelder/openfast

max_issues_repo_name: ndevelder/openfast

max_forks_repo_name: ndevelder/openfast


text:
      SUBROUTINE FILE_PROMPT( TEMPLATE_STR, PROMPT_STR, DIR_OUT,
     & NO_PROMPT, HIT, NEW_HIT, FILE_NAME, EOF )

!***********************************************************************
!* Prompt/Process File Name/Template/Hitlist
!*
!* Language: Fortran
!*
!* Author: Stuart G. Mentzer
!*
!* Date: 1999/10/18
!***********************************************************************

      ! Headers
      INCLUDE 'file.fi'
      INCLUDE 'uds_fxn.fi'

      ! Arguments ______________________________________________________

      CHARACTER*(*) TEMPLATE_STR ! Input file name or template
      ! Ignored if HIT.OPEN / Blanked after usage

      CHARACTER*(*) PROMPT_STR ! Initial prompt string

      CHARACTER*(*) DIR_OUT ! Output directory

      LOGICAL NO_PROMPT ! Specifies not to prompt for inputs

      RECORD /FILE/  HIT ! Hit file

      LOGICAL NEW_HIT ! Set to indicate new hit file entry or list

      CHARACTER*(*) FILE_NAME ! File name returned from top of list

      LOGICAL EOF ! Indicates user-entered EOF (Ctrl-Z)

      ! Variables ______________________________________________________

      LOGICAL PROMPT, DONE, NO_MATCH

      INTEGER IOS, IOS_W, I, LN, N_REC_OLD

      ! Functions ______________________________________________________

      CHARACTER UPCASE*50

      EXTERNAL UPCASE

      ! Initializations
      NEW_HIT = .FALSE.
      EOF = .FALSE.
      PROMPT = ( .NOT. NO_PROMPT )

      ! Check/read existing hitlist
      IF ( HIT.OPEN ) THEN ! Read from existing hitlist
        READ( HIT.UNIT, '(A)', IOSTAT=IOS ) FILE_NAME
        IF ( IOS .EQ. 0 ) THEN ! Return with the file name
          HIT.N_REC = HIT.N_REC - 1
          INQUIRE( FILE=FILE_NAME, NAME=FILE_NAME, IOSTAT=IOS )
          RETURN
        ELSE ! Hitlist file end or error
          CALL FILE_DEL( HIT )
        END IF
      ELSE ! Initialize HIT file structure
        CALL FILE_INIT( HIT )
      END IF

      ! Check/process input template / Prompt for template
  101 IF ( BLANK( TEMPLATE_STR ) ) THEN ! No input template
        IF ( .NOT. PROMPT ) THEN ! Return blank file name
          FILE_NAME = ' '
          RETURN
        END IF
        ! Prompt for file name/template/hitlist
        WRITE( *, '(/1X,A/'' >> '',$)' )
     &     PROMPT_STR(:L_TRIM(PROMPT_STR))
        READ( *, '(A)', END=199 ) FILE_NAME
        IF ( BLANK( FILE_NAME ) ) RETURN
      ELSE ! Use input template
        IF ( FILE_NAME .NE. TEMPLATE_STR ) PROMPT = .FALSE.
        FILE_NAME = TEMPLATE_STR
        TEMPLATE_STR = ' '
      END IF

      ! Check/process template or hitlist
      N_REC_OLD = 0
      CALL FILE_LIST( FILE_NAME, DIR_OUT, HIT, NO_MATCH )
      IF ( NO_MATCH ) THEN ! No matching file - Loop
        GO TO 101
      ELSE IF ( ( HIT.OPEN ) .AND. ( PROMPT ) ) THEN
        DONE = .FALSE.
        DO WHILE ( .NOT. DONE )
          WRITE( *, '(/2A,2X,A/'' >> '',$)' )
     &       ' Next file name/template/hitlist?   ',
     &       '( S - Show, U - Undo, E - Erase )',
     &       '[no more]'
          READ( *, '(A)', IOSTAT=IOS ) FILE_NAME
          IF ( IOS .EQ. -1 ) THEN ! Ctrl-Z entered - Start over
            CALL FILE_DEL( HIT )
            GO TO 101
          ELSE IF ( ( IOS .NE. 0 ) .OR. ( BLANK( FILE_NAME ) ) )
     &       THEN ! Done
            DONE = .TRUE.
          ELSE IF ( UPCASE( FILE_NAME ) .EQ. 'S' ) THEN ! Show list
            REWIND( HIT.UNIT )
            DO I = 1, HIT.N_REC
              READ( HIT.UNIT, '(A)', IOSTAT=IOS ) FILE_NAME
              IF ( IOS .EQ. 0 ) THEN
                WRITE( *, *, IOSTAT=IOS_W )
     &             FILE_NAME(:L_TRIM(FILE_NAME))
              END IF
            END DO
          ELSE IF ( UPCASE( FILE_NAME ) .EQ. 'U' ) THEN ! Undo last
            REWIND( HIT.UNIT )
            IF ( N_REC_OLD .GT. 0 ) THEN
              DO I = 1, N_REC_OLD - 1 ! Skip to last old record
                READ( HIT.UNIT, '(A)', IOSTAT=IOS )
              END DO
              READ( HIT.UNIT, '(A)', IOSTAT=IOS ) FILE_NAME
              IF ( IOS .NE. 0 ) THEN ! List corrupted
                WRITE( *, * ) '*** List corrupted - Deleting'
                CALL FILE_DEL( HIT )
                GO TO 101
              END IF
              BACKSPACE( HIT.UNIT )
              LN = LEN_TRIM( FILE_NAME )
              IF ( LN .GT. 0 ) THEN
                WRITE( HIT.UNIT, '(A)', IOSTAT=IOS_W )
     &             FILE_NAME(:LN)
              ELSE
                WRITE( HIT.UNIT, '()', IOSTAT=IOS_W )
              END IF
              HIT.N_REC = N_REC_OLD
            ELSE ! Delete list file
              CALL FILE_DEL( HIT )
              GO TO 101
            END IF
          ELSE IF ( UPCASE( FILE_NAME ) .EQ. 'E' ) THEN ! Erase list
            CALL FILE_DEL( HIT )
            GO TO 101
          ELSE ! Process new input
            N_REC_OLD = HIT.N_REC
            CALL FILE_LIST( FILE_NAME, DIR_OUT, HIT, NO_MATCH )
          END IF
        END DO
      END IF

      ! Set up to return a file name
      IF ( HIT.OPEN ) THEN ! Read from new hitlist
        REWIND( HIT.UNIT )
        READ( HIT.UNIT, '(A)', IOSTAT=IOS ) FILE_NAME
        IF ( IOS .NE. 0 ) THEN ! Hitlist end or error
          CALL FILE_DEL( HIT )
          FILE_NAME = ' '
          RETURN
        END IF
        HIT.N_REC = HIT.N_REC - 1
        IF ( HIT.N_REC .LE. 0 ) CALL FILE_DEL( HIT ) ! Return 1 file
      END IF
      INQUIRE( FILE=FILE_NAME, NAME=FILE_NAME, IOSTAT=IOS )
      NEW_HIT = .TRUE.

      RETURN

      ! Ctrl-Z handler
  199 CALL FILE_DEL( HIT )
      FILE_NAME = ' '
      EOF = .TRUE.
      RETURN

      END

ext: for

max_stars_repo_name: DeadParrot/NHTSA-Tools

max_issues_repo_name: DeadParrot/NHTSA-Tools

max_forks_repo_name: DeadParrot/NHTSA-Tools


text:
	SUBROUTINE DSDASH  ( szdsh, size, iret )
C************************************************************************
C* DSDASH								*
C* 									*
C* This subroutine sets the line dashing scale.  If the parameter is	*
C* not positive, no change is made.					*
C* 									*
C* DSDASH  ( SZDSH, SIZE, IRET )					*
C*                                                                    	*
C* Input parameters:							*
C* 	SZDSH		REAL		Line dashing scale		*
C*                                                                    	*
C* Output parameters:							*
C*	SIZE		REAL		Line dashing scale		*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Jacobs/NCEP	 3/98						*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVACT.CMN'
	INCLUDE		'DEVCHR.CMN'
C-----------------------------------------------------------------------
	iret = NORMAL
C*
	IF  ( szdsh .gt. 0. )  tszdsh = szdsh
	size = tszdsh
C
C*	If the driver is VG, send the attributes to the device
C*	for output.
C
	IF  ( ddev .eq. 'VG' )  THEN
	    CALL HSDASH ( size, iret )
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! { dg-do run }
!
! CO_SUM/CO_MIN/CO_MAX
!
program test
  implicit none
  intrinsic co_max
  intrinsic co_min
  intrinsic co_sum
  call test_min
  call test_max
  call test_sum
contains
  subroutine test_max
    integer :: val
    val = this_image ()
    call co_max (val, result_image=1)
    if (this_image() == 1) then
      !write(*,*) "Maximal value", val
      if (val /= num_images()) STOP 1
    end if
  end subroutine test_max

  subroutine test_min
    integer :: val
    val = this_image ()
    call co_min (val, result_image=1)
    if (this_image() == 1) then
      !write(*,*) "Minimal value", val
      if (val /= 1) STOP 2
    end if
  end subroutine test_min

  subroutine test_sum
    integer :: val, n
    val = this_image ()
    call co_sum (val, result_image=1)
    if (this_image() == 1) then
      !write(*,*) "The sum is ", val
      n = num_images()
      if (val /= (n**2 + n)/2) STOP 3
    end if
  end subroutine test_sum
end program test


text:
	INCLUDE 'VICMAIN_FOR'

	SUBROUTINE MAIN44
	IMPLICIT NONE
	INTEGER	COUNT, DEF
	INTEGER	INPCOUNT, OUTCOUNT
	INTEGER	DIM, OUTDIM
	INTEGER	I, GRFILE
	REAL	FIRST, SECOND, THIRD(40)
	REAL	ZEROS(40)/40*0.0/
	REAL	X, Y, X1, Y1, X2, Y2, DX, DY, DIST, INTERVAL, DELTA, T
	LOGICAL	EOF, EOL,   ALLZERO
	CHARACTER*16  OPTION

	CALL XVPARM ('OPTION', OPTION, COUNT, DEF, 0)

	CALL XVPARM ('DIM', DIM, COUNT, DEF, 0)

	IF (OPTION(1:6) .EQ. 'APPEND') THEN

	    CALL wrgr_grutil (1, 1, DIM)
	    CALL XVPCNT ('INP', INPCOUNT)
	    DO I = 1, INPCOUNT
		CALL rdgr_grutil (I, I+1, DIM)
		EOF = .FALSE.
		DO WHILE (.NOT. EOF)
		    CALL NEXTGR (I+1, EOF, FIRST, SECOND, THIRD)
		    IF (.NOT. EOF) THEN
			EOL = .FALSE.
			DO WHILE (.NOT. EOL)
			    CALL PUTGR (1, FIRST, SECOND, THIRD)
			    CALL GETGR (I+1, EOL, EOF, FIRST, SECOND, THIRD)
			ENDDO
			CALL PUTGR (1, 0.0, 0.0, ZEROS)
		    ENDIF
		ENDDO
		CALL CLGR (I+1)
	    ENDDO
	    CALL CLGR (1)

	ELSE IF (OPTION(1:4) .EQ. 'SWAP') THEN

	    CALL rdgr_grutil (1, 1, DIM)
	    CALL wrgr_grutil (1, 2, DIM)
	    EOF = .FALSE.
	    DO WHILE (.NOT. EOF)
		CALL NEXTGR (1, EOF, FIRST, SECOND, THIRD)
		IF (.NOT. EOF) THEN
		    EOL = .FALSE.
		    DO WHILE (.NOT. EOL)
			CALL PUTGR (2, SECOND, FIRST, THIRD)
			CALL GETGR (1, EOL, EOF, FIRST, SECOND, THIRD)
		    ENDDO
		    CALL PUTGR (2, 0.0, 0.0, ZEROS)
		ENDIF
	    ENDDO
	    CALL CLGR (1)
	    CALL CLGR (2)

	ELSE IF (OPTION(1:6) .EQ. '2DTO3D') THEN

	    CALL wrgr_grutil (1, 1, 3)
	    CALL XVPCNT ('INP', INPCOUNT)
	    DO I = 1, INPCOUNT
		CALL rdgr_grutil (I, I+1, 2)
		THIRD(1) = FLOAT(I)
		EOF = .FALSE.
		DO WHILE (.NOT. EOF)
		    CALL NEXTGR (I+1, EOF, FIRST, SECOND)
		    IF (.NOT. EOF) THEN
			EOL = .FALSE.
			DO WHILE (.NOT. EOL)
			    CALL PUTGR (1, FIRST, SECOND, THIRD)
			    CALL GETGR (I+1, EOL, EOF, FIRST, SECOND)
			ENDDO
			CALL PUTGR (1, 0.0, 0.0, 0.0)
		    ENDIF
		ENDDO
		CALL CLGR (I+1)
	    ENDDO
	    CALL CLGR (1)

	ELSE IF (OPTION(1:6) .EQ. '3DTO2D') THEN

	    CALL rdgr_grutil (1, 1, 3)
	    CALL XVPCNT ('OUT', OUTCOUNT)
	    DO I = 1, OUTCOUNT
		CALL wrgr_grutil (I, I+1, 2)
	    ENDDO
	    EOF = .FALSE.
	    DO WHILE (.NOT. EOF)
		CALL NEXTGR (1, EOF, FIRST, SECOND, THIRD)
		GRFILE = INT(THIRD(1))
		IF (.NOT. EOF .AND.
     +			GRFILE .GE. 1 .AND. GRFILE .LE. OUTCOUNT) THEN
		    EOL = .FALSE.
		    DO WHILE (.NOT. EOL)
			CALL PUTGR (GRFILE+1, FIRST, SECOND)
			CALL GETGR (1, EOL, EOF, FIRST, SECOND, THIRD)
		    ENDDO
		    CALL PUTGR (GRFILE+1, 0.0, 0.0)
		ENDIF
	    ENDDO
	    CALL CLGR (1)
	    DO I = 1, OUTCOUNT
		CALL CLGR (I+1)
	    ENDDO

	ELSE IF (OPTION(1:5) .EQ. 'REDIM') THEN

	    CALL XVPARM ('OUTDIM', OUTDIM, COUNT, DEF, 0)
	    CALL rdgr_grutil (1, 1, DIM)
	    CALL wrgr_grutil (1, 2, OUTDIM)
	    SECOND = 0.0
	    DO I = 3, OUTDIM
		THIRD(I-2) = 0.0
	    ENDDO
	    EOF = .FALSE.
	    DO WHILE (.NOT. EOF)
		CALL NEXTGR (1, EOF, FIRST, SECOND, THIRD)
		IF (.NOT. EOF) THEN
		    EOL = .FALSE.
		    DO WHILE (.NOT. EOL)
			CALL PUTGR (2, FIRST, SECOND, THIRD)
			CALL GETGR (1, EOL, EOF, FIRST, SECOND, THIRD)
		    ENDDO
		    CALL PUTGR (2, 0.0, 0.0, ZEROS)
		ENDIF
	    ENDDO
	    CALL CLGR (1)
	    CALL CLGR (2)

	ELSE IF (OPTION(1:4) .EQ. 'FILL') THEN

	    CALL XVP ('INTERVAL', INTERVAL, COUNT)
	    CALL rdgr_grutil (1, 1, DIM)
	    CALL wrgr_grutil (1, 2, DIM)

	    DO WHILE (.TRUE.)
C			Scan for the beginning of a line string
		X1 = 0.0
		Y1 = 0.0
		DO WHILE (X1 .EQ. 0.0 .AND. Y1 .EQ. 0.0)
		    CALL GETGR (1, ALLZERO, EOF, Y1, X1, THIRD)
		    IF ( EOF ) GO TO 30
		ENDDO
C			Put out the first point in the linestring
		CALL PUTGR (2, Y1, X1, THIRD)
						! Get the next pair
		CALL GETGR (1, ALLZERO, EOF, Y2, X2, THIRD)
		IF ( EOF ) GO TO 30
		EOL = (X2 .EQ. 0.0 .AND. Y2 .EQ. 0.0)

		IF (.NOT. EOL ) THEN
		    DO WHILE (.NOT. EOL .AND. .NOT. EOF)
			DY = Y2 - Y1
			DX = X2 - X1
			DIST = SQRT( DY**2 + DX**2 )
			IF (DIST .GT. 0.0) THEN
			    DELTA = INTERVAL/DIST
			ELSE
			    DELTA = 2.0
			ENDIF
				! Fill in the points in between
			T = DELTA
			DO WHILE (T .LT. 1.0)
			    Y = Y1 + T*DY
			    X = X1 + T*DX
			    CALL PUTGR (2, Y, X, THIRD)
			    T = T + DELTA
			ENDDO
			CALL PUTGR (2, Y2, X2, THIRD)
				! Get the next point in linestring
			Y1 = Y2
			X1 = X2
			CALL GETGR (1, ALLZERO, EOF, Y2, X2, THIRD)
			EOL = (X2 .EQ. 0.0 .AND. Y2 .EQ. 0.0)
			IF ( EOF ) GO TO 30
		    ENDDO
		ENDIF
		CALL PUTGR (2, 0.0, 0.0, ZEROS)
	    ENDDO
   30	    CONTINUE
	    CALL CLGR (1)
	    CALL CLGR (2)

	ELSE IF (OPTION(1:6) .EQ. 'SMOOTH') THEN

	    CALL XVP ('INTERVAL', INTERVAL, COUNT)
	    INTERVAL = INTERVAL**2
	    CALL rdgr_grutil (1, 1, DIM)
	    CALL wrgr_grutil (1, 2, DIM)

	    DO WHILE (.TRUE.)
C			Scan for the beginning of a line string
		X1 = 0.0
		Y1 = 0.0
		DO WHILE (X1 .EQ. 0.0 .AND. Y1 .EQ. 0.0)
		    CALL GETGR (1, ALLZERO, EOF, Y1, X1, THIRD)
		    IF ( EOF ) GO TO 40
		ENDDO
C			Put out the first point in the linestring
		CALL PUTGR (2, Y1, X1, THIRD)

		EOL = .FALSE.
		DO WHILE (.NOT. EOL)
				! Get the next point in linestring
		    CALL GETGR (1, ALLZERO, EOF, Y2, X2, THIRD )
		    IF ( EOF ) GO TO 40
		    EOL = (X2 .EQ. 0.0 .AND. Y2 .EQ. 0.0)
		    IF (.NOT. EOL) THEN
C		       Output if distance from last output point is large enough
			DIST = (Y2-Y1)**2 + (X2-X1)**2
			IF (DIST .GE. INTERVAL) THEN
			    CALL PUTGR (2, Y2, X2, THIRD)
			    Y1 = Y2
			    X1 = X2
			ENDIF
		    ENDIF
		ENDDO
		CALL PUTGR (2, 0.0, 0.0, ZEROS)
	    ENDDO
 40	    CONTINUE
	    CALL CLGR (1)
	    CALL CLGR (2)

	ENDIF

	RETURN
	END

c
c	Code added to check for missing files - fr 088285
c

      subroutine wrgr_grutil(instance,vicarfile,dimension)

      integer instance, vicarfile, dimension
      integer status, wrgr

      status = wrgr(instance,vicarfile,dimension)
      if ( status .ne. 1 ) call signalgr(vicarfile,status,1)

      return
      end

      subroutine rdgr_grutil(instance,vicarfile,dimension)

      integer instance, vicarfile, dimension
      integer status, rdgr

      status = rdgr(instance,vicarfile,dimension)
      if ( status .ne. 1 ) call signalgr(vicarfile,status,1)

      return
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download

!> This module supports both the obs_impact_tool and reading in the
!> table at runtime for use during the assimilation phase, to
!> alter the impact of observations on the state vector.
!>
!> these routines build a table(ntypes, nqtys) which can be
!> indexed quickly at runtime.  this happens in the assimilation
!> loop so performance matters.

module obs_impact_mod

use      types_mod, only : r8, obstypelength, missing_r8
use  utilities_mod, only : error_handler, E_ERR, E_MSG,       &
                           open_file, close_file, get_next_filename, to_upper
use  obs_kind_mod        ! all qtys/types, so impossible to enumerate them here
use parse_args_mod, only : get_args_from_string

implicit none
private

public :: create_impact_table,   &
          allocate_impact_table, &
          read_impact_table,     &
          free_impact_table

character(len=*), parameter :: source = 'obs_impact_mod.f90'

integer :: ios, i

! make all strings the same as the max length of the parameter limit
integer, parameter :: string_length = obstypelength

integer, parameter :: MAXLINELEN = 512
character(len=MAXLINELEN) :: readbuf
character(len=6) :: readformat = '(A512)'

! global line number for error messages
integer :: linecount
character(len=MAXLINELEN+32) :: errline

integer, parameter :: MAXWORDS = 1024
integer :: wordcount
character(len=string_length+1) :: wordarray(MAXWORDS)
integer :: itemlist(MAXWORDS)

! table of contents entry types
integer, parameter :: ENTRY_UNKNOWN  = 0
integer, parameter :: ENTRY_DARTQTY  = 1
integer, parameter :: ENTRY_DARTTYPE = 2
integer, parameter :: ENTRY_GROUP    = 3
integer, parameter :: ENTRY_KEYWORD  = 4

! predefined keywords, group names, and dart type/qty names
type type_entry
  character(len=string_length) :: ename
  integer :: etype
  integer :: evalue
end type type_entry

type type_toc
   integer :: toc_count
   integer :: type_count
   integer :: qty_count
   type(type_entry), allocatable :: toc_entries(:)
end type

! this is the single, global table of contents.
! FIXME: toc should maybe be dictionary (dict?)
! but for modularity, we include it as an argument
! in all the calls, generally last except for optional args.

! user-defined groups
type type_group_entry
  integer :: member_count
  integer, allocatable :: group_members(:)
end type type_group_entry

type type_group
  integer :: group_count
  type(type_group_entry), allocatable :: groups(:)
end type type_group

! this is the single, global table of groups.
! but for modularity, we include it as an argument
! in all the calls, generally last except for optional args.

integer :: this_state
integer, parameter :: STATE_UNDEFINED = 0
integer, parameter :: STATE_DEFGROUP = 1
integer, parameter :: STATE_DEFNOTGROUP = 2
integer, parameter :: STATE_DEFIMPACT = 3

! are we processing ALL, ALLTYPES, or ALLQTYS
integer :: this_category

! define a comment character (#, %, !) which makes the rest
! of the input line a comment
character, parameter :: commentchar = '#'

character(len=string_length) :: current_group_name = ''
integer :: current_group_id = 0
logical :: otherwork

logical, save :: module_initialized = .false.

character(len=512) :: msgstring, msgstring2, msgstring3

! .true. gives more output
logical :: debug = .false.

contains

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

! TOOL:
subroutine create_impact_table(input_filename, output_filename, debug_flag)
character(len=*),  intent(in) :: input_filename
character(len=*),  intent(in) :: output_filename
logical,           intent(in), optional :: debug_flag

! this is the routine that reads in the config file
! and creates an output file that's suitable for reading
! at runtime.  so this is the 'create tool' part.

type(type_toc)        :: toc
type(type_group)      :: group_toc
real(r8), allocatable :: table(:,:)
integer               :: funit

! initialization and setup
call initialize_module(debug)

! set up space for the output table
call allocate_impact_table(table)

! read in types and qty strings, compute sizes,
! allocate space for the toc and group types

call init_strings_and_toc(toc, group_toc)

this_state = STATE_UNDEFINED

! read in ascii control file with and build table
call error_handler(E_MSG, 'obs_impact_tool', ' reading impacts from: '//trim(input_filename))

linecount = 0
!print *, 'file contents: '

! loop until read returns nothing
funit = open_file(input_filename, action='read')
readloop: do
   ! get the next line
   read(funit, readformat, iostat=ios) readbuf
   if (ios /= 0) exit readloop
   linecount = linecount + 1

!print *, 'read input line, was: '//trim(readbuf)

   ! prepare msg in case of error
   write(errline, '(A,I6,A)') 'error on line ', linecount, '; line contents: '

   ! parse it into white-space separated words
   call get_args_from_string(readbuf, wordcount, wordarray)
   if (wordcount < 1) cycle readloop

   ! count actual number of words (e.g. stop at comment char)
   ! and note which ones are already in the TOC.  do not add
   ! new ones if they do not yet exist.
   call search_toc(wordcount, wordarray, itemlist, toc)

   ! change state based on keywords, add new group names to toc
   call state_change(wordcount, wordarray, itemlist, toc, group_toc, this_state, otherwork)

   if (otherwork) then
      if (this_state == STATE_DEFGROUP) then
         ! add to groups if defining groups.
         call add_to_group(wordcount, itemlist, current_group_id, group_toc, toc)
      else if (this_state == STATE_DEFNOTGROUP) then
         ! remove from groups if defining exclusion groups.
         call sub_from_group(wordcount, itemlist, current_group_id, this_category, group_toc, toc)
      else if (this_state == STATE_DEFIMPACT) then
         ! build impact table based on what is going on now
         call update_impact(table, wordcount, itemlist, wordarray, toc, group_toc)
      endif
   endif

enddo readloop
call close_file(funit)

call error_handler(E_MSG, 'obs_impact_tool', ' end of input')

! anything still missing will be set to 1.0
where (table(:,:) == missing_r8) table(:,:) = 1.0_r8

call output_impact_table(output_filename, table, toc)

if (.true.) then  ! was debug
   write(msgstring, *) 'closing file ',  trim(output_filename)
   call error_handler(E_MSG, 'obs_impact_tool', msgstring)
endif

call finalize_tables(toc, group_toc, table)

end subroutine create_impact_table

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

! TOOL & RUNTIME:
subroutine allocate_impact_table(table, ntypes, nqtys)

real(r8), allocatable, intent(out) :: table(:,:)
integer, optional,     intent(out) :: ntypes
integer, optional,     intent(out) :: nqtys

integer :: qty_count, type_count

! output table is dimensioned (numtypes, 0:numqtys)
! space for results, and initial values
! default to 'unset'.  at the end, anything unset will be
! changed to 1.0, which means impact like usual with no change.

qty_count = get_num_quantities()
type_count = get_num_types_of_obs()

allocate(table(type_count, 0:qty_count))
table(:,:) = missing_r8

if (present(ntypes)) ntypes = type_count
if (present(nqtys)) nqtys = qty_count

end subroutine allocate_impact_table

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

! this was more like the create_impact_table routine but has since
! been replaced with a routine that reads 3 words per line, errors
! out if not 3, and then calls the set routine immediately.  input
! lines have to be:
!  type_name qty_name real_value
! at this stage in the process.  it identifies words which aren't
! type or qty names without a dictionary or state machine.

! RUNTIME:
subroutine read_impact_table(sourcefile, table, allow_any_values, anyvals_string)

character(len=*), intent(in)    :: sourcefile
real(r8),         intent(inout) :: table(:,0:)
logical,          intent(in)    :: allow_any_values
character(len=*), intent(in), optional :: anyvals_string

integer :: i, j
integer :: funit
character(len=obstypelength) :: typename, qtyname
real(r8) :: rvalue

! assumes the table is already allocated the proper size.
! could check here for that?

call error_handler(E_MSG, 'obs_impact_table', ' reading impacts from: '//trim(sourcefile))

linecount = 0
!print *, 'file contents: '

! loop until read returns nothing
funit = open_file(sourcefile, action='read')
readloop: do
   ! get the next line
   read(funit, readformat, iostat=ios) readbuf
   if (ios /= 0) exit readloop
   linecount = linecount + 1

!print *, 'read input line, was: '//trim(readbuf)

   ! prepare msg in case of error
   write(errline, '(A,I6,A)') 'error on line ', linecount, '; line contents: '

   ! parse it into 3 white-space separated words, two strings and one real
   read(readbuf, *, iostat=ios) typename, qtyname, rvalue
   if (ios /= 0) then
      call error_handler(E_ERR, 'read_impact_table', &
              'lines must contain a specific type, a generic qty, and a real value', &
              source, text2=errline, text3=readbuf)
   endif
!print *, trim(typename)//' '//trim(qtyname)//' ', rvalue

   call set_impact(table, typename, qtyname, rvalue, allow_any_values, anyvals_string)

enddo readloop
call close_file(funit)

! almost ready for use - set missing entries to 1.0_r8
where (table(:,:) == missing_r8) table(:,:) = 1.0_r8

! debug only?
do j=0, ubound(table, 2)
   do i=1, ubound(table, 1)
      if (table(i, j) /= 1.0_r8) then
         write(errline, '(2A33,F12.6)') get_name_for_type_of_obs(i), get_name_for_quantity(j), table(i,j)
         call error_handler(E_MSG, 'read_impact_table', errline)
      endif
      ! if (table(i, j) /= 1.0_r8) print *, i, j, table(i,j)
   enddo
enddo

end subroutine read_impact_table

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

! RUNTIME:
subroutine set_impact(table, typename, qtyname, rvalue, allow_any_values, anyvals_string)

real(r8),         intent(inout) :: table(:,0:)
character(len=*), intent(in)    :: typename
character(len=*), intent(in)    :: qtyname
real(r8),         intent(in)    :: rvalue
logical,          intent(in)    :: allow_any_values
character(len=*), intent(in), optional :: anyvals_string

! change this to true to force values to be 0 or 1
logical, save :: fully_on_or_off = .false.

integer :: index1, index2

! expect exactly 3 tokens on these lines:
!  type  qty  real(r8)

index1 = get_index_for_type_of_obs(typename)
index2 = get_index_for_quantity(qtyname)
!print *, 'in set_impact, type, qty index values are: ', index1, index2

if (index1 < 0) then
   call error_handler(E_ERR, 'obs_impact', &
                     'first word on line is unrecognized specific obs TYPE', &
                     source, text2=errline, text3=readbuf)
endif
if (index2 < 0) then
   call error_handler(E_ERR, 'obs_impact', &
                     'second word on line is unrecognized generic QTY', &
                     source, text2=errline, text3=readbuf)
endif

! options for the actual impact value:
!   1. anything goes
!   2. values restricted
!       a. must be 0.0 or 1.0 only
!       b. must be between 0.0 and 1.0 inclusive

if (.not. allow_any_values) then
   if (fully_on_or_off) then
      if (rvalue /= 0.0_r8 .and. rvalue /= 1.0_r8) then
         call error_handler(E_ERR, 'obs_impact', &
                           'impact values must be 0 or 1', &
                            source, text2=readbuf)
      endif
   else
      if (present(anyvals_string)) then
         msgstring3='set "'//trim(anyvals_string)//'=.true." in namelist to allow'
      else
         msgstring3=""
      endif
      if (rvalue < 0.0_r8 .or. rvalue > 1.0_r8) then
         call error_handler(E_ERR, 'obs_impact', &
                           'impact values must be between 0 and 1, inclusive', &
                           source, text2=readbuf, text3=msgstring3)
      endif
   endif
endif

! ok, we have the first, second, and value.  update the table.

if ((table(index1, index2) /= missing_r8) .and. &
    (table(index1, index2) /= rvalue)) then
   write(msgstring, *) 'error; different impact factor already set for this pair'
   write(msgstring2, *) 'previous value was', table(index1, index2), ' new value is ', rvalue
   write(msgstring3, *) trim(typename), ' and ', trim(qtyname)
   call error_handler(E_ERR, 'obs_impact', msgstring, &
                     source, text2=msgstring2, text3=msgstring3)

endif

table(index1, index2) = rvalue

end subroutine set_impact

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

! TOOL & RUNTIME:
subroutine free_impact_table(table)

real(r8), intent(inout), allocatable :: table(:,:)

if (allocated(table)) deallocate(table)

end subroutine free_impact_table

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

! TOOL & RUNTIME:
subroutine initialize_module(debug_flag)
logical, intent(in), optional :: debug_flag

if (module_initialized) return

module_initialized = .true.

if (present(debug_flag)) debug = debug_flag

end subroutine initialize_module

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

! TOOL:
subroutine init_strings_and_toc(toc, group_toc)

type(type_toc),        intent(inout) :: toc
type(type_group),      intent(inout) :: group_toc

integer :: i, type_count, qty_count, ntoc
integer :: groupsize, membersize
character(len=obstypelength), allocatable :: knownqtys(:)
character(len=obstypelength), allocatable :: knowntypes(:)

! first, get all possible generic qtys from obs_kind_mod
! BEWARE!! the first qty number is 0!!!!!
qty_count = get_num_quantities()
allocate(knownqtys(0:qty_count))
do i=0, qty_count
   knownqtys(i) = get_name_for_quantity(i)
   call to_upper(knownqtys(i))
enddo

! then, get all possible specific types from obs_kind_mod
type_count = get_num_types_of_obs()
allocate(knowntypes(type_count))
do i=1, type_count
   knowntypes(i) = get_name_for_type_of_obs(i)
   call to_upper(knowntypes(i))
enddo

! toc size  - arbitrary size, could be larger
ntoc = (type_count + qty_count) * 4

! space for groups - both sizes completely arbitrary
groupsize = qty_count * 2
membersize = qty_count * 4
group_toc%group_count = groupsize
allocate(group_toc%groups(groupsize))
do i=1, groupsize
   group_toc%groups(i)%member_count = 0
   allocate(group_toc%groups(i)%group_members(membersize))
enddo

allocate(toc%toc_entries(ntoc))
toc%toc_count = 0

call build_toc(qty_count, knownqtys, type_count, knowntypes, toc)

deallocate(knownqtys, knowntypes)

end subroutine init_strings_and_toc

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

! TOOL:
subroutine finalize_tables(toc, group_toc, table)

type(type_toc),        intent(inout) :: toc
type(type_group),      intent(inout) :: group_toc
real(r8), allocatable, intent(inout), optional :: table(:,:)

integer :: i

do i=1, group_toc%group_count
   deallocate(group_toc%groups(i)%group_members)
enddo
deallocate(group_toc%groups)
group_toc%group_count = 0

deallocate(toc%toc_entries)
toc%toc_count = 0
toc%qty_count = 0
toc%type_count = 0

if (present(table)) call free_impact_table(table)

end subroutine finalize_tables

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

subroutine search_toc(wordcount, wordarray, itemlist, toc)

integer,          intent(inout) :: wordcount
character(len=*), intent(in)    :: wordarray(:)
integer,          intent(out)   :: itemlist(:)
type(type_toc),   intent(inout) :: toc

integer :: wid
character(len=string_length) :: nextword
logical :: exists

! initialize the return values so we can just return
! when we reach the end of the list
itemlist(:) = -1

! make sure there is at least one word
if (wordcount < 1) return

do wid=1, wordcount

   nextword = wordarray(wid)
   call to_upper(nextword)

   if (is_comment(nextword)) then
      wordcount = wid-1
      return
   endif

   ! look up id numbers for keywords which are already in the toc
   ! but do not add new ones yet
   exists = in_toc(nextword, toc, itemlist(wid))

enddo

end subroutine search_toc

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

subroutine build_toc(nqtys, qtynames, ntypes, typenames, toc)

integer,          intent(in)    :: nqtys
character(len=*), intent(in)    :: qtynames(0:)
integer,          intent(in)    :: ntypes
character(len=*), intent(in)    :: typenames(:)
type(type_toc),   intent(inout) :: toc

integer :: i

! add keyword entries to toc table
call add_keywords(toc)

! a bit specialized but types and qty counts seem to be
! needed everywhere, so put them where we can access them.
call add_limits(toc, ntypes, nqtys)

! at some point we should make types and qtys disjoint
! sets of integer parameter values, because then there is
! no confusion about using one when you meant the other.
! (e.g. do not use 1->n at all; start one at 1000 and the
! other at 5000, and subtract the offset each time before
! using.)

! the first qty has a value of 0
do i=0, nqtys
   if (qtynames(i) /= 'UNKNOWN') then
      call add_toc(qtynames(i), toc, ENTRY_DARTQTY, i)
   endif
enddo

do i=1, ntypes
   if (typenames(i) /= 'UNKNOWN') then
      call add_toc(typenames(i), toc, ENTRY_DARTTYPE, i)
   endif
enddo

end subroutine build_toc

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

function in_toc(nextword, toc, t_index)

character(len=*),  intent(in)  :: nextword
type(type_toc),    intent(in)  :: toc
integer, optional, intent(out) :: t_index
logical                        :: in_toc

integer :: i

if (present(t_index)) t_index = -1
in_toc = .false.

do i=1, toc%toc_count
   if (toc%toc_entries(i)%ename == nextword) then
      if (present(t_index)) t_index = i
      in_toc = .true.
      return
   endif
enddo

! fall through intentionally and return false (and index -1) here
!print *, 'did not find '//trim(nextword)//' in toc'

end function in_toc

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

subroutine add_toc(nextword, toc, etype, evalue)

character(len=*),  intent(in)    :: nextword
type(type_toc),    intent(inout) :: toc
integer,           intent(in)    :: etype
integer, optional, intent(in)    :: evalue

integer :: i

do i=1, toc%toc_count
   if (toc%toc_entries(i)%ename == nextword) then
      call error_handler(E_ERR, 'add_toc:', 'word already in toc: '//trim(nextword), &
                       source, text2='cannot have duplicate entries')
   endif
enddo

toc%toc_entries(i)%etype = etype
toc%toc_entries(i)%ename = nextword
if (present(evalue)) then
   toc%toc_entries(i)%evalue = evalue
else
   toc%toc_entries(i)%evalue = -1
endif

toc%toc_count = toc%toc_count + 1

end subroutine add_toc

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

function get_name_by_value(evalue, etype, toc)

integer,          intent(in) :: evalue
integer,          intent(in) :: etype
type(type_toc),   intent(in) :: toc
character(len=string_length) :: get_name_by_value

integer :: i

get_name_by_value = 'UNKNOWN'

ENTRIES: do i=1, toc%toc_count
   if (toc%toc_entries(i)%etype /= etype) cycle ENTRIES

   if (toc%toc_entries(i)%evalue == evalue) then
      get_name_by_value = toc%toc_entries(i)%ename
      return
   endif
enddo ENTRIES

end function get_name_by_value

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

subroutine print_toc(toc, entry)

type(type_toc),    intent(in) :: toc
integer, optional, intent(in) :: entry

integer :: i

if (present(entry)) then
   write(*,'(A34,I5)') 'table of contents: item : ', entry
   write(*, '(A32,I3,I5)') trim(toc%toc_entries(entry)%ename), &
                              toc%toc_entries(entry)%etype, toc%toc_entries(entry)%evalue
else
   write(*,'(A34,I5)') 'table of contents: item count: ', toc%toc_count
   do i=1, toc%toc_count
      write(*, '(I5,A32,I3,I5)') i, trim(toc%toc_entries(i)%ename), &
                              toc%toc_entries(i)%etype, toc%toc_entries(i)%evalue
   enddo
endif

end subroutine print_toc

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

subroutine add_keywords(toc)

type(type_toc),   intent(inout) :: toc

! if you need to add any new keywords here is
! where you do that.  you also have to add the
! code to process what they do some place else.

call add_toc("GROUP",    toc, ENTRY_KEYWORD)
call add_toc("IMPACT",   toc, ENTRY_KEYWORD)

call add_toc("END",      toc, ENTRY_KEYWORD)
call add_toc("EXCEPT",   toc, ENTRY_KEYWORD)

call add_toc("UNKNOWN",  toc, ENTRY_KEYWORD)
call add_toc("ALL",      toc, ENTRY_KEYWORD)

call add_toc("ALLQTYS",  toc, ENTRY_KEYWORD)
call add_toc("ALLTYPES", toc, ENTRY_KEYWORD)

end subroutine add_keywords

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

function is_keyword(item, name, toc)

integer,            intent(in) :: item
character(len=*),   intent(in) :: name
type(type_toc),     intent(in) :: toc
logical                        :: is_keyword

is_keyword = .false.

if ((get_toc_name(item, toc) == name) .and. &
    (get_toc_type(item, toc) == ENTRY_KEYWORD)) is_keyword = .true.

! returns here if not keyword

end function is_keyword

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

subroutine add_limits(toc, ntypes, nqtys)

type(type_toc),   intent(inout) :: toc
integer,          intent(in)    :: ntypes
integer,          intent(in)    :: nqtys

toc%qty_count = nqtys
toc%type_count = ntypes

end subroutine add_limits

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

subroutine state_change(wordcount, wordarray, itemlist, toc, group_toc, this_state, otherwork)

integer,          intent(inout) :: wordcount
character(len=*), intent(in)    :: wordarray(:)
integer,          intent(in)    :: itemlist(:)
type(type_toc),   intent(inout) :: toc
type(type_group), intent(inout) :: group_toc
integer,          intent(inout) :: this_state
logical,          intent(out)   :: otherwork

! wordcount has been updated to take into account comment chars
! words have been uppercased already, and itemlist is -1 if the
! word is not yet in the toc.

! if the current line is changing state, there is no other work
! to be done.  if it is not a state line, it is adding items to
! a block or defining a relationship and needs to be processed further.
otherwork = .false.

! blank lines, comment lines, etc cannot change state
if (wordcount < 1) return

if (itemlist(1) < 0) then
   call error_handler(E_ERR, 'state_change',  &
             'first word on line is unrecognized keyword, QTY, TYPE, or group name', &
             source, text2=errline, text3=readbuf)
endif

select case (this_state)

   case (STATE_UNDEFINED)
      if (is_keyword(itemlist(1), 'GROUP', toc)) then
         !GROUP - start definition
         if (wordcount < 2) then
            call error_handler(E_ERR, 'state_change', &
                      'GROUP keyword must be followed by a group name', &
                      source, text2=errline, text3=readbuf)
         endif
         if (itemlist(2) > 0) then
            call error_handler(E_ERR, 'state_change', &
                      'group name '//trim(wordarray(2))//' already in use', &
                      source, text2=errline, text3=readbuf)
         endif

         this_state = STATE_DEFGROUP
         current_group_name = wordarray(2)
         call to_upper(current_group_name)
         current_group_id = current_group_id + 1
         call add_toc(current_group_name, toc, ENTRY_GROUP, current_group_id)
!print *, 'adding group name to toc, name,id:', trim(wordarray(2)), current_group_id

      else if (is_keyword(itemlist(1), 'IMPACT', toc)) then
         !IMPACT - start definition
         if (wordcount > 1) then
            call error_handler(E_ERR, 'state_change', &
                      'IMPACT keyword should not be followed by additional text', &
                      source, text2=errline, text3=readbuf)
         endif

         this_state = STATE_DEFIMPACT

      else if (is_keyword(itemlist(1), 'END', toc)) then
         !END - error if not inside another block
         call error_handler(E_ERR, 'state_change', &
                   'END keyword found while not inside a GROUP or IMPACT block', &
                   source, text2=errline, text3=readbuf)

      else
         ! if we are already outside a block, there should not be any other text
         call error_handler(E_ERR, 'state_change', &
                   'text found while not inside a GROUP or IMPACT block', &
                   source, text2=errline, text3=readbuf)
      endif

   case (STATE_DEFGROUP, STATE_DEFNOTGROUP)
      if (is_keyword(itemlist(1), 'GROUP', toc)) then
         !GROUP while already in a group defn - error
         call error_handler(E_ERR, 'state_change', 'GROUP keyword found inside a GROUP block', &
                            source, text2=errline, text3=readbuf)
      else if (is_keyword(itemlist(1), 'IMPACT', toc)) then
         !IMPACT while already in a group defn - error
         call error_handler(E_ERR, 'state_change', 'IMPACT keyword found inside a GROUP block', &
                            source, text2=errline, text3=readbuf)
      else if (is_keyword(itemlist(1), 'END', toc)) then
         !END - change state if word 2 matches, else error
         if (wordcount < 2) then
            call error_handler(E_ERR, 'state_change', &
                      'END keyword found but missing second word GROUP', &
                      source, text2=errline, text3=readbuf)
         endif
         if (.not. is_keyword(itemlist(2), 'GROUP', toc)) then
            call error_handler(E_ERR, 'state_change', &
                      'END keyword found but second word not GROUP', &
                      source, text2=errline, text3=readbuf)
         endif

         this_state = STATE_UNDEFINED

      else if ((is_keyword(itemlist(1), 'ALL',      toc)) .or. &
               (is_keyword(itemlist(1), 'ALLQTYS', toc)) .or. &
               (is_keyword(itemlist(1), 'ALLTYPES', toc))) then

         !ALL and friends - add all, and then if word 2 is EXCEPT, take away
         if (this_state == STATE_DEFNOTGROUP) then
            call error_handler(E_ERR, 'state_change', &
                      'ALL keyword found but already excluding items', &
                      source, text2=errline, text3=readbuf)
         endif

         this_category = itemlist(1)

         ! require group be empty to start.
         if (get_group_count(current_group_id, group_toc) /= 0) then
            call error_handler(E_ERR, 'state_change', &
                      'ALL EXCEPT must be first lines inside a group', &
                      source, text2=errline, text3=readbuf)
         endif

         ! add all qtys, types, or types+qtys to group.
         call add_all_to_group(toc, current_group_id, itemlist(1), group_toc)

         if (wordcount >= 2) then
            if (.not. is_keyword(itemlist(2), 'EXCEPT', toc)) then
               call error_handler(E_ERR, 'state_change', &
                         'ALL keyword found but second word not EXCEPT', &
                         source, text2=errline, text3=readbuf)
            endif

!print *, 'changing state to defnotgroup'
            ! finish any qtys/groups on rest of this line.
            if (wordcount > 2) &
               call sub_from_group(wordcount-2, itemlist(3:wordcount), current_group_id, itemlist(1), group_toc, toc)

         endif

         ! subsequent lines remove qtys and groups from group.
         this_state = STATE_DEFNOTGROUP

      else
         ! anything else is fine; no state change
         otherwork = .true.

      endif

   case (STATE_DEFIMPACT)
      if (is_keyword(itemlist(1), 'GROUP', toc)) then
            !GROUP while already in a impact defn - error
            call error_handler(E_ERR, 'state_change', 'GROUP keyword found inside an IMPACT block', &
                               source, text2=errline, text3=readbuf)
      else if (is_keyword(itemlist(1), 'IMPACT', toc)) then
            !IMPACT while already in a impact defn - error
            call error_handler(E_ERR, 'state_change', 'IMPACT keyword found inside an IMPACT block', &
                                source, text2=errline, text3=readbuf)
      else if (is_keyword(itemlist(1), 'END', toc)) then
            !END - change state if word 2 matches, else error
            if (wordcount < 2) then
               call error_handler(E_ERR, 'state_change', &
                         'END keyword found but missing second word IMPACT', &
                         source, text2=errline, text3=readbuf)
            endif
            if (.not. is_keyword(itemlist(2), 'IMPACT', toc)) then
               call error_handler(E_ERR, 'state_change', &
                         'END keyword found but second word not IMPACT', &
                         source, text2=errline, text3=readbuf)
            endif

            this_state = STATE_UNDEFINED

      else
            ! anything else is fine; no state change
            otherwork = .true.

      endif

  case default
    call error_handler(E_ERR,'state_change:','internal error: bad state value',source)
end select

end subroutine state_change

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

function has_nextword(wstring, nextword)

character(len=*), intent(in)    :: wstring
character(len=*), intent(out)   :: nextword
logical                         :: has_nextword

! assume there is a next word on the line
has_nextword = .true.

! upper case the string
nextword = wstring
call to_upper(nextword)

! return false if there is no next word because of comment char
if (is_comment(nextword)) has_nextword = .false.

end function has_nextword

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

function is_comment(nexttoken)

character(len=*), intent(in) :: nexttoken
logical                      :: is_comment

if (nexttoken(1:1) == commentchar) then
   is_comment = .true.
else
   is_comment = .false.
endif

end function is_comment

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

subroutine add_to_group(wordcount, itemlist, groupid, group_toc, toc)

integer,          intent(in)    :: wordcount
integer,          intent(in)    :: itemlist(:)
integer,          intent(in)    :: groupid
type(type_group), intent(inout) :: group_toc
type(type_toc),   intent(inout) :: toc

integer :: i, j, num_mem

call check_groupindex(groupid, 'add_to_group', group_toc)

do i = 1, wordcount
   if (itemlist(i) < 0) then
      call error_handler(E_ERR, 'add_to_group', &
                         'unrecognized item cannot be added to group', &
                         source, text2=errline, text3=readbuf)
   endif

   ! check for dups
   do j=1, group_toc%groups(groupid)%member_count
      if (group_toc%groups(groupid)%group_members(j) == itemlist(i)) then
         call error_handler(E_ERR, 'add_to_group', &
                            'duplicate member found; '//trim(get_toc_name(itemlist(i),toc))//' already a member of group', &
                            source, text2=errline, text3=readbuf)
      endif
   enddo

   ! make sure they are qtys or types, do not allow nested groups
   if (get_toc_type(itemlist(i), toc) == ENTRY_GROUP) then
      call error_handler(E_ERR, 'add_to_group', &
                         'group members must be DART QTYS or TYPES (nested groups not allowed)', &
                         source, text2=errline, text3=readbuf)
   endif

   num_mem = group_toc%groups(groupid)%member_count + 1

   group_toc%groups(groupid)%member_count = num_mem
   group_toc%groups(groupid)%group_members(num_mem) = itemlist(i)

enddo

end subroutine add_to_group

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

subroutine add_all_to_group(toc, groupid, category, group_toc)

type(type_toc),   intent(in)    :: toc
integer,          intent(in)    :: groupid
integer,          intent(in)    :: category
type(type_group), intent(inout) :: group_toc

integer :: i, itemtype, num_mem
logical :: doqtys, dotypes

doqtys = .false.
dotypes = .false.

if (is_keyword(category, 'ALL', toc)) then
   doqtys = .true.
   dotypes = .true.
else if (is_keyword(category, 'ALLQTYS',  toc)) then
   doqtys = .true.
else if (is_keyword(category, 'ALLTYPES', toc)) then
   dotypes = .true.
else
   call error_handler(E_ERR, 'add_all_to_group', &
                      'internal error: ALLxxx keyword unrecognized', source)
endif

call check_groupindex(groupid, 'add_all_to_group', group_toc)

! for all qtys and/or types in toc, add them as requested
toc_loop: do i = 1, toc%toc_count

   itemtype = get_toc_type(i, toc)
   if ((itemtype == ENTRY_DARTQTY  .and. doqtys) .or. &
       (itemtype == ENTRY_DARTTYPE .and. dotypes)) then

      num_mem = group_toc%groups(groupid)%member_count + 1

      !print *, 'adding index i, with name,val ', trim(get_toc_name(i, toc)), get_toc_value(i, toc)
      group_toc%groups(groupid)%member_count = num_mem
      group_toc%groups(groupid)%group_members(num_mem) = i
   endif

enddo toc_loop

end subroutine add_all_to_group

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

subroutine sub_from_group(wordcount, itemlist, groupid, category, group_toc, toc)

integer,          intent(in)    :: wordcount
integer,          intent(in)    :: itemlist(:)
integer,          intent(in)    :: groupid
integer,          intent(in)    :: category
type(type_group), intent(inout) :: group_toc
type(type_toc),   intent(inout) :: toc

integer :: i, j, itemtype
logical :: isgroup
integer :: itemval, itemcount, subgroupid
logical :: doqtys, dotypes

doqtys = .false.
dotypes = .false.

if (is_keyword(category, 'ALL', toc)) then
   doqtys = .true.
   dotypes = .true.
else if (is_keyword(category, 'ALLQTYS',  toc)) then
   doqtys = .true.
else if (is_keyword(category, 'ALLTYPES', toc)) then
   dotypes = .true.
else
   call error_handler(E_ERR, 'sub_from_group', &
             'internal error: ALLxxx keyword unrecognized', source)
endif

call check_groupindex(groupid, 'sub_from_group', group_toc)

do i = 1, wordcount
   if (itemlist(i) < 0) then
      call error_handler(E_ERR, 'sub_from_group', &
                'EXCEPT items must be known QTYS, TYPES, or groups', &
                source, text2=errline, text3=readbuf)

   endif

   ! see if itemlist(i) is an item or a group.  if item, remove it.
   ! if group, iterate through the group and remove the member items.
   ! for now, do NOT allow nested groups, although at some point
   ! i guess we could start recursing.  but it makes the
   ! error checking a nightmare, and hard to detect logical loops.
   ! for now, just say no.

   call iteminfo(itemlist(i), itemval, isgroup, itemcount, toc, group_toc, qtyonly = .false.)

   if (isgroup) subgroupid = itemval

   do j=1, itemcount

      if (isgroup) then

         itemval = get_group_member_by_index(subgroupid, j, group_toc)

         itemtype = get_toc_type(itemval, toc)
         if (.not. (itemtype == ENTRY_DARTQTY  .and. doqtys) .and. &
             .not. (itemtype == ENTRY_DARTTYPE .and. dotypes)) then
            call error_handler(E_ERR, 'sub_from_group', &
                      'error using EXCEPT; item '//trim(get_toc_name(itemval, toc))//' not a group member', &
                      source, text2=errline, text3=readbuf)
         endif

      else
         itemval = itemlist(i)
         if (.not. is_group_member(groupid, itemval, group_toc)) then
            call error_handler(E_ERR, 'sub_from_group', &
                      'error using EXCEPT; item '//trim(get_toc_name(itemval, toc))//' not a member of group', &
                      source, text2=errline, text3=readbuf)
         endif
      endif

      call del_group_member(groupid, itemval, group_toc)

   enddo

enddo

end subroutine sub_from_group

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

function is_group_member(groupid, itemnum, group_toc)

integer,          intent(in) :: groupid
integer,          intent(in) :: itemnum
type(type_group), intent(in) :: group_toc
logical                      :: is_group_member

is_group_member = .false.

call check_groupindex(groupid, 'is_group_member', group_toc)

if (group_toc%groups(groupid)%member_count < 0) then
   call error_handler(E_ERR, 'is_group_member', &
                     'no members in requested group', source)
endif

do i=1, group_toc%groups(groupid)%member_count
   if (group_toc%groups(groupid)%group_members(i) == itemnum) then
      is_group_member = .true.
      return
   endif
enddo

end function is_group_member

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

function get_group_member_by_index(groupid, membernum, group_toc)

integer,          intent(in) :: groupid
integer,          intent(in) :: membernum
type(type_group), intent(in) :: group_toc
integer                      :: get_group_member_by_index

call check_groupindex(groupid, 'get_group_member_by_index', group_toc)

if (group_toc%groups(groupid)%member_count <= 0) then
   call error_handler(E_ERR, 'get_group_member_by_index', &
             'no members in requested group', source)
endif

if (membernum > group_toc%groups(groupid)%member_count) then
   call error_handler(E_ERR, 'get_group_member_by_index', &
             'membernum larger than number of members in this group', source)
endif

get_group_member_by_index = group_toc%groups(groupid)%group_members(membernum)

end function get_group_member_by_index

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

subroutine del_group_member(groupid, itemnum, group_toc)

integer,          intent(in)    :: groupid
integer,          intent(in)    :: itemnum
type(type_group), intent(inout) :: group_toc

integer :: i, membernum
integer :: num_mem

call check_groupindex(groupid, 'del_group_member', group_toc)

membernum = -1
FINDME: do i=1, group_toc%groups(groupid)%member_count
   if (group_toc%groups(groupid)%group_members(i) == itemnum) then
      membernum = i
      exit FINDME
   endif
enddo FINDME

if (membernum < 0 .or. membernum > group_toc%groups(groupid)%member_count) then
   call error_handler(E_ERR, 'del_group_member', &
             'membernum not a valid member number', source)
endif

num_mem = group_toc%groups(groupid)%member_count - 1

do i=membernum, num_mem
   group_toc%groups(groupid)%group_members(i) = group_toc%groups(groupid)%group_members(i+1)
enddo

group_toc%groups(groupid)%member_count = num_mem

end subroutine del_group_member

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

function get_group_count(groupid, group_toc)

integer,          intent(in) :: groupid
type(type_group), intent(in) :: group_toc
integer                      :: get_group_count

call check_groupindex(groupid, 'get_group_count', group_toc)

if (group_toc%groups(groupid)%member_count < 0) then
   call error_handler(E_ERR, 'get_group_count', &
             'no members in requested group', source)
endif

get_group_count = group_toc%groups(groupid)%member_count

end function get_group_count

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

subroutine print_group_contents(groupid, group_toc)

integer,          intent(in) :: groupid
type(type_group), intent(in) :: group_toc

integer :: i, num

call check_groupindex(groupid, 'print_group_contents', group_toc)

num = group_toc%groups(groupid)%member_count

write(*,*) 'group id number ', groupid
write(*,*) 'member count ', num

do i=1, num
   write(*, '(A14,I5,I5)') 'member, id ', i, group_toc%groups(groupid)%group_members(i)
enddo

end subroutine print_group_contents

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

subroutine check_groupindex(groupid, caller, group_toc)

integer,          intent(in) :: groupid
character(len=*), intent(in) :: caller
type(type_group), intent(in) :: group_toc

if (groupid <= 0) then
   call error_handler(E_ERR, caller, &
             'internal error: group index must be 1 or larger', source)
endif

if (groupid > group_toc%group_count) then
   call error_handler(E_ERR, caller, &
             'internal error: group index larger than number of groups', source)
endif

end subroutine check_groupindex

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

subroutine find_all_types_for_qty(givenqty, typecount, typesfound, typeidlist)

integer, intent(in)  :: givenqty
integer, intent(in)  :: typecount
integer, intent(out) :: typesfound
integer, intent(out) :: typeidlist(:)

integer :: i, qtyindex

! given a generic qty, find all the specific types which
! map to this qty.  i believe this has to be done by iterating
! through every type to see what the qty is.

! this routine is dealing with the raw indices into the qtys
! and types array - it no longer is using the table of contents values.

typesfound = 0
typeidlist(:) = -1

do i=1, typecount
   qtyindex = get_quantity_for_type_of_obs(i)
   if (qtyindex == givenqty) then
      typesfound = typesfound + 1
      typeidlist(typesfound) = i
   endif
enddo

end subroutine find_all_types_for_qty

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

subroutine update_impact(table, wordcount, itemlist, wordarray, toc, group_toc)

real(r8),         intent(inout) :: table(:,0:)
integer,          intent(in)    :: wordcount
integer,          intent(in)    :: itemlist(:)
character(len=*), intent(in)    :: wordarray(:)
type(type_toc),   intent(in)    :: toc
type(type_group), intent(in)    :: group_toc

real(r8) :: rvalue
integer :: itemval1, itemval2, groupid1, groupid2
integer :: groupitem1, groupitem2
integer :: itemcount1, itemcount2
logical :: isgroup1, isgroup2
integer :: i, j

! expect exactly 3 tokens on these lines:
!  qty_or_type_or_group   qty_or_group   real(r8)
! (where 2nd group cannot contain types)

call check_impact_line(wordcount, itemlist, wordarray, rvalue)

! ok, we have the first, second, and value.  update the table.
! any combination of qty names or group names is allowed.

call iteminfo(itemlist(1), itemval1, isgroup1, itemcount1, toc, group_toc, qtyonly=.false.)
call iteminfo(itemlist(2), itemval2, isgroup2, itemcount2, toc, group_toc, qtyonly=.true.)

! both single items
if (.not. isgroup1 .and. .not. isgroup2) then

   call fill_impact_table_values(table, itemval1, itemval2, rvalue, toc, itemlist(1), 'k/k')

! single item to group
else if (.not. isgroup1 .and. isgroup2) then
   groupid2 = itemval2
   do i=1, itemcount2
      groupitem2 = get_group_member_by_index(groupid2, i, group_toc)
      call require_only_qtys_in_group(groupitem2, groupid2, toc)
      itemval2 = get_toc_value(groupitem2, toc)
      call fill_impact_table_values(table, itemval1, itemval2, rvalue, toc, itemlist(1), 'k/g')
   enddo

! group to single item
else if (isgroup1 .and. .not. isgroup2) then
   groupid1 = itemval1
   itemval2 = get_toc_value(itemlist(2), toc)
   do i=1, itemcount1
      groupitem1 = get_group_member_by_index(groupid1, i, group_toc)
      itemval1 = get_toc_value(groupitem1, toc)
      call fill_impact_table_values(table, itemval1, itemval2, rvalue, toc, groupitem1, 'g/k')
   enddo

! both groups
else if (isgroup1 .and. isgroup2) then
   groupid1 = itemval1
   groupid2 = itemval2
   do i=1, itemcount1
      groupitem1 = get_group_member_by_index(groupid1, i, group_toc)
      itemval1 = get_toc_value(groupitem1, toc)
      do j=1, itemcount2
         groupitem2 = get_group_member_by_index(groupid2, j, group_toc)
         call require_only_qtys_in_group(groupitem2, groupid2, toc)
         itemval2 = get_toc_value(groupitem2, toc)
         call fill_impact_table_values(table, itemval1, itemval2, rvalue, toc, groupitem1, 'g/g')
      enddo
   enddo

else
   call error_handler(E_ERR, 'update_impact', &
             'first and second entries on line must be either a QTY or a group name', &
             source, text2=errline, text3=readbuf)
endif

end subroutine update_impact

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

subroutine check_impact_line(wordcount, itemlist, wordarray, rvalue)

integer,          intent(in)  :: wordcount
integer,          intent(in)  :: itemlist(:)
character(len=*), intent(in)  :: wordarray(:)
real(r8),         intent(out) :: rvalue

if (wordcount /= 3) then
   call error_handler(E_ERR, 'check_impact_line', &
                     'impact lines require exactly 3 words per line', &
                     source, text2=errline, text3=readbuf)
endif
if (itemlist(1) < 0) then
   call error_handler(E_ERR, 'check_impact_line', &
                     'first word on line is unrecognized QTY, TYPE, or group name', &
                     source, text2=errline, text3=readbuf)
endif
if (itemlist(2) < 0) then
   call error_handler(E_ERR, 'check_impact_line', &
                     'second word on line is unrecognized QTY or group name', &
                     source, text2=errline, text3=readbuf)
endif

call extract_value(wordarray(3), rvalue)

end subroutine check_impact_line

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

! TOOL:
subroutine fill_impact_table_values(table, tableindex1, tableindex2, rvalue, toc, tocindex1, label)

real(r8),         intent(inout)        :: table(:,0:)
integer,          intent(in)           :: tableindex1
integer,          intent(in)           :: tableindex2
real(r8),         intent(in)           :: rvalue
type(type_toc),   intent(in)           :: toc
integer,          intent(in)           :: tocindex1
character(len=*), intent(in), optional :: label

integer :: i, typeid(toc%type_count), typesfound

! expand tableindex1 if type here and loop over list
if (get_toc_type(tocindex1, toc) == ENTRY_DARTQTY) then

   ! this routine is now dealing with the absolute indices into the types
   ! and qtys array; not the indices in the table of contents.  the values
   ! returned can be used directly to set the table entries.

   call find_all_types_for_qty(tableindex1, toc%type_count, typesfound, typeid)

   do i=1, typesfound
      call set_impact_table_values(table, typeid(i), tableindex2, rvalue, toc, label)
   enddo

else

   call set_impact_table_values(table, tableindex1, tableindex2, rvalue, toc, label)

endif

end subroutine fill_impact_table_values

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

! TOOL:
subroutine set_impact_table_values(table, tableindex1, tableindex2, rvalue, toc, label)

real(r8),         intent(inout)        :: table(:,0:)
integer,          intent(in)           :: tableindex1
integer,          intent(in)           :: tableindex2
real(r8),         intent(in)           :: rvalue
type(type_toc),   intent(in)           :: toc
character(len=*), intent(in), optional :: label

if (present(label) .and. debug) then
   print *, trim(label), ' set_impact_table_values: ', tableindex1, tableindex2
endif

if ((table(tableindex1, tableindex2) /= missing_r8) .and. &
    (table(tableindex1, tableindex2) /= rvalue)) then
   write(msgstring,  *) 'error; different impact factor already set for this pair'
   write(msgstring2, *) 'previous value was', table(tableindex1, tableindex2), ' new value is ', rvalue
   write(msgstring3, *) trim(get_name_by_value(tableindex1, ENTRY_DARTTYPE, toc)), ' and ',  &
                        trim(get_name_by_value(tableindex2, ENTRY_DARTQTY,  toc))
   call error_handler(E_ERR, 'update_impact', msgstring, &
                     source, text2=msgstring2, text3=msgstring3)

endif

table(tableindex1, tableindex2) = rvalue

end subroutine set_impact_table_values

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

subroutine iteminfo(item, itemid, isgroup, itemcount, toc, group_toc, qtyonly)

integer,          intent(in)  :: item
integer,          intent(out) :: itemid
logical,          intent(out) :: isgroup
integer,          intent(out) :: itemcount
type(type_toc),   intent(in)  :: toc
type(type_group), intent(in)  :: group_toc
logical,          intent(in)  :: qtyonly

integer :: entryinfo

entryinfo = get_toc_type(item, toc)

! our usage of this is asymmetric, but it is not clear how
! much the tool should know about this.
if(entryinfo == ENTRY_GROUP) then
   itemid = get_toc_value(item, toc)
   isgroup = .true.
   itemcount = get_group_count(itemid, group_toc)
   ! FIXME: same for group members for item 2 - not symmetric? check here?

else if (qtyonly .and. entryinfo == ENTRY_DARTTYPE) then
   call error_handler(E_ERR, 'iteminfo', &
             'entry on right, '//trim(get_toc_name(item,toc))//' can only be a DART QTY', &
             source, text2=errline, text3=readbuf)

else if (qtyonly .and. entryinfo == ENTRY_DARTQTY) then
   itemid = get_toc_value(item, toc)
   isgroup = .false.
   itemcount = 1

else if (.not. qtyonly .and. &
        (entryinfo == ENTRY_DARTQTY .or. entryinfo == ENTRY_DARTTYPE)) then
   itemid = get_toc_value(item, toc)
   isgroup = .false.
   itemcount = 1

else

   call error_handler(E_ERR, 'iteminfo', &
             'unrecognized input; expecting a DART QTY, TYPE, or GROUP name', &
             source, text2=errline, text3=readbuf)
endif

end subroutine iteminfo

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

subroutine require_only_qtys_in_group(item, groupid, toc)

integer,        intent(in) :: item
integer,        intent(in) :: groupid
type(type_toc), intent(in) :: toc

character(len=string_length) :: gname

if (get_toc_type(item, toc) /= ENTRY_DARTQTY) then

   gname = get_name_by_value(groupid, ENTRY_GROUP, toc)
   call error_handler(E_ERR, 'require_only_qtys_in_group', &
             'if in righthand column, group '//trim(gname)//' must include only DART QTYS', &
             source, text2=errline, text3=readbuf)
endif

end subroutine require_only_qtys_in_group

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

subroutine must_be_specific_type(item, toc)

integer,        intent(in) :: item
type(type_toc), intent(in) :: toc

character(len=string_length) :: ename

if (get_toc_type(item, toc) /= ENTRY_DARTTYPE) then

   ename = get_toc_name(item, toc)
   call error_handler(E_ERR, 'must_be_specific_type', &
                     'lefthand column must be a DART TYPE, not '//trim(ename), &
                     source, text2=errline, text3=readbuf)
endif

end subroutine must_be_specific_type

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

subroutine must_be_generic_qty(item, toc)

integer,        intent(in) :: item
type(type_toc), intent(in) :: toc

character(len=string_length) :: ename

if (get_toc_type(item, toc) /= ENTRY_DARTQTY) then

   ename = get_toc_name(item, toc)
   call error_handler(E_ERR, 'must_be_generic_qty', &
                     'righthand column must be a DART QTY, not '//trim(ename), &
                     source, text2=errline, text3=readbuf)
endif

end subroutine must_be_generic_qty

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

function get_toc_type(item, toc)

integer,        intent(in) :: item
type(type_toc), intent(in) :: toc
integer                    :: get_toc_type

if (item > toc%toc_count) then
   call error_handler(E_ERR, 'get_toc_type', &
             'called with an index larger than the number of items in toc', source)
endif

get_toc_type = toc%toc_entries(item)%etype

end function get_toc_type

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

function get_toc_value(item, toc)

integer,        intent(in) :: item
type(type_toc), intent(in) :: toc
integer                    :: get_toc_value

if (item > toc%toc_count) then
   call error_handler(E_ERR, 'get_toc_value', &
             'called with an index larger than the number of items in toc', source)
endif

get_toc_value = toc%toc_entries(item)%evalue

end function get_toc_value

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

function get_toc_name(item, toc)

integer,        intent(in)   :: item
type(type_toc), intent(in)   :: toc
character(len=string_length) :: get_toc_name

if (item > toc%toc_count) then
   call error_handler(E_ERR, 'get_toc_name', &
             'called with an index larger than the number of items in toc', source)
endif

get_toc_name = toc%toc_entries(item)%ename

end function get_toc_name

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

subroutine extract_value(wordarray, rvalue)

character(len=*), intent(in)  :: wordarray
real(r8),         intent(out) :: rvalue

! how to validate this?  look for chars other than 0-9 . - e d ?
read(wordarray, *, iostat=ios) rvalue
if (ios /= 0) then
   call error_handler(E_ERR, 'extract_value', &
             'bad numeric value, expecting third item on line to be a number', &
             source, text2=errline, text3=readbuf)
endif

end subroutine extract_value

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

subroutine output_impact_table(oname, table, toc)

character(len=*), intent(in) :: oname
real(r8),         intent(in) :: table(:,0:)
type(type_toc),   intent(in) :: toc

integer :: i, j, funit
integer :: ntypes, nqtys

funit = open_file(oname, 'formatted', action='write')

ntypes = ubound(table, 1)
nqtys = ubound(table, 2)

do j=0, nqtys
   do i=1, ntypes
      if (table(i, j) /= 1.0_r8) then
         write(funit, '(A34,A34,F18.6)') &
                        trim(get_name_by_value(i, ENTRY_DARTTYPE, toc)), &
                        trim(get_name_by_value(j, ENTRY_DARTQTY,  toc)), &
                        table(i, j)
      endif
   enddo
enddo

call close_file(funit)

end subroutine output_impact_table

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

end module obs_impact_mod


text:
program ieee_test
use iso_fortran_env
use ieee_arithmetic
real(kind=REAL64) :: x
if (ieee_support_datatype(x)) then
    stop 1
else
    stop 0
endif
end program ieee_test

max_stars_repo_name: dmitrygx/HPCInfo

max_issues_repo_name: dmitrygx/HPCInfo

max_forks_repo_name: dmitrygx/HPCInfo