text:
!
! SRF_Utility
!
! Module containing SRF application routines
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 18-Mar-2003
!                       paul.vandelst@ssec.wisc.edu
!
!       Modified by:    Yong Chen, CIRA/CSU/JCSDA 22-Aug-2008
!                       Yong.Chen@noaa.gov
!

MODULE SRF_Utility

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds,            ONLY: fp
  USE Message_Handler,       ONLY: SUCCESS, FAILURE, WARNING, Display_Message
  USE Interpolate_Utility,   ONLY: Polynomial_Interpolate
  USE Integrate_Utility,     ONLY: Integral
  USE SRF_Define,            ONLY: SRF_type, &
                                   Associated_SRF, &
                                   Destroy_SRF, &
                                   Allocate_SRF, &
                                   Assign_SRF, &
                                   Frequency_SRF, &
                                   Integrate_SRF
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: Interpolate_SRF
  PUBLIC :: Convolve_with_SRF
  PUBLIC :: Calc_SRF_First_Mom

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_RCS_ID = &
  ! Maximum msg length
  INTEGER, PARAMETER :: ML = 256
  ! Keyword set value
  INTEGER, PARAMETER :: SET = 1
  ! Literal constants
  REAL(fp), PARAMETER :: ZERO = 0.0_fp
  REAL(fp), PARAMETER :: ONE  = 1.0_fp
  REAL(fp), PARAMETER :: ONEpointFIVE  = 1.5_fp

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                         ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                            ##
!################################################################################
!################################################################################

!------------------------------------------------------------------------------
!:sdoc:+
!
! NAME:
!       Interpolate_SRF
!
! PURPOSE:
!       Function to interpolate an input SRF to another frequency grid.
!
! CALLING SEQUENCE:
!       Error_STatus = Interpolate_SRF( SRF                    , &  ! Input
!                                       iSRF                   , &  ! In/Output
!                                       Order      =Order      , &  ! Optional input
!                                       RCS_Id     =RCS_Id     , &  ! Revision control
!                                       Message_Log=Message_Log  )  ! Error messaging
!
! INPUT ARGUMENTS:
!       SRF:          SRF structure containing the instrument channel
!                     response to be interpolated.
!                     UNITS:      N/A
!                     TYPE:       SRF_type
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN)
!
!       iSRF:         On input contains the components used to generate the
!                     interpolation frequency grid:
!                       iSRF%f1_Band(:)
!                       iSRF%f2_Band(:)
!                       iSRF%npts_Band(:)
!                     and the interpolation frequency grid itself:
!                       iSRF%Frequency(:)
!                     UNITS:      N/A
!                     TYPE:       SRF_type
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUT ARGUMENTS:
!       iSRF:         On output, contains the interpolated SRF response
!                     in the component,
!                       iSRF%Response(:)
!                     UNITS:      N/A
!                     TYPE:       SRF_type
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN OUT)
!
! OPTIONAL INPUT ARGUMENTS:
!       Order:        The order of the interpolating polynomial. This
!                     argument, if supplied, is passed to the interpolation
!                     function where, if not specified, linear interpolation
!                     (order = 1) is the default. If specified, must be an
!                     odd number > 0.
!                     UNITS:      None
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Message_Log:  Character string specifying a filename in which any
!                     msgs will be logged. If not specified, or if an
!                     error occurs opening the log file, the default action
!                     is to output msgs to standard output.
!                     UNITS:      None
!                     TYPE:       Character
!                     DIMENSION:  Scalar, LEN = *
!                     ATTRIBUTES: INTENT(IN), OPTIONAL
!
! OPTIONAL OUTPUT ARGUMENTS:
!       RCS_Id:       Character string containing the Revision Control
!                     System Id field for the module.
!                     UNITS:      None
!                     TYPE:       CHARACTER
!                     DIMENSION:  Scalar, LEN = *
!                     ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status: The return value is an integer defining the error status.
!                     The error codes are defined in the Message_Handler module.
!                     If == SUCCESS the SRF integration was successful
!                        == FAILURE an error occurred
!                     UNITS:      N/A
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!
! COMMENTS:
!       Note the INTENT on the output SRF argument is IN OUT rather
!       than just OUT. This is necessary because the argument may be defined on
!       input. To prevent memory leaks, the IN OUT INTENT is a must.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 18-Mar-2003
!                       paul.vandelst@ssec.wisc.edu
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION Interpolate_SRF( SRF         , &  ! Input
                            iSRF        , &  ! Output
                            Order       , &  ! Optional input
                            RCS_Id      , &  ! Revision control
                            Message_Log ) &  ! Error messaging
                          RESULT( Error_Status )
    ! Arguments
    TYPE(SRF_type),         INTENT(IN)     :: SRF
    TYPE(SRF_type),         INTENT(IN OUT) :: iSRF
    INTEGER,      OPTIONAL, INTENT(IN)     :: Order
    CHARACTER(*), OPTIONAL, INTENT(OUT)    :: RCS_Id
    CHARACTER(*), OPTIONAL, INTENT(IN)     :: Message_Log
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Interpolate_SRF'
    REAL(fp),     PARAMETER :: EPS = EPSILON(ONE)
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: i1, i2, j1, j2, m, ni, nj
    REAL(fp) :: df

    ! Setup
    ! -----
    Error_Status = SUCCESS
    IF ( PRESENT(RCS_Id) ) RCS_Id = MODULE_RCS_ID

    ! Check structure association
    IF ( .NOT. Associated_SRF( SRF ) ) THEN
      Error_Status = FAILURE
      msg = 'Some or all INPUT SRF pointer members are NOT associated.'
      CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
      RETURN
    END IF

    IF ( .NOT. Associated_SRF( iSRF ) ) THEN
      Error_Status = FAILURE
      msg = 'Some or all OUTPUT iSRF pointer members are NOT associated.'
      CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
      RETURN
    END IF

    ! Set the output SRF info
    iSRF%Sensor_Id        = SRF%Sensor_Id
    iSRF%WMO_Satellite_ID = SRF%WMO_Satellite_ID
    iSRF%WMO_Sensor_ID    = SRF%WMO_Sensor_ID
    iSRF%Sensor_Type      = SRF%Sensor_Type
    iSRF%Channel          = SRF%Channel

    ! Begin loop over bands for interpolation
    ! ---------------------------------------
    i2 = 0; j2 = 0
    Band_Loop: DO m = 1, SRF%n_Bands

      ! Determine where the old SRF slots into its frequency grid for this band
      i1 = i2 + 1
      i2 = SRF%npts_Band(m) + i1 - 1

      ! Determine where the new SRF slots into its frequency grid for this band
      j1 = j2 + 1
      j2 = iSRF%npts_Band(m) + j1 - 1

      ! Perform the interpolation
      Error_Status = Polynomial_Interpolate(  SRF%Frequency(i1:i2)  , &  ! Input,  X
                                              SRF%Response(i1:i2)   , &  ! Input,  Y
                                             iSRF%Frequency(j1:j2)  , &  ! Input,  Xint
                                             iSRF%Response(j1:j2)   , &  ! Output, Yint
                                             Order      =Order      , &  ! Optional input
                                             Message_Log=Message_Log  )
      IF ( Error_Status /= SUCCESS ) THEN
        WRITE( msg,'("Error interpolating band ",i0," SRF response")' ) m
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF

    END DO Band_Loop

    ! Integrate the interpolated SRF
    ! ------------------------------
    Error_Status = Integrate_SRF( iSRF,Message_Log=Message_Log )
    IF ( Error_Status /= SUCCESS ) THEN
      CALL Display_Message( ROUTINE_NAME, &
                            'Error occurred integrating interpolated SRF', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

  END FUNCTION Interpolate_SRF

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Convolve_with_SRF
!
! PURPOSE:
!       Procedure to convolve a spectrum with an instrument channel SRF.
!
! CALLING SEQUENCE:
!       Error_Status = Convolve_with_SRF( Frequency              , &  ! Input
!                                         Spectrum               , &  ! Input
!                                         SRF                    , &  ! Input
!                                         Convolved_Spectrum     , &  ! Output
!                                         Interpolate=Interpolate, &  ! Optional input
!                                         Integrate  =Integrate  , &  ! Optional input
!                                         Order      =Order      , &  ! Optional input
!                                         RCS_Id     =RCS_Id     , &  ! Revision control
!                                         Message_Log=Message_Log  )  ! Error messaging
!
! INPUT ARGUMENTS:
!       Frequency:           Frequency grid of the input spectrum.
!                            Must be evenly spaced in frequency.
!                            UNITS:      Variable (cm^-1 or GHz)
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Rank-1
!                            ATTRIBUTES: INTENT(IN)
!
!       Spectrum:            Spectrum to be convolved with the SRF response.
!                            UNITS:      Variable.
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Rank-1
!                            ATTRIBUTES: INTENT(IN)
!
!       SRF:                 SRF structure containing the instrument channel
!                            response to be used in the convolution.
!                            UNITS:      N/A
!                            TYPE:       SRF_type
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       Convolved_Spectrum:  The input spectrum convolved with the input SRF
!                            response.
!                            UNITS:      Same as input Spectrum argument
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUT ARGUMENTS:
!       Interpolate:         Set this argument to force interpolation of the SRF
!                            response to the input frequency grid. If not set,
!                            the default is to assume that the SRF and spectrum
!                            frequency grids match up for the required span.
!                            If = 0; do not interpolate the SRF response (DEFAULT)
!                               = 1; interpolate input SRF response to the input
!                                    frequency grid.
!                            UNITS:      N/A
!                            TYPE:       INTEGER
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Integrate:           Set this argument to perform numerical integration
!                            when doing the convolutions. If not set, the default
!                            action is to simply sum the results.
!                            If = 0; use summation to compute convolved value. (DEFAULT)
!                               = 1; perform numerical integration to compute the
!                                    convolved value.
!                            UNITS:      N/A
!                            TYPE:       INTEGER
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Order:               The order of the interpolating polynomial. This
!                            argument, if supplied, is passed to the interpolation
!                            function where, if not specified, linear interpolation
!                            (order = 1) is the default. If specified, must be an
!                            odd number > 0.
!                            UNITS:      None
!                            TYPE:       INTEGER
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Message_Log:         Character string specifying a filename in which any
!                            msgs will be logged. If not specified, or if an
!                            error occurs opening the log file, the default action
!                            is to output msgs to standard output.
!                            UNITS:      None
!                            TYPE:       Character
!                            DIMENSION:  Scalar, LEN = *
!                            ATTRIBUTES: INTENT(IN), OPTIONAL
!
! OPTIONAL OUTPUT ARGUMENTS:
!       RCS_Id:              Character string containing the Revision Control
!                            System Id field for the module.
!                            UNITS:      None
!                            TYPE:       CHARACTER
!                            DIMENSION:  Scalar, LEN = *
!                            ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:        The return value is an integer defining the error status.
!                            The error codes are defined in the Message_Handler module.
!                            If == SUCCESS the spectrum convolution was successful
!                               == FAILURE an error occurred
!                            UNITS:      N/A
!                            TYPE:       INTEGER
!                            DIMENSION:  Scalar
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 18-Mar-2003
!                       paul.vandelst@ssec.wisc.edu
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION Convolve_with_SRF( Frequency         , &  ! Input
                              Spectrum          , &  ! Input
                              SRF               , &  ! Input
                              Convolved_Spectrum, &  ! Output
                              Interpolate       , &  ! Optional input
                              Integrate         , &  ! Optional input
                              Order             , &  ! Optional input
                              RCS_Id            , &  ! Revision control
                              Message_Log       ) &  ! Error messaging
                            RESULT( Error_Status )
    ! Arguments
    REAL(fp),               INTENT(IN)  :: Frequency(:)
    REAL(fp),               INTENT(IN)  :: Spectrum(:)
    TYPE(SRF_type),         INTENT(IN)  :: SRF
    REAL(fp),               INTENT(OUT) :: Convolved_Spectrum
    INTEGER,      OPTIONAL, INTENT(IN)  :: Interpolate
    INTEGER,      OPTIONAL, INTENT(IN)  :: Integrate
    INTEGER,      OPTIONAL, INTENT(IN)  :: Order
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: RCS_Id
    CHARACTER(*), OPTIONAL, INTENT(IN)  :: Message_Log
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Convolve_with_SRF'
    INTEGER, PARAMETER :: MAX_N_BANDS = 4
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: Destroy_Status
    LOGICAL :: Frequency_Match
    LOGICAL :: Summation
    TYPE(SRF_type) :: uSRF
    INTEGER :: n_Frequencies, n_Points
    INTEGER :: i1, i2, j1, j2, m
    REAL(fp) :: df, cspc
    REAL(fp) :: f1_Band(MAX_N_BANDS), f2_Band(MAX_N_BANDS)
    INTEGER  :: npts_Band(MAX_N_BANDS)

    ! Set up
    ! ------
    Error_Status = SUCCESS
    IF ( PRESENT(RCS_Id) ) RCS_Id = MODULE_RCS_ID

    ! Check for consistent array sizes
    n_Frequencies = SIZE(Frequency)
    IF ( SIZE(Spectrum) /= n_Frequencies ) THEN
      Error_Status = FAILURE
      msg = 'Size of input frequency and spectrum arrays are different'
      CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
      RETURN
    END IF

    ! Default is no interpolation as SRF and spectrum frequency match is assumed...
    Frequency_Match = .TRUE.
    ! ...unless the Interpolate argument is set
    IF ( PRESENT( Interpolate ) ) THEN
      IF ( Interpolate == SET ) Frequency_Match = .FALSE.
    END IF

    ! Default is to sum the convolved SRF and spectrum...
    Summation = .TRUE.
    ! ...unless the Integrate argument is set
    IF ( PRESENT( Integrate ) ) THEN
      IF ( Integrate == SET ) Summation = .FALSE.
    END IF

    ! Compute the input spectrum frequency interval
    dF = Frequency(2) - Frequency(1)

    ! Copy or interpolate SRF as required
    ! -----------------------------------
    IF ( Frequency_Match ) THEN
      ! Copy the input SRF to a local structure
      Error_Status = Assign_SRF( SRF,uSRF,Message_Log=Message_Log )
      IF ( Error_Status /= SUCCESS ) THEN
        msg = 'Error copying SRF'
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF

    ELSE
      ! Find the frequency points that correspond to the input SRF band edge(s)
      DO m = 1, SRF%n_Bands
        j1 = f_Index( SRF%f1_Band(m),Frequency(1),dF,lowf_edge=.true. )
        j2 = f_Index( SRF%f2_Band(m),Frequency(1),dF )
        f1_Band(m)   = Frequency(j1)
        f2_Band(m)   = Frequency(j2)
        npts_Band(m) = j2-j1+1
      END DO
      ! Allocate local SRF structure
      n_Points = SUM(npts_Band(1:SRF%n_Bands))
      Error_Status = Allocate_SRF( n_Points,uSRF,n_Bands=SRF%n_Bands,Message_Log=Message_Log)
      IF ( Error_Status /= SUCCESS ) THEN
        msg = 'Error allocating SRF'
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF
      ! Fill the band components
      uSRF%f1_Band   = f1_Band(1:SRF%n_Bands)
      uSRF%f2_Band   = f2_Band(1:SRF%n_Bands)
      uSRF%npts_Band = npts_Band(1:SRF%n_Bands)
      ! Compute the frequency grid
      Error_Status = Frequency_SRF( uSRF,Message_Log=Message_Log)
      IF ( Error_Status /= SUCCESS ) THEN
        msg = 'Error computing SRF frequency grid'
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF
      ! Interpolate the SRF
      Error_Status = Interpolate_SRF( SRF,uSRF,Order=Order,Message_Log=Message_Log )
      IF ( Error_Status /= SUCCESS ) THEN
        msg = 'Error interpolating SRF'
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF
      ! Integrate the SRF
      Error_Status = Integrate_SRF( uSRF,Message_Log=Message_Log)
      IF ( Error_Status /= SUCCESS ) THEN
        msg = 'Error integrating SRF'
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF
    END IF

    ! Perform the convolution
    ! -----------------------
    ! Initialise the point offsets and sum
    Convolved_Spectrum = ZERO
    i2 = 0
    ! Begin loop over bands
    Band_Loop: DO m = 1, uSRF%n_Bands

      ! Determine where the current band SRF slots into its frequency grid
      i1 = i2 + 1
      i2 = uSRF%npts_Band(m) + i1 - 1

      ! Determine where the current band SRF slots into the spectrum frequency grid
      j1 = f_Index( uSRF%f1_Band(m),Frequency(1),dF )
      j2 = f_Index( uSRF%f2_Band(m),Frequency(1),dF )

      ! Perform the convolution
      IF ( Summation ) THEN
        ! Just sum it all up
        cspc = SUM(uSRF%Response(i1:i2)*Spectrum(j1:j2))*dF
      ELSE
        ! Integrate it
        Error_Status = Integral( Frequency(j1:j2), &
                                 uSRF%Response(i1:i2)*Spectrum(j1:j2), &
                                 cspc )
        IF ( Error_Status /= SUCCESS ) THEN
          WRITE( msg,'("Error integrating SRF*SPC for band ",i0)' ) m
          CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
          RETURN
        END IF
      END IF

      ! Accumulate the result
      Convolved_Spectrum = Convolved_Spectrum + cspc
    END DO Band_Loop

    ! Normalise the result
    IF ( Summation ) THEN
      Convolved_Spectrum = Convolved_Spectrum/uSRF%Summation_SRF
    ELSE
      Convolved_Spectrum = Convolved_Spectrum/uSRF%Integrated_SRF
    END IF

    ! Clean up
    ! --------
    Destroy_Status = Destroy_SRF( uSRF,Message_Log=Message_Log )
    IF ( Destroy_Status /= SUCCESS ) THEN
      CALL Display_Message( ROUTINE_NAME, &
                            'Error destroying local uSRF structure', &
                            WARNING, &
                            Message_Log=Message_Log )
    END IF

  END FUNCTION Convolve_with_SRF

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Calc_SRF_First_Mom
!
! PURPOSE:
!       Procedure to calculate the first moment of an SRF.
!
! CALLING SEQUENCE:
!       Error_Status = Calc_SRF_First_Mom( SRF                    , &  ! Input
!                                          First_Moment           , &  ! Output
!                                          Message_Log=Message_Log  )  ! Error messaging
!
! INPUT ARGUMENTS:
!
!          SRF:              SRF structure containing the instrument channel
!                            response to be used in the convolution.
!                            UNITS:      N/A
!                            TYPE:       SRF_type
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!
!          First_Moment:     Calculated first moment of the SRF.
!                            UNITS:      cm^-1
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUT ARGUMENTS:
!
!  Message_Log:              Character string specifying a filename in which any
!                            msgs will be logged. If not specified, or if an
!                            error occurs opening the log file, the default action
!                            is to output msgs to standard output.
!                            UNITS:      None
!                            TYPE:       Character
!                            DIMENSION:  Scalar, LEN = *
!                            ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!
!       Error_Status:        The return value is an integer defining the error status.
!                            The error codes are defined in the Message_Handler module.
!                            If == SUCCESS the spectrum convolution was successful
!                               == FAILURE an error occurred
!                            UNITS:      N/A
!                            TYPE:       INTEGER
!                            DIMENSION:  Scalar
!
! CREATION HISTORY:
!       Written by:    David Groff, SAIC 8-June-2009
!                      david.groff@noaa.gov
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION Calc_SRF_First_Mom( SRF               , &  ! Input
                               First_Moment      , &  ! Output
                               Message_Log) &  ! Error messaging
                               RESULT( Error_Status )
    ! Arguments
    TYPE(SRF_type),           INTENT(IN OUT) :: SRF
    CHARACTER(*),   OPTIONAL, INTENT(IN)     :: Message_Log
    REAL(fp),                 INTENT(OUT)    :: First_Moment
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Calc_SRF_First_Mom'
    INTEGER :: CUBIC_ORDER = 3
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: i1, i2, m, n
    REAL(fp):: Band_Weight

    ! Setup
    ! -----
    Error_Status = SUCCESS

    ! ALL pointers must be associated
    IF ( .NOT. Associated_SRF( SRF ) ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Some or all INPUT SRF pointer members are NOT associated.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Check the number of bands
    IF ( SRF%n_Bands < 1 ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'SRF structure must contain at least 1 band', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Check the number of points
    IF ( SRF%n_Points < 2 ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'SRF structure must contain at least 2 points.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Check the number of points in each band
    IF ( ANY(SRF%npts_Band < 2) ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'SRF must contain at least 2 points for each band.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Check the total points
    IF ( SUM(SRF%npts_Band) /= SRF%n_Points ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'SRF must have consistent data points.', &
                            Error_Status, &
                            Message_Log=Message_Log )
    END IF

    ! Compute the first moment of the SRF
    ! -----------------------------------
    ! Initialize the first moment
    First_Moment = ZERO
    Band_Weight = ZERO
    ! Initialise the offset counter
    n = 0

    ! Loop over the bands
    DO m = 1, SRF%n_Bands

      ! The point limits for this band
      i1 = n + 1
      i2 = SRF%npts_Band(m) + n

      ! Calculate Weight for the band
      Error_Status = Integral( SRF%Frequency(i1:i2),                      &
                               SRF%Response(i1:i2)*SRF%Frequency(i1:i2),  &
                               Band_Weight                                )
      IF ( Error_Status /= SUCCESS ) THEN
        WRITE( msg,'("Error occurred integrating channel ",i0," SRF")' ) SRF%Channel
        CALL Display_Message( ROUTINE_NAME,TRIM(msg),Error_Status,Message_Log=Message_Log )
        RETURN
      END IF

      ! Accumulate the band weights
      First_Moment = First_Moment + Band_Weight

      ! Update the offset counter
      n = n + SRF%npts_Band(m)
    END DO

    First_Moment = First_Moment/SRF%Integrated_SRF

  END FUNCTION Calc_SRF_First_Mom

!##################################################################################
!##################################################################################
!##                                                                              ##
!##                          ## PRIVATE MODULE ROUTINES ##                       ##
!##                                                                              ##
!##################################################################################
!##################################################################################

!------------------------------------------------------------------------------------
!
! NAME:
!       f_Interval
!
! PURPOSE:
!       Function to compute the interval of a frequency grid from its parameters
!
! CALLING SEQUENCE:
!       df = f_Interval(f1,f2,n)
!
! INPUT ARGUMENTS:
!       f1:  Begin frequency of the grid.
!            UNITS:      Variable (cm^-1 or GHz)
!            TYPE:       REAL(fp)
!            DIMENSION:  Scalar
!            ATTRIBUTES: INTENT(IN)
!
!       f2:  End frequency of the grid.
!            UNITS:      Same as input f1 argument.
!            TYPE:       REAL(fp)
!            DIMENSION:  Scalar
!            ATTRIBUTES: INTENT(IN)
!
!       n:   Number of points in the frequency grid.
!            UNITS:      N/A
!            TYPE:       INTEGER
!            DIMENSION:  Scalar
!            ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       df:  Interval of the frequency grid.
!            UNITS:      Same as input f1 argument.
!            TYPE:       REAL(fp)
!            DIMENSION:  Scalar
!
!------------------------------------------------------------------------------------

  FUNCTION f_Interval(f1,f2,n) RESULT(dF)
    REAL(fp), INTENT(IN) :: f1, f2
    INTEGER , INTENT(IN) :: n
    REAL(fp) :: dF
    dF = (f2-f1)/REAL(n-1,fp)
  END FUNCTION f_Interval

!------------------------------------------------------------------------------------
!
! NAME:
!       f_Index
!
! PURPOSE:
!       Subroutine to determine the index value of a frequency grid that
!       correspond to a supplied frequency values.
!
! CALLING SEQUENCE:
!       idx = f_Index( f, f1_grid, df_grid, lowf_edge )
!
! INPUT ARGUMENTS:
!       f:          Frequency for which an index value is required.
!                   UNITS:      Variable (cm^-1 or GHz)
!                   TYPE:       REAL(fp)
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN)
!
!       f1_grid:    Frequency of the first point of a frequency grid.
!                   Should be <= f.
!                   UNITS:      Same as input f argument.
!                   TYPE:       REAL(fp)
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN)
!
!       df_grid:    Interval of the frequency grid.
!                   UNITS:      Same as input f argument.
!                   TYPE:       REAL(fp)
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN)
!
! INPUT ARGUMENTS:
!       lowf_edge:  Logical flag to indicate what SRF band edge we are at. This
!                   flag affects how the index value is computed.
!                   If .TRUE.  then idx is adjusted so that f1(idx) >= f
!                      .FALSE. then idx is adjusted so that f1(idx) <= f [DEFAULT]
!                   UNITS:      N/A
!                   TYPE:       LOGICAL
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!       idx:        Index of the frequency grid, f1, that most closely corresponds
!                   to the passed SRF band edge frequency.
!                   UNITS:      N/A
!                   TYPE:       INTEGER
!                   DIMENSION:  Scalar
!
!------------------------------------------------------------------------------------

  FUNCTION f_Index(f,f1_grid,df_grid,lowf_edge) RESULT(idx)
    REAL(fp),          INTENT(IN) :: f, f1_grid, df_grid
    LOGICAL, OPTIONAL, INTENT(IN) :: lowf_edge
    INTEGER :: idx
    LOGICAL :: highf_edge
    REAL(fp) :: f1

    ! Assume we are at the high frequency edge of a band...
    highf_edge = .TRUE.
    ! ...unless the low frequency edge keyword is set.
    IF ( PRESENT(lowf_edge) ) THEN
      IF ( lowf_edge ) highf_edge = .FALSE.
    END IF

    ! Compute the index
    idx = INT(ONEpointFIVE + ( (f-f1_grid)/df_grid ))

    ! Adjust the index based on whether we are
    ! at a low or high frequency band edge.
    f1 = f1_grid + REAL(idx-1)*df_grid
    IF ( highf_edge ) THEN
      IF ( f1 > f ) idx = idx-1
    ELSE
      IF ( f1 < f ) idx = idx+1
    END IF

  END FUNCTION f_Index

END MODULE SRF_Utility

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
FUNCTION CSNDZEQ(IOPT,SNDDIA,GPDIASED,TAUR,TAUB,SNDDMX,DEP,SSG,WS)

  ! *** CALCULATES NEAR BED REFERENCE CONCENTRATION REFERENCE HEIGHT (DIMENSIONLESS)
  !
  ! CHANGE RECORD
  ! DATE MODIFIED     BY               DESCRIPTION
  !----------------------------------------------------------------------!
  ! 2011-03           PAUL M. CRAIG    REWRITTEN TO F90
  ! 2011-03-02        PAUL M. CRAIG    FIXED CASE WHEN WS <= 0

  IMPLICIT NONE

  INTEGER :: IOPT
  REAL    :: CSNDZEQ,SNDDIA,GPDIASED,TAUR,TAUB,SNDDMX,DEP,SSG,WS,TMPVAL,REY,TAURS,VAL,VAL1,ZEQ1

  IF( IOPT == 1 )THEN
    ! *** IOPT=1  BASED ON
    ! ***
    ! *** GARCIA, M., AND G. PARKER, 1991: ENTRAINMENT OF BED SEDIMENT
    ! *** INTO SUSPENSION, J. HYDRAULIC ENGINEERING, 117, 414-435.

    CSNDZEQ=0.05

  ELSEIF( IOPT == 2 )THEN
    !
    ! *** IOPT=2  BASED ON
    ! ***
    ! *** SMITH, J. D., AND S. R. MCLEAN, 1977: SPATIALLY AVERAGED FLOW
    ! *** OVER A WAVY SURFACE, J. GEOPHYSICAL RESEARCH, 82, 1735-1746.

    TMPVAL=26.3*SNDDMX*(TAUB-TAUR)/GPDIASED
    TMPVAL=TMPVAL*SNDDIA/SNDDMX
    TMPVAL=TMPVAL/DEP
    CSNDZEQ=MAX(TMPVAL,0.01)

  ELSEIF( IOPT == 3 )THEN
    ! *** IOPT=3  BASED ON
    ! ***
    ! *** VAN RIJN, L. C., 1984: SEDIMENT TRANSPORT, PART II: SUSPENDED
    ! *** LOAD TRANSPORT, J. HYDRAULIC ENGINEERING, 110, 1623-1641.
    IF( WS > 0. )THEN
      REY=1.E4*SNDDIA*( (9.8*(SSG-1.))**0.333 )
      IF( REY <= 10.) TAURS=(4.*WS/REY)**2
      IF( REY > 10.) TAURS=0.016*WS*WS
      VAL=(TAUB/TAURS)-1.
      VAL=MIN(MAX(VAL,0.),100.)
      VAL1=1.-EXP(-0.5*VAL)
      VAL1=0.11*VAL1*(25.-VAL)
      ZEQ1=0.5*VAL1*(DEP**0.7)*(SNDDMX**0.3)
      ZEQ1=ZEQ1/DEP
      CSNDZEQ=MAX(ZEQ1,0.01)
    ELSE
      CSNDZEQ=0.01
    ENDIF

  ELSEIF( IOPT  ==  4 .OR. IOPT == 5 )THEN
    ! *** HAMRICK'S SEDFLUME OPTION
    CSNDZEQ=0.01

  ELSE
    ! *** BAD OPTION
    STOP ' BAD CSNDZEQ OPTION'

  ENDIF

  RETURN

END

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
!
! bugs.eclipse.org bug # 210980
!
! T_CONTAINS was inlined for a contains-stmt.  Thus a contains statement can't
! have a label as it currently stands.  The contains-stmt needs to be put back as
! a statement rather than just an inlined token.
!

module foo

10 contains

subroutine boo
20 continue
end subroutine boo

end module foo

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
Address(Denali Drive) bears the Alaskan Native name for Mount McKinley, North Americas highest mountain.

 It runs from Shasta Drive up to West Covell Boulevard in West Davis.  There is a medical building, apartments and homes next to Denali drive, but there do not appear to be any addresses on Denali.

 Intersecting Streets
  Shasta Drive
  Muir Woods Place
  Shenandoah Place
  Lassen Place
  Bryce Lane
  Isle Royale Lane
  Glacier Drive
  West Covell Boulevard

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! sdl2_pixels.f90
!
! Fortran 2008 interfaces to SDL_pixels.h.
!
! Author:  Philipp Engel
! GitHub:  https://github.com/interkosmos/fortran-sdl2/
! Licence: ISC
module sdl2_pixels
    use, intrinsic :: iso_c_binding
    use :: c_util
    implicit none

    ! SDL_PixelFormatEnum
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_UNKNOWN     = int(z'0')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_INDEX1LSB   = int(z'11100100')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_INDEX1MSB   = int(z'11200100')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_INDEX4LSB   = int(z'12100400')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_INDEX4MSB   = int(z'12200400')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_INDEX8      = int(z'13000801')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGB332      = int(z'14110801')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGB444      = int(z'15120c02')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGB555      = int(z'15130f02')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGR555      = int(z'15530f02')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ARGB4444    = int(z'15321002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGBA4444    = int(z'15421002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ABGR4444    = int(z'15721002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGRA4444    = int(z'15821002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ARGB1555    = int(z'15331002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGBA5551    = int(z'15441002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ABGR1555    = int(z'15731002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGRA5551    = int(z'15841002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGB565      = int(z'15151002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGR565      = int(z'15551002')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGB24       = int(z'17101803')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGR24       = int(z'17401803')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGB888      = int(z'16161804')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGBX8888    = int(z'16261804')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGR888      = int(z'16561804')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGRX8888    = int(z'16661804')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ARGB8888    = int(z'16362004')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_RGBA8888    = int(z'16462004')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ABGR8888    = int(z'16762004')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_BGRA8888    = int(z'16862004')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_ARGB2101010 = int(z'16372004')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_YV12        = int(z'32315659')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_IYUV        = int(z'56555949')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_YUY2        = int(z'32595559')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_UYVY        = int(z'59565955')
    integer(kind=c_int), parameter :: SDL_PIXELFORMAT_YVYU        = int(z'55595659')

    integer(kind=c_int), parameter :: SDL_ALPHA_OPAQUE      = 255
    integer(kind=c_int), parameter :: SDL_ALPHA_TRANSPARENT = 0

    ! SDL_Color
    type, bind(c) :: sdl_color
        integer(kind=c_uint8_t) :: r
        integer(kind=c_uint8_t) :: g
        integer(kind=c_uint8_t) :: b
        integer(kind=c_uint8_t) :: a
    end type sdl_color

    ! SDL_Palette
    type, bind(c) :: sdl_palette
        integer(kind=c_int)      :: ncolors
        type(c_ptr)              :: colors
        integer(kind=c_uint32_t) :: version
        integer(kind=c_int)      :: ref_count
    end type sdl_palette

    ! SDL_PixelFormat
    type, bind(c) :: sdl_pixel_format
        integer(kind=c_uint32_t) :: format
        type(c_ptr)              :: palette
        integer(kind=c_uint8_t)  :: bits_per_pixel
        integer(kind=c_uint8_t)  :: bytes_per_pixel
        integer(kind=c_uint8_t)  :: padding(2)
        integer(kind=c_uint32_t) :: r_mask
        integer(kind=c_uint32_t) :: g_mask
        integer(kind=c_uint32_t) :: b_mask
        integer(kind=c_uint32_t) :: a_mask
        integer(kind=c_uint8_t)  :: r_loss
        integer(kind=c_uint8_t)  :: g_loss
        integer(kind=c_uint8_t)  :: b_loss
        integer(kind=c_uint8_t)  :: a_loss
        integer(kind=c_uint8_t)  :: r_shift
        integer(kind=c_uint8_t)  :: g_shift
        integer(kind=c_uint8_t)  :: b_shift
        integer(kind=c_uint8_t)  :: a_shift
        integer(kind=c_int)      :: ref_count
        type(c_ptr)              :: next
    end type sdl_pixel_format

    public :: sdl_alloc_format_
    public :: sdl_free_format
    public :: sdl_get_rgb
    public :: sdl_get_rgba
    public :: sdl_map_rgb_
    public :: sdl_map_rgba_

    interface
        ! SDL_PixelFormat* SDL_AllocFormat(Uint32 pixel_format)
        function sdl_alloc_format_(pixel_format) bind(c, name='SDL_AllocFormat')
            import :: c_ptr, c_uint32_t
            integer(kind=c_uint32_t), intent(in), value :: pixel_format
            type(c_ptr)                                 :: sdl_alloc_format_
        end function sdl_alloc_format_

        ! Uint32 SDL_MapRGB(const SDL_PixelFormat *format, Uint8 r, Uint8 g, Uint8 b)
        function sdl_map_rgb_(format, r, g, b) bind(c, name='SDL_MapRGB')
            import :: c_uint8_t, c_uint32_t, sdl_pixel_format
            type(sdl_pixel_format) , intent(in)        :: format
            integer(kind=c_uint8_t), intent(in), value :: r
            integer(kind=c_uint8_t), intent(in), value :: g
            integer(kind=c_uint8_t), intent(in), value :: b
            integer(kind=c_uint32_t)                    :: sdl_map_rgb_
        end function sdl_map_rgb_

        ! Uint32 SDL_MapRGBA(const SDL_PixelFormat *format, Uint8 r, Uint8 g, Uint8 b, Uint8 a)
        function sdl_map_rgba_(format, r, g, b, a) bind(c, name='SDL_MapRGBA')
            import :: c_uint8_t, c_uint32_t, sdl_pixel_format
            type(sdl_pixel_format) , intent(in)        :: format
            integer(kind=c_uint8_t), intent(in), value :: r
            integer(kind=c_uint8_t), intent(in), value :: g
            integer(kind=c_uint8_t), intent(in), value :: b
            integer(kind=c_uint8_t), intent(in), value :: a
            integer(kind=c_uint32_t)                    :: sdl_map_rgba_
        end function sdl_map_rgba_

        ! void SDL_FreeFormat(SDL_PixelFormat* format)
        subroutine sdl_free_format(format) bind(c, name='SDL_FreeFormat')
            import :: sdl_pixel_format
            type(sdl_pixel_format), intent(in) :: format
        end subroutine sdl_free_format

        ! void SDL_GetRGB(Uint32 pixel, const SDL_PixelFormat *format, Uint8 *r, Uint8 *g, Uint8 *b)
        subroutine sdl_get_rgb(pixel, format, r, g, b) bind(c, name='SDL_GetRGB')
            import :: c_ptr, c_uint16_t, c_uint32_t, sdl_pixel_format
            integer(kind=c_uint32_t), intent(in), value :: pixel
            type(sdl_pixel_format),   intent(in)        :: format
            integer(kind=c_uint16_t), intent(inout)     :: r
            integer(kind=c_uint16_t), intent(inout)     :: g
            integer(kind=c_uint16_t), intent(inout)     :: b
        end subroutine sdl_get_rgb

        ! void SDL_GetRGBA(Uint32 pixel, const SDL_PixelFormat *format, Uint8 *r, Uint8 *g, Uint8 *b, Uint8 *a)
        subroutine sdl_get_rgba(pixel, format, r, g, b, a) bind(c, name='SDL_GetRGBA')
            import :: c_ptr, c_uint16_t, c_uint32_t, sdl_pixel_format
            integer(kind=c_uint32_t), intent(in), value :: pixel
            type(sdl_pixel_format),   intent(in)        :: format
            integer(kind=c_uint16_t), intent(inout)     :: r
            integer(kind=c_uint16_t), intent(inout)     :: g
            integer(kind=c_uint16_t), intent(inout)     :: b
            integer(kind=c_uint16_t), intent(inout)     :: a
        end subroutine sdl_get_rgba
    end interface
end module sdl2_pixels


text:
program pgm2mask

  use iso_varying_string
  use cmdline_arguments, only: get_options, bad_options, have_args, option_exists, &
      has_value, get_value, assignment(=), next_arg, num_args
  use string_functions, only: join, real
  use file_functions, only: exists, freeunit, stderr, stdout, open
  use pnm_class
  use precision

  implicit none

  integer :: i, j, num, error, nx, ny, unit, maskvalue

  character(len=2000) :: fname, outname

  type (pnm_object) :: pnm
  type (varying_string) :: myoptions(3)

  integer(i8_kind), dimension(:,:), allocatable :: total

  integer, dimension(:,:), allocatable :: data, indices

  ! These are our accepted command line options (see subroutine usage for
  ! an explanation)

  myoptions(1) = 'help'
  myoptions(2) = 'outfile'
  myoptions(3) = 'missing'

  ! This call parses the command line arguments for command line options
  call get_options(myoptions, error)

  ! Check we weren't passed duff options -- spit the dummy if we were
  if (error > 0) then
     write(stderr,*) 'ERROR! Unknown options: ',join(bad_options()," ")
     call usage
     STOP
  end if

  ! Check if we just want to print the usage
  if (option_exists('help')) then
     call usage
     STOP
  end if

  if (num_args() < 1) then
     write(stderr,*) 'ERROR! Must supply pgm file(s) as command-line arguments'
     call usage
     STOP
  end if

  maskvalue = 0

  do while (have_args())

     ! The files are on the command line
     fname = next_arg()

     call read(pnm,trim(fname))

     nx = size(pnm,1)
     ny = size(pnm,2)

     ! Allocate the space required for the data
     allocate(data(nx,ny),indices(nx,ny))

     data = pnm

     ! Reverse the order of the y-axis
     data = data(:,ny:1:-1)

     num = 0
     do j = 1, ny
        do i = 1, nx
           num = num+1
           indices(i,j) = num
        end do
     end do

     ! Use the last input filename as a template for the output filename.
     ! Search for a ".pgm" suffix .. if we find one than
     ! replace it with "_ave.pgm", otherwise just slap "_ave.pgm" on the end
     i = index(fname,".pgm",back=.TRUE.)
     if (i == 0) then
        ! Couldn't find a matching suffix, so look for the end of the
        ! filename string in the character variable
        i = verify(fname," ",back=.TRUE.) + 1
     end if
     outname = fname
     outname(i:) = ".mask"
     print *,i,trim(outname)

     i = rle_encode(trim(outname), pack(indices, data == maskvalue))

     deallocate(data,indices)

  end do

contains

  subroutine usage

    write(stderr,*)
    write(stderr,*) 'Convert pgm files to mask file format'
    write(stderr,*)
    write(stderr,*) 'Usage: pgm2mask [--help] pgmfile(s)'
    write(stderr,*)
    write(stderr,*) '  --help    - print this message'
    write(stderr,*)

  end subroutine usage

  integer function rle_encode (filename, indices) result(nrecords)

    implicit none

    character(len=*), intent(in)  :: filename
    integer, intent(in)           :: indices(:)

    ! Local variables
    integer :: i, num, index, unitnum, masklength

    logical :: masking

    unitnum = freeunit()
    num = size(indices)

    open (unit=unitnum, file=filename, form='unformatted', status='replace')

    masking = .false.
    masklength = 1
    index = indices(1)

    write(unitnum) 'RLE'
    write(unitnum) num

    ! Return value
    nrecords = 0

    if (num /= 0) then
       do i = 2, num
          if ((indices(i) - indices(i-1)) == 1) then
             masklength=masklength+1
          else
             write(unitnum)index,masklength
             nrecords = nrecords + 1
             index = indices(i)
             masklength = 1
          end if
       end do

       nrecords = nrecords + 1
       ! Must write out the last record
       write(unitnum)index,masklength

    end if

    close(unitnum)

  end function rle_encode

end program pgm2mask

max_stars_repo_name: Welberry/pgmutils

max_issues_repo_name: Welberry/pgmutils

max_forks_repo_name: Welberry/pgmutils


text:
  module Mod_nsm_ComputeVorticity
   use typre
   use Mod_PointerSetter
   use Mod_nsm_BaseElmope
   use Mod_nsm_InterpolateGradients
   implicit none
   private
   public SetPointersComputeVorticity

   type, extends(PointerSetter) :: SPComputeVorticity
contains
      procedure :: SpecificSet => SpecificSetComputeVorticity
   end type
   type(SPComputeVorticity) :: SetPointersComputeVorticity

   !Vorticity
   real(rp), allocatable    :: elvort(:,:), gpvort(:)

contains

   !----------------------------------------------------------------------------
   !Setting Pointers
   subroutine SpecificSetComputeVorticity(d)
      implicit none
      class(SPComputeVorticity) :: d
      logical :: aux_logic

      !Logics for deciding if we need to compute vorticity
      aux_logic = .false.
      if (a%npp_stepi(10) /= 0) then
         if (mod(a%istep,a%npp_stepi(10))==0) then
            aux_logic = .true.
         endif
      endif

      !If vorticity needs to be computed
      if (aux_logic .eqv. .true.) then
         call SetPointersInterpolateGradients%Set
         call ConcatenateProcedures(ProcHook_Initializations,AllocVorti)
         call ConcatenateProcedures(ProcHook_InGaussElmatsAssembly,ComputeVorticity)
         call ConcatenateProcedures(ProcHook_AssemblyEndite,AssemblyVorticity)
         call ConcatenateProcedures(ProcHook_Finalizations,DeallocVorti)
         call ConcatenateProcedures(ProcHook_PostLoop,SmoothVorticity)
      endif
   end subroutine

   !Vorticity
   subroutine AllocVorti
      implicit none

      call a%Memor%alloc(3,e%mnode,elvort,'elvort','nsm_EndElmope')
      call a%Memor%alloc(3,gpvort,'gpvort','nsm_EndElmope')
      a%vorti(:,:) = 0_rp
   end subroutine

   subroutine ComputeVorticity
      implicit none
      integer(ip)  :: inode

      gpvort(3) = grvel(2,1)-grvel(1,2)

      if (e%ndime == 3) then
         gpvort(1) = grvel(3,2)-grvel(2,3)
         gpvort(2) = grvel(1,3)-grvel(3,1)
      endif
      do inode=1,e%pnode
         elvort(:,inode)=elvort(:,inode)+e%shape(inode,e%igaus)*gpvort*dvol
      end do
   end subroutine

   subroutine AssemblyVorticity
      implicit none
      integer(ip)  :: inode

!       do inode=1,e%pnode
!          a%vorti(:,e%lnods(inode)) = a%vorti(:,e%lnods(inode)) + elvort(:,inode)
!       end do
      call a%Mesh%AssemblyToArray(e,size(elvort,1),elvort,a%vorti)
      elvort(:,:)=0_rp
   end subroutine

   subroutine SmoothVorticity
      implicit none

      call a%Mesh%Smooth(3,a%vorti)

   end subroutine

   subroutine DeallocVorti
      implicit none

      call a%Memor%dealloc(3,e%mnode,elvort,'elvort','nsm_EndElmope')
      call a%Memor%dealloc(3,gpvort,'gpvort','nsm_EndElmope')
   end subroutine
end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
module buffered_ttgt_tc_driver_module
    use, intrinsic :: iso_fortran_env, only : int64
    use :: data_api, only : &
            stream, &
            scratch_buffer

    use :: tensor_api, only : &
            tensor, &
            scalar

    use :: tc_ttgt_module, only : tc_ttgt
    use :: ttgt_descriptor_module, only : ttgt_descriptor
    use :: tensor_permute_module, only : tensor_permute
    use :: matrix_multiplication_module, only : matrix_multiplication

    implicit none
    private

    public :: buffered_ttgt_tc_driver

    type, extends(tc_ttgt) :: buffered_ttgt_tc_driver
        type(scratch_buffer) :: scratch
    contains
        procedure :: synchronize => synchronize
        procedure :: cleanup => cleanup
        procedure :: setup => setup
        procedure :: set_scratch_buffer => set_scratch_buffer
        procedure :: allocate_tensor => allocate_tensor
        procedure :: destroy_permuted_tensor => destroy_permuted_tensor
        procedure :: destroy_permuted_tensors => destroy_permuted_tensors
        procedure :: clear => clear
    end type buffered_ttgt_tc_driver

    interface buffered_ttgt_tc_driver
        module procedure constructor_empty
        module procedure constructor
    end interface buffered_ttgt_tc_driver
contains
    function constructor_empty() result(this)
        type(buffered_ttgt_tc_driver) :: this

        call this%clear()
    end function constructor_empty

    function constructor(descr, permuter, mm, scratch) result(this)
        type(ttgt_descriptor), intent(in) :: descr
        class(tensor_permute), intent(in) :: permuter
        class(matrix_multiplication), intent(in) :: mm
        type(scratch_buffer), intent(in) :: scratch
        type(buffered_ttgt_tc_driver) :: this

        this = buffered_ttgt_tc_driver()

        call this%setup(descr, permuter, mm, scratch)
    end function constructor

    subroutine destroy_permuted_tensors(this, c, a, b, astream)
        class(buffered_ttgt_tc_driver), intent(inout) :: this
        class(tensor), intent(inout) :: c, a, b
        type(stream), intent(in), optional :: astream

        call this%destroy_permuted_tensor(b, this%descr%permute_b)
        call this%destroy_permuted_tensor(a, this%descr%permute_a)
        call this%destroy_permuted_tensor(c, this%descr%permute_c)
        call this%scratch%checkpoint(astream)
    end subroutine destroy_permuted_tensors

    subroutine synchronize(this, astream)
        class(buffered_ttgt_tc_driver), intent(in) :: this
        type(stream), intent(in), optional :: astream

        call this%permuter%synchronize(astream)
    end subroutine synchronize

    subroutine cleanup(this)
        class(buffered_ttgt_tc_driver), intent(inout) :: this

        if (allocated(this%permuter)) then
            call this%permuter%cleanup()
            deallocate(this%permuter)
        end if

        if (allocated(this%mm)) then
            call this%mm%cleanup()
            deallocate(this%mm)
        end if

        call this%builder%cleanup()
        call this%scratch%cleanup()

        call this%clear()
    end subroutine cleanup

    subroutine setup(this, descr, permuter, mm, scratch)
        class(buffered_ttgt_tc_driver), intent(inout) :: this
        type(ttgt_descriptor), intent(in) :: descr
        class(tensor_permute), intent(in) :: permuter
        class(matrix_multiplication), intent(in) :: mm
        type(scratch_buffer), intent(in) :: scratch

        call this%set_descriptor(descr)
        call this%set_permuter(permuter)
        call this%set_matrix_multiplication(mm)
        call this%set_scratch_buffer(scratch)
    end subroutine setup

    subroutine set_scratch_buffer(this, scratch)
        class(buffered_ttgt_tc_driver), intent(inout) :: this
        type(scratch_buffer), intent(in) ::scratch

        this%scratch = scratch
        call this%scratch%initialize()
    end subroutine set_scratch_buffer

    subroutine allocate_tensor(this, src, dst, order, initialize)
        class(buffered_ttgt_tc_driver), intent(inout) :: this
        class(tensor), intent(in) :: src
        class(tensor), allocatable, intent(inout) :: dst
        integer, dimension(:), intent(in) :: order
        logical, intent(in) :: initialize

        call this%builder%allocate_and_create_in_scratch(dst, this%scratch, src%datatype, src%dims(order), initialize)
    end subroutine allocate_tensor

    subroutine destroy_permuted_tensor(this, dst, permute)
        class(buffered_ttgt_tc_driver), intent(inout) :: this
        class(tensor), intent(inout) :: dst
        logical, intent(in) :: permute

        if ( permute ) then
            call this%scratch%destroy(dst%storage)
        end if
        call dst%release()
    end subroutine destroy_permuted_tensor

    subroutine clear(this)
        class(buffered_ttgt_tc_driver), intent(inout) :: this

        this%descr = ttgt_descriptor()
    end subroutine clear
end module buffered_ttgt_tc_driver_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


text: Anne Richardson MSW, CHT appears to be a hypnotherapy hypnotherapist.  Possibly the only one in Davis?


!
! Copyright (C) by Argonne National Laboratory
!     See COPYRIGHT in top-level directory
!

subroutine MPI_Accumulate_f08ts(origin_addr, origin_count, origin_datatype, target_rank, &
    target_disp, target_count, target_datatype, op, win, ierror)
    use, intrinsic :: iso_c_binding, only : c_int
    use :: mpi_f08, only : MPI_Datatype, MPI_Op, MPI_Win
    use :: mpi_f08, only : MPI_ADDRESS_KIND
    use :: mpi_c_interface, only : c_Datatype, c_Op, c_Win
    use :: mpi_c_interface, only : MPIR_Accumulate_cdesc

    implicit none

    type(*), dimension(..), intent(in), asynchronous :: origin_addr
    integer, intent(in) :: origin_count
    integer, intent(in) :: target_rank
    integer, intent(in) :: target_count
    type(MPI_Datatype), intent(in) :: origin_datatype
    integer(MPI_ADDRESS_KIND), intent(in) :: target_disp
    type(MPI_Datatype), intent(in) :: target_datatype
    type(MPI_Op), intent(in) :: op
    type(MPI_Win), intent(in) :: win
    integer, optional, intent(out) :: ierror

    integer(c_int) :: origin_count_c
    integer(c_int) :: target_rank_c
    integer(c_int) :: target_count_c
    integer(c_Datatype) :: origin_datatype_c
    integer(MPI_ADDRESS_KIND) :: target_disp_c
    integer(c_Datatype) :: target_datatype_c
    integer(c_Op) :: op_c
    integer(c_Win) :: win_c
    integer(c_int) :: ierror_c

    if (c_int == kind(0)) then
        ierror_c = MPIR_Accumulate_cdesc(origin_addr, origin_count, origin_datatype%MPI_VAL, target_rank, &
            target_disp, target_count, target_datatype%MPI_VAL, op%MPI_VAL, win%MPI_VAL)
    else
        origin_count_c = origin_count
        origin_datatype_c = origin_datatype%MPI_VAL
        target_rank_c = target_rank
        target_disp_c = target_disp
        target_count_c = target_count
        target_datatype_c = target_datatype%MPI_VAL
        op_c = op%MPI_VAL
        win_c = win%MPI_VAL
        ierror_c = MPIR_Accumulate_cdesc(origin_addr, origin_count_c, origin_datatype_c, target_rank_c, &
            target_disp_c, target_count_c, target_datatype_c, op_c, win_c)
    end if

    if (present(ierror)) ierror = ierror_c

end subroutine MPI_Accumulate_f08ts


text:
C
C     file gnbnaux.f
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 1999 by UCAR                   *
C     *                                                               *
C     *       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                      FISHPACK version 4.1                     *
C     *                                                               *
C     *     A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF      *
C     *                                                               *
C     *      SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS        *
C     *                                                               *
C     *                             BY                                *
C     *                                                               *
C     *        JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET         *
C     *                                                               *
C     *                             OF                                *
C     *                                                               *
C     *         THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH          *
C     *                                                               *
C     *                BOULDER, COLORADO  (80307)  U.S.A.             *
C     *                                                               *
C     *                   WHICH IS SPONSORED BY                       *
C     *                                                               *
C     *              THE NATIONAL SCIENCE FOUNDATION                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C PACKAGE GNBNAUX
C
C LATEST REVISION        NOVEMBER 1988
C
C PURPOSE                TO PROVIDE AUXILIARY ROUTINES FOR FISHPACK
C                        ENTRIES GENBUN AND POISTG.
C
C USAGE                  THERE ARE NO USER ENTRIES IN THIS PACKAGE.
C                        THE ROUTINES IN THIS PACKAGE ARE NOT INTENDED
C                        TO BE CALLED BY USERS, BUT RATHER BY ROUTINES
C                        IN PACKAGES GENBUN AND POISTG.
C
C SPECIAL CONDITIONS     NONE
C
C I/O                    NONE
C
C PRECISION              SINGLE
C
C REQUIRED LIBRARY       COMF FROM FISHPACK
C FILES
C
C LANGUAGE               FORTRAN
C
C HISTORY                WRITTEN IN 1979 BY ROLAND SWEET OF NCAR'S
C                        SCIENTIFIC COMPUTING DIVISION.  MADE AVAILABLE
C                        ON NCAR'S PUBLIC LIBRARIES IN JANUARY, 1980.
C
C PORTABILITY            FORTRAN 77
C ********************************************************************
      SUBROUTINE COSGEN (N,IJUMP,FNUM,FDEN,A)
      DIMENSION       A(*)
C
C
C     THIS SUBROUTINE COMPUTES REQUIRED COSINE VALUES IN ASCENDING
C     ORDER.  WHEN IJUMP .GT. 1 THE ROUTINE COMPUTES VALUES
C
C        2*COS(J*PI/L) , J=1,2,...,L AND J .NE. 0(MOD N/IJUMP+1)
C
C     WHERE L = IJUMP*(N/IJUMP+1).
C
C
C     WHEN IJUMP = 1 IT COMPUTES
C
C            2*COS((J-FNUM)*PI/(N+FDEN)) ,  J=1, 2, ... ,N
C
C     WHERE
C        FNUM = 0.5, FDEN = 0.0,  FOR REGULAR REDUCTION VALUES
C        FNUM = 0.0, FDEN = 1.0, FOR B-R AND C-R WHEN ISTAG = 1
C        FNUM = 0.0, FDEN = 0.5, FOR B-R AND C-R WHEN ISTAG = 2
C        FNUM = 0.5, FDEN = 0.5, FOR B-R AND C-R WHEN ISTAG = 2
C                                IN POISN2 ONLY.
C
C
      PI = PIMACH(DUM)
      IF (N .EQ. 0) GO TO 105
      IF (IJUMP .EQ. 1) GO TO 103
      K3 = N/IJUMP+1
      K4 = K3-1
      PIBYN = PI/FLOAT(N+IJUMP)
      DO 102 K=1,IJUMP
         K1 = (K-1)*K3
         K5 = (K-1)*K4
         DO 101 I=1,K4
            X = K1+I
            K2 = K5+I
            A(K2) = -2.*COS(X*PIBYN)
  101    CONTINUE
  102 CONTINUE
      GO TO 105
  103 CONTINUE
      NP1 = N+1
      Y = PI/(FLOAT(N)+FDEN)
      DO 104 I=1,N
         X = FLOAT(NP1-I)-FNUM
         A(I) = 2.*COS(X*Y)
  104 CONTINUE
  105 CONTINUE
      RETURN
      END
      SUBROUTINE MERGE (TCOS,I1,M1,I2,M2,I3)
      DIMENSION       TCOS(*)
C
C     THIS SUBROUTINE MERGES TWO ASCENDING STRINGS OF NUMBERS IN THE
C     ARRAY TCOS.  THE FIRST STRING IS OF LENGTH M1 AND STARTS AT
C     TCOS(I1+1).  THE SECOND STRING IS OF LENGTH M2 AND STARTS AT
C     TCOS(I2+1).  THE MERGED STRING GOES INTO TCOS(I3+1).
C
C
      J1 = 1
      J2 = 1
      J = I3
      IF (M1 .EQ. 0) GO TO 107
      IF (M2 .EQ. 0) GO TO 104
  101 J = J+1
      L = J1+I1
      X = TCOS(L)
      L = J2+I2
      Y = TCOS(L)
      IF (X-Y) 102,102,103
  102 TCOS(J) = X
      J1 = J1+1
      IF (J1 .GT. M1) GO TO 106
      GO TO 101
  103 TCOS(J) = Y
      J2 = J2+1
      IF (J2 .LE. M2) GO TO 101
      IF (J1 .GT. M1) GO TO 109
  104 K = J-J1+1
      DO 105 J=J1,M1
         M = K+J
         L = J+I1
         TCOS(M) = TCOS(L)
  105 CONTINUE
      GO TO 109
  106 CONTINUE
      IF (J2 .GT. M2) GO TO 109
  107 K = J-J2+1
      DO 108 J=J2,M2
         M = K+J
         L = J+I2
         TCOS(M) = TCOS(L)
  108 CONTINUE
  109 CONTINUE
      RETURN
      END
      SUBROUTINE TRIX (IDEGBR,IDEGCR,M,A,B,C,Y,TCOS,D,W)
C
C     SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS WHERE THE
C     COEFFICIENT MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY
C     TRIDIAGONAL  ( . . . , A(I), B(I), C(I), . . . ).
C
      DIMENSION       A(*)       ,B(*)       ,C(*)       ,Y(*)       ,
     1                TCOS(*)    ,D(*)       ,W(*)
      MM1 = M-1
      IFB = IDEGBR+1
      IFC = IDEGCR+1
      L = IFB/IFC
      LINT = 1
      DO 108 K=1,IDEGBR
         X = TCOS(K)
         IF (K .NE. L) GO TO 102
         I = IDEGBR+LINT
         XX = X-TCOS(I)
         DO 101 I=1,M
            W(I) = Y(I)
            Y(I) = XX*Y(I)
  101    CONTINUE
  102    CONTINUE
         Z = 1./(B(1)-X)
         D(1) = C(1)*Z
         Y(1) = Y(1)*Z
         DO 103 I=2,MM1
            Z = 1./(B(I)-X-A(I)*D(I-1))
            D(I) = C(I)*Z
            Y(I) = (Y(I)-A(I)*Y(I-1))*Z
  103    CONTINUE
         Z = B(M)-X-A(M)*D(MM1)
         IF (Z .NE. 0.) GO TO 104
         Y(M) = 0.
         GO TO 105
  104    Y(M) = (Y(M)-A(M)*Y(MM1))/Z
  105    CONTINUE
         DO 106 IP=1,MM1
            I = M-IP
            Y(I) = Y(I)-D(I)*Y(I+1)
  106    CONTINUE
         IF (K .NE. L) GO TO 108
         DO 107 I=1,M
            Y(I) = Y(I)+W(I)
  107    CONTINUE
         LINT = LINT+1
         L = (LINT*IFB)/IFC
  108 CONTINUE
      RETURN
      END
      SUBROUTINE TRI3 (M,A,B,C,K,Y1,Y2,Y3,TCOS,D,W1,W2,W3)
      DIMENSION       A(*)       ,B(*)       ,C(*)       ,K(4)       ,
     1                TCOS(*)    ,Y1(*)      ,Y2(*)      ,Y3(*)      ,
     2                D(*)       ,W1(*)      ,W2(*)      ,W3(*)
C
C     SUBROUTINE TO SOLVE THREE LINEAR SYSTEMS WHOSE COMMON COEFFICIENT
C     MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY
C
C                  TRIDIAGONAL (...,A(I),B(I),C(I),...)
C
      MM1 = M-1
      K1 = K(1)
      K2 = K(2)
      K3 = K(3)
      K4 = K(4)
      IF1 = K1+1
      IF2 = K2+1
      IF3 = K3+1
      IF4 = K4+1
      K2K3K4 = K2+K3+K4
      IF (K2K3K4 .EQ. 0) GO TO 101
      L1 = IF1/IF2
      L2 = IF1/IF3
      L3 = IF1/IF4
      LINT1 = 1
      LINT2 = 1
      LINT3 = 1
      KINT1 = K1
      KINT2 = KINT1+K2
      KINT3 = KINT2+K3
  101 CONTINUE
      DO 115 N=1,K1
         X = TCOS(N)
         IF (K2K3K4 .EQ. 0) GO TO 107
         IF (N .NE. L1) GO TO 103
         DO 102 I=1,M
            W1(I) = Y1(I)
  102    CONTINUE
  103    IF (N .NE. L2) GO TO 105
         DO 104 I=1,M
            W2(I) = Y2(I)
  104    CONTINUE
  105    IF (N .NE. L3) GO TO 107
         DO 106 I=1,M
            W3(I) = Y3(I)
  106    CONTINUE
  107    CONTINUE
         Z = 1./(B(1)-X)
         D(1) = C(1)*Z
         Y1(1) = Y1(1)*Z
         Y2(1) = Y2(1)*Z
         Y3(1) = Y3(1)*Z
         DO 108 I=2,M
            Z = 1./(B(I)-X-A(I)*D(I-1))
            D(I) = C(I)*Z
            Y1(I) = (Y1(I)-A(I)*Y1(I-1))*Z
            Y2(I) = (Y2(I)-A(I)*Y2(I-1))*Z
            Y3(I) = (Y3(I)-A(I)*Y3(I-1))*Z
  108    CONTINUE
         DO 109 IP=1,MM1
            I = M-IP
            Y1(I) = Y1(I)-D(I)*Y1(I+1)
            Y2(I) = Y2(I)-D(I)*Y2(I+1)
            Y3(I) = Y3(I)-D(I)*Y3(I+1)
  109    CONTINUE
         IF (K2K3K4 .EQ. 0) GO TO 115
         IF (N .NE. L1) GO TO 111
         I = LINT1+KINT1
         XX = X-TCOS(I)
         DO 110 I=1,M
            Y1(I) = XX*Y1(I)+W1(I)
  110    CONTINUE
         LINT1 = LINT1+1
         L1 = (LINT1*IF1)/IF2
  111    IF (N .NE. L2) GO TO 113
         I = LINT2+KINT2
         XX = X-TCOS(I)
         DO 112 I=1,M
            Y2(I) = XX*Y2(I)+W2(I)
  112    CONTINUE
         LINT2 = LINT2+1
         L2 = (LINT2*IF1)/IF3
  113    IF (N .NE. L3) GO TO 115
         I = LINT3+KINT3
         XX = X-TCOS(I)
         DO 114 I=1,M
            Y3(I) = XX*Y3(I)+W3(I)
  114    CONTINUE
         LINT3 = LINT3+1
         L3 = (LINT3*IF1)/IF4
  115 CONTINUE
      RETURN
C
C REVISION HISTORY---
C
C SEPTEMBER 1973    VERSION 1
C APRIL     1976    VERSION 2
C JANUARY   1978    VERSION 3
C DECEMBER  1979    VERSION 3.1
C OCTOBER   1980    CHANGED SEVERAL DIVIDES OF FLOATING INTEGERS
C                   TO INTEGER DIVIDES TO ACCOMODATE CRAY-1 ARITHMETIC.
C FEBRUARY  1985    DOCUMENTATION UPGRADE
C NOVEMBER  1988    VERSION 3.2, FORTRAN 77 CHANGES
C-----------------------------------------------------------------------
      END

max_stars_repo_name: ipelupessy/omuse

max_issues_repo_name: ipelupessy/omuse

max_forks_repo_name: ipelupessy/omuse