CRTM_GeometryInfo.f90

Go to the documentation of this file.

!
! CRTM_GeometryInfo
!
! Application module for the GeometryInfo structure.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 19-May-2004
!                       paul.vandelst@noaa.gov
!
!                       Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!

MODULE crtm_geometryinfo

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module use
  USE type_kinds     , ONLY: fp
  USE message_handler, ONLY: success, warning, failure, display_message
  USE date_utility   , ONLY: dayofyear
  USE crtm_parameters, ONLY: zero, one, two  , &
                             twopi           , &
                             earth_radius    , &
                             satellite_height, &
                             degrees_to_radians, &
                             max_trans_zenith_angle
  USE crtm_geometryinfo_define, ONLY: crtm_geometry_type    , &
                                      crtm_geometryinfo_type, &
                                      crtm_geometryinfo_isvalid
  ! Disable implicit typing
  IMPLICIT NONE


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: crtm_geometryinfo_compute
  PUBLIC :: saconv, saconv_tl, saconv_ad 
  PUBLIC :: vaconv, vaconv_tl, vaconv_ad 
  PUBLIC :: crtm_geometryinfo_version 
  
  
  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module
  CHARACTER(*),  PARAMETER :: module_version_id = &
  '$Id: CRTM_GeometryInfo.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
  ! Metres->kilometres conversion factor
  REAL(fp), PARAMETER :: m_to_km = 1.0e-03_fp


CONTAINS


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

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_GeometryInfo_Compute
! 
! PURPOSE:
!       Elemental subroutine to compute the derived geometry from the user
!       specified components of the CRTM GeometryInfo structure.
!
! CALLING SEQUENCE:
!       CALL CRTM_GeometryInfo_Compute( GeometryInfo )
!
! INPUTS:
!       GeometryInfo:  The GeometryInfo object containing the user
!                      defined inputs
!                      UNITS:      N/A
!                      TYPE:       CRTM_GeometryInfo_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       GeometryInfo:  The GeometryInfo structure with the non-user
!                      components filled.
!                      UNITS:      N/A
!                      TYPE:       CRTM_GeometryInfo_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN OUT)
!
! SIDE EFFECTS:
!       This function changes the values of the non-user components of the
!       GeometryInfo object.
!
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE crtm_geometryinfo_compute( gInfo )
    TYPE(crtm_geometryinfo_type), INTENT(IN OUT) :: ginfo
    
    REAL(fp) :: cosv

    ! Compute the derived values
    ! ...Derived sensor angles
    ginfo%Sensor_Scan_Radian    = degrees_to_radians * ginfo%user%Sensor_Scan_Angle
    ginfo%Sensor_Zenith_Radian  = degrees_to_radians * ginfo%user%Sensor_Zenith_Angle
    ginfo%Sensor_Azimuth_Radian = degrees_to_radians * ginfo%user%Sensor_Azimuth_Angle
    ginfo%Cosine_Sensor_Zenith  = cos(ginfo%Sensor_Zenith_Radian)
    ginfo%Secant_Sensor_Zenith  = one / ginfo%Cosine_Sensor_Zenith
    ! ...Check user zenith angle. If it is larger than the transmittance algorithms' limit
    !    then use the limiting value in the algorithms.  The optical depths will be scaled 
    !    back to the those at the user zenith angle.
    IF( ginfo%user%Sensor_Zenith_Angle > max_trans_zenith_angle )THEN
      ginfo%Trans_Zenith_Radian = degrees_to_radians * max_trans_zenith_angle
      ginfo%Secant_Trans_Zenith = one / cos(ginfo%Trans_Zenith_Radian)
    ELSE
      ginfo%Trans_Zenith_Radian = ginfo%Sensor_Zenith_Radian
      ginfo%Secant_Trans_Zenith = ginfo%Secant_Sensor_Zenith
    END IF
    ! ...Distance ratio, but only if zenith angle is large enough
    IF ( abs(ginfo%user%Sensor_Zenith_Angle) > one ) THEN
      ginfo%Distance_Ratio = abs(sin(ginfo%Sensor_Scan_Radian)/sin(ginfo%Sensor_Zenith_Radian))
    END IF
    ! ...Derived source angles
    ginfo%Source_Zenith_Radian  = degrees_to_radians * ginfo%user%Source_Zenith_Angle
    ginfo%Source_Azimuth_Radian = degrees_to_radians * ginfo%user%Source_Azimuth_Angle
    cosv = cos(ginfo%Source_Zenith_Radian)
    IF( cosv /= zero )THEN
      ginfo%Secant_Source_Zenith = one / cosv
    ELSE
      ! ...Set to a large number with the sign of the original value
      ginfo%Secant_Source_Zenith = huge(cosv) * cosv/abs(cosv)  ! SIGN(HUGE(cosv),cosv) ?
    ENDIF
      
    ! ...Derived flux angles
    ginfo%Flux_Zenith_Radian = degrees_to_radians * ginfo%user%Flux_Zenith_Angle
    ginfo%Secant_Flux_Zenith = one / cos(ginfo%Flux_Zenith_Radian)
    ! ...AU ratio term
    ginfo%AU_ratio2 = compute_au_ratio2( ginfo%user%Year, ginfo%user%Month, ginfo%user%Day )
                    
  END SUBROUTINE crtm_geometryinfo_compute


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!      SACONV  
!
! PURPOSE:
!      Elemental subroutine to compute the sensor zenith angle at a given altitude
!      for a given surface sensor zenith angle.
!
! CALLING SEQUENCE:
!        CALL SACONV( Sensor_Zenith_Radian, &  ! Input
!                     Altitude            , &  ! Input
!                     Local_Zenith_Radian   )  ! Output
!
! INPUTS:
!       Sensor_Zenith_Radian: Sensor zenith angle at the Earth's surface
!                             UNITS:      radians
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN) 
!
!       Altitude:             The altitude at which the local sensor zenith angle
!                             is required 
!                             UNITS:      metres
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!  
! OUTPUTS:
!       Local_Zenith_Radian:  The sensor zenith angle at the supplied altitude.
!                             UNITS:      radians
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input
!                             ATTRIBUTES: INTENT(OUT)
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!
!                       Paul van Delst, 18-Nov-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE saconv( &
    Sensor_Zenith_Radian, &  ! Input
    Altitude            , &  ! Input
    Local_Zenith_Radian   )  ! Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: sensor_zenith_radian
    REAL(fp), INTENT(IN)  :: altitude 
    REAL(fp), INTENT(OUT) :: local_zenith_radian
    ! Local variables
    REAL(fp) :: ra

    ! Compute the radius, in km, of the point at which to calc the angle
    ra = earth_radius + (m_to_km * altitude)

    ! Compute the angle
    local_zenith_radian = asin((earth_radius / ra) * sin(sensor_zenith_radian))
  
  END SUBROUTINE saconv


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!      SACONV_TL  
!
! PURPOSE:
!      Tangent-linear form of elemental subroutine to compute the sensor zenith
!      angle at a given altitude for a given surface sensor zenith angle.
!
! CALLING SEQUENCE:
!        CALL SACONV_TL( Sensor_Zenith_Radian , &  ! FWD Input
!                        Altitude             , &  ! FWD Input
!                        Altitude_TL          , &  ! TL  Input
!                        Local_Zenith_Radian_TL )  ! TL  Output
!
! INPUTS:
!       Sensor_Zenith_Radian:    Sensor zenith angle at the Earth's surface
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN) 
!
!       Altitude:                The altitude at which the local sensor zenith
!                                angle is required 
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!  
!       Altitude_TL:             Tangent-linear altitude. 
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!  
! OUTPUTS:
!       Local_Zenith_Radian_TL:  The tangent-linear sensor zenith angle at the
!                                supplied altitude.
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Same as input
!                                ATTRIBUTES: INTENT(OUT)
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!
!                       Paul van Delst, 20-Nov-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE saconv_tl( &
    Sensor_Zenith_Radian,  &  ! FWD Input
    Altitude,              &  ! FWD Input
    Altitude_TL,           &  ! TL  Input
    Local_Zenith_Radian_TL )  ! TL  Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: sensor_zenith_radian 
    REAL(fp), INTENT(IN)  :: altitude 
    REAL(fp), INTENT(IN)  :: altitude_tl 
    REAL(fp), INTENT(OUT) :: local_zenith_radian_tl
    ! Local variables
    REAL(fp) :: ra, ra_tl   
    REAL(fp) :: local_zenith_radian

    ! Compute the radius, in km, of the point at which to calc the angle
    ra    = earth_radius + (m_to_km * altitude)
    ra_tl = m_to_km * altitude_tl

    ! Compute the tangent-linear angle
    CALL saconv( sensor_zenith_radian, altitude, local_zenith_radian )
    local_zenith_radian_tl = -tan(local_zenith_radian) * ra_tl / ra
  
  END SUBROUTINE saconv_tl


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!      SACONV_AD  
!
! PURPOSE:
!      Adjoint form of elemental subroutine to compute the sensor zenith
!      angle at a given altitude for a given surface sensor zenith angle.
!
! CALLING SEQUENCE:
!        CALL SACONV_AD( Sensor_Zenith_Radian  , &  ! FWD Input
!                        Altitude              , &  ! FWD Input
!                        Local_Zenith_Radian_AD, &  ! AD  Input
!                        Altitude_AD             )  ! AD  Output
!
! INPUTS:
!       Sensor_Zenith_Radian:    Sensor zenith angle at the Earth's surface
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN) 
!
!       Altitude:                The altitude at which the local sensor zenith
!                                angle is required 
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!  
!       Local_Zenith_Radian_AD:  The adjoint sensor zenith angle at the
!                                supplied altitude.
!                                *** SET TO ZERO ON EXIT ***
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN OUT)
!  
! OUTPUTS:
!       Altitude_AD:             Adjoint altitude. 
!                                *** MUST HAVE VALUE ON ENTRY ***
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Same as input
!                                ATTRIBUTES: INTENT(IN OUT)
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!
!                       Paul van Delst, 23-Nov-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE saconv_ad( &
    Sensor_Zenith_Radian,   &  ! FWD Input
    Altitude,               &  ! FWD Input
    Local_Zenith_Radian_AD, &  ! AD  Input
    Altitude_AD             )  ! AD  Output
    ! Arguments
    REAL(fp), INTENT(IN)     :: sensor_zenith_radian 
    REAL(fp), INTENT(IN)     :: altitude 
    REAL(fp), INTENT(IN OUT) :: local_zenith_radian_ad  
    REAL(fp), INTENT(IN OUT) :: altitude_ad
    ! Local variables
    REAL(fp) :: local_zenith_radian  
    REAL(fp) :: ra, ra_ad   

    ! Forward model calculations
    ra = earth_radius + (m_to_km * altitude)
    CALL saconv( sensor_zenith_radian, altitude, local_zenith_radian )
           
    ! Compute the angle adjoint
    ra_ad = -tan(local_zenith_radian) * local_zenith_radian_ad / ra
    local_zenith_radian_ad = zero
    altitude_ad = altitude_ad + (m_to_km * ra_ad)
   
  END SUBROUTINE saconv_ad

  
!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!      VACONV  
!
! PURPOSE:
!      Elemental subroutine to compute the sensor zenith angle at a given altitude
!      for a given sensor scan angle.
!
! CALLING SEQUENCE:
!        CALL VACONV( Sensor_Scan_Radian, & ! Input
!                     Satellite_Altitude, & ! Input
!                     Altitude,           & ! Input
!                     Local_Zenith_Radian ) ! Output
!
! INPUTS:
!       Sensor_Scan_Radian:   Sensor scan angle.
!                             UNITS:      radians
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN) 
!
!       Satellite_Altitude:   The satellite altitude 
!                             UNITS:      kilometres
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!  
!       Altitude:             The altitude at which the local sensor zenith angle
!                             is required 
!                             UNITS:      metres
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Local_Zenith_Radian:  The sensor zenith angle at the supplied altitude.
!                             UNITS:      radians
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input
!                             ATTRIBUTES: INTENT(OUT)
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!
!                       Paul van Delst, 23-Nov-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE vaconv( &
    Sensor_Scan_Radian, & ! Input
    Satellite_Altitude, & ! Input
    Altitude,           & ! Input
    Local_Zenith_Radian ) ! Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: sensor_scan_radian
    REAL(fp), INTENT(IN)  :: satellite_altitude 
    REAL(fp), INTENT(IN)  :: altitude 
    REAL(fp), INTENT(OUT) :: local_zenith_radian
    ! Local variables
    REAL(fp) :: ra
    REAL(fp) :: rs

    ! Radius calculations
    ! ...Compute the radius, in km, of the point at which to calc the angle
    ra = earth_radius + (m_to_km * altitude)
    ! ...The radius of the satellite orbit, in km.
    rs = earth_radius + satellite_altitude

    ! Compute the angle
    local_zenith_radian = asin((rs / ra) * sin(sensor_scan_radian))
  
  END SUBROUTINE vaconv


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!      VACONV_TL
!
! PURPOSE:
!      Tangent-linear form of elemental subroutine to compute the sensor zenith
!      angle at a given altitude for a given sensor scan angle.
!
! CALLING SEQUENCE:
!        CALL VACONV_TL( Sensor_Scan_Radian,    &  ! FWD Input
!                        Satellite_Altitude,    &  ! FWD Input
!                        Altitude,              &  ! FWD Input
!                        Altitude_TL,           &  ! TL  Input
!                        Local_Zenith_Radian_TL )  ! TL  Output
!
! INPUTS:
!       Sensor_Scan_Radian:      Sensor scan angle.
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN) 
!
!       Satellite_Altitude:      The satellite altitude 
!                                UNITS:      kilometres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!  
!       Altitude:                The altitude at which the local sensor zenith
!                                angle is required 
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!
!       Altitude_TL:             Tangent-linear altitude. 
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!  
! OUTPUTS:
!       Local_Zenith_Radian_TL:  The tangent-linear sensor zenith angle at the
!                                supplied altitude.
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Same as input
!                                ATTRIBUTES: INTENT(OUT)
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!
!                       Paul van Delst, 23-Nov-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE vaconv_tl( &
    Sensor_Scan_Radian,    &  ! FWD Input
    Satellite_Altitude,    &  ! FWD Input
    Altitude,              &  ! FWD Input
    Altitude_TL,           &  ! TL  Input
    Local_Zenith_Radian_TL )  ! TL  Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: sensor_scan_radian 
    REAL(fp), INTENT(IN)  :: satellite_altitude 
    REAL(fp), INTENT(IN)  :: altitude 
    REAL(fp), INTENT(IN)  :: altitude_tl 
    REAL(fp), INTENT(OUT) :: local_zenith_radian_tl
    ! Local variables
    REAL(fp) :: ra, ra_tl   
    REAL(fp) :: rs   
    REAL(fp) :: local_zenith_radian

    ! Radius calculations
    ! ...Compute the radius, in km, of the point at which to calc the angle
    ra    = earth_radius + (m_to_km * altitude)
    ra_tl = m_to_km * altitude_tl
    ! ...The radius of the satellite orbit, in km.
    rs = earth_radius + satellite_altitude

    ! Compute the tangent-linear angle
    CALL vaconv( sensor_scan_radian, satellite_altitude, altitude, local_zenith_radian )
    local_zenith_radian_tl = -tan(local_zenith_radian) * ra_tl / ra
  
  END SUBROUTINE vaconv_tl


!----------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!      VACONV_AD
!
! PURPOSE:
!      Adjoint form of elemental subroutine to compute the sensor zenith
!      angle at a given altitude for a given sensor scan angle.
!
! CALLING SEQUENCE:
!        CALL VACONV_AD( Sensor_Scan_Radian,    &  ! FWD Input
!                        Satellite_Altitude,    &  ! FWD Input
!                        Altitude,              &  ! FWD Input
!                        Local_Zenith_Radian_AD )  ! TL  Output
!                        Altitude_AD,           &  ! TL  Input
!
! INPUTS:
!       Sensor_Scan_Radian:      Sensor scan angle.
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN) 
!
!       Satellite_Altitude:      The satellite altitude 
!                                UNITS:      kilometres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!  
!       Altitude:                The altitude at which the local sensor zenith
!                                angle is required 
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN)
!
!       Local_Zenith_Radian_AD:  The adjoint sensor zenith angle at the
!                                supplied altitude.
!                                *** SET TO ZERO ON EXIT ***
!                                UNITS:      radians
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Scalar or any rank
!                                ATTRIBUTES: INTENT(IN OUT)
!  
! OUTPUTS:
!       Altitude_AD:             Adjoint altitude. 
!                                *** MUST HAVE VALUE ON ENTRY ***
!                                UNITS:      metres
!                                TYPE:       REAL(fp)
!                                DIMENSION:  Same as input
!                                ATTRIBUTES: INTENT(IN OUT)
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CSU/CIRA 10-May-2006
!                       Yong.Chen@noaa.gov
!
!                       Paul van Delst, 23-Nov-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!--------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE vaconv_ad( &
    Sensor_Scan_Radian,     &  ! FWD Input
    Satellite_Altitude,     &  ! FWD Input
    Altitude,               &  ! FWD Input
    Local_Zenith_Radian_AD, &  ! AD  Input
    Altitude_AD             )  ! AD  Output
    ! Arguments
    REAL(fp), INTENT(IN)     :: sensor_scan_radian 
    REAL(fp), INTENT(IN)     :: satellite_altitude 
    REAL(fp), INTENT(IN)     :: altitude 
    REAL(fp), INTENT(IN OUT) :: local_zenith_radian_ad 
    REAL(fp), INTENT(IN OUT) :: altitude_ad
    ! Local variables
    REAL(fp) :: ra, ra_ad   
    REAL(fp) :: rs   
    REAL(fp) :: local_zenith_radian

    ! Forward model calcuations
    ra = earth_radius + (m_to_km * altitude)
    rs = earth_radius + satellite_altitude
    CALL vaconv( sensor_scan_radian, satellite_altitude, altitude, local_zenith_radian )

    ! Compute the angle adjoint
    ra_ad = -tan(local_zenith_radian) * local_zenith_radian_ad / ra
    local_zenith_radian_ad = zero
    altitude_ad = altitude_ad + (m_to_km * ra_ad)
  
  END SUBROUTINE vaconv_ad


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

  SUBROUTINE crtm_geometryinfo_version( Id )
    CHARACTER(*), INTENT(OUT) :: id
    id = module_version_id
  END SUBROUTINE crtm_geometryinfo_version


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

!--------------------------------------------------------------------------------
!
! NAME:
!       Compute_AU_ratio2
!
! PURPOSE:
!       Elemental function to compute the square of the ratio between the
!       semi-major axis of the Earth's elliptical orbit about the Sun, and
!       the actual Earth-Sun distance.
!
! CALLING SEQUENCE:
!       AU_ratio2 = Compute_AU_ratio2( Year, Month, Day )
!
! INPUTS:
!       Year:          The year in 4-digit format, e.g. 1997.
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar or same as geo input
!                      ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Month:         The month of the year (1-12).
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar or same as geo input
!                      ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Day:           The day of the month (1-28/29/30/31).
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar or same as geo input
!                      ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!       AU_Ratio2:     The square of the ratio between the semi-major axis
!                      of the Earth's elliptical orbit about the Sun, and
!                      the actual Earth-Sun distance.
!                      UNITS:      N/A
!                      TYPE:       REAL(fp)
!                      DIMENSION:  Same as inputs.
!
! PROCEDURE:
!       Adapted from eqn. 2.2.9 in 
!
!       Liou, K.N., 2002, "An Introduction to Atmospheric Radiation",
!         2nd Ed., Academic Press, San Diego, California
!
!                  __ N 
!       ( a )^2   \
!       (---)   =  >  a(i).cos(i.t) + b(i).sin(i.t)
!       ( r )     /__
!                     i=0
!
!       where
!       
!         a = semi-major axis of Earth's orbit about the Sun
!         
!         r = Earth-Sun distance at time t
!         
!             2.PI.(DoY-1)
!         t = ------------
!               Max_DoY
!
!         DoY = day of the year
!         
!         Max_DoY = number of days in the year (365 or 366 for leap year)
!                        
!       and the coefficients are taken from table 2.2 in Liou(2002):
!       
!           i     a(i)       b(i)
!         ---------------------------
!           0   1.000110       0
!           1   0.034221   0.001280
!           2   0.000719   0.000077
!
!--------------------------------------------------------------------------------

  ELEMENTAL FUNCTION compute_au_ratio2( Year, Month, Day ) RESULT( AU_ratio2 )
    ! Arguments
    INTEGER, INTENT(IN) :: year 
    INTEGER, INTENT(IN) :: month
    INTEGER, INTENT(IN) :: day  
    ! Function result
    REAL(fp) :: au_ratio2
    ! Local parameters
    INTEGER , PARAMETER :: n = 2
    REAL(fp), PARAMETER :: a(0:n) = (/ 1.000110_fp, 0.034221_fp, 0.000719_fp /)
    REAL(fp), PARAMETER :: b(0:n) = (/ zero       , 0.001280_fp, 0.000077_fp /)
    ! Local variables
    INTEGER  :: i
    REAL(fp) :: doy, max_doy, t, it

    ! Determine the days-of-year
    doy     = REAL(DayOfYear( Day, Month, Year ), fp)
    max_doy = REAL(DayOfYear( 31, 12, Year ), fp)
    
    ! Location of Earth in orbit relative to perihelion
    t = twopi * (doy-one)/max_doy
    
    ! Compute the ratio term
    au_ratio2 = zero
    DO i = 0, n
      it = REAL(i,fp) * t
      au_ratio2 = au_ratio2 + (a(i) * cos(it)) + (b(i) * sin(it))
    END DO
    
  END FUNCTION compute_au_ratio2

END MODULE crtm_geometryinfo