text:
!
! SEcategory_Define
!
! Module defining the SEcategory object.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 17-Aug-2011
!                       paul.vandelst@noaa.gov

MODULE SEcategory_Define

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds           , ONLY: fp, Long, Double
  USE Message_Handler      , ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  USE Compare_Float_Numbers, ONLY: OPERATOR(.EqualTo.)
  USE File_Utility         , ONLY: File_Open, File_Exists
  USE Binary_File_Utility  , ONLY: Open_Binary_File        , &
                                   WriteGAtts_Binary_File  , &
                                   ReadGAtts_Binary_File   , &
                                   WriteLogical_Binary_File, &
                                   ReadLogical_Binary_File
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Parameters
  PUBLIC :: SECATEGORY_DATATYPE
  ! Datatypes
  PUBLIC :: SEcategory_type
  ! Operators
  PUBLIC :: OPERATOR(==)
  ! Procedures
  PUBLIC :: SEcategory_Associated
  PUBLIC :: SEcategory_Destroy
  PUBLIC :: SEcategory_Create
  PUBLIC :: SEcategory_Inspect
  PUBLIC :: SEcategory_ValidRelease
  PUBLIC :: SEcategory_Info
  PUBLIC :: SEcategory_Name
  PUBLIC :: SEcategory_Index
  PUBLIC :: SEcategory_DefineVersion
  PUBLIC :: SEcategory_SetValue
  PUBLIC :: SEcategory_GetValue
  PUBLIC :: SEcategory_InquireFile
  PUBLIC :: SEcategory_ReadFile
  PUBLIC :: SEcategory_WriteFile

  ! ---------------------
  ! Procedure overloading
  ! ---------------------
  INTERFACE OPERATOR(==)
    MODULE PROCEDURE SEcategory_Equal
  END INTERFACE OPERATOR(==)

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
    '$Id: SEcategory_Define.f90,v 1.1 2013/04/07 22:56:27 jguo Exp $'
  ! Datatype information
  CHARACTER(*), PARAMETER :: SECATEGORY_DATATYPE = 'SEcategory'
  ! Release and version
  INTEGER, PARAMETER :: SECATEGORY_RELEASE = 3  ! This determines structure and file formats.
  INTEGER, PARAMETER :: SECATEGORY_VERSION = 1  ! This is just the default data version.
  ! Close status for write errors
  CHARACTER(*), PARAMETER :: WRITE_ERROR_STATUS = 'DELETE'
  ! Literal constants
  REAL(fp), PARAMETER :: ZERO = 0.0_fp
  REAL(fp), PARAMETER :: ONE  = 1.0_fp
  ! String lengths
  INTEGER,  PARAMETER :: ML = 256 ! Message length
  INTEGER,  PARAMETER :: SL =  80 ! String length

  ! ----------------------------------
  ! SEcategory data type definitions
  ! ----------------------------------
  !:tdoc+:
  TYPE :: SEcategory_type
    ! Allocation indicator
    LOGICAL :: Is_Allocated = .FALSE.
    ! Datatype information
    CHARACTER(SL) :: Datatype_Name = SECATEGORY_DATATYPE
    ! Release and version information
    INTEGER(Long) :: Release = SECATEGORY_RELEASE
    INTEGER(Long) :: Version = SECATEGORY_VERSION
    ! Surface classification name
    CHARACTER(SL) :: Classification_Name = ''
    ! Dimensions
    INTEGER(Long) :: String_Length   = SL
    INTEGER(Long) :: n_Frequencies   = 0  ! L dim.
    INTEGER(Long) :: n_Surface_Types = 0  ! N dim.
    ! Dimensional vectors
    REAL(Double),  ALLOCATABLE :: Frequency(:)             ! Lx1
    CHARACTER(SL), ALLOCATABLE :: Surface_Type(:)          ! Nx1
    ! Surface type validity array
    LOGICAL,       ALLOCATABLE :: Surface_Type_IsValid(:)  ! Nx1
    ! Reflectance LUT data
    REAL(Double),  ALLOCATABLE :: Reflectance(:,:)         ! LxN
  END TYPE SEcategory_type
  !:tdoc-:

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                           ## PUBLIC PROCEDURES ##                          ##
!##                                                                            ##
!################################################################################
!################################################################################

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Associated
!
! PURPOSE:
!       Elemental function to test the status of the allocatable components
!       of the SEcategory structure.
!
! CALLING SEQUENCE:
!       Status = SEcategory_Associated( SEcategory )
!
! OBJECTS:
!       SEcategory: Structure which is to have its member's
!                   status tested.
!                   UNITS:      N/A
!                   TYPE:       SEcategory_type
!                   DIMENSION:  Scalar or any rank
!                   ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Status:     The return value is a logical value indicating the
!                   status of the NLTE members.
!                    .TRUE.  - if ANY of the SEcategory allocatable members
!                              are in use.
!                    .FALSE. - if ALL of the SEcategory allocatable members
!                              are not in use.
!                   UNITS:      N/A
!                   TYPE:       LOGICAL
!                   DIMENSION:  Same as input
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION SEcategory_Associated( self ) RESULT( Status )
    TYPE(SEcategory_type), INTENT(IN) :: self
    LOGICAL :: Status
    Status = self%Is_Allocated
  END FUNCTION SEcategory_Associated

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Destroy
!
! PURPOSE:
!       Elemental subroutine to re-initialize SEcategory objects.
!
! CALLING SEQUENCE:
!       CALL SEcategory_Destroy( SEcategory )
!
! OBJECTS:
!       SEcategory:   Re-initialized SEcategory structure.
!                     UNITS:      N/A
!                     TYPE:       SEcategory_type
!                     DIMENSION:  Scalar or any rank
!                     ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE SEcategory_Destroy( self )
    TYPE(SEcategory_type), INTENT(OUT) :: self
    self%Is_Allocated = .FALSE.
    self%n_Frequencies   = 0
    self%n_Surface_Types = 0
  END SUBROUTINE SEcategory_Destroy

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Create
!
! PURPOSE:
!       Elemental subroutine to create an instance of an SEcategory object.
!
! CALLING SEQUENCE:
!       CALL SEcategory_Create( SEcategory     , &
!                               n_Frequencies  , &
!                               n_Surface_Types  )
!
! OBJECTS:
!       SEcategory:         SEcategory object structure.
!                           UNITS:      N/A
!                           TYPE:       SEcategory_type
!                           DIMENSION:  Scalar or any rank
!                           ATTRIBUTES: INTENT(OUT)
!
! INPUTS:
!       n_Frequencies:      Number of spectral frequencies for which there are
!                           reflectance data.
!                           Must be > 0.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Conformable with the SEcategory object
!                           ATTRIBUTES: INTENT(IN)
!
!       n_Surface_Types:    Number of land surface types for which is are
!                           reflectance data.
!                           Must be > 0.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Conformable with the SEcategory object
!                           ATTRIBUTES: INTENT(IN)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE SEcategory_Create( &
    self           , &  ! Output
    n_Frequencies  , &  ! Input
    n_Surface_Types  )  ! Input
    ! Arguments
    TYPE(SEcategory_type), INTENT(OUT) :: self
    INTEGER              , INTENT(IN)  :: n_Frequencies
    INTEGER              , INTENT(IN)  :: n_Surface_Types
    ! Local variables
    INTEGER :: alloc_stat

    ! Check input
    IF ( n_Frequencies   < 1 .OR. &
         n_Surface_Types < 1 ) RETURN

    ! Perform the allocation
    ALLOCATE( self%Surface_Type_IsValid( n_Surface_Types ), &
              self%Frequency( n_Frequencies ), &
              self%Surface_Type( n_Surface_Types ), &
              self%Reflectance( n_Frequencies, n_Surface_Types ), &
              STAT = alloc_stat )
    IF ( alloc_stat /= 0 ) RETURN

    ! Initialise
    ! ...Dimensions
    self%n_Frequencies   = n_Frequencies
    self%n_Surface_Types = n_Surface_Types
    ! ...Arrays
    self%Surface_Type_IsValid = .TRUE.
    self%Frequency            = ZERO
    self%Surface_Type         = ''
    self%Reflectance          = ZERO

    ! Set allocation indicator
    self%Is_Allocated = .TRUE.

  END SUBROUTINE SEcategory_Create

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Inspect
!
! PURPOSE:
!       Subroutine to print the contents of a SEcategory object to stdout.
!
! CALLING SEQUENCE:
!       CALL SEcategory_Inspect( SEcategory )
!
! OBJECTS:
!       SEcategory:    SEcategory object to display.
!                      UNITS:      N/A
!                      TYPE:       SEcategory_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE SEcategory_Inspect( self)
    TYPE(SEcategory_type), INTENT(IN) :: self
    INTEGER :: n
    WRITE(*,'(1x,"SEcategory OBJECT")')
    ! Release/version info
    WRITE(*,'(3x,"Release.Version     :",1x,i0,".",i0)') self%Release, self%Version
    ! Surface classification name
    WRITE(*,'(3x,"Classification_Name :",1x,a)') TRIM(self%Classification_Name)
    ! Dimensions
    WRITE(*,'(3x,"n_Frequencies       :",1x,i0)') self%n_Frequencies
    WRITE(*,'(3x,"n_Surface_Types     :",1x,i0)') self%n_Surface_Types
    IF ( .NOT. SEcategory_Associated(self) ) RETURN
    ! Dimension arrays
    WRITE(*,'(3x,"Frequency :")')
    WRITE(*,'(5(1x,es13.6,:))') self%Frequency
    WRITE(*,'(3x,"Surface_Type - (IsValid) :")')
    DO n = 1, self%n_Surface_Types
      WRITE(*,'(5x,a," - (",l1,")")') TRIM(self%Surface_Type(n)), self%Surface_Type_IsValid(n)
    END DO
    ! Reflectance array
    WRITE(*,'(3x,"Reflectance :")')
    DO n = 1, self%n_Surface_Types
      WRITE(*,'(5x,a)') self%Surface_Type(n)
      WRITE(*,'(5(1x,es13.6,:))') self%Reflectance(:,n)
    END DO
  END SUBROUTINE SEcategory_Inspect

!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_ValidRelease
!
! PURPOSE:
!       Function to check the SEcategory Release value.
!
! CALLING SEQUENCE:
!       IsValid = SEcategory_ValidRelease( SEcategory )
!
! INPUTS:
!       SEcategory:    SEcategory object for which the Release component
!                      is to be checked.
!                      UNITS:      N/A
!                      TYPE:       SEcategory_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       IsValid:       Logical value defining the release validity.
!                      UNITS:      N/A
!                      TYPE:       LOGICAL
!                      DIMENSION:  Scalar
!
!:sdoc-:
!----------------------------------------------------------------------------------

  FUNCTION SEcategory_ValidRelease( self ) RESULT( IsValid )
    ! Arguments
    TYPE(SEcategory_type), INTENT(IN) :: self
    ! Function result
    LOGICAL :: IsValid
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'SEcategory_ValidRelease'
    ! Local variables
    CHARACTER(ML) :: msg

    ! Set up
    IsValid = .TRUE.

    ! Check release is not too old
    IF ( self%Release < SECATEGORY_RELEASE ) THEN
      IsValid = .FALSE.
      WRITE( msg,'("An SEcategory data update is needed. ", &
                  &"SEcategory release is ",i0,". Valid release is ",i0,"." )' ) &
                  self%Release, SECATEGORY_RELEASE
      CALL Display_Message( ROUTINE_NAME, msg, INFORMATION ); RETURN
    END IF

    ! Check release is not too new
    IF ( self%Release > SECATEGORY_RELEASE ) THEN
      IsValid = .FALSE.
      WRITE( msg,'("An SEcategory software update is needed. ", &
                  &"SEcategory release is ",i0,". Valid release is ",i0,"." )' ) &
                  self%Release, SECATEGORY_RELEASE
      CALL Display_Message( ROUTINE_NAME, msg, INFORMATION ); RETURN
    END IF

  END FUNCTION SEcategory_ValidRelease

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Info
!
! PURPOSE:
!       Subroutine to return a string containing version and dimension
!       information about a SEcategory object.
!
! CALLING SEQUENCE:
!       CALL SEcategory_Info( SEcategory, Info )
!
! OBJECTS:
!       SEcategory:    SEcategory object about which info is required.
!                      UNITS:      N/A
!                      TYPE:       SEcategory_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Info:          String containing version and dimension information
!                      about the SEcategory object.
!                      UNITS:      N/A
!                      TYPE:       CHARACTER(*)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE SEcategory_Info( self, Info )
    ! Arguments
    TYPE(SEcategory_type), INTENT(IN)  :: self
    CHARACTER(*),          INTENT(OUT) :: Info
    ! Parameters
    INTEGER, PARAMETER :: CARRIAGE_RETURN = 13
    INTEGER, PARAMETER :: LINEFEED = 10
    ! Local variables
    CHARACTER(2000) :: Long_String

    ! Write the required data to the local string
    WRITE( Long_String, &
           '(a,1x,"SEcategory RELEASE.VERSION: ",i2,".",i2.2,a,3x, &
           &"CLASSIFICATION: ",a,",",2x,&
           &"N_FREQUENCIES=",i0,2x,&
           &"N_SURFACE_TYPES=",i0 )' ) &
           ACHAR(CARRIAGE_RETURN)//ACHAR(LINEFEED), &
           self%Release, self%Version, &
           ACHAR(CARRIAGE_RETURN)//ACHAR(LINEFEED), &
           TRIM(self%Classification_Name), &
           self%n_Frequencies, &
           self%n_Surface_Types

    ! Trim the output based on the
    ! dummy argument string length
    Info = Long_String(1:MIN(LEN(Info), LEN_TRIM(Long_String)))

  END SUBROUTINE SEcategory_Info

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Name
!
! PURPOSE:
!       Pure function to return the datatype name of an SEcategory object.
!
! CALLING SEQUENCE:
!       datatype_name = SEcategory_Name( SEcategory )
!
! OBJECTS:
!       SEcategory:   SEcategory object structure.
!                     UNITS:      N/A
!                     TYPE:       SEcategory_type
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Status:       The return value is a the character string containing
!                     the datatype name of the structure.
!                     UNITS:      N/A
!                     TYPE:       CHARACTER
!                     DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE FUNCTION SEcategory_Name( self ) RESULT( datatype_name )
    ! Arguments
    TYPE(SEcategory_type), INTENT(IN) :: self
    ! Function result
    CHARACTER(LEN(self%Datatype_Name)) :: datatype_name

    datatype_name = self%Datatype_Name

  END FUNCTION SEcategory_Name

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_Index
!
! PURPOSE:
!       Pure function to return the index of a particular surface type
!       in the SEcategory object.
!
! CALLING SEQUENCE:
!       idx = SEcategory_Index( SEcategory, Surface_Type )
!
! OBJECTS:
!       SEcategory:    Valid, allocated SEcategory object from which
!                      the surface type index is to be retrieved.
!                      UNITS:      N/A
!                      TYPE:       SEcategory_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! INPUTS:
!       Surface_Type:  Character string containing the name of the
!                      surface type for which the index is required.
!                      UNITS:      N/A
!                      TYPE:       CHARACTER(*)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       idx:           The index along the surface type dimension corresponding
!                      to the requested surface type name.
!                      If no surface type match is found, the value 0 is returned.
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE FUNCTION SEcategory_Index( &
    self         , &  ! Input
    Surface_Type ) &  ! Input
  RESULT( idx )
    ! Arguments
    TYPE(SEcategory_type), INTENT(IN) :: self
    CHARACTER(*)         , INTENT(IN) :: Surface_Type
    ! Function result
    INTEGER :: idx
    ! Local variables
    INTEGER :: i

    ! Setup
    idx = 0
    IF ( .NOT. SEcategory_Associated(self) ) RETURN

    ! Match surface type and assign
    DO i = 1, self%n_Surface_Types
      IF ( self%Surface_Type(i) == Surface_Type ) THEN
        idx = i
        RETURN
      END IF
    END DO

  END FUNCTION SEcategory_Index

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_DefineVersion
!
! PURPOSE:
!       Subroutine to return the module version information.
!
! CALLING SEQUENCE:
!       CALL SEcategory_DefineVersion( Id )
!
! OUTPUTS:
!       Id:    Character string containing the version Id information
!              for the module.
!              UNITS:      N/A
!              TYPE:       CHARACTER(*)
!              DIMENSION:  Scalar
!              ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE SEcategory_DefineVersion( Id )
    CHARACTER(*), INTENT(OUT) :: Id
    Id = MODULE_VERSION_ID
  END SUBROUTINE SEcategory_DefineVersion

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_SetValue
!
! PURPOSE:
!       Subroutine to set the contents of a valid SEcategory object.
!
! CALLING SEQUENCE:
!       CALL SEcategory_SetValue( SEcategory, &
!                                 Version              = Version             , &
!                                 Classification_Name  = Classification_Name , &
!                                 Frequency            = Frequency           , &
!                                 Surface_Type         = Surface_Type        , &
!                                 Surface_Type_IsValid = Surface_Type_IsValid, &
!                                 Reflectance          = Reflectance           )
!
! OBJECTS:
!       SEcategory:           Valid, allocated SEcategory object for which
!                             values are to be set.
!                             UNITS:      N/A
!                             TYPE:       SEcategory_type
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: INTENT(IN OUT)
!
! OPTIONAL INPUTS:
!       Version:              Integer indicating the data version. If not specified
!                             the value of the module parameter SECATEGORY_VERSION
!                             is used.
!                             UNITS:      N/A
!                             TYPE:       INTEGER
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Classification_Name:  String identifying the classification system used
!                             for the surface types.
!                             UNITS:      N/A
!                             TYPE:       INTEGER
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Frequency:            Real array to which the Frequency component of the
!                             SEcategory object is to be set. The size of the
!                             input must match the allocated size of the component,
!                             otherwise all the component values are set to zero.
!                             UNITS:      N/A
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Rank-1 (L)
!                             ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Surface_Type:         Character array to which the Surface_Type component
!                             of the SEcategory object is to be set. The size of the
!                             input must match the allocated size of the component,
!                             otherwise all the component values are set to a blank string.
!                             UNITS:      N/A
!                             TYPE:       CHARACTER(*)
!                             DIMENSION:  Rank-1 (N)
!                             ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Surface_Type_IsValid: Logical array to specify if a particular surface type is
!                             valid for the context in which it is to be used.
!                             UNITS:      N/A
!                             TYPE:       LOGICAL
!                             DIMENSION:  Rank-1 (N)
!                             ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Reflectance:          Real array to which the Reflectance component of the
!                             SEcategory object is to be set. The size of the
!                             input must match the allocated size of the component,
!                             otherwise all the component values are set to zero.
!                             UNITS:      N/A
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Rank-2 (L x N)
!                             ATTRIBUTES: INTENT(IN), OPTIONAL
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE SEcategory_SetValue( &
    self                , &  ! Input
    Version             , &  ! Optional input
    Classification_Name , &  ! Optional input
    Frequency           , &  ! Optional input
    Surface_Type        , &  ! Optional input
    Surface_Type_IsValid, &  ! Optional input
    Reflectance           )  ! Optional input
    ! Arguments
    TYPE(SEcategory_type) , INTENT(IN OUT) :: self
    INTEGER     , OPTIONAL, INTENT(IN)     :: Version
    CHARACTER(*), OPTIONAL, INTENT(IN)     :: Classification_Name
    REAL(fp)    , OPTIONAL, INTENT(IN)     :: Frequency(:)
    CHARACTER(*), OPTIONAL, INTENT(IN)     :: Surface_Type(:)
    LOGICAL     , OPTIONAL, INTENT(IN)     :: Surface_Type_IsValid(:)
    REAL(fp)    , OPTIONAL, INTENT(IN)     :: Reflectance(:,:)

    IF ( .NOT. SEcategory_Associated(self) ) RETURN

    IF ( PRESENT(Version) ) self%Version = Version
    IF ( PRESENT(Classification_Name) ) self%Classification_Name = Classification_Name

    IF ( PRESENT(Frequency) ) THEN
      IF ( SIZE(Frequency) == self%n_Frequencies ) THEN
        self%Frequency = Frequency
      ELSE
        self%Frequency = ZERO
      END IF
    END IF

    IF ( PRESENT(Surface_Type) ) THEN
      IF ( SIZE(Surface_Type) == self%n_Surface_Types ) THEN
        self%Surface_Type = Surface_Type
      ELSE
        self%Surface_Type = ''
      END IF
    END IF

    IF ( PRESENT(Surface_Type_IsValid) ) THEN
      IF ( SIZE(Surface_Type_IsValid) == self%n_Surface_Types ) THEN
        self%Surface_Type_IsValid = Surface_Type_IsValid
      ELSE
        self%Surface_Type_IsValid = .FALSE.
      END IF
    END IF

    IF ( PRESENT(Reflectance) ) THEN
      IF ( SIZE(Reflectance,DIM=1) == self%n_Frequencies .AND. &
           SIZE(Reflectance,DIM=2) == self%n_Surface_Types ) THEN
        self%Reflectance = Reflectance
      ELSE
        self%Reflectance = ZERO
      END IF
    END IF

  END SUBROUTINE SEcategory_SetValue

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_GetValue
!
! PURPOSE:
!       Subroutine to get the contents of a valid SEcategory object.
!
! CALLING SEQUENCE:
!       CALL SEcategory_GetValue( SEcategory, &
!                                 Surface_Type_ToGet   = Surface_Type_ToGet  , &
!                                 Version              = Version             , &
!                                 Classification_Name  = Classification_Name , &
!                                 n_Frequencies        = n_Frequencies       , &
!                                 n_Surface_Types      = n_Surface_Types     , &
!                                 Frequency            = Frequency           , &
!                                 Surface_Type         = Surface_Type        , &
!                                 Surface_Type_IsValid = Surface_Type_IsValid, &
!                                 Reflectance          = Reflectance         , &
!                                 Surface_Reflectance  = Surface_Reflectance   )
!
! OBJECTS:
!       SEcategory:              Valid, allocated SEcategory object from which
!                                values are to be retrieved.
!                                UNITS:      N/A
!                                TYPE:       SEcategory_type
!                                DIMENSION:  Scalar
!                                ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       Surface_Type_ToGet:      Character string containing a valid surface type
!                                name in the SEcategory object.
!                                NOTE: - This argument is used in conjuction with
!                                        the Surface_Reflectance dummy output
!                                        argument to retrieve the reflectance of a
!                                        particular surface type.
!                                      - This argument is ignored if the optional
!                                        Surface_Reflectance argument is not also
!                                        provided.
!                                UNITS:      N/A
!                                TYPE:       CHARACTER(*)
!                                DIMENSION:  Scalar
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! OPTIONAL OUTPUTS:
!       Version:                 Integer indicating the data version of the object.
!                                UNITS:      N/A
!                                TYPE:       INTEGER
!                                DIMENSION:  Scalar
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Classification_Name:     String identifying the classification system used
!                                for the surface types.
!                                UNITS:      N/A
!                                TYPE:       INTEGER
!                                DIMENSION:  Scalar
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Frequency:               Real array to which the Frequency component of the
!                                SEcategory object will be assigned. The actual
!                                argument must be declared as allocatable.
!                                UNITS:      N/A
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Rank-1 (L)
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL, ALLOCATABLE
!
!       Surface_Type:            Character array to which the Surface_Type component
!                                of the SEcategory object will be assigned. The actual
!                                argument must be declared as allocatable.
!                                UNITS:      N/A
!                                TYPE:       CHARACTER(*)
!                                DIMENSION:  Rank-1 (N)
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL, ALLOCATABLE
!
!       Surface_Type_IsValid:    Logical array that specifies if a particular surface
!                                type is valid for the context in which it is to be used.
!                                UNITS:      N/A
!                                TYPE:       LOGICAL
!                                DIMENSION:  Rank-1 (N)
!                                ATTRIBUTES: INTENT(IN), OPTIONAL, ALLOCATABLE
!
!       Reflectance:             Real array to which the Reflectance component of the
!                                SEcategory object will be assigned. The actual
!                                argument must be declared as allocatable.
!                                UNITS:      N/A
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Rank-2 (L x N)
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL, ALLOCATABLE
!
!       Surface_Reflectance:     Real array to which the Reflectance component for a
!                                given surface type in SEcategory object will be
!                                assigned. The actual argument must be declared as
!                                allocatable.
!                                NOTE: - This argument is used in conjuction with
!                                        the Surface_Type_ToGet dummy input
!                                        argument to retrieve the reflectance of a
!                                        particular surface type.
!                                      - This argument is ignored if the optional
!                                        Surface_Type_ToGet argument is not
!                                        also provided.
!                                UNITS:      N/A
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Rank-1 (L)
!                                ATTRIBUTES: INTENT(OUT), OPTIONAL, ALLOCATABLE
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE SEcategory_GetValue( &
    self                , &  ! Input
    Surface_Type_ToGet  , &  ! Optional input
    Version             , &  ! Optional output
    Classification_Name , &  ! Optional output
    n_Frequencies       , &  ! Optional output
    n_Surface_Types     , &  ! Optional output
    Frequency           , &  ! Optional output
    Surface_Type        , &  ! Optional output
    Surface_Type_IsValid, &  ! Optional output
    Reflectance         , &  ! Optional output
    Surface_Reflectance   )  ! Optional output
    ! Arguments
    TYPE(SEcategory_type),               INTENT(IN)  :: self
    CHARACTER(*),              OPTIONAL, INTENT(IN)  :: Surface_Type_ToGet
    INTEGER     ,              OPTIONAL, INTENT(OUT) :: Version
    CHARACTER(*),              OPTIONAL, INTENT(OUT) :: Classification_Name
    INTEGER     ,              OPTIONAL, INTENT(OUT) :: n_Frequencies
    INTEGER     ,              OPTIONAL, INTENT(OUT) :: n_Surface_Types
    REAL(fp)    , ALLOCATABLE, OPTIONAL, INTENT(OUT) :: Frequency(:)
    CHARACTER(*), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: Surface_Type(:)
    LOGICAL     , ALLOCATABLE, OPTIONAL, INTENT(OUT) :: Surface_Type_IsValid(:)
    REAL(fp)    , ALLOCATABLE, OPTIONAL, INTENT(OUT) :: Reflectance(:,:)
    REAL(fp)    , ALLOCATABLE, OPTIONAL, INTENT(OUT) :: Surface_Reflectance(:)
    ! Local variables
    INTEGER :: i

    IF ( .NOT. SEcategory_Associated(self) ) RETURN

    IF ( PRESENT(Version            ) ) Version             = self%Version
    IF ( PRESENT(Classification_Name) ) Classification_Name = self%Classification_Name
    IF ( PRESENT(n_Frequencies      ) ) n_Frequencies       = self%n_Frequencies
    IF ( PRESENT(n_Surface_Types    ) ) n_Surface_Types     = self%n_Surface_Types

    IF ( PRESENT(Frequency) ) THEN
      ALLOCATE(Frequency(self%n_Frequencies))
      Frequency = self%Frequency
    END IF

    IF ( PRESENT(Surface_Type) ) THEN
      ALLOCATE(Surface_Type(self%n_Surface_Types))
      Surface_Type = self%Surface_Type
    END IF

    IF ( PRESENT(Surface_Type_IsValid) ) THEN
      ALLOCATE(Surface_Type_IsValid(self%n_Surface_Types))
      Surface_Type_IsValid = self%Surface_Type_IsValid
    END IF

    IF ( PRESENT(Reflectance) ) THEN
      ALLOCATE(Reflectance(self%n_Frequencies, self%n_Surface_Types))
      Reflectance = self%Reflectance
    END IF

    IF ( PRESENT(Surface_Type_ToGet) .AND. PRESENT(Surface_Reflectance) ) THEN
      ! Match surface type and assign
      i = SEcategory_Index(self,Surface_Type_ToGet)
      IF ( i > 0 ) THEN
        ALLOCATE(Surface_Reflectance(self%n_Frequencies))
        Surface_Reflectance = self%Reflectance(:,i)
      END IF
    END IF

  END SUBROUTINE SEcategory_GetValue

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_InquireFile
!
! PURPOSE:
!       Function to inquire SEcategory object files.
!
! CALLING SEQUENCE:
!       Error_Status = SEcategory_InquireFile( &
!                        Filename                         , &
!                        n_Frequencies   = n_Frequencies  , &
!                        n_Surface_Types = n_Surface_Types, &
!                        Release         = Release        , &
!                        Version         = Version          )
!
! INPUTS:
!       Filename:           Character string specifying the name of the
!                           data file to inquire.
!                           UNITS:      N/A
!                           TYPE:       CHARACTER(*)
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(IN)
!
! OPTIONAL OUTPUTS:
!       n_Frequencies:      Number of spectral frequencies for which there are
!                           reflectance data.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       n_Surface_Types:    Number of land surface types for which is are
!                           reflectance data.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Release:            The data/file release number. Used to check
!                           for data/software mismatch.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Version:            The data/file version number. Used for
!                           purposes only in identifying the dataset for
!                           a particular release.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           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 file inquire was successful
!                              == FAILURE an unrecoverable error occurred.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION SEcategory_InquireFile( &
    Filename       , &  ! Input
    n_Frequencies  , &  ! Optional output
    n_Surface_Types, &  ! Optional output
    Release        , &  ! Optional output
    Version        , &  ! Optional output
    Title          , &  ! Optional output
    History        , &  ! Optional output
    Comment        ) &  ! Optional output
  RESULT( err_stat )
    ! Arguments
    CHARACTER(*),           INTENT(IN)  :: Filename
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Frequencies
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Surface_Types
    INTEGER     , OPTIONAL, INTENT(OUT) :: Release
    INTEGER     , OPTIONAL, INTENT(OUT) :: Version
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Title
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: History
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Comment
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'SEcategory_InquireFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat
    INTEGER :: fid
    TYPE(SEcategory_type) :: SEcategory

    ! Setup
    err_stat = SUCCESS
    ! ...Check that the file exists
    IF ( .NOT. File_Exists( Filename ) ) THEN
      msg = 'File '//TRIM(Filename)//' not found.'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Open the file
    err_stat = Open_Binary_File( Filename, fid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error opening '//TRIM(Filename)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read and check the datatype name
    err_stat = Read_Datatype( fid, SEcategory%Datatype_name )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading Datatype_Name'
      CALL Inquire_Cleanup(); RETURN
    END IF
    IF ( TRIM(SEcategory%Datatype_Name) /= SECATEGORY_DATATYPE ) THEN
      msg = SECATEGORY_DATATYPE//' datatype name check failed.'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the release and version
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Release, &
      SEcategory%Version
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading Release/Version - '//TRIM(io_msg)
      CALL Inquire_Cleanup(); RETURN
    END IF
    IF ( .NOT. SEcategory_ValidRelease( SEcategory ) ) THEN
      msg = 'SEcategory Release check failed.'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the dimensions
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%n_Frequencies  , &
      SEcategory%n_Surface_Types
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimension values from '//TRIM(Filename)//' - '//TRIM(io_msg)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the global attributes
    err_stat = ReadGAtts_Binary_File( &
                 fid, &
                 Title   = Title  , &
                 History = History, &
                 Comment = Comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading global attributes'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Close the file
    CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
    IF ( io_stat /= 0 ) THEN
      msg = 'Error closing '//TRIM(Filename)//' - '//TRIM(io_msg)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Assign the return arguments
    IF ( PRESENT(n_Frequencies  ) ) n_Frequencies   = SEcategory%n_Frequencies
    IF ( PRESENT(n_Surface_Types) ) n_Surface_Types = SEcategory%n_Surface_Types
    IF ( PRESENT(Release        ) ) Release         = SEcategory%Release
    IF ( PRESENT(Version        ) ) Version         = SEcategory%Version

  CONTAINS

    SUBROUTINE Inquire_CleanUp()
      ! Close file if necessary
      IF ( File_Open(fid) ) THEN
        CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
        IF ( io_stat /= 0 ) &
          msg = TRIM(msg)//'; Error closing input file during error cleanup - '//TRIM(io_msg)
      END IF
      ! Set error status and print error message
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
    END SUBROUTINE Inquire_CleanUp

  END FUNCTION SEcategory_InquireFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_ReadFile
!
! PURPOSE:
!       Function to read SEcategory object files.
!
! CALLING SEQUENCE:
!       Error_Status = SEcategory_ReadFile( &
!                        SEcategory         , &
!                        Filename           , &
!                        No_Close = No_Close, &
!                        Quiet    = Quiet     )
!
! OBJECTS:
!       SEcategory:     SEcategory object containing the data read from file.
!                       UNITS:      N/A
!                       TYPE:       SEcategory_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(OUT)
!
! INPUTS:
!       Filename:       Character string specifying the name of a
!                       SEcategory data file to read.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       No_Close:       Set this logical argument to *NOT* close the datafile
!                       upon exiting this routine. This option is required if
!                       the SEcategory data is embedded within another file.
!                       If == .FALSE., File is closed upon function exit [DEFAULT].
!                          == .TRUE.,  File is NOT closed upon function exit
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Quiet:          Set this logical argument to suppress INFORMATION
!                       messages being printed to stdout
!                       If == .FALSE., INFORMATION messages are OUTPUT [DEFAULT].
!                          == .TRUE.,  INFORMATION messages are SUPPRESSED.
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       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 file read was successful
!                          == FAILURE, an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION SEcategory_ReadFile( &
    SEcategory, &  ! Output
    Filename  , &  ! Input
    No_Close  , &  ! Optional input
    Quiet     , &  ! Optional input
    Title     , &  ! Optional output
    History   , &  ! Optional output
    Comment   , &  ! Optional output
    Debug     ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(SEcategory_type),  INTENT(OUT) :: SEcategory
    CHARACTER(*),           INTENT(IN)  :: Filename
    LOGICAL,      OPTIONAL, INTENT(IN)  :: No_Close
    LOGICAL,      OPTIONAL, INTENT(IN)  :: Quiet
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Title
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: History
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Comment
    LOGICAL,      OPTIONAL, INTENT(IN)  :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'SEcategory_ReadFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    LOGICAL :: close_file
    LOGICAL :: noisy
    INTEGER :: io_stat
    INTEGER :: fid
    TYPE(SEcategory_type) :: dummy

    ! Setup
    err_stat = SUCCESS
    ! ...Check No_Close argument
    close_file = .TRUE.
    IF ( PRESENT(No_Close) ) close_file = .NOT. No_Close
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(Debug) ) THEN
      IF ( Debug ) noisy = .TRUE.
    END IF

    ! Check if the file is open.
    IF ( File_Open( Filename ) ) THEN
      ! ...Inquire for the logical unit number
      INQUIRE( FILE=Filename, NUMBER=fid )
      ! ...Ensure it's valid
      IF ( fid < 0 ) THEN
        msg = 'Error inquiring '//TRIM(Filename)//' for its FileID'
        CALL Read_CleanUp(); RETURN
      END IF
    ELSE
      ! ...Open the file if it exists
      IF ( File_Exists( Filename ) ) THEN
        err_stat = Open_Binary_File( Filename, fid )
        IF ( err_Stat /= SUCCESS ) THEN
          msg = 'Error opening '//TRIM(Filename)
          CALL Read_CleanUp(); RETURN
        END IF
      ELSE
        msg = 'File '//TRIM(Filename)//' not found.'
        CALL Read_CleanUp(); RETURN
      END IF
    END IF

    ! Read and check the datatype name
    err_stat = Read_Datatype( fid, dummy%Datatype_name )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading Datatype_Name'
      CALL Read_Cleanup(); RETURN
    END IF
    IF ( TRIM(dummy%Datatype_Name) /= SECATEGORY_DATATYPE ) THEN
      msg = SECATEGORY_DATATYPE//' datatype name check failed.'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read and check the release and version
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      dummy%Release, &
      dummy%Version
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading Release/Version - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    IF ( .NOT. SEcategory_ValidRelease( dummy ) ) THEN
      msg = 'SEcategory Release check failed.'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the dimensions
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      dummy%n_Frequencies  , &
      dummy%n_Surface_Types
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading data dimensions - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Allocate the object
    CALL SEcategory_Create( &
           SEcategory          , &
           dummy%n_Frequencies  , &
           dummy%n_Surface_Types  )
    IF ( .NOT. SEcategory_Associated( SEcategory ) ) THEN
      msg = 'SEcategory object allocation failed.'
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Explicitly assign the version number
    SEcategory%Version = dummy%Version

    ! Read the global attributes
    err_stat = ReadGAtts_Binary_File( &
                 fid, &
                 Title   = Title  , &
                 History = History, &
                 Comment = Comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading global attributes'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the surface classification name
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Classification_Name
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading classification name - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the coefficient data
    ! ...Read the surface type names
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Surface_Type
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading surface type names - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Read the surface type validity array
    err_stat = ReadLogical_Binary_file( fid, SEcategory%Surface_Type_IsValid )
    IF ( err_stat /= SUCCESs ) THEN
      msg = 'Error reading surface type validity array'
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Read the dimensional vectors
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Frequency
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimensional vectors - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Read the reflectance data
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Reflectance
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading reflectance data - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF

    ! Close the file
    IF ( close_file ) THEN
      CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
      IF ( io_stat /= 0 ) THEN
        msg = 'Error closing '//TRIM(Filename)//' - '//TRIM(io_msg)
        CALL Read_Cleanup(); RETURN
      END IF
    END IF

    ! Output an info message
     IF ( noisy ) THEN
       CALL SEcategory_Info( SEcategory, msg )
       CALL Display_Message( ROUTINE_NAME, 'FILE: '//TRIM(Filename)//'; '//TRIM(msg), INFORMATION )
     END IF

   CONTAINS

     SUBROUTINE Read_CleanUp()
       IF ( File_Open(Filename) ) THEN
         CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
         IF ( io_stat /= 0 ) &
           msg = TRIM(msg)//'; Error closing input file during error cleanup - '//TRIM(io_msg)
       END IF
       CALL SEcategory_Destroy( SEcategory )
       err_stat = FAILURE
       CALL Display_Message( ROUTINE_NAME, msg, err_stat )
     END SUBROUTINE Read_CleanUp

  END FUNCTION SEcategory_ReadFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       SEcategory_WriteFile
!
! PURPOSE:
!       Function to write SEcategory object files.
!
! CALLING SEQUENCE:
!       Error_Status = SEcategory_WriteFile( &
!                        SEcategory         , &
!                        Filename           , &
!                        No_Close = No_Close, &
!                        Quiet    = Quiet     )
!
! OBJECTS:
!       SEcategory:     SEcategory object containing the data to write to file.
!                       UNITS:      N/A
!                       TYPE:       SEcategory_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! INPUTS:
!       Filename:       Character string specifying the name of a
!                       SEcategory format data file to write.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       No_Close:       Set this logical argument to *NOT* close the datafile
!                       upon exiting this routine. This option is required if
!                       the SEcategory data is to be embedded within another file.
!                       If == .FALSE., File is closed upon function exit [DEFAULT].
!                          == .TRUE.,  File is NOT closed upon function exit
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Quiet:          Set this logical argument to suppress INFORMATION
!                       messages being printed to stdout
!                       If == .FALSE., INFORMATION messages are OUTPUT [DEFAULT].
!                          == .TRUE.,  INFORMATION messages are SUPPRESSED.
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       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 file write was successful
!                          == FAILURE, an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION SEcategory_WriteFile( &
    SEcategory, &  ! Input
    Filename  , &  ! Input
    No_Close  , &  ! Optional input
    Quiet     , &  ! Optional input
    Title     , &  ! Optional input
    History   , &  ! Optional input
    Comment   , &  ! Optional input
    Debug     ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(SEcategory_type),  INTENT(IN) :: SEcategory
    CHARACTER(*),           INTENT(IN) :: Filename
    LOGICAL,      OPTIONAL, INTENT(IN) :: No_Close
    LOGICAL,      OPTIONAL, INTENT(IN) :: Quiet
    CHARACTER(*), OPTIONAL, INTENT(IN) :: Title
    CHARACTER(*), OPTIONAL, INTENT(IN) :: History
    CHARACTER(*), OPTIONAL, INTENT(IN) :: Comment
    LOGICAL,      OPTIONAL, INTENT(IN) :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'SEcategory_WriteFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    LOGICAL :: close_file
    LOGICAL :: noisy
    INTEGER :: io_stat
    INTEGER :: fid

    ! Setup
    err_stat = SUCCESS
    ! ...Check No_Close argument
    close_file = .TRUE.
    IF ( PRESENT(No_Close) ) close_file = .NOT. No_Close
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(Debug) ) THEN
      IF ( Debug ) noisy = .TRUE.
    END IF
    ! ...Check there is data to write
    IF ( .NOT. SEcategory_Associated( SEcategory ) ) THEN
      msg = 'SEcategory object is empty.'
      CALL Write_Cleanup(); RETURN
    END IF

    ! Check if the file is open.
    IF ( File_Open( FileName ) ) THEN
      ! ...Inquire for the logical unit number
      INQUIRE( FILE=Filename, NUMBER=fid )
      ! ...Ensure it's valid
      IF ( fid < 0 ) THEN
        msg = 'Error inquiring '//TRIM(Filename)//' for its FileID'
        CALL Write_CleanUp(); RETURN
      END IF
    ELSE
      ! ...Open the file for output
      err_stat = Open_Binary_File( Filename, fid, For_Output=.TRUE. )
      IF ( err_Stat /= SUCCESS ) THEN
        msg = 'Error opening '//TRIM(Filename)
        CALL Write_CleanUp(); RETURN
      END IF
    END IF

    ! Write the datatype name
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      LEN(SEcategory%Datatype_Name)
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing Datatype_Name length - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Datatype_Name
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing Datatype_Name - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the release and version
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Release, &
      SEcategory%Version
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing Release/Version - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the dimensions
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%n_Frequencies  , &
      SEcategory%n_Surface_Types
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing data dimensions - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the global attributes
    err_stat = WriteGAtts_Binary_File( &
                 fid, &
                 Write_Module = MODULE_VERSION_ID, &
                 Title        = Title  , &
                 History      = History, &
                 Comment      = Comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing global attributes'
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the surface classification name
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Classification_Name
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing classification name - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the coefficient data
    ! ...Write the surface type names
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Surface_Type
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing surface type names - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF
    ! ...Write the surface type validity array
    err_stat = WriteLogical_Binary_file( fid, SEcategory%Surface_Type_IsValid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing surface type validity array'
      CALL Write_Cleanup(); RETURN
    END IF
    ! ...Write the dimensional vectors
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Frequency
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing dimensional vectors - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF
    ! ...Write the reflectance data
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      SEcategory%Reflectance
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing reflectance data - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Close the file
    IF ( close_file ) THEN
      CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
      IF ( io_stat /= 0 ) THEN
        msg = 'Error closing '//TRIM(Filename)//' - '//TRIM(io_msg)
        CALL Write_Cleanup(); RETURN
      END IF
    END IF

    ! Output an info message
     IF ( noisy ) THEN
       CALL SEcategory_Info( SEcategory, msg )
       CALL Display_Message( ROUTINE_NAME, 'FILE: '//TRIM(Filename)//'; '//TRIM(msg), INFORMATION )
     END IF

   CONTAINS

     SUBROUTINE Write_Cleanup()
       IF ( File_Open(Filename) ) THEN
         CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
         IF ( io_stat /= 0 ) &
           msg = TRIM(msg)//'; Error closing output file during error cleanup - '//TRIM(io_msg)
       END IF
       err_stat = FAILURE
       CALL Display_Message( ROUTINE_NAME, msg, err_stat )
     END SUBROUTINE Write_Cleanup

  END FUNCTION SEcategory_WriteFile

!################################################################################
!################################################################################
!##                                                                            ##
!##                          ## PRIVATE PROCEDURES ##                          ##
!##                                                                            ##
!################################################################################
!################################################################################

!--------------------------------------------------------------------------------
!
! NAME:
!       SEcategory_Equal
!
! PURPOSE:
!       Elemental function to test the equality of two SEcategory objects.
!       Used in OPERATOR(==) interface block.
!
! CALLING SEQUENCE:
!       is_equal = SEcategory_Equal( x, y )
!
!         or
!
!       IF ( x == y ) THEN
!         ...
!       END IF
!
! OBJECTS:
!       x, y:          Two SEcategory objects to be compared.
!                      UNITS:      N/A
!                      TYPE:       SEcategory_type
!                      DIMENSION:  Scalar or any rank
!                      ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       is_equal:      Logical value indicating whether the inputs are equal.
!                      UNITS:      N/A
!                      TYPE:       LOGICAL
!                      DIMENSION:  Same as inputs.
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION SEcategory_Equal( x, y ) RESULT( is_equal )
    TYPE(SEcategory_type), INTENT(IN) :: x, y
    LOGICAL :: is_equal

    ! Set up
    is_equal = .FALSE.

    ! Check the object association status
    IF ( (.NOT. SEcategory_Associated(x)) .OR. &
         (.NOT. SEcategory_Associated(y))      ) RETURN

   ! Check contents
    ! ...Release/version info
    IF ( (x%Release /= y%Release) .OR. &
         (x%Version /= y%Version) ) RETURN
    ! ...Classification name
    IF ( (x%Classification_Name /= y%Classification_Name) ) RETURN
    ! ...Dimensions
    IF ( (x%n_Frequencies   /= y%n_Frequencies   ) .OR. &
         (x%n_Surface_Types /= y%n_Surface_Types ) ) RETURN
    ! ...Arrays
    IF ( ALL(x%Frequency    .EqualTo. y%Frequency    ) .AND. &
         ALL(x%Surface_Type     ==    y%Surface_Type ) .AND. &
         ALL(x%Reflectance  .EqualTo. y%Reflectance  ) ) &
      is_equal = .TRUE.

  END FUNCTION SEcategory_Equal

  ! Function to read the datatype name from file

  FUNCTION Read_Datatype( fid, datatype_name ) RESULT( err_stat )
    ! Arguments
    INTEGER     , INTENT(IN)  :: fid
    CHARACTER(*), INTENT(OUT) :: datatype_name
    ! Function result
    INTEGER :: err_stat
    ! Local variables
    CHARACTER(1), ALLOCATABLE :: dummy(:)
    INTEGER :: i, strlen
    INTEGER :: io_stat
    INTEGER :: alloc_stat

    ! Set up
    err_stat = FAILURE
    datatype_name = ''

    ! Get the string length
    READ( fid, IOSTAT=io_stat ) strlen
    IF ( io_stat /= 0 ) RETURN

    ! Allocate dummy string array
    ALLOCATE( dummy(strlen), STAT=alloc_stat )
    IF ( alloc_stat /= 0 ) RETURN

    ! Read the string into the dummy array
    READ( fid, IOSTAT=io_stat ) dummy
    IF ( io_stat /= 0 ) RETURN

    ! Transfer array into string
    DO i = 1, MIN(strlen,LEN(datatype_name))
      datatype_name(i:i) = dummy(i)
    END DO

    ! Done
    err_stat = SUCCESS
  END FUNCTION Read_Datatype

END MODULE SEcategory_Define

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


! chk_writealloc.f90 --
!     Check: Does the compiler allow you to write a derived type with an allocatable component
!     via list-directed output?
!
!     Note: this is not allowed according to gfortran and Intel Fortran, but PG Fortran
!     does allow it. Consider this an extension!
!
program chk_writealloc
    implicit none

    type mydata
        character(len=:), allocatable :: string
    end type mydata

    type(mydata) :: data

    data%string = 'A simple string'

    write(*,*) data

    write(*,*) 'Apparently the compiler allows this as an extension'
end program chk_writealloc


text:
      subroutine gg_ww_int(p,msq)
      implicit none
c--- Author: R. K. Ellis, May 2011
c--- For now, work in the approximation of two massless isodoublets
c--- Box contributions are then complete
c--- Triangle (vector) pieces always vanish
c--- Triangle (axial) pieces cancel for massless isodoublets

      include 'constants.f'
      include 'ewcouple.f'
      include 'qcdcouple.f'
      include 'masses.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      include 'Higgsint.f'
      include 'scale.f'
      include 'noglue.f'
      include 'anom_higgs.f'
      include 'process.f'
      include 'docheck.f'
      include 'first.f'
      integer h1,h2,nu,i,j,k,om,del1,del2,k12h,k34h,k56h11,k34h11,e
      double precision p(mxpart,4),msq(fn:nf,fn:nf),msqgg,fac
      double precision mfsq,tau,tauinv,rt,pttwo,rescale
      double complex Avec(2,2),Ahiggs(2,2),Agen3(2,2),Atot(2,2),
     & faccont,fachiggs,amphiggs,f,e3De4,sum(2,2,-2:0)

      double precision phi,muk,rho,ssig,csig,theta,
     & p1true(4),p2true(4),p3true(4),p4true(4),p5true(4),p6true(4),
     & dot,s12,s34,s56,dot1256,afac,bfac,gden,delta,
     & dot1234,dot3456,pmax,ptWsafetycut_massive,ptWsafetycut_massless
      double complex fvs,fvf,box(2,2,-2:0),triang(2,2,-2:0),
     & bub(2,2,-2:0)
      character*9 pp,pm
      logical includegens1and2,includegen3,
     & caseggWW4l,caseHWWHpI,caseHWWint
      parameter(pp='q+qb-g+g+',pm='q+qb-g+g-')
      parameter(del1=7,del2=8)
      parameter(k12h=9,k34h=10,k56h11=11,k34h11=12)
      save caseggWW4l,caseHWWHpI,caseHWWint
!$omp threadprivate(caseggWW4l,caseHWWHpI,caseHWWint)

c--- set this to true to include generations 1 and 2 of (light) quarks
      includegens1and2=.true.
c--- set this to true to include 3rd generation of quarks (t,b)
      includegen3=.true.

c--- omit massive loops for pt(W) < "ptWsafetycut_massive" (for num. stability)
      ptWsafetycut_massive=2d0

c--- omit massless loops for pt(W) < "ptWsafetycut_massless" (for num. stability)
      ptWsafetycut_massless=0.05d0

      if (first) then
        write(6,*)'****************************************************'
        write(6,*)'*                                                  *'
        if (includegens1and2) then
        write(6,*)'*  gg->WW box loop includes gens. 1 and 2          *'
        else
        write(6,*)'*  gg->WW box loop does not include gens. 1 and 2  *'
        endif
        if (includegen3) then
        write(6,*)'*  gg->WW box loop includes 3rd generation         *'
        else
        write(6,*)'*  gg->WW box loop does not include 3rd gen.       *'
        endif
        write(6,*)'*                                                  *'
        write(6,*)'****************************************************'
        first=.false.
        if (includegen3) call qlinit
        caseggWW4l=.false.
        caseHWWHpI=.false.
        caseHWWint=.false.
        if ((case .eq. 'ggWW4l')
     & .or. (case .eq. 'WWqqbr')) caseggWW4l=.true.
        if (case .eq. 'HWWH+i') caseHWWHpI=.true.
        if (case .eq. 'HWWint') caseHWWint=.true.
      endif

c--- if neither contribution is included print warning message and stop
      if ((includegens1and2 .eqv. .false.) .and.
     &    (includegen3      .eqv. .false.)) then
         write(6,*) 'Box loop is set to zero, please edit gg_WW_int.f'
         stop
      endif
c--- if noglue print warning message and stop
      if (noglue) then
         write(6,*) 'Please set noglue .false. in input file'
         stop
      endif

c--- logical variable "docheck"
c---   .true .  --> print out coefficients of integrals at special point
c---   .false.  --> run as normal
      docheck=.false.

c--- set flag to signal calculation of Higgs related contributions only,
c---  e.g. entering interference
      if     (caseggWW4l) then
        Higgsint=.false.
      elseif ((caseHWWHpi) .or. (caseHWWint)) then
        Higgsint=.true.
      else
        write(6,*) 'Unexpected case in gg_WW_int: ',case
        stop
      endif

      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      enddo
      enddo

c--- rescale
      pmax=abs(p(1,4))
      do i=2,6
      if (abs(p(i,4)) .gt. pmax) pmax=abs(p(i,4))
      enddo
      pmax=1d0
      do i=1,6
      do j=1,4
      p(i,j)=p(i,j)/pmax
      enddo
      enddo
      mt=mt/pmax
      musq=musq/pmax**2

c--- These lines set up the point for the numerical check
       if (docheck) then
         Higgsint=.false.
         includegen3=.true.
         include 'kinpoint.f'
         mt=0.5d0
         hmass=5.6d0
         hwidth=3.21d0
         zmass=0.53342d0
         zwidth=0.3414523523d0
         do nu=1,4
         om=nu-1
         if (nu.eq.1) om=4
         p(1,om)=p1true(nu)
         p(2,om)=p2true(nu)
         p(3,om)=p3true(nu)
         p(4,om)=p4true(nu)
         p(5,om)=p5true(nu)
         p(6,om)=p6true(nu)
         enddo
      endif

c--- if computing 3rd generation, set up extra flattened vectors
      if (includegen3) then
C--- define flattened vectors (k12 and k56)
      s12=2d0*dot(p,1,2)
      s56=2d0*dot(p,5,6)
      s34=2d0*dot(p,3,4)
      dot1256=0.5d0*(s34-s12-s56)
      delta=dot1256**2-s12*s56
      gden=dot1256+sqrt(delta)
      afac=s12/gden
      bfac=s56/gden

      do nu=1,4
      p(del2,nu)=one/(one-afac*bfac)
     & *(p(1,nu)+p(2,nu)-afac*(p(5,nu)+p(6,nu)))
      p(del1,nu)=one/(one-afac*bfac)
     & *(p(5,nu)+p(6,nu)-bfac*(p(1,nu)+p(2,nu)))
      enddo

C--- define flattened vectors (k12 and k34)
      dot1234=0.5d0*(s56-s12-s34)
      delta=dot1234**2-s12*s34
      gden=dot1234+sqrt(delta)
      afac=s12/gden
      bfac=s34/gden

      do nu=1,4
      p(k12h,nu)=one/(one-afac*bfac)
     & *(p(1,nu)+p(2,nu)-afac*(p(3,nu)+p(4,nu)))
      p(k34h,nu)=one/(one-afac*bfac)
     & *(p(3,nu)+p(4,nu)-bfac*(p(1,nu)+p(2,nu)))
      enddo

C--- define flattened vectors (k56 and k34)
      dot3456=0.5d0*(s12-s56-s34)
      delta=dot3456**2-s56*s34
      gden=dot3456+sqrt(delta)
      afac=s56/gden
      bfac=s34/gden

      do nu=1,4
      p(k56h11,nu)=one/(one-afac*bfac)
     & *(p(5,nu)+p(6,nu)-afac*(p(3,nu)+p(4,nu)))
      p(k34h11,nu)=one/(one-afac*bfac)
     & *(p(3,nu)+p(4,nu)-bfac*(p(5,nu)+p(6,nu)))
      enddo

      endif
c--- end of 3rd generation initialization

c--- set up spinor products (including for flat vectors, for 3rd gen)
      if (docheck) then
        call spinorz(12,p,za,zb)        ! Use spinorz for numerical check
      else
        if (includegen3) then
          call spinoru(12,p,za,zb)
        else
          call spinoru(6,p,za,zb)
        endif
      endif

c--- fill amplitudes with contributions of continuum W pairs
c--- note: 1 = minus, 2=plus

c      Avec(2,2)=a64v(mm,3,4,1,2,6,5,zb,za)*(-half*im)
c      Avec(2,1)=a64v(mp,3,4,1,2,6,5,zb,za)*(-half*im)
c      Avec(1,2)=a64v(pm,3,4,1,2,6,5,zb,za)*(-half*im)
c      Avec(1,1)=a64v(pp,3,4,1,2,6,5,zb,za)*(-half*im)

c      write(6,*) 'Avec(1,1)',Avec(1,1)*2d0
c      write(6,*) 'Avec(1,2)',Avec(1,2)*2d0
c      write(6,*) 'Avec(2,1)',Avec(2,1)*2d0
c      write(6,*) 'Avec(2,2)',Avec(2,2)*2d0
c      write(6,*)

C A simpler way of doing the same thing

c--- fill amplitudes used for generations 1 and 2
      if (includegens1and2) then
      Avec(2,2)=im*(fvs(pp,4,3,1,2,5,6,za,zb)+fvf(pp,4,3,1,2,5,6,za,zb))
      Avec(2,1)=im*(fvs(pm,4,3,1,2,5,6,za,zb)+fvf(pm,4,3,1,2,5,6,za,zb))
      Avec(1,2)=im*(fvs(pm,3,4,1,2,6,5,zb,za)+fvf(pm,3,4,1,2,6,5,zb,za))
      Avec(1,1)=im*(fvs(pp,3,4,1,2,6,5,zb,za)+fvf(pp,3,4,1,2,6,5,zb,za))
      else
        do h1=1,2
        do h2=1,2
        Avec(h1,h2)=czip
        enddo
        enddo
      endif

c      write(6,*) 'Avec(2,1)',Avec(1,1)
c      write(6,*) 'Avec(2,2)',Avec(1,2)
c      write(6,*) 'Avec(2,1)',Avec(2,1)
c      write(6,*) 'Avec(2,2)',Avec(2,2)

c--- factor two for complete massless isodoublets
      faccont=dcmplx(2d0)

c--- fill amplitudes used for 3rd generation
      if (includegen3) then
        do h1=1,2
        do h2=1,2
        do e=-2,0
        box(h1,h2,e)=czip
        triang(h1,h2,e)=czip
        bub(h1,h2,e)=czip
        enddo
        enddo
        enddo

c--- compute integrals and their coefficients
        if (docheck) then
          call massivebox(1,2,3,4,5,6,za,zb,box)
          call massivetri(1,2,3,4,5,6,za,zb,triang)
          call massivebub(1,2,3,4,5,6,za,zb,bub)
        else
          call massivebox6(1,2,3,4,5,6,za,zb,box)
          call massivetri6(1,2,3,4,5,6,za,zb,triang)
          call massivebub(1,2,3,4,5,6,za,zb,bub)
        endif

c--- This contribution is finite so we only retain "0" piece
        e=0
        do h1=1,2
        do h2=1,2
          sum(h1,h2,e)=box(h1,h2,e)+triang(h1,h2,e)+bub(h1,h2,e)
          Agen3(h1,h2)=sum(h1,h2,0)
        enddo
        enddo
      else
        do h1=1,2
        do h2=1,2
        Agen3(h1,h2)=czip
        enddo
        enddo
      endif

c--- rescale back
      do i=1,6
      do j=1,4
      p(i,j)=p(i,j)*pmax
      enddo
      enddo
      mt=mt*pmax
      musq=musq*pmax**2

c--- fill amplitudes with contributions of Higgs: top loop
      mfsq=mt**2
      tau=s(1,2)/(4d0*mfsq)
      tauinv=1d0/tau
      fachiggs=cone/dcmplx(s(1,2)-hmass**2,hmass*hwidth)
c--- shat-dependent width
c      fachiggs=cone/dcmplx(s(1,2)-hmass**2,hwidth*s(1,2)/hmass)
c--- piece proportional to hmass*hwidth
c      fachiggs=im*imag(fachiggs)
c--- piece proportional to (s(1,2)-hmass**2)
c      fachiggs=real(fachiggs)

      if (tau .le. 1d0) then
         f=dcmplx(dasin(sqrt(tau))**2)
      elseif (tau .gt. 1d0) then
         rt=sqrt(1d0-tauinv)
         f=-0.25d0*(dcmplx(log((1d0+rt)/(1d0-rt)))-im*pi)**2
      else
         f=czip
      endif
      e3De4=2d0*za(3,5)*zb(6,4)/(s(3,4)*s(5,6))
      amphiggs=mfsq*(cone+(cone-dcmplx(tauinv))*f)*im*e3De4
      Ahiggs(1,1)=fachiggs*amphiggs*za(1,2)/zb(2,1)
      Ahiggs(1,2)=czip
      Ahiggs(2,1)=czip
      Ahiggs(2,2)=fachiggs*amphiggs*zb(1,2)/za(2,1)

c--- fill amplitudes with contributions of Higgs: bottom loop
      mfsq=mb**2
      tau=s(1,2)/(4d0*mfsq)
      tauinv=1d0/tau

      if (tau .le. 1d0) then
         f=dcmplx(dasin(sqrt(tau))**2)
      elseif (tau .gt. 1d0) then
         rt=sqrt(1d0-tauinv)
         f=-0.25d0*(dcmplx(log((1d0+rt)/(1d0-rt)))-im*pi)**2
      else
         f=czip
      endif
      amphiggs=mfsq*(cone+(cone-dcmplx(tauinv))*f)*im*e3De4

      Ahiggs(1,1)=Ahiggs(1,1)+fachiggs*amphiggs*za(1,2)/zb(2,1)
      Ahiggs(2,2)=Ahiggs(2,2)+fachiggs*amphiggs*zb(1,2)/za(2,1)

      if (docheck) then
c--- numerical check includes top loops only
        Ahiggs(1,1)=Ahiggs(1,1)-fachiggs*amphiggs*za(1,2)/zb(2,1)
        Ahiggs(2,2)=Ahiggs(2,2)-fachiggs*amphiggs*zb(1,2)/za(2,1)
       do h1=1,2
        do h2=1,2
        write(6,*) 'h1,h2, massless',h1,h2,faccont*Avec(h1,h2)
        write(6,*) 'h1,h2, massive ',h1,h2,Agen3(h1,h2)
        write(6,*) 'h1,h2, Higgs   ',h1,h2,Ahiggs(h1,h2)
        write(6,*) 'h1,h2, ---- SUM',h1,h2,
     &               faccont*Avec(h1,h2)+Agen3(h1,h2)+Ahiggs(h1,h2)
        enddo
        enddo
c        pause
      endif

c--- ensure numerical stability: set massive loops to zero
c--- for pt(W) < "ptWsafetycut_massive" GeV
      if (pttwo(3,4,p) .lt. ptWsafetycut_massive) then
        do h1=1,2
        do h2=1,2
          Agen3(h1,h2)=czip
        enddo
        enddo
      endif

c--- ensure numerical stability: set massless loops to zero
c--- for pt(W) < "ptWsafetycut_massless" GeV
      if (pttwo(3,4,p) .lt. ptWsafetycut_massless) then
        do h1=1,2
        do h2=1,2
          Avec(h1,h2)=czip
        enddo
        enddo
      endif

c--- Rescale for width study
      if((keep_smhiggs_norm).and.(anom_higgs)) then
         rescale=chi_higgs**2
         Ahiggs(:,:)=Ahiggs(:,:)*rescale
      endif

      msqgg=0d0
      do h1=1,2
      do h2=1,2
      Atot(h1,h2)=faccont*Avec(h1,h2)+Agen3(h1,h2)+Ahiggs(h1,h2)

      if     (caseggWW4l) then
c--- This accumulates total contributions
        msqgg=msqgg+cdabs(Atot(h1,h2))**2
      elseif (caseHWWHpi) then
c--- This only accumulates contributions containing the Higgs diagram,
c---  i.e. the Higgs diagrams squared and the interference
        msqgg=msqgg+cdabs(Atot(h1,h2))**2
     &             -cdabs(faccont*Avec(h1,h2)+Agen3(h1,h2))**2
      elseif (caseHWWint) then
c--- This only accumulates the interference
        msqgg=msqgg+cdabs(Atot(h1,h2))**2
     &             -cdabs(faccont*Avec(h1,h2)+Agen3(h1,h2))**2
     &             -cdabs(Ahiggs(h1,h2))**2
      else
        write(6,*) 'Unexpected case in gg_WW_int: ',case
        stop
      endif

      enddo
      enddo

c--- overall factor from diagrams
      fac=avegg*V*(2d0*gwsq*gsq/(16d0*pisq)*gwsq/2d0)**2
     & *s(3,4)**2/((s(3,4)-wmass**2)**2+(wwidth*wmass)**2)
     & *s(5,6)**2/((s(5,6)-wmass**2)**2+(wwidth*wmass)**2)
      fac=fac/pmax**4

      msq(0,0)=msqgg*fac
      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


! -----------------------------------------------------------------
! Programmer(s): David J. Gardner @ LLNL
!                Daniel R. Reynolds @ SMU
! -----------------------------------------------------------------
! LLNS Copyright Start
! Copyright (c) 2014, Lawrence Livermore National Security
! This work was performed under the auspices of the U.S. Department
! of Energy by Lawrence Livermore National Laboratory in part under
! Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
! Produced at the Lawrence Livermore National Laboratory.
! All rights reserved.
! For details, see the LICENSE file.
! LLNS Copyright End
! -----------------------------------------------------------------
! This file contains a Fortran module for interfacing directly with
! the SUNDIALS serial NVector using the ISO_C_BINDING module.
! -----------------------------------------------------------------

module fnvector_serial_mod

  !======= Interfaces =========
  interface

     ! =================================================================
     ! Constructors
     ! =================================================================

     ! -----------------------------------------------------------------
     ! N_VNew_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VNew_Serial(vec_length) &
          bind(C,name='N_VNew_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       integer(c_long), value :: vec_length
     end function FN_VNew_Serial

     ! -----------------------------------------------------------------
     ! N_VNewEmpty_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VNewEmpty_Serial(vec_length) &
          bind(C,name='N_VNewEmpty_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       integer(c_long), value :: vec_length
     end function FN_VNewEmpty_Serial

     ! -----------------------------------------------------------------
     ! N_VMake_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VMake_Serial(length, v_data) &
          bind(C,name='N_VMake_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       integer(c_long), value :: length
       real(c_double)         :: v_data(length)
     end function FN_VMake_Serial

     ! -----------------------------------------------------------------
     ! N_VCloneVectorArray_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VCloneVectorArray_Serial(count, w) &
          bind(C,name='N_VCloneVectorArray_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       integer(c_int), value :: count
       type(c_ptr),    value :: w
     end function FN_VCloneVectorArray_Serial

     ! -----------------------------------------------------------------
     ! N_VCloneVectorArrayEmpty_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VCloneVectorArrayEmpty_Serial(count, w) &
          bind(C,name='N_VCloneVectorArrayEmpty_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       integer(c_int), value :: count
       type(c_ptr),    value :: w
     end function FN_VCloneVectorArrayEmpty_Serial

     ! =================================================================
     ! Destructors
     ! =================================================================

     ! -----------------------------------------------------------------
     ! N_VDestroy_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VDestroy_Serial(v) &
          bind(C,name='N_VDestroy_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: v
     end subroutine FN_VDestroy_Serial

     ! -----------------------------------------------------------------
     ! N_VDestroyVectorArray_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VDestroyVectorArray_Serial(vs, count) &
          bind(C,name='N_VDestroyVectorArray_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr)    :: vs
       integer(c_int) :: count
     end subroutine FN_VDestroyVectorArray_Serial

     ! =================================================================
     ! Other routines
     ! =================================================================

     ! -----------------------------------------------------------------
     ! N_VGetLength_Serial
     ! -----------------------------------------------------------------

     integer(c_long) function FN_VGetLength_Serial(v) &
          bind(C,name='N_VGetLength_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: v
     end function FN_VGetLength_Serial

     ! -----------------------------------------------------------------
     ! N_VPrint_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VPrint_Serial(v) &
          bind(C,name='N_VPrint_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: v
     end subroutine FN_VPrint_Serial

     ! -----------------------------------------------------------------
     ! NOT INTERFACED: N_VPrintFile_Serial
     ! -----------------------------------------------------------------

     ! =================================================================
     ! Operations
     ! =================================================================

     ! -----------------------------------------------------------------
     ! N_VGetVectorID: NOT INTERFACED
     ! -----------------------------------------------------------------

     ! -----------------------------------------------------------------
     ! N_VCloneEmpty_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VCloneEmpty_Serial(w) &
          bind(C,name='N_VCloneEmpty_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: w
     end function FN_VCloneEmpty_Serial

     ! -----------------------------------------------------------------
     ! N_VClone_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VClone_Serial(w) &
          bind(C,name='N_VClone_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: w
     end function FN_VClone_Serial

     ! -----------------------------------------------------------------
     ! N_VSpace_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VSpace_Serial(v, lrw, liw) &
          bind(C,name='N_VSpace_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: v
       integer(c_long)    :: lrw
       integer(c_long)    :: liw
     end subroutine FN_VSpace_Serial

     ! -----------------------------------------------------------------
     ! N_VGetArrayPointer_Serial
     ! -----------------------------------------------------------------

     type(c_ptr) function FN_VGetArrayPointer_Serial(vec) &
          bind(C,name='N_VGetArrayPointer_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: vec
     end function FN_VGetArrayPointer_Serial

     ! -----------------------------------------------------------------
     ! N_VSetArrayPointer_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VSetArrayPointer_Serial(v_data, v) &
          bind(C,name='N_VSetArrayPointer_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       real(c_double)     :: v_data
       type(c_ptr), value :: v
     end subroutine FN_VSetArrayPointer_Serial

     ! -----------------------------------------------------------------
     ! N_VLinearSum_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VLinearSum_Serial(a, x, b, y, z) &
          bind(C,name='N_VLinearSum_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       real(c_double), value :: a
       type(c_ptr),    value :: x
       real(c_double), value :: b
       type(c_ptr),    value :: y
       type(c_ptr),    value :: z
     end subroutine FN_VLinearSum_Serial

     ! -----------------------------------------------------------------
     ! N_VConst_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VConst_Serial(c, z) &
          bind(C,name='N_VConst_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       real(c_double), value :: c
       type(c_ptr),    value :: z
     end subroutine FN_VConst_Serial

     ! -----------------------------------------------------------------
     ! N_VProd_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VProd_Serial(x, y, z) &
          bind(C,name='N_VProd_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: y
       type(c_ptr), value :: z
     end subroutine FN_VProd_Serial

     ! -----------------------------------------------------------------
     ! N_VDiv_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VDiv_Serial(x, y, z) &
          bind(C,name='N_VDiv_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: y
       type(c_ptr), value :: z
     end subroutine FN_VDiv_Serial

     ! -----------------------------------------------------------------
     ! N_VScale_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VScale_Serial(c, x, z) &
          bind(C,name='N_VScale_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       real(c_double), value :: c
       type(c_ptr),    value :: x
       type(c_ptr),    value :: z
     end subroutine FN_VScale_Serial

     ! -----------------------------------------------------------------
     ! N_VAbs_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VAbs_Serial(x, z) &
          bind(C,name='N_VAbs_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: z
     end subroutine FN_VAbs_Serial

     ! -----------------------------------------------------------------
     ! N_VInv_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VInv_Serial(x, z) &
          bind(C,name='N_VInv_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: z
     end subroutine FN_VInv_Serial

     ! -----------------------------------------------------------------
     ! N_VAddConst
     ! -----------------------------------------------------------------

     subroutine FN_VAddConst_Serial(x, b, z) &
          bind(C,name='N_VAddConst_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: x
       real(c_double), value :: b
       type(c_ptr),    value :: z
     end subroutine FN_VAddConst_Serial

     ! -----------------------------------------------------------------
     ! N_VDotProd_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VDotProd_Serial(x, y) &
          bind(C,name='N_VDotProd_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: y
     end function FN_VDotProd_Serial

     ! -----------------------------------------------------------------
     ! N_VMaxNorm_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VMaxNorm_Serial(x) &
          bind(C,name='N_VMaxNorm_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
     end function FN_VMaxNorm_Serial

     ! -----------------------------------------------------------------
     ! N_VWrmsNorm_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VWrmsNorm_Serial(x, w) &
          bind(C,name='N_VWrmsNorm_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: w
     end function FN_VWrmsNorm_Serial

     ! -----------------------------------------------------------------
     ! N_VWrmsNormMask_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VWrmsNormMask_Serial(x, w, id) &
          bind(C,name='N_VWrmsNormMask_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: w
       type(c_ptr), value :: id
     end function FN_VWrmsNormMask_Serial

     ! -----------------------------------------------------------------
     ! N_VMin_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VMin_Serial(x) &
          bind(C,name='N_VMin_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
     end function FN_VMin_Serial

     ! -----------------------------------------------------------------
     ! N_VWL2Norm_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VWL2Norm_Serial(x, w) &
          bind(C,name='N_VWL2Norm_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: w
     end function FN_VWL2Norm_Serial

     ! -----------------------------------------------------------------
     ! N_VL1Norm_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VL1Norm_Serial(x) &
          bind(C,name='N_VL1Norm_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
     end function FN_VL1Norm_Serial

     ! -----------------------------------------------------------------
     ! N_VCompare_Serial
     ! -----------------------------------------------------------------

     subroutine FN_VCompare_Serial(c, x, z) &
          bind(C,name='N_VCompare_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       real(c_double), value :: c
       type(c_ptr),    value :: x
       type(c_ptr),    value :: z
     end subroutine FN_VCompare_Serial

     ! -----------------------------------------------------------------
     ! N_VInvTest_Serial
     ! -----------------------------------------------------------------

     integer(c_int) function FN_VInvTest_Serial(x, z) &
          bind(C,name='N_VInvTest_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x
       type(c_ptr), value :: z
     end function FN_VInvTest_Serial

     ! -----------------------------------------------------------------
     ! N_VConstrMask_Serial
     ! -----------------------------------------------------------------

     integer(c_int) function FN_VConstrMask_Serial(c, x, m) &
          bind(C,name='N_VConstrMask_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: c
       type(c_ptr), value :: x
       type(c_ptr), value :: m
     end function FN_VConstrMask_Serial

     ! -----------------------------------------------------------------
     ! N_VMinQuotient_Serial
     ! -----------------------------------------------------------------

     real(c_double) function FN_VMinQuotient_Serial(num, denom) &
          bind(C,name='N_VMinQuotient_Serial')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: num
       type(c_ptr), value :: denom
     end function FN_VMinQuotient_Serial

  end interface

contains

  ! ================================================================
  ! Helpful routines
  ! ================================================================

  ! ----------------------------------------------------------------
  ! FN_VGetData_Serial
  !
  ! Extracts data array from a serial SUNDIALS N_Vector
  ! ----------------------------------------------------------------

  subroutine FN_VGetData_Serial(vec, f_array)

    !======= Inclusions ===========
    use, intrinsic :: iso_c_binding

    !======= Declarations =========
    implicit none

    ! calling variables
    type(c_ptr)             :: vec
    integer(c_long)         :: length
    real(c_double), pointer :: f_array(:)

    ! C pointer for N_Vector interal data array
    type(c_ptr) :: c_array

    !======= Internals ============

    ! get data pointer from N_Vector
    c_array = FN_VGetArrayPointer_Serial(vec)

    ! get vector length
    length = FN_VGetLength_Serial(vec)

    ! convert c pointer to f pointer
    call c_f_pointer(c_array, f_array, (/length/))

  end subroutine FN_VGetData_Serial

end module fnvector_serial_mod


text:
	SUBROUTINE AW_PIDB ( bdata, irept, lenbyt, iff, iret )
C************************************************************************
C* AW_PIDB                                                              *
C*                                                                      *
C* This subroutine is the product identification block.  It identifies  *
C* the origin of the product, the classification, retention time,       *
C* product identifier, and file time.  This subroutine represents       *
C* mode 1, submode 1.                                                   *
C* 								        *
C* AW_PIDB ( BDATA, IREPT, LENBYT, IFF, IRET )                        	*
C*								        *
C* Input parameters: 						        *
C*      BDATA(*) BYTE		Byte data array                         *
C*      IREPT   INTEGER         Beginning byte position for the block   *
C*      LENBYT  INTEGER         Length of the block in bytes            *
C*      IFF     INTEGER         Length and checksum indicator           *
C*								        *
C* Output parameters:						        *
C*	IRET	INTEGER		Return code                             *
C**								        *
C* Log:									*
C* A. Hardy/GSC          7/98                                           *
C* T. Piper/GSC		11/98	Updated	prolog				*
C* A. Hardy/GSC         11/98   Added assignment do loop                *
C* A. Hardy/GSC		 3/99	Updated	prolog				*
C* A. Hardy/GSC          2/00   Increase the file byte size 20k -> 80k  *
C* A. Hardy/GSC         10/00   Fixed display;Put file number and dbug  *
C*         			in awcmn.cmn				*
C* T. Piper/SAIC	 2/03	Changed b to bdata and added INT	*
C************************************************************************
 	INCLUDE         'awcmn.cmn'
C*
	BYTE		bdata(*)
        CHARACTER       clas
	CHARACTER*40    orgid, orgcnt, prodid, flid
        INTEGER         icnt, irtm, imon, iday, ihour,
     +                  imin, itotln, iyear
        INTEGER         ibte(80000)
C------------------------------------------------------------------------
        IF ( dbug .eq. 'y' ) THEN
            write(flun,*)'In AW_PIDB'
	END IF
 	iret = 0
	icnt = irept + 3
	itotln = irept + lenbyt
	IF ( iff .eq. 0 ) itotln = itotln - 1
C
C*      Reassign byte array to an integer array.
C
        DO i = 1, 80000
            ibte(i) = bdata(i)
        END DO
C
C*      Get the 4 character originator identification, classification,
C*      and retention time.
C
	DO i = 1, 4
	    orgid( i:i) = CHAR(INT(bdata(icnt+i)))
	END DO
	icnt = icnt + 5
	clas =  CHAR(INT(bdata(icnt)))
        irtm =  ibits ( ibte(icnt+1), 0, 8 )
        IF ( dbug .eq. 'y' ) THEN
	    write(flun,11)'ORGINATING ID','CLASS','RETENTION TIME'
	    write(flun,10)orgid, clas, irtm
 11         format(a13,5x,a5,5x,a14)
 10         format(a6,25x,a,5x,i5)
	END IF
C
C*      Get the file indicator and product identifier.
C
        flid = CHAR(INT(bdata(icnt+2)))
	icnt = icnt + 2
	DO i = 1, 9
	    prodid(i:i) = CHAR(INT(bdata(icnt+i)))
	END DO
	icnt = icnt + 10
        IF ( dbug .eq. 'y' ) THEN
	    write(flun,29)'FILE INDICATOR','PRODUCT ID'
 29         format(a15,15x,a11)
	    write(flun,30)flid, prodid
 30         format(5x,a1,25x,a9)
	END IF
C
C*      Get the product file time: year, month, day, hour and min.
C
	CALL AW_ADBT (bdata, icnt, iyear, iret )
C
	imon   = ibits ( ibte(icnt+2), 0, 8 )
	iday   = ibits ( ibte(icnt+3), 0, 8 )
	ihour  = ibits ( ibte(icnt+4), 0, 8 )
	imin   = ibits ( ibte(icnt+5), 0, 8 )
	icnt = icnt + 5
        IF ( dbug .eq. 'y' ) THEN
            write(flun,31)'YEAR','MONTH','DAY','HOUR',' MIN'
 31         format(4(a5,1x),a5)
            write(flun,32)iyear, imon, iday, ihour, imin
 32         format(4(i5,1x),i5)
	END IF
C
C*      Get the rest of the product identifier, if there is more.
C
        IF ( icnt .ne. itotln ) THEN
	    DO i = 1, 6
	        orgcnt(i:i) = CHAR(INT(bdata( icnt+i)))
	    END DO
            IF ( dbug .eq. 'y' ) THEN
	        write(flun,*)'ORIGNATOR INDICATOR CONTINUED'
	        write(flun,33)orgcnt
 33             format(5x,a6)
	    END IF
        END IF
C
C*      Check for checksum.
C
        IF ( dbug .eq. 'y' )  THEN
            IF ( iff .eq. 0 ) write(flun,*)'Checksum found.'
        END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak