CRTM_RTSolution_Define.f90

Go to the documentation of this file.

!
! CRTM_RTSolution_Define
!
! Module defining the CRTM RTSolution structure and containing routines
! to manipulate it.
!
!
! CREATION HISTORY:
!   Written by:     Paul van Delst, 13-May-2004
!                   paul.vandelst@noaa.gov
!

MODULE crtm_rtsolution_define


  ! ------------------
  ! Environment set up
  ! ------------------
  ! Intrinsic modules
  USE iso_fortran_env      , ONLY: output_unit
  ! Module use statements
  USE type_kinds           , ONLY: fp
  USE message_handler      , ONLY: success, failure, warning, information, display_message
  USE compare_float_numbers, ONLY: default_n_sigfig, &
                                   OPERATOR(.equalto.), &
                                   compares_within_tolerance
  USE file_utility         , ONLY: file_open, file_exists
  USE binary_file_utility  , ONLY: open_binary_file      , &
                                   writegatts_binary_file, &
                                   readgatts_binary_file
  USE sensorinfo_parameters, ONLY: invalid_sensor, &
                                   invalid_wmo_satellite_id, &
                                   invalid_wmo_sensor_id
  USE crtm_parameters      , ONLY: strlen
  ! Disable all implicit typing
  IMPLICIT NONE


  ! --------------------
  ! Default visibilities
  ! --------------------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: crtm_rtsolution_type
  ! Operators
  PUBLIC :: OPERATOR(==)
  PUBLIC :: OPERATOR(+)
  PUBLIC :: OPERATOR(-)
  ! Public procedures
  PUBLIC :: crtm_rtsolution_associated
  PUBLIC :: crtm_rtsolution_destroy
  PUBLIC :: crtm_rtsolution_create
  PUBLIC :: crtm_rtsolution_zero
  PUBLIC :: crtm_rtsolution_inspect
  PUBLIC :: crtm_rtsolution_defineversion
  PUBLIC :: crtm_rtsolution_compare
  PUBLIC :: crtm_rtsolution_statistics
  PUBLIC :: crtm_rtsolution_inquirefile
  PUBLIC :: crtm_rtsolution_readfile
  PUBLIC :: crtm_rtsolution_writefile


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

  INTERFACE OPERATOR(+)
    MODULE PROCEDURE crtm_rtsolution_add
  END INTERFACE OPERATOR(+)

  INTERFACE OPERATOR(-)
    MODULE PROCEDURE crtm_rtsolution_subtract
  END INTERFACE OPERATOR(-)

  INTERFACE OPERATOR(**)
    MODULE PROCEDURE crtm_rtsolution_exponent
  END INTERFACE OPERATOR(**)

  INTERFACE OPERATOR(/)
    MODULE PROCEDURE crtm_rtsolution_normalise
  END INTERFACE OPERATOR(/)

  INTERFACE sqrt
    MODULE PROCEDURE crtm_rtsolution_sqrt
  END INTERFACE sqrt

  INTERFACE crtm_rtsolution_inspect
    MODULE PROCEDURE scalar_inspect
    MODULE PROCEDURE rank2_inspect
  END INTERFACE crtm_rtsolution_inspect


  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: module_version_id = &
  '$Id: CRTM_RTSolution_Define.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
  ! Literal constants
  REAL(fp), PARAMETER :: zero = 0.0_fp
  ! Message string length
  INTEGER, PARAMETER :: ml = 256
  ! File status on close after write error
  CHARACTER(*), PARAMETER :: write_error_status = 'DELETE'


  ! -------------------------------
  ! RTSolution data type definition
  ! -------------------------------
  !:tdoc+:
  TYPE :: crtm_rtsolution_type
    ! Allocation indicator
    LOGICAL :: is_allocated = .false.
    ! Dimensions
    INTEGER :: n_layers = 0  ! K
    ! Sensor information
    CHARACTER(STRLEN) :: sensor_id        = ''
    INTEGER           :: wmo_satellite_id = invalid_wmo_satellite_id
    INTEGER           :: wmo_sensor_id    = invalid_wmo_sensor_id
    INTEGER           :: sensor_channel   = 0
    ! RT algorithm information
    CHARACTER(STRLEN*5) :: rt_algorithm_name = ''
    ! Internal variables. Users do not need to worry about these.
    LOGICAL :: scattering_flag = .true.
    INTEGER :: n_full_streams  = 0
    INTEGER :: n_stokes        = 0
    ! Forward radiative transfer intermediate results for a single channel
    !    These components are not defined when they are used as TL, AD
    !    and K variables
    REAL(fp) :: sod                     = zero  ! Scattering Optical Depth
    REAL(fp) :: surface_emissivity      = zero
    REAL(fp) :: surface_reflectivity    = zero
    REAL(fp) :: up_radiance             = zero
    REAL(fp) :: down_radiance           = zero
    REAL(fp) :: down_solar_radiance     = zero
    REAL(fp) :: surface_planck_radiance = zero
    REAL(fp), ALLOCATABLE :: upwelling_overcast_radiance(:)   ! K
    REAL(fp), ALLOCATABLE :: upwelling_radiance(:)   ! K
    REAL(fp), ALLOCATABLE :: layer_optical_depth(:)  ! K
    ! Radiative transfer results for a single channel/node
    REAL(fp) :: radiance               = zero
    REAL(fp) :: brightness_temperature = zero
  END TYPE crtm_rtsolution_type
  !:tdoc-:


CONTAINS


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

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_Associated
!
! PURPOSE:
!   Elemental function to test the status of the allocatable components
!   of a CRTM RTSolution object.
!
! CALLING SEQUENCE:
!   Status = CRTM_RTSolution_Associated( RTSolution )
!
! OBJECTS:
!   RTSolution:   RTSolution structure which is to have its member's
!                 status tested.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar or any rank
!                 ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!   Status:       The return value is a logical value indicating the
!                 status of the RTSolution members.
!                   .TRUE.  - if the array components are allocated.
!                   .FALSE. - if the array components are not allocated.
!                 UNITS:      N/A
!                 TYPE:       LOGICAL
!                 DIMENSION:  Same as input RTSolution argument
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_associated( RTSolution ) RESULT( Status )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rtsolution
    LOGICAL :: status
    status = rtsolution%Is_Allocated
  END FUNCTION crtm_rtsolution_associated


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_Destroy
!
! PURPOSE:
!   Elemental subroutine to re-initialize CRTM RTSolution objects.
!
! CALLING SEQUENCE:
!   CALL CRTM_RTSolution_Destroy( RTSolution )
!
! OBJECTS:
!   RTSolution:   Re-initialized RTSolution structure.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar OR any rank
!                 ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE crtm_rtsolution_destroy( RTSolution )
    TYPE(crtm_rtsolution_type), INTENT(OUT) :: rtsolution
    rtsolution%Is_Allocated = .false.
    rtsolution%n_Layers = 0
  END SUBROUTINE crtm_rtsolution_destroy


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_Create
!
! PURPOSE:
!   Elemental subroutine to create an instance of the CRTM RTSolution object.
!
! CALLING SEQUENCE:
!   CALL CRTM_RTSolution_Create( RTSolution, n_Layers )
!
! OBJECTS:
!   RTSolution:   RTSolution structure.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar or any rank
!                 ATTRIBUTES: INTENT(OUT)
!
! INPUTS:
!   n_Layers:     Number of layers for which there is RTSolution data.
!                 Must be > 0.
!                 UNITS:      N/A
!                 TYPE:       INTEGER
!                 DIMENSION:  Same as RTSolution object
!                 ATTRIBUTES: INTENT(IN)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE crtm_rtsolution_create( RTSolution, n_Layers )
    ! Arguments
    TYPE(crtm_rtsolution_type), INTENT(OUT) :: rtsolution
    INTEGER,                    INTENT(IN)  :: n_layers
    ! Local variables
    INTEGER :: alloc_stat

    ! Check input
    IF ( n_layers < 1 ) RETURN

    ! Perform the allocation
    ALLOCATE( rtsolution%Upwelling_Radiance(n_layers), &
              rtsolution%Upwelling_Overcast_Radiance(n_layers), &
              rtsolution%Layer_Optical_Depth(n_layers), &
              stat = alloc_stat )
    IF ( alloc_stat /= 0 ) RETURN

    ! Initialise
    ! ...Dimensions
    rtsolution%n_Layers = n_layers
    ! ...Arrays
    rtsolution%Upwelling_Radiance  = zero
    rtsolution%Upwelling_Overcast_Radiance  = zero
    rtsolution%Layer_Optical_Depth = zero

    ! Set allocation indicator
    rtsolution%Is_Allocated = .true.

  END SUBROUTINE crtm_rtsolution_create


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_Zero
!
! PURPOSE:
!   Elemental subroutine to zero out the data components
!   in a CRTM RTSolution object.
!
! CALLING SEQUENCE:
!   CALL CRTM_RTSolution_Zero( rts )
!
! OUTPUTS:
!   rts:          CRTM RTSolution structure in which the data components
!                 are to be zeroed out.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar or any rank
!                 ATTRIBUTES: INTENT(IN OUT)
!
! COMMENTS:
!   - The dimension components of the structure are *NOT* set to zero.
!   - The sensor infomration and RT algorithm components are
!     *NOT* reset in this routine.
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE crtm_rtsolution_zero( RTSolution )
    TYPE(crtm_rtsolution_type), INTENT(IN OUT) :: rtsolution

    ! Zero out the scalar data components
    rtsolution%SOD                     = zero
    rtsolution%Surface_Emissivity      = zero
    rtsolution%Surface_Reflectivity    = zero
    rtsolution%Up_Radiance             = zero
    rtsolution%Down_Radiance           = zero
    rtsolution%Down_Solar_Radiance     = zero
    rtsolution%Surface_Planck_Radiance = zero
    rtsolution%Radiance                = zero
    rtsolution%Brightness_Temperature  = zero

    ! Zero out the array data components
    IF ( crtm_rtsolution_associated(rtsolution) ) THEN
      rtsolution%Upwelling_Radiance  = zero
      rtsolution%Upwelling_Overcast_Radiance  = zero
      rtsolution%Layer_Optical_Depth = zero
    END IF

  END SUBROUTINE crtm_rtsolution_zero


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_Inspect
!
! PURPOSE:
!   Subroutine to print the contents of a CRTM RTSolution object to stdout.
!
! CALLING SEQUENCE:
!   CALL CRTM_RTSolution_Inspect( RTSolution, Unit=unit )
!
! INPUTS:
!   RTSolution:  CRTM RTSolution object to display.
!                UNITS:      N/A
!                TYPE:       CRTM_RTSolution_type
!                DIMENSION:  Scalar or Rank-2 (n_channels x n_profiles)
!                ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!   Unit:        Unit number for an already open file to which the output
!                will be written.
!                If the argument is specified and the file unit is not
!                connected, the output goes to stdout.
!                UNITS:      N/A
!                TYPE:       INTEGER
!                DIMENSION:  Scalar
!                ATTRIBUTES: INTENT(IN), OPTIONAL
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE scalar_inspect( RTSolution, Unit )
    ! Arguments
    TYPE(CRTM_RTSolution_type), INTENT(IN) :: RTSolution
    INTEGER,          OPTIONAL, INTENT(IN) :: Unit
    ! Local variables
    INTEGER :: fid

    ! Setup
    fid = output_unit
    IF ( PRESENT(unit) ) THEN
      IF ( file_open(unit) ) fid = unit
    END IF


    WRITE(fid,'(1x,"RTSolution OBJECT")')
    ! Display components
    WRITE(fid,'(3x,"Sensor Id                : ",a )') trim(rtsolution%Sensor_ID)
    WRITE(fid,'(3x,"WMO Satellite Id         : ",i0)') rtsolution%WMO_Satellite_ID
    WRITE(fid,'(3x,"WMO Sensor Id            : ",i0)') rtsolution%WMO_Sensor_ID
    WRITE(fid,'(3x,"Channel                  : ",i0)') rtsolution%Sensor_Channel
    WRITE(fid,'(3x,"RT Algorithm Name        : ",a )') rtsolution%RT_Algorithm_Name
    WRITE(fid,'(3x,"Scattering Optical Depth : ",es13.6)') rtsolution%SOD
    WRITE(fid,'(3x,"Surface Emissivity       : ",es13.6)') rtsolution%Surface_Emissivity
    WRITE(fid,'(3x,"Surface Reflectivity     : ",es13.6)') rtsolution%Surface_Reflectivity
    WRITE(fid,'(3x,"Up Radiance              : ",es13.6)') rtsolution%Up_Radiance
    WRITE(fid,'(3x,"Down Radiance            : ",es13.6)') rtsolution%Down_Radiance
    WRITE(fid,'(3x,"Down Solar Radiance      : ",es13.6)') rtsolution%Down_Solar_Radiance
    WRITE(fid,'(3x,"Surface Planck Radiance  : ",es13.6)') rtsolution%Surface_Planck_Radiance
    WRITE(fid,'(3x,"Radiance                 : ",es13.6)') rtsolution%Radiance
    WRITE(fid,'(3x,"Brightness Temperature   : ",es13.6)') rtsolution%Brightness_Temperature
    IF ( .NOT. crtm_rtsolution_associated(rtsolution) ) RETURN
    WRITE(fid,'(3x,"n_Layers : ",i0)') rtsolution%n_Layers
    WRITE(fid,'(3x,"Upwelling Radiance       :")')
    WRITE(fid,'(5(1x,es13.6,:))') rtsolution%Upwelling_Radiance
    WRITE(fid,'(3x,"Layer Optical Depth      :")')
    WRITE(fid,'(5(1x,es13.6,:))') rtsolution%Layer_Optical_Depth
  END SUBROUTINE scalar_inspect


  SUBROUTINE rank2_inspect( RTSolution, Unit )
    TYPE(CRTM_RTSolution_type), INTENT(IN) :: RTSolution(:,:)
    INTEGER,          OPTIONAL, INTENT(IN) :: Unit
    INTEGER :: fid
    INTEGER :: i, n_channels
    INTEGER :: j, n_profiles

    fid = output_unit
    IF ( PRESENT(unit) ) THEN
      IF ( file_open(unit) ) fid = unit
    END IF

    n_channels = SIZE(rtsolution,1)
    n_profiles = SIZE(rtsolution,2)
    DO j = 1, n_profiles
      DO i = 1, n_channels
        WRITE(fid, fmt='(1x,"PROFILE INDEX:",i0,", CHANNEL INDEX:",i0," - ")', advance='NO') j,i
        CALL scalar_inspect(rtsolution(i,j), unit=unit)
      END DO
    END DO
  END SUBROUTINE rank2_inspect



!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_DefineVersion
!
! PURPOSE:
!   Subroutine to return the module version information.
!
! CALLING SEQUENCE:
!   CALL CRTM_RTSolution_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 crtm_rtsolution_defineversion( Id )
    CHARACTER(*), INTENT(OUT) :: id
    id = module_version_id
  END SUBROUTINE crtm_rtsolution_defineversion


!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!   CRTM_RTSolution_Compare
!
! PURPOSE:
!   Elemental function to compare two CRTM_RTSolution objects to within
!   a user specified number of significant figures.
!
! CALLING SEQUENCE:
!   is_comparable = CRTM_RTSolution_Compare( x, y, n_SigFig=n_SigFig )
!
! OBJECTS:
!   x, y:          Two CRTM RTSolution objects to be compared.
!                  UNITS:      N/A
!                  TYPE:       CRTM_RTSolution_type
!                  DIMENSION:  Scalar or any rank
!                  ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!   n_SigFig:      Number of significant figure to compare floating point
!                  components.
!                  UNITS:      N/A
!                  TYPE:       INTEGER
!                  DIMENSION:  Conformable with inputs
!                  ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!   is_comparable: Logical value indicating whether the inputs are
!                  comparable.
!                  UNITS:      N/A
!                  TYPE:       LOGICAL
!                  DIMENSION:  Same as inputs.
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_compare( &
    x, &
    y, &
    n_SigFig ) &
  result( is_comparable )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: x, y
    INTEGER,          OPTIONAL, INTENT(IN) :: n_sigfig
    LOGICAL :: is_comparable
    ! Variables
    INTEGER :: n

    ! Set up
    is_comparable = .false.
    IF ( PRESENT(n_sigfig) ) THEN
      n = abs(n_sigfig)
    ELSE
      n = default_n_sigfig
    END IF

    ! Check the structure association status
    IF ( crtm_rtsolution_associated(x) .NEQV. crtm_rtsolution_associated(y) ) RETURN

    ! Check the sensor information
    IF ( (x%Sensor_ID        /= y%Sensor_ID       ) .OR. &
         (x%WMO_Satellite_ID /= y%WMO_Satellite_ID) .OR. &
         (x%WMO_Sensor_ID    /= y%WMO_Sensor_ID   ) .OR. &
         (x%Sensor_Channel   /= y%Sensor_Channel  ) ) RETURN

    ! Check the RT algorithm name
    IF ( x%RT_Algorithm_Name /= y%RT_Algorithm_Name ) RETURN

    ! Check the scalar components
    IF ( .NOT. compares_within_tolerance(x%SOD                    , y%SOD                    , n) .OR. &
         .NOT. compares_within_tolerance(x%Surface_Emissivity     , y%Surface_Emissivity     , n) .OR. &
         .NOT. compares_within_tolerance(x%Surface_Reflectivity   , y%Surface_Reflectivity   , n) .OR. &
         .NOT. compares_within_tolerance(x%Up_Radiance            , y%Up_Radiance            , n) .OR. &
         .NOT. compares_within_tolerance(x%Down_Radiance          , y%Down_Radiance          , n) .OR. &
         .NOT. compares_within_tolerance(x%Down_Solar_Radiance    , y%Down_Solar_Radiance    , n) .OR. &
         .NOT. compares_within_tolerance(x%Surface_Planck_Radiance, y%Surface_Planck_Radiance, n) .OR. &
         .NOT. compares_within_tolerance(x%Radiance               , y%Radiance               , n) .OR. &
         .NOT. compares_within_tolerance(x%Brightness_Temperature , y%Brightness_Temperature , n) ) RETURN

    ! Check the array components
    IF ( crtm_rtsolution_associated(x) .AND. crtm_rtsolution_associated(y) ) THEN
      IF ( (.NOT. all(compares_within_tolerance(x%Upwelling_Overcast_Radiance, y%Upwelling_Overcast_Radiance, n))) .OR. &
           (.NOT. all(compares_within_tolerance(x%Upwelling_Radiance         , y%Upwelling_Radiance         , n))) .OR. &
           (.NOT. all(compares_within_tolerance(x%Layer_Optical_Depth        , y%Layer_Optical_Depth        , n))) ) RETURN
    END IF

    ! If we get here, the structures are comparable
    is_comparable = .true.

  END FUNCTION crtm_rtsolution_compare


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_Statistics
!
! PURPOSE:
!   Function to compute the statistics of an array of CRTM RTSolution objects.
!
! CALLING SEQUENCE:
!   Error_Status = CRTM_RTSolution_Statistics( rts, rts_stats )
!
! INPUTS:
!   rts:          RTSolution object array for which statistics are required.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Rank-2 (n_channels x n_profiles)
!                 ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!   rts_stats:    Allocatable RTSolution object array containing the statistics
!                 for each channel.
!                   rts_stats(:,1) contains the profile average
!                   rts_stats(:,2) contains the profile standard deviation
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Rank-2 (n_channels x 2)
!                 ATTRIBUTES: INTENT(IN), ALLOCATABLE
!
! 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 crtm_rtsolution_statistics(rts, rts_stats) RESULT( err_stat )
    ! Arguments
    TYPE(crtm_rtsolution_type),              INTENT(IN)  :: rts(:,:)
    TYPE(crtm_rtsolution_type), ALLOCATABLE, INTENT(OUT) :: rts_stats(:,:)
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_RTSolution_Statistics'
    ! Function variables
    CHARACTER(ML) :: err_msg
    CHARACTER(ML) :: alloc_msg
    INTEGER :: alloc_stat
    INTEGER :: n_channels, l
    INTEGER :: n_profiles, m
    REAL(fp) :: factor

    ! Setup
    err_stat = success
    n_channels = SIZE(rts, dim=1)
    n_profiles = SIZE(rts, dim=2)
    factor = REAL(n_profiles,fp)
    
    
    ! Allocate the output stats object array
    ALLOCATE( rts_stats(n_channels, 2), &
              stat = alloc_stat, errmsg = alloc_msg )
    IF ( alloc_stat /= 0 ) THEN
      err_msg = 'Error allocating output RTSolution structure - '//trim(alloc_msg)
      err_stat = failure
      CALL display_message( routine_name, err_msg, err_stat ); RETURN
    END IF
    rts_stats(:,1) = rts(:,1)
    rts_stats(:,2) = rts(:,1)
    CALL crtm_rtsolution_zero(rts_stats)


    ! Compute the average
    DO m = 1, n_profiles
      DO l = 1, n_channels
        rts_stats(l,1) = rts_stats(l,1) + rts(l,m)
      END DO
    END DO
    rts_stats(:,1) = rts_stats(:,1)/factor


    ! Compute the standard deviation
    DO m = 1, n_profiles
      DO l = 1, n_channels
        rts_stats(l,2) = rts_stats(l,2) + (rts(l,m) - rts_stats(l,1))**2
      END DO
    END DO
    rts_stats(:,2) = sqrt(rts_stats(:,2)/factor)


    ! Replace the algorithm identifier
    rts_stats(:,1)%RT_Algorithm_Name = 'Object average'
    rts_stats(:,2)%RT_Algorithm_Name = 'Object standard deviation'
    
  END FUNCTION crtm_rtsolution_statistics


!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_InquireFile
!
! PURPOSE:
!   Function to inquire CRTM RTSolution object files.
!
! CALLING SEQUENCE:
!   Error_Status = CRTM_RTSolution_InquireFile( Filename               , &
!                                               n_Channels = n_Channels, &
!                                               n_Profiles = n_Profiles  )
!
! INPUTS:
!   Filename:       Character string specifying the name of a
!                   CRTM RTSolution data file to read.
!                   UNITS:      N/A
!                   TYPE:       CHARACTER(*)
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN)
!
! OPTIONAL OUTPUTS:
!   n_Channels:     The number of spectral channels for which there is
!                   data in the file.
!                   UNITS:      N/A
!                   TYPE:       INTEGER
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
!   n_Profiles:     The number of profiles in the data file.
!                   UNITS:      N/A
!                   TYPE:       INTEGER
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
! 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 crtm_rtsolution_inquirefile( &
    Filename   , &  ! Input
    n_Channels , &  ! Optional output
    n_Profiles ) &  ! Optional output
  result( err_stat )
    ! Arguments
    CHARACTER(*),           INTENT(IN)  :: filename
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_channels
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_profiles
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_RTSolution_InquireFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat
    INTEGER :: fid
    INTEGER :: l, m

    ! Set up
    err_stat = success
    ! Check that the file exists
    IF ( .NOT. file_exists( trim(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 the number of channels,profiles
    READ( fid,iostat=io_stat,iomsg=io_msg ) l, m
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimensions from '//trim(filename)//' - '//trim(io_msg)
      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

    ! Set the return arguments
    IF ( PRESENT(n_channels) ) n_channels = l
    IF ( PRESENT(n_profiles) ) n_profiles = m

  CONTAINS

    SUBROUTINE inquire_cleanup()
      IF ( file_open( filename ) ) THEN
        CLOSE( fid,iostat=io_stat,iomsg=io_msg )
        IF ( io_stat /= success ) &
          msg = trim(msg)//'; Error closing input file during error cleanup - '//trim(io_msg)
      END IF
      err_stat = failure
      CALL display_message( routine_name, msg, err_stat )
    END SUBROUTINE inquire_cleanup

  END FUNCTION crtm_rtsolution_inquirefile


!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_ReadFile
!
! PURPOSE:
!   Function to read CRTM RTSolution object files.
!
! CALLING SEQUENCE:
!   Error_Status = CRTM_RTSolution_ReadFile( Filename                , &
!                                            RTSolution              , &
!                                            Quiet      = Quiet      , &
!                                            n_Channels = n_Channels , &
!                                            n_Profiles = n_Profiles , &
!
! INPUTS:
!   Filename:     Character string specifying the name of an
!                 RTSolution format data file to read.
!                 UNITS:      N/A
!                 TYPE:       CHARACTER(*)
!                 DIMENSION:  Scalar
!                 ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!   RTSolution:   CRTM RTSolution object array containing the RTSolution
!                 data.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Rank-2 (n_Channels x n_Profiles)
!                 ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!   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
!
! OPTIONAL OUTPUTS:
!   n_Channels:   The number of channels for which data was read.
!                 UNITS:      N/A
!                 TYPE:       INTEGER
!                 DIMENSION:  Scalar
!                 ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
!   n_Profiles:   The number of profiles for which data was read.
!                 UNITS:      N/A
!                 TYPE:       INTEGER
!                 DIMENSION:  Scalar
!                 ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
!
! 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 crtm_rtsolution_readfile( &
    Filename   , &  ! Input
    RTSolution , &  ! Output
    Quiet      , &  ! Optional input
    n_Channels , &  ! Optional output
    n_Profiles , &  ! Optional output
    Old_Version, &  ! Optional input (Allow reading of previous version files)
    debug      ) &  ! Optional input (Debug output control)
  result( err_stat )
    ! Arguments
    CHARACTER(*),               INTENT(IN)  :: filename
    TYPE(crtm_rtsolution_type), INTENT(OUT) :: rtsolution(:,:)
    LOGICAL,          OPTIONAL, INTENT(IN)  :: quiet
    INTEGER,          OPTIONAL, INTENT(OUT) :: n_channels
    INTEGER,          OPTIONAL, INTENT(OUT) :: n_profiles
    LOGICAL,          OPTIONAL, INTENT(IN)  :: old_version
    LOGICAL,          OPTIONAL, INTENT(IN)  :: debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_RTSolution_ReadFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat
    LOGICAL :: noisy
    INTEGER :: fid
    INTEGER :: l, n_file_channels, n_input_channels
    INTEGER :: m, n_file_profiles, n_input_profiles


    ! Set up
    err_stat = success
    ! ...Check Quiet argument
    noisy = .true.
    IF ( PRESENT(quiet) ) noisy = .NOT. quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(debug) ) noisy = debug


    ! Open the file
    err_stat = open_binary_file( filename, fid )
    IF ( err_stat /= success ) THEN
      msg = 'Error opening '//trim(filename)
      CALL read_cleanup(); RETURN
    END IF


    ! Read the dimensions
    READ( fid,iostat=io_stat,iomsg=io_msg ) n_file_channels, n_file_profiles
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimensions from '//trim(filename)//' - '//trim(io_msg)
      CALL read_cleanup(); RETURN
    END IF
    ! ...Check if n_Channels in file is > size of output array
    n_input_channels = SIZE(rtsolution,dim=1)
    IF ( n_file_channels > n_input_channels ) THEN
      WRITE( msg,'("Number of channels, ",i0," > size of the output RTSolution", &
                  &" array dimension, ",i0,". Only the first ",i0, &
                  &" channels will be read.")' ) &
                  n_file_channels, n_input_channels, n_input_channels
      CALL display_message( routine_name, msg, warning )
    END IF
    n_input_channels = min(n_input_channels, n_file_channels)
    ! ...Check if n_Profiles in file is > size of output array
    n_input_profiles = SIZE(rtsolution,dim=2)
    IF ( n_file_profiles > n_input_profiles ) THEN
      WRITE( msg,'( "Number of profiles, ",i0," > size of the output RTSolution", &
                   &" array dimension, ",i0,". Only the first ",i0, &
                   &" profiles will be read.")' ) &
                   n_file_profiles, n_input_profiles, n_input_profiles
      CALL display_message( routine_name, msg, warning )
    END IF
    n_input_profiles = min(n_input_profiles, n_file_profiles)


    ! Loop over all the profiles and channels
    profile_loop: DO m = 1, n_input_profiles
      channel_loop: DO l = 1, n_input_channels
        err_stat = read_record( fid, rtsolution(l,m), old_version=old_version )
        IF ( err_stat /= success ) THEN
          WRITE( msg,'("Error reading RTSolution element (",i0,",",i0,") from ",a)' ) &
                 l, m, trim(filename)
          CALL read_cleanup(); RETURN
        END IF
      END DO channel_loop
    END DO profile_loop


    ! 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 read_cleanup(); RETURN
    END IF


    ! Set the return values
    IF ( PRESENT(n_channels) ) n_channels = n_input_channels
    IF ( PRESENT(n_profiles) ) n_profiles = n_input_profiles


    ! Output an info message
    IF ( noisy ) THEN
      WRITE( msg,'("Number of channels and profiles read from ",a,": ",i0,1x,i0)' ) &
             trim(filename), n_input_channels, n_input_profiles
      CALL display_message( routine_name, 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 crtm_rtsolution_destroy( rtsolution )
      err_stat = failure
      CALL display_message( routine_name, msg, err_stat )
    END SUBROUTINE read_cleanup

  END FUNCTION crtm_rtsolution_readfile


!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_RTSolution_WriteFile
!
! PURPOSE:
!   Function to write CRTM RTSolution object files.
!
! CALLING SEQUENCE:
!   Error_Status = CRTM_RTSolution_WriteFile( Filename     , &
!                                             RTSolution   , &
!                                             Quiet = Quiet  )
!
! INPUTS:
!   Filename:     Character string specifying the name of the
!                 RTSolution format data file to write.
!                 UNITS:      N/A
!                 TYPE:       CHARACTER(*)
!                 DIMENSION:  Scalar
!                 ATTRIBUTES: INTENT(IN)
!
!   RTSolution:   CRTM RTSolution object array containing the RTSolution
!                 data.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Rank-2 (n_Channels x n_Profiles)
!                 ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!   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
!
! SIDE EFFECTS:
!   - If the output file already exists, it is overwritten.
!   - If an error occurs during *writing*, the output file is deleted before
!     returning to the calling routine.
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION crtm_rtsolution_writefile( &
    Filename   , &  ! Input
    RTSolution , &  ! Input
    Quiet      , &  ! Optional input
    Debug      ) &  ! Optional input (Debug output control)
  result( err_stat )
    ! Arguments
    CHARACTER(*),               INTENT(IN) :: filename
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rtsolution(:,:)
    LOGICAL,          OPTIONAL, INTENT(IN) :: quiet
    LOGICAL,          OPTIONAL, INTENT(IN) :: debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_RTSolution_WriteFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat
    LOGICAL :: noisy
    INTEGER :: fid
    INTEGER :: l, n_output_channels
    INTEGER :: m, n_output_profiles

    ! Set up
    err_stat = success
    ! ...Check Quiet argument
    noisy = .true.
    IF ( PRESENT(quiet) ) noisy = .NOT. quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(debug) ) noisy = debug
    n_output_channels = SIZE(rtsolution,dim=1)
    n_output_profiles = SIZE(rtsolution,dim=2)


    ! Open the file
    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


    ! Write the dimensions
    WRITE( fid,iostat=io_stat,iomsg=io_msg ) n_output_channels, n_output_profiles
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing dimensions to '//trim(filename)//' - '//trim(io_msg)
      CALL write_cleanup(); RETURN
    END IF


    ! Write the data
    profile_loop: DO m = 1, n_output_profiles
      channel_loop: DO l = 1, n_output_channels
        err_stat = write_record( fid, rtsolution(l,m) )
        IF ( err_stat /= success ) THEN
          WRITE( msg,'("Error writing RTSolution element (",i0,",",i0,") to ",a)' ) &
                 l, m, trim(filename)
          CALL write_cleanup(); RETURN
        END IF
      END DO channel_loop
    END DO profile_loop


    ! Close the file (if error, no delete)
    CLOSE( fid,status='KEEP',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


    ! Output an info message
    IF ( noisy ) THEN
      WRITE( msg,'("Number of channels and profiles written to ",a,": ",i0,1x,i0 )' ) &
             trim(filename), n_output_channels, n_output_profiles
      CALL display_message( routine_name, msg, information )
    END IF

  CONTAINS

    SUBROUTINE write_cleanup()
      IF ( file_open( filename ) ) THEN
        CLOSE( fid,status=write_error_status,iostat=io_stat,iomsg=io_msg )
        IF ( io_stat /= 0 ) &
          msg = trim(msg)//'; Error deleting 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 crtm_rtsolution_writefile



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

!------------------------------------------------------------------------------
!
! NAME:
!   CRTM_RTSolution_Equal
!
! PURPOSE:
!   Elemental function to test the equality of two CRTM_RTSolution objects.
!   Used in OPERATOR(==) interface block.
!
! CALLING SEQUENCE:
!   is_equal = CRTM_RTSolution_Equal( x, y )
!
!     or
!
!   IF ( x == y ) THEN
!     ...
!   END IF
!
! OBJECTS:
!   x, y:          Two CRTM RTSolution objects to be compared.
!                  UNITS:      N/A
!                  TYPE:       CRTM_RTSolution_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 crtm_rtsolution_equal( x, y ) RESULT( is_equal )
    TYPE(crtm_rtsolution_type) , INTENT(IN)  :: x, y
    LOGICAL :: is_equal

    ! Setup
    is_equal = .false.

    ! Check the structure association status
    IF ( crtm_rtsolution_associated(x) .NEQV. crtm_rtsolution_associated(y) ) RETURN

    ! Check scalars
    IF ( (x%n_Layers == y%n_Layers) .AND. &
         (x%Sensor_ID         == y%Sensor_ID        ) .AND. &
         (x%WMO_Satellite_ID  == y%WMO_Satellite_ID ) .AND. &
         (x%WMO_Sensor_ID     == y%WMO_Sensor_ID    ) .AND. &
         (x%Sensor_Channel    == y%Sensor_Channel   ) .AND. &
         (x%RT_Algorithm_Name == y%RT_Algorithm_Name) .AND. &
         (x%SOD                     .equalto. y%SOD                    ) .AND. &
         (x%Surface_Emissivity      .equalto. y%Surface_Emissivity     ) .AND. &
         (x%Surface_Reflectivity    .equalto. y%Surface_Reflectivity   ) .AND. &
         (x%Up_Radiance             .equalto. y%Up_Radiance            ) .AND. &
         (x%Down_Radiance           .equalto. y%Down_Radiance          ) .AND. &
         (x%Down_Solar_Radiance     .equalto. y%Down_Solar_Radiance    ) .AND. &
         (x%Surface_Planck_Radiance .equalto. y%Surface_Planck_Radiance) .AND. &
         (x%Radiance                .equalto. y%Radiance               ) .AND. &
         (x%Brightness_Temperature  .equalto. y%Brightness_Temperature ) ) &
      is_equal = .true.

    ! Check arrays (which may or may not be allocated)
    IF ( crtm_rtsolution_associated(x) .AND. crtm_rtsolution_associated(y) ) THEN
      is_equal = is_equal .AND. &
                 all(x%Upwelling_Overcast_Radiance .equalto. y%Upwelling_Overcast_Radiance ) .AND. &
                 all(x%Upwelling_Radiance          .equalto. y%Upwelling_Radiance          ) .AND. &
                 all(x%Layer_Optical_Depth         .equalto. y%Layer_Optical_Depth         )
    END IF

  END FUNCTION crtm_rtsolution_equal


!--------------------------------------------------------------------------------
!
! NAME:
!   CRTM_RTSolution_Add
!
! PURPOSE:
!   Pure function to add two CRTM RTSolution objects.
!   Used in OPERATOR(+) interface block.
!
! CALLING SEQUENCE:
!   rtssum = CRTM_RTSolution_Add( rts1, rts2 )
!
!     or
!
!   rtssum = rts1 + rts2
!
!
! INPUTS:
!   rts1, rts2: The RTSolution objects to add.
!               UNITS:      N/A
!               TYPE:       CRTM_RTSolution_type
!               DIMENSION:  Scalar
!               ATTRIBUTES: INTENT(IN OUT)
!
! RESULT:
!   rtssum:     RTSolution object containing the summed components.
!               UNITS:      N/A
!               TYPE:       CRTM_RTSolution_type
!               DIMENSION:  Scalar
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_add( rts1, rts2 ) RESULT( rtssum )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rts1, rts2
    TYPE(crtm_rtsolution_type) :: rtssum
    INTEGER :: k

    ! Check input
    ! ...If input structure association status different, do nothing
    IF ( crtm_rtsolution_associated(rts1) .NEQV. crtm_rtsolution_associated(rts2) ) RETURN
    ! ...If input structure for different sensors, do nothing
    IF ( (rts1%Sensor_ID         /= rts2%Sensor_ID        ) .AND. &
         (rts1%WMO_Satellite_ID  /= rts2%WMO_Satellite_ID ) .AND. &
         (rts1%WMO_Sensor_ID     /= rts2%WMO_Sensor_ID    ) .AND. &
         (rts1%Sensor_Channel    /= rts2%Sensor_Channel   ) ) RETURN

    ! Copy the first structure
    rtssum = rts1

    ! And add the second one's components to it
    ! ...Handle RT_Algorithm_Name
    rtssum%RT_Algorithm_Name = 'Object add'
    ! ...The scalar values
    rtssum%SOD                     = rtssum%SOD                     + rts2%SOD
    rtssum%Surface_Emissivity      = rtssum%Surface_Emissivity      + rts2%Surface_Emissivity
    rtssum%Surface_Reflectivity    = rtssum%Surface_Reflectivity    + rts2%Surface_Reflectivity
    rtssum%Up_Radiance             = rtssum%Up_Radiance             + rts2%Up_Radiance
    rtssum%Down_Radiance           = rtssum%Down_Radiance           + rts2%Down_Radiance
    rtssum%Down_Solar_Radiance     = rtssum%Down_Solar_Radiance     + rts2%Down_Solar_Radiance
    rtssum%Surface_Planck_Radiance = rtssum%Surface_Planck_Radiance + rts2%Surface_Planck_Radiance
    rtssum%Radiance                = rtssum%Radiance                + rts2%Radiance
    rtssum%Brightness_Temperature  = rtssum%Brightness_Temperature  + rts2%Brightness_Temperature
    ! ...The arrays (which may or may not be allocated)
    IF ( crtm_rtsolution_associated(rts1) .AND. crtm_rtsolution_associated(rts2) ) THEN
      k = rts1%n_Layers
      rtssum%Upwelling_Overcast_Radiance(1:k)  = rtssum%Upwelling_Overcast_Radiance(1:k)  + &
                                                   rts2%Upwelling_Overcast_Radiance(1:k)

      rtssum%Upwelling_Radiance(1:k) = rtssum%Upwelling_Radiance(1:k) + &
                                         rts2%Upwelling_Radiance(1:k)

      rtssum%Layer_Optical_Depth(1:k) = rtssum%Layer_Optical_Depth(1:k) + &
                                          rts2%Layer_Optical_Depth(1:k)
    END IF

  END FUNCTION crtm_rtsolution_add


!--------------------------------------------------------------------------------
!
! NAME:
!   CRTM_RTSolution_Subtract
!
! PURPOSE:
!   Pure function to subtract two CRTM RTSolution objects.
!   Used in OPERATOR(-) interface block.
!
! CALLING SEQUENCE:
!   rtsdiff = CRTM_RTSolution_Subtract( rts1, rts2 )
!
!     or
!
!   rtsdiff = rts1 - rts2
!
!
! INPUTS:
!   rts1, rts2: The RTSolution objects to difference.
!               UNITS:      N/A
!               TYPE:       CRTM_RTSolution_type
!               DIMENSION:  Scalar
!               ATTRIBUTES: INTENT(IN OUT)
!
! RESULT:
!   rtsdiff:    RTSolution object containing the differenced components.
!               UNITS:      N/A
!               TYPE:       CRTM_RTSolution_type
!               DIMENSION:  Scalar
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_subtract( rts1, rts2 ) RESULT( rtsdiff )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rts1, rts2
    TYPE(crtm_rtsolution_type) :: rtsdiff
    INTEGER :: k

    ! Check input
    ! ...If input structure association status different, do nothing
    IF ( crtm_rtsolution_associated(rts1) .NEQV. crtm_rtsolution_associated(rts2) ) RETURN
    ! ...If input structure for different sensors, do nothing
    IF ( (rts1%Sensor_ID         /= rts2%Sensor_ID        ) .AND. &
         (rts1%WMO_Satellite_ID  /= rts2%WMO_Satellite_ID ) .AND. &
         (rts1%WMO_Sensor_ID     /= rts2%WMO_Sensor_ID    ) .AND. &
         (rts1%Sensor_Channel    /= rts2%Sensor_Channel   ) ) RETURN

    ! Copy the first structure
    rtsdiff = rts1

    ! And subtract the second one's components from it
    ! ...Handle RT_Algorithm_Name
    rtsdiff%RT_Algorithm_Name = 'Object subtract'
    ! ...The scalar values
    rtsdiff%SOD                     = rtsdiff%SOD                     - rts2%SOD
    rtsdiff%Surface_Emissivity      = rtsdiff%Surface_Emissivity      - rts2%Surface_Emissivity
    rtsdiff%Surface_Reflectivity    = rtsdiff%Surface_Reflectivity    - rts2%Surface_Reflectivity
    rtsdiff%Up_Radiance             = rtsdiff%Up_Radiance             - rts2%Up_Radiance
    rtsdiff%Down_Radiance           = rtsdiff%Down_Radiance           - rts2%Down_Radiance
    rtsdiff%Down_Solar_Radiance     = rtsdiff%Down_Solar_Radiance     - rts2%Down_Solar_Radiance
    rtsdiff%Surface_Planck_Radiance = rtsdiff%Surface_Planck_Radiance - rts2%Surface_Planck_Radiance
    rtsdiff%Radiance                = rtsdiff%Radiance                - rts2%Radiance
    rtsdiff%Brightness_Temperature  = rtsdiff%Brightness_Temperature  - rts2%Brightness_Temperature
    ! ...The arrays (which may or may not be allocated)
    IF ( crtm_rtsolution_associated(rts1) .AND. crtm_rtsolution_associated(rts2) ) THEN
      k = rts1%n_Layers
      rtsdiff%Upwelling_Overcast_Radiance(1:k)  = rtsdiff%Upwelling_Overcast_Radiance(1:k)  - &
                                                     rts2%Upwelling_Overcast_Radiance(1:k)

      rtsdiff%Upwelling_Radiance(1:k) = rtsdiff%Upwelling_Radiance(1:k) - &
                                           rts2%Upwelling_Radiance(1:k)

      rtsdiff%Layer_Optical_Depth(1:k) = rtsdiff%Layer_Optical_Depth(1:k) - &
                                            rts2%Layer_Optical_Depth(1:k)
    END IF

  END FUNCTION crtm_rtsolution_subtract


!--------------------------------------------------------------------------------
!
! NAME:
!   CRTM_RTSolution_Exponent
!
! PURPOSE:
!   Elemental utility function to raise the components of a CRTM RTSolution
!   object to the supplied integer power
!
!   Used to compute RTSolution statistics
!
! CALLING SEQUENCE:
!   rts_power = CRTM_RTSolution_Exponent( rts, power )
!
! OBJECTS:
!   rts:          RTSolution structure which is to have its members values
!                 raised to an integer power.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar or any rank
!                 ATTRIBUTES: INTENT(IN)
!
! INPUTS:
!   power:        Exponent power to be used.
!                 UNITS:      N/A
!                 TYPE:       INRTGER
!                 DIMENSION:  Same as input rts object
!                 ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!   rts_power:    RTSolution structure containing the exponent result.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Same as input rts object
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_exponent( rts, power ) RESULT( rts_power )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rts
    INTEGER                   , INTENT(IN) :: power
    TYPE(crtm_rtsolution_type) :: rts_power
    INTEGER :: k

    ! Copy the structure
    rts_power = rts

    ! Raise the components to the supplied power
    ! ...Handle RT_Algorithm_Name
    rts_power%RT_Algorithm_Name = 'Object exponent'
    ! ...The scalar values
    rts_power%SOD                     = (rts_power%SOD)**power
    rts_power%Surface_Emissivity      = (rts_power%Surface_Emissivity)**power
    rts_power%Surface_Reflectivity    = (rts_power%Surface_Reflectivity)**power
    rts_power%Up_Radiance             = (rts_power%Up_Radiance)**power
    rts_power%Down_Radiance           = (rts_power%Down_Radiance)**power
    rts_power%Down_Solar_Radiance     = (rts_power%Down_Solar_Radiance)**power
    rts_power%Surface_Planck_Radiance = (rts_power%Surface_Planck_Radiance)**power
    rts_power%Radiance                = (rts_power%Radiance)**power
    rts_power%Brightness_Temperature  = (rts_power%Brightness_Temperature)**power
    ! ...The arrays (which may or may not be allocated)
    IF ( crtm_rtsolution_associated(rts) ) THEN
      k = rts%n_Layers
      rts_power%Upwelling_Overcast_Radiance(1:k) = (rts_power%Upwelling_Overcast_Radiance(1:k))**power
      rts_power%Upwelling_Radiance(1:k)          = (rts_power%Upwelling_Radiance(1:k))**power
      rts_power%Layer_Optical_Depth(1:k)         = (rts_power%Layer_Optical_Depth(1:k))**power
    END IF

  END FUNCTION crtm_rtsolution_exponent


!--------------------------------------------------------------------------------
!
! NAME:
!   CRTM_RTSolution_Normalise
!
! PURPOSE:
!   Elemental utility function to normalise the components of
!   a CRTM RTSolution object.
!
!   Used to compute RTSolution statistics
!
! CALLING SEQUENCE:
!   rts_normal = CRTM_RTSolution_Normalise( rts, factor )
!
! OBJECTS:
!   rts:          RTSolution structure which is to have its members values
!                 normalised.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar or any rank
!                 ATTRIBUTES: INTENT(IN)
!
! INPUTS:
!   factor:       Normalisation factor to be used.
!                 UNITS:      N/A
!                 TYPE:       REAL(fp)
!                 DIMENSION:  Same as input rts object
!                 ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!   rts_normal:   RTSolution structure containing the normalised input
!                 object.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Same as input rts object
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_normalise( rts, factor ) RESULT( rts_normal )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rts
    REAL(fp)                  , INTENT(IN) :: factor
    TYPE(crtm_rtsolution_type) :: rts_normal
    INTEGER :: k

    ! Copy the structure
    rts_normal = rts

    ! Raise the components to the supplied normal
    ! ...Handle RT_Algorithm_Name
    rts_normal%RT_Algorithm_Name = 'Object normalise'
    ! ...The scalar values
    rts_normal%SOD                     = rts_normal%SOD/factor
    rts_normal%Surface_Emissivity      = rts_normal%Surface_Emissivity/factor
    rts_normal%Surface_Reflectivity    = rts_normal%Surface_Reflectivity/factor
    rts_normal%Up_Radiance             = rts_normal%Up_Radiance/factor
    rts_normal%Down_Radiance           = rts_normal%Down_Radiance/factor
    rts_normal%Down_Solar_Radiance     = rts_normal%Down_Solar_Radiance/factor
    rts_normal%Surface_Planck_Radiance = rts_normal%Surface_Planck_Radiance/factor
    rts_normal%Radiance                = rts_normal%Radiance/factor
    rts_normal%Brightness_Temperature  = rts_normal%Brightness_Temperature/factor
    ! ...The arrays (which may or may not be allocated)
    IF ( crtm_rtsolution_associated(rts) ) THEN
      k = rts%n_Layers
      rts_normal%Upwelling_Overcast_Radiance(1:k) = rts_normal%Upwelling_Overcast_Radiance(1:k)/factor
      rts_normal%Upwelling_Radiance(1:k)          = rts_normal%Upwelling_Radiance(1:k)/factor
      rts_normal%Layer_Optical_Depth(1:k)         = rts_normal%Layer_Optical_Depth(1:k)/factor
    END IF

  END FUNCTION crtm_rtsolution_normalise


!--------------------------------------------------------------------------------
!
! NAME:
!   CRTM_RTSolution_Sqrt
!
! PURPOSE:
!   Elemental utility function to compute the SQRT() of the components of
!   a CRTM RTSolution object.
!
!   Used to compute RTSolution statistics
!
! CALLING SEQUENCE:
!   rts_sqrt = CRTM_RTSolution_Sqrt( rts )
!
! OBJECTS:
!   rts:          RTSolution structure which is to have the square root
!                 taken of its member's vales.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Scalar or any rank
!                 ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!   rts_sqrt:     RTSolution structure containing the square root of the
!                 input object.
!                 UNITS:      N/A
!                 TYPE:       CRTM_RTSolution_type
!                 DIMENSION:  Same as input rts object
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION crtm_rtsolution_sqrt( rts ) RESULT( rts_sqrt )
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rts
    TYPE(crtm_rtsolution_type) :: rts_sqrt
    INTEGER :: k

    ! Copy the structure
    rts_sqrt = rts

    ! Raise the components to the supplied normal
    ! ...Handle RT_Algorithm_Name
    rts_sqrt%RT_Algorithm_Name = 'Object SQRT()'
    ! ...The scalar values
    rts_sqrt%SOD                     = sqrt(rts_sqrt%SOD)
    rts_sqrt%Surface_Emissivity      = sqrt(rts_sqrt%Surface_Emissivity)
    rts_sqrt%Surface_Reflectivity    = sqrt(rts_sqrt%Surface_Reflectivity)
    rts_sqrt%Up_Radiance             = sqrt(rts_sqrt%Up_Radiance)
    rts_sqrt%Down_Radiance           = sqrt(rts_sqrt%Down_Radiance)
    rts_sqrt%Down_Solar_Radiance     = sqrt(rts_sqrt%Down_Solar_Radiance)
    rts_sqrt%Surface_Planck_Radiance = sqrt(rts_sqrt%Surface_Planck_Radiance)
    rts_sqrt%Radiance                = sqrt(rts_sqrt%Radiance)
    rts_sqrt%Brightness_Temperature  = sqrt(rts_sqrt%Brightness_Temperature)
    ! ...The arrays (which may or may not be allocated)
    IF ( crtm_rtsolution_associated(rts) ) THEN
      k = rts%n_Layers
      rts_sqrt%Upwelling_Overcast_Radiance(1:k) = sqrt(rts_sqrt%Upwelling_Overcast_Radiance(1:k))
      rts_sqrt%Upwelling_Radiance(1:k)          = sqrt(rts_sqrt%Upwelling_Radiance(1:k))
      rts_sqrt%Layer_Optical_Depth(1:k)         = sqrt(rts_sqrt%Layer_Optical_Depth(1:k))
    END IF

  END FUNCTION crtm_rtsolution_sqrt


!
! NAME:
!   Read_Record
!
! PURPOSE:
!   Utility function to read a single RTSolution data record
!

  FUNCTION read_record( &
    fid, &  ! Input
    rts, &  ! Output
    old_version) &  ! Optional input
  result( err_stat )
    ! Arguments
    INTEGER,                    INTENT(IN)  :: fid
    TYPE(crtm_rtsolution_type), INTENT(OUT) :: rts
    LOGICAL,          OPTIONAL, INTENT(IN)  :: old_version
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_RTSolution_ReadFile(Record)'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat
    INTEGER :: n_layers
    LOGICAL :: current_version

    ! Set up
    err_stat = success
    ! ...PRocess optional arguments
    current_version = .true.
    IF ( PRESENT(old_version) ) current_version = .NOT. old_version

    ! Read the dimensions
    READ( fid,iostat=io_stat,iomsg=io_msg ) n_layers
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimensions - '//trim(io_msg)
      CALL read_record_cleanup(); RETURN
    END IF


    ! Allocate the RTSolution structure if necessary
    IF ( n_layers > 0 ) THEN
      CALL crtm_rtsolution_create( rts, n_layers )
      IF ( .NOT. crtm_rtsolution_associated( rts ) ) THEN
        msg = 'Error creating output object.'
        CALL read_record_cleanup(); RETURN
      END IF
    END IF


    ! Read the sensor info
    READ( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%Sensor_Id       , &
      rts%WMO_Satellite_Id, &
      rts%WMO_Sensor_Id   , &
      rts%Sensor_Channel
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading sensor information - '//trim(io_msg)
      CALL read_record_cleanup(); RETURN
    END IF


    ! Read the RT algorithm name
    READ( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%RT_Algorithm_Name
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading RT Algorithm Name'//trim(io_msg)
      CALL read_record_cleanup(); RETURN
    END IF


    ! Read the forward radiative transfer intermediate results
    READ( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%SOD                    , &
      rts%Surface_Emissivity     , &
      rts%Surface_Reflectivity   , &
      rts%Up_Radiance            , &
      rts%Down_Radiance          , &
      rts%Down_Solar_Radiance    , &
      rts%Surface_Planck_Radiance
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading scalar intermediate results - '//trim(io_msg)
      CALL read_record_cleanup(); RETURN
    END IF
    IF ( n_layers > 0 ) THEN
      IF ( current_version ) THEN
        READ( fid,iostat=io_stat,iomsg=io_msg ) &
          rts%Upwelling_Overcast_Radiance , &
          rts%Upwelling_Radiance, &
          rts%Layer_Optical_Depth
      ELSE
        READ( fid,iostat=io_stat,iomsg=io_msg ) &
          rts%Upwelling_Radiance , &
          rts%Layer_Optical_Depth
      END IF
      IF ( io_stat /= 0 ) THEN
        msg = 'Error reading array intermediate results - '//trim(io_msg)
        CALL read_record_cleanup(); RETURN
      END IF
    END IF


    ! Read the radiative transfer results
    READ( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%Radiance              , &
      rts%Brightness_Temperature
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading result data - '//trim(io_msg)
      CALL read_record_cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE read_record_cleanup()
      CALL crtm_rtsolution_destroy( rts )
      CLOSE( fid,iostat=io_stat,iomsg=io_msg )
      IF ( io_stat /= success ) &
        msg = trim(msg)//'; Error closing file during error cleanup - '//trim(io_msg)
      err_stat = failure
      CALL display_message( routine_name, msg, err_stat )
    END SUBROUTINE read_record_cleanup

  END FUNCTION read_record


!
! NAME:
!   Write_Record
!
! PURPOSE:
!   Function to write a single RTSolution data record
!

  FUNCTION write_record( &
    fid, &  ! Input
    rts) &  ! Input
  result( err_stat )
    ! Arguments
    INTEGER,                    INTENT(IN) :: fid
    TYPE(crtm_rtsolution_type), INTENT(IN) :: rts
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_RTSolution_WriteFile(Record)'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Set up
    err_stat = success


    ! Write the data dimensions
    WRITE( fid,iostat=io_stat,iomsg=io_msg ) rts%n_Layers
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing dimensions - '//trim(io_msg)
      CALL write_record_cleanup(); RETURN
    END IF


    ! Write the sensor info
    WRITE( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%Sensor_Id       , &
      rts%WMO_Satellite_Id, &
      rts%WMO_Sensor_Id   , &
      rts%Sensor_Channel
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing sensor information - '//trim(io_msg)
      CALL write_record_cleanup(); RETURN
    END IF


    ! Write the sensor info
    WRITE( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%RT_Algorithm_Name
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing RT Algorithm Name'//trim(io_msg)
      CALL write_record_cleanup(); RETURN
    END IF


    ! Write the forward radiative transfer intermediate results
    WRITE( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%SOD                    , &
      rts%Surface_Emissivity     , &
      rts%Surface_Reflectivity   , &
      rts%Up_Radiance            , &
      rts%Down_Radiance          , &
      rts%Down_Solar_Radiance    , &
      rts%Surface_Planck_Radiance
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing scalar intermediate results - '//trim(io_msg)
      CALL write_record_cleanup(); RETURN
    END IF
    IF ( rts%n_Layers > 0 ) THEN
      WRITE( fid,iostat=io_stat,iomsg=io_msg ) &
        rts%Upwelling_Overcast_Radiance , &
        rts%Upwelling_Radiance, &
        rts%Layer_Optical_Depth
      IF ( io_stat /= 0 ) THEN
        msg = 'Error writing array intermediate results - '//trim(io_msg)
        CALL write_record_cleanup(); RETURN
      END IF
    END IF


    ! Write the radiative transfer results
    WRITE( fid,iostat=io_stat,iomsg=io_msg ) &
      rts%Radiance              , &
      rts%Brightness_Temperature
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing result data - '//trim(io_msg)
      CALL write_record_cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE write_record_cleanup()
      CLOSE( fid,status=write_error_status,iostat=io_stat,iomsg=io_msg )
      IF ( io_stat /= 0 ) &
        msg = trim(msg)//'; Error closing file during error cleanup - '//trim(io_msg)
      err_stat = failure
      CALL display_message( routine_name, msg, err_stat )
    END SUBROUTINE write_record_cleanup

  END FUNCTION write_record

END MODULE crtm_rtsolution_define