!
! CRTM_IR_Water_SfcOptics
!
! Module to compute the surface optical properties for WATER surfaces at
! infrared frequencies required for determining the WATER surface
! contribution to the radiative transfer.
!
! This module is provided to allow developers to "wrap" their existing
! codes inside the provided functions to simplify integration into
! the main CRTM_SfcOptics module.
!
!
! CREATION HISTORY:
! Written by: Paul van Delst, 25-Jun-2005
! paul.vandelst@noaa.gov
!
MODULE CRTM_IR_Water_SfcOptics
! -----------------
! Environment setup
! -----------------
! Module use
USE Type_Kinds, ONLY: fp
USE Message_Handler, ONLY: SUCCESS, Display_Message
USE CRTM_Parameters, ONLY: ZERO, ONE, TWO, FOUR, PI, DEGREES_TO_RADIANS
USE CRTM_SpcCoeff, ONLY: SC, SpcCoeff_IsSolar
USE CRTM_Surface_Define, ONLY: CRTM_Surface_type
USE CRTM_GeometryInfo_Define, ONLY: CRTM_GeometryInfo_type
USE CRTM_SfcOptics_Define, ONLY: CRTM_SfcOptics_type
USE CRTM_IRSSEM, ONLY: IRSSEM_type=>iVar_type, &
CRTM_Compute_IRSSEM, &
CRTM_Compute_IRSSEM_TL, &
CRTM_Compute_IRSSEM_AD
USE CRTM_IRwaterCoeff, ONLY: IRwaterC
! Disable implicit typing
IMPLICIT NONE
! ------------
! Visibilities
! ------------
! Everything private by default
PRIVATE
! Data types
PUBLIC :: iVar_type
! Science routines
PUBLIC :: Compute_IR_Water_SfcOptics
PUBLIC :: Compute_IR_Water_SfcOptics_TL
PUBLIC :: Compute_IR_Water_SfcOptics_AD
! -----------------
! Module parameters
! -----------------
CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
'$Id: CRTM_IR_Water_SfcOptics.f90 99117 2017-11-27 18:37:14Z tong.zhu@noaa.gov $'
! Coefficients for Sigma**2 in the Cox & Munk slope probability density function
REAL(fp), PARAMETER :: CM_1 = 0.003_fp, CM_2 = 5.12e-3_fp
! --------------------------------------
! Structure definition to hold forward
! variables across FWD, TL, and AD calls
! --------------------------------------
TYPE :: iVar_type
PRIVATE
! Variables in routines rough sea BRDF
REAL(fp) :: pdf ! slope distribution function
REAL(fp) :: W ! BRDF = W*pdf
REAL(fp) :: tan2_theta_f ! tan(theta_f)**2
! IRSSEM data structure
TYPE(IRSSEM_type) :: IRSSEM
END TYPE iVar_type
CONTAINS
!################################################################################
!################################################################################
!## ##
!## ## PUBLIC MODULE ROUTINES ## ##
!## ##
!################################################################################
!################################################################################
!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! Compute_IR_Water_SfcOptics
!
! PURPOSE:
! Function to compute the surface emissivity and reflectivity at infrared
! frequencies over a water surface.
!
! This function is a wrapper for third party code.
!
! CALLING SEQUENCE:
! Error_Status = Compute_IR_Water_SfcOptics( &
! Surface , & ! Input
! GeometryInfo, & ! Input
! SensorIndex , & ! Input
! ChannelIndex, & ! Output
! SfcOptics , & ! Output
! iVar ) ! Internal variable output
!
! INPUTS:
! Surface: CRTM_Surface structure containing the surface state
! data.
! UNITS: N/A
! TYPE: CRTM_Surface_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! GeometryInfo: CRTM_GeometryInfo structure containing the
! view geometry information.
! UNITS: N/A
! TYPE: CRTM_GeometryInfo_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SensorIndex: Sensor index id. This is a unique index associated
! with a (supported) sensor used to access the
! shared coefficient data for a particular sensor.
! See the ChannelIndex argument.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ChannelIndex: Channel index id. This is a unique index associated
! with a (supported) sensor channel used to access the
! shared coefficient data for a particular sensor's
! channel.
! See the SensorIndex argument.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
! SfcOptics: CRTM_SfcOptics structure containing the surface
! optical properties required for the radiative
! transfer calculation. On input the Angle component
! is assumed to contain data.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN OUT)
!
! iVar: Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
! The contents of this structure are NOT accessible
! outside of the CRTM_IR_Water_SfcOptics module.
! UNITS: N/A
! TYPE: iVar_type
! DIMENSION: Scalar
! ATTRIBUTES: 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 computation was sucessful
! == FAILURE an unrecoverable error occurred
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
!
! COMMENTS:
! Note the INTENT on the output SfcOptics argument is IN OUT rather
! than just OUT as it is assumed to contain some data upon input.
!
!:sdoc-:
!----------------------------------------------------------------------------------
FUNCTION Compute_IR_Water_SfcOptics( &
Surface , & ! Input
GeometryInfo, & ! Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
SfcOptics , & ! Output
iVar ) & ! Internal variable output
RESULT( Error_Status )
! Arguments
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface
TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo
INTEGER, INTENT(IN) :: SensorIndex
INTEGER, INTENT(IN) :: ChannelIndex
TYPE(CRTM_SfcOptics_type), INTENT(IN OUT) :: SfcOptics
TYPE(iVar_type), INTENT(IN OUT) :: iVar
! Function result
INTEGER :: Error_Status
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_IR_Water_SfcOptics'
! Local variables
INTEGER :: j, nZ, iZ
REAL(fp) :: Frequency
REAL(fp) :: Relative_Azimuth_Radian, brdf
! Set up
Error_Status = SUCCESS
! ...Short name for angle dimensions
nZ = SfcOptics%n_Angles
iZ = SfcOptics%Index_Sat_Ang
! ...Retrieve data from structures
Frequency = SC(SensorIndex)%Wavenumber(ChannelIndex)
! Compute IR sea surface emissivity
Error_Status = CRTM_Compute_IRSSEM( &
IRwaterC , & ! Input model coefficients
Surface%Wind_Speed , & ! Input
Frequency , & ! Input
SfcOptics%Angle(1:nZ) , & ! Input
iVar%IRSSEM , & ! Internal variable output
SfcOptics%Emissivity(1:nZ,1) ) ! Output
IF ( Error_Status /= SUCCESS ) THEN
CALL Display_Message( ROUTINE_NAME, &
'Error computing IR sea surface emissivity', &
Error_Status )
RETURN
END IF
! Compute the solar direct BRDF
IF ( SpcCoeff_IsSolar(SC(SensorIndex), ChannelIndex=ChannelIndex) ) THEN
IF( GeometryInfo%Source_Zenith_Radian < PI/TWO ) THEN
Relative_Azimuth_Radian = GeometryInfo%Sensor_Azimuth_Radian - &
GeometryInfo%Source_Azimuth_Radian
CALL BRDF_Rough_Sea(SC(SensorIndex)%Wavenumber(ChannelIndex), &
GeometryInfo%Source_Zenith_Radian, &
Relative_Azimuth_Radian, &
GeometryInfo%Sensor_Zenith_Radian, &
Surface%Wind_Speed, &
brdf, &
iVar)
SfcOptics%Direct_Reflectivity(1:nZ,1) = brdf
ELSE
SfcOptics%Direct_Reflectivity(1:nZ,1) = ZERO
END IF
END IF
! Surface reflectance (currently assumed to be specular ALWAYS)
DO j = 1, nZ
SfcOptics%Reflectivity(j,1,j,1) = ONE-SfcOptics%Emissivity(j,1)
END DO
END FUNCTION Compute_IR_Water_SfcOptics
!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! Compute_IR_Water_SfcOptics_TL
!
! PURPOSE:
! Function to compute the tangent-linear surface emissivity and
! reflectivity at infrared frequencies over a water surface.
!
! This function is a wrapper for third party code.
!
! CALLING SEQUENCE:
! Error_Status = Compute_IR_Water_SfcOptics_TL( &
! Surface , &
! SfcOptics , &
! Surface_TL , &
! GeometryInfo, &
! SensorIndex , &
! ChannelIndex, &
! SfcOptics_TL, &
! iVar )
!
! INPUTS:
! Surface: CRTM_Surface structure containing the surface state
! data.
! UNITS: N/A
! TYPE: CRTM_Surface_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! Surface_TL: CRTM_Surface structure containing the tangent-linear
! surface state data.
! UNITS: N/A
! TYPE: CRTM_Surface_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SfcOptics: CRTM_SfcOptics structure containing the surface
! optical properties required for the radiative
! transfer calculation.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! GeometryInfo: CRTM_GeometryInfo structure containing the
! view geometry information.
! UNITS: N/A
! TYPE: CRTM_GeometryInfo_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SensorIndex: Sensor index id. This is a unique index associated
! with a (supported) sensor used to access the
! shared coefficient data for a particular sensor.
! See the ChannelIndex argument.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ChannelIndex: Channel index id. This is a unique index associated
! with a (supported) sensor channel used to access the
! shared coefficient data for a particular sensor's
! channel.
! See the SensorIndex argument.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! iVar: Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
! The contents of this structure are NOT accessible
! outside of the CRTM_IR_Water_SfcOptics module.
! UNITS: N/A
! TYPE: iVar_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
! SfcOptics_TL: CRTM_SfcOptics structure containing the tangent-linear
! surface optical properties required for the tangent-
! linear radiative transfer calculation.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN 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 computation was sucessful
! == FAILURE an unrecoverable error occurred
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
!
! COMMENTS:
! Note the INTENT on the output SfcOptics_TL argument is IN OUT rather
! than just OUT. This is necessary because the argument may be defined
! upon input.
!
!:sdoc-:
!----------------------------------------------------------------------------------
FUNCTION Compute_IR_Water_SfcOptics_TL( &
Surface , & ! Input
SfcOptics , & ! Input
Surface_TL , & ! Input
GeometryInfo, & ! Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
SfcOptics_TL, & ! Output
iVar ) & ! Internal variable input
RESULT ( Error_Status )
! Arguments
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface_TL
TYPE(CRTM_SfcOptics_type), INTENT(IN) :: SfcOptics
TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo
INTEGER, INTENT(IN) :: SensorIndex
INTEGER, INTENT(IN) :: ChannelIndex
TYPE(CRTM_SfcOptics_type), INTENT(IN OUT) :: SfcOptics_TL
TYPE(iVar_type), INTENT(IN) :: iVar
! Function result
INTEGER :: Error_Status
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_IR_Water_SfcOptics_TL'
! Local variables
INTEGER :: j, nZ, iZ
REAL(fp) :: Relative_Azimuth_Radian, brdf_TL
! Set up
Error_Status = SUCCESS
! ...Short name for angle dimensions
nZ = SfcOptics%n_Angles
iZ = SfcOptics%Index_Sat_Ang
! Compute tangent-linear IR sea surface emissivity
Error_Status = CRTM_Compute_IRSSEM_TL( &
IRwaterC , & ! Input model coefficients
Surface_TL%Wind_Speed , & ! Input
iVar%IRSSEM , & ! Internal variable input
SfcOptics_TL%Emissivity(1:nZ,1) ) ! Output
IF ( Error_Status /= SUCCESS ) THEN
CALL Display_Message( ROUTINE_NAME, &
'Error computing Tangent_linear IR sea surface emissivity', &
Error_Status )
RETURN
END IF
! Compute the tangent-linear solar direct BRDF
IF ( SpcCoeff_IsSolar(SC(SensorIndex), ChannelIndex=ChannelIndex) ) THEN
IF( GeometryInfo%Source_Zenith_Radian < PI/TWO ) THEN
Relative_Azimuth_Radian = GeometryInfo%Sensor_Azimuth_Radian - &
GeometryInfo%Source_Azimuth_Radian
CALL BRDF_Rough_Sea_TL(Surface%Wind_Speed, &
Surface_TL%Wind_Speed, &
brdf_TL, &
iVar)
SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = brdf_TL
ELSE
SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = ZERO
END IF
END IF
! Surface reflectance (currently assumed to be specular ALWAYS)
DO j = 1, nZ
SfcOptics_TL%Reflectivity(j,1,j,1) = -SfcOptics_TL%Emissivity(j,1)
END DO
END FUNCTION Compute_IR_Water_SfcOptics_TL
!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! Compute_IR_Water_SfcOptics_AD
!
! PURPOSE:
! Function to compute the adjoint surface emissivity and
! reflectivity at infrared frequencies over a water surface.
!
! This function is a wrapper for third party code.
!
! CALLING SEQUENCE:
! Error_Status = Compute_IR_Water_SfcOptics_AD( &
! Surface , &
! SfcOptics , &
! SfcOptics_AD, &
! GeometryInfo, &
! SensorIndex , &
! ChannelIndex, &
! Surface_AD , &
! iVar )
!
! INPUTS:
! Surface: CRTM_Surface structure containing the surface state
! data.
! UNITS: N/A
! TYPE: CRTM_Surface_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SfcOptics: CRTM_SfcOptics structure containing the surface
! optical properties required for the radiative
! transfer calculation.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SfcOptics_AD: CRTM_SfcOptics structure containing the adjoint
! surface optical properties required for the adjoint
! radiative transfer calculation.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN OUT)
!
! GeometryInfo: CRTM_GeometryInfo structure containing the
! view geometry information.
! UNITS: N/A
! TYPE: CRTM_GeometryInfo_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SensorIndex: Sensor index id. This is a unique index associated
! with a (supported) sensor used to access the
! shared coefficient data for a particular sensor.
! See the ChannelIndex argument.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ChannelIndex: Channel index id. This is a unique index associated
! with a (supported) sensor channel used to access the
! shared coefficient data for a particular sensor's
! channel.
! See the SensorIndex argument.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! iVar: Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
! The contents of this structure are NOT accessible
! outside of the CRTM_IR_Water_SfcOptics module.
! UNITS: N/A
! TYPE: iVar_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
! Surface_AD: CRTM_Surface structure containing the adjoint
! surface state data.
! UNITS: N/A
! TYPE: CRTM_Surface_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN 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 computation was sucessful
! == FAILURE an unrecoverable error occurred
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
!
! COMMENTS:
! Note the INTENT on the input SfcOptics_AD argument is IN OUT rather
! than just OUT. This is necessary because components of this argument
! may need to be zeroed out upon output.
!
! Note the INTENT on the output Surface_AD argument is IN OUT rather
! than just OUT. This is necessary because the argument may be defined
! upon input.
!
!:sdoc-:
!----------------------------------------------------------------------------------
FUNCTION Compute_IR_Water_SfcOptics_AD( &
Surface , & ! Input
SfcOptics , & ! Input
SfcOptics_AD, & ! Input
GeometryInfo, & ! Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
Surface_AD , & ! Output
iVar ) & ! Internal variable input
RESULT ( Error_Status )
! Arguments
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface
TYPE(CRTM_SfcOptics_type), INTENT(IN) :: SfcOptics
TYPE(CRTM_SfcOptics_type), INTENT(IN OUT) :: SfcOptics_AD
TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo
INTEGER, INTENT(IN) :: SensorIndex
INTEGER, INTENT(IN) :: ChannelIndex
TYPE(CRTM_Surface_type), INTENT(IN OUT) :: Surface_AD
TYPE(iVar_type), INTENT(IN) :: iVar
! Function result
INTEGER :: Error_Status
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_IR_Water_SfcOptics_AD'
! Local variables
INTEGER :: j, nZ, iZ
REAL(fp) :: Relative_Azimuth_Radian, brdf_AD
! Set up
Error_Status = SUCCESS
! ...Short name for angle dimensions
nZ = SfcOptics%n_Angles
iZ = SfcOptics%Index_Sat_Ang
! Surface reflectance (currently assumed to be specular ALWAYS)
DO j = nZ, 1, -1
SfcOptics_AD%Emissivity(j,1) = SfcOptics_AD%Emissivity(j,1) - &
SfcOptics_AD%Reflectivity(j,1,j,1)
SfcOptics_AD%Reflectivity(j,1,j,1) = ZERO
END DO
! Solar direct BRDF
IF ( SpcCoeff_IsSolar(SC(SensorIndex), ChannelIndex=ChannelIndex) ) THEN
IF( GeometryInfo%Source_Zenith_Radian < PI/TWO ) THEN
Relative_Azimuth_Radian = GeometryInfo%Sensor_Azimuth_Radian - &
GeometryInfo%Source_Azimuth_Radian
brdf_AD = SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1))
SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = ZERO
CALL BRDF_Rough_Sea_AD(Surface%Wind_Speed, &
brdf_AD, &
Surface_AD%Wind_Speed, &
iVar)
END IF
SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = ZERO
END IF
! Compute sdjoint IRSSEM sea surface emissivity
Error_Status = CRTM_Compute_IRSSEM_AD( &
IRwaterC , & ! Input model coefficients
SfcOptics_AD%Emissivity(1:nZ,1), & ! Input
iVar%IRSSEM , & ! Internal Variable Input
Surface_AD%Wind_Speed ) ! Output
IF ( Error_Status /= SUCCESS ) THEN
CALL Display_Message( ROUTINE_NAME, &
'Error computing Adjoint IR sea surface emissivity', &
Error_Status )
RETURN
END IF
END FUNCTION Compute_IR_Water_SfcOptics_AD
!##################################################################################
!##################################################################################
!## ##
!## ## PRIVATE MODULE ROUTINES ## ##
!## ##
!##################################################################################
!##################################################################################
!--------------------------------------------------------------------
! Compute rough sea Bi-directional Reflectance Distribution Function (BRDF)
! for IR solar reflection
! Inputs:
! Frequency - Frequency (cm-1)
! theta_s - sun zenith angle (Radian)
! dphi - relative sun azimuth agnle, relative to the senosr's azimuth angle (Radian)
! theta_r - senor zenith angle (Radian)
! Wind_Speed - wind speed (m/s)
!
! note: dphi is such defined that if the observation direction and the sun direction
! ara in the same vertical plane, then dphi = 180 degree
!
! sun Zenith sensor
! \ | /
! theta_s | Theta_r
! \ | /
! \ | /
! \ | /
! \ | /
! \ |/
! -------------------------
!
! output:
! brdf - value of the Bi-directional Reflectance Distribution Function at
! the given condition.
! In/out:
! iVar - Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
!
! Written by Y. Han, May 5, 2009
!--------------------------------------------------------------------
SUBROUTINE BRDF_Rough_Sea(Frequency, theta_s, dphi, theta_r, Wind_Speed, &
brdf, iVar)
REAL(fp), INTENT(IN) :: Frequency
REAL(fp), INTENT(IN) :: theta_s
REAL(fp), INTENT(IN) :: dphi
REAL(fp), INTENT(IN) :: theta_r
REAL(fp), INTENT(IN) :: Wind_Speed
REAL(fp), INTENT(OUT) :: brdf
TYPE(iVar_type), INTENT(IN OUT) :: iVar
! LOCAL
REAL(fp), PARAMETER :: MIN_THETA = 1.0e-15_fp
REAL(fp) :: sin_theta_s, cos_theta_s, sin_theta_r, cos_theta_r, sin2_theta_s, &
sin2_theta_r, cos_dphi, CosSum2, sec4_theta_f, &
cos2_alpha, cos_alpha, alpha, rho
! various intermidiate variables
sin_theta_s = SIN(theta_s)
cos_theta_s = MAX(COS(theta_s), MIN_THETA) ! make sure COS(theta_s) > 0
sin_theta_r = SIN(theta_r)
cos_theta_r = MAX(COS(theta_r), MIN_THETA) ! make sure COS(theta_r) > 0
sin2_theta_s = sin_theta_s*sin_theta_s
sin2_theta_r = sin_theta_r*sin_theta_r
cos_dphi = COS(dphi)
CosSum2 = (cos_theta_s + cos_theta_r)**2
! Compute specular reflection angle alpha
cos2_alpha = (ONE + sin_theta_s*sin_theta_r*cos_dphi + cos_theta_s*cos_theta_r)/TWO
cos_alpha = SQRT( cos2_alpha )
alpha = ACOS(MIN(cos_alpha, ONE))
! Compute Fresnel_reflectance
rho = Fresnel_Reflectance(Frequency, alpha)
! Compute Tan(theta_f)**2, where theta_f is the zenith angle of the normal
! of the facet at the point of reflection.
iVar%tan2_theta_f = (sin2_theta_s + sin2_theta_r &
+ TWO*sin_theta_s*sin_theta_r*cos_dphi) &
/ CosSum2
! Compute splope probability density function
CALL Slope_pdf(iVar%tan2_theta_f, Wind_Speed, iVar%pdf)
! Compute BRDF
sec4_theta_f = (iVar%tan2_theta_f + ONE)**2
iVar%W = PI * rho * sec4_theta_f / (FOUR*cos_theta_r*cos_theta_s)
brdf = iVar%W * iVar%pdf
END SUBROUTINE BRDF_Rough_Sea
SUBROUTINE BRDF_Rough_Sea_TL(Wind_Speed, &
Wind_Speed_TL, brdf_TL, iVar)
REAL(fp), INTENT(IN) :: Wind_Speed
REAL(fp), INTENT(IN) :: Wind_Speed_TL
REAL(fp), INTENT(OUT) :: brdf_TL
TYPE(iVar_type), INTENT(IN) :: iVar
! LOCAL
REAL(fp) :: pdf_TL
CALL Slope_pdf_TL(iVar%tan2_theta_f, Wind_Speed, iVar%pdf, Wind_Speed_TL, pdf_TL)
brdf_TL = iVar%W * pdf_TL
END SUBROUTINE BRDF_Rough_Sea_TL
SUBROUTINE BRDF_Rough_Sea_AD(Wind_Speed, &
brdf_AD, Wind_Speed_AD, iVar)
REAL(fp), INTENT(IN) :: Wind_Speed
REAL(fp), INTENT(INOUT) :: brdf_AD
REAL(fp), INTENT(INOUT) :: Wind_Speed_AD
TYPE(iVar_type), INTENT(IN) :: iVar
! LOCAL
REAL(fp) :: pdf_AD
pdf_AD = iVar%W * brdf_AD
brdf_AD = ZERO
CALL Slope_pdf_AD(iVar%tan2_theta_f, Wind_Speed, iVar%pdf, pdf_AD, Wind_Speed_AD)
END SUBROUTINE BRDF_Rough_Sea_AD
! Compute facet slope distribution function (pdf)
! Inputs: tan2_theta_f - tan(theta_f)**2
! Wind_Speed (m/s)
! Outputs: pdf
SUBROUTINE Slope_pdf(tan2_theta_f, Wind_Speed, pdf)
REAL(fp), INTENT(IN) :: tan2_theta_f
REAL(fp), INTENT(IN) :: Wind_Speed
REAL(fp), INTENT(OUT) :: pdf
! Local
REAL(fp) :: Sigma2
! Cox & Munk slope probability density function
Sigma2 = CM_1 + CM_2*Wind_Speed
pdf = EXP(-tan2_theta_f / Sigma2) / (PI*Sigma2)
END SUBROUTINE Slope_pdf
SUBROUTINE Slope_pdf_TL(tan2_theta_f, Wind_Speed, pdf, Wind_Speed_TL, pdf_TL)
REAL(fp), INTENT(IN) :: tan2_theta_f
REAL(fp), INTENT(IN) :: Wind_Speed
REAL(fp), INTENT(IN) :: pdf
REAL(fp), INTENT(IN) :: Wind_Speed_TL
REAL(fp), INTENT(OUT) :: pdf_TL
! LOCAL
REAL(fp) :: Sigma2, Sigma2_TL
Sigma2 = CM_1 + CM_2*Wind_Speed
Sigma2_TL = CM_2*Wind_Speed_TL
pdf_TL = ( pdf*(tan2_theta_f/Sigma2 - ONE)/Sigma2 )*Sigma2_TL
END SUBROUTINE Slope_pdf_TL
SUBROUTINE Slope_pdf_AD(tan2_theta_f, Wind_Speed, pdf, pdf_AD, Wind_Speed_AD)
REAL(fp), INTENT(IN) :: tan2_theta_f
REAL(fp), INTENT(IN) :: Wind_Speed
REAL(fp), INTENT(IN) :: pdf
REAL(fp), INTENT(IN OUT) :: pdf_AD
REAL(fp), INTENT(IN OUT) :: Wind_Speed_AD
! LOCAL
REAL(fp) :: Sigma2, Sigma2_AD
Sigma2 = CM_1 + CM_2*Wind_Speed
Sigma2_AD = ( pdf*(tan2_theta_f/Sigma2 - ONE)/Sigma2 )*pdf_AD
pdf_AD = ZERO
Wind_Speed_AD = Wind_Speed_AD + 5.12e-3_fp*Sigma2_AD
END SUBROUTINE Slope_pdf_AD
!---------------------------------------------------------
! Compute Fresnel sea surface reflectivity
! Inputs:
! Frequency - Frequency cm-1
! Ang_i - incident angle (Radian)
! output (as a function return):
! r - Fresnel reflectivity
! Written by Y. Han, May 5, 2009
!---------------------------------------------------------
FUNCTION Fresnel_Reflectance(Frequency, Ang_i) RESULT( r )
REAL(fp), INTENT(IN) :: Frequency
REAL(fp), INTENT(IN) :: Ang_i
REAL(fp) :: r
! LOCAL
REAL(fp) :: rh, rv
COMPLEX(fp) :: CCos_i, CCos_t, n, z
! call function to compute complex refractive index
n = Ref_Index(Frequency)
! Fresnel reflectivity
z = CMPLX(SIN(Ang_i), ZERO, fp)/n
CCos_t = SQRT(CMPLX(ONE, ZERO, fp) - z*z)
CCos_i = CMPLX(COS(Ang_i), ZERO, fp)
rv = ( ABS( (n*CCos_i - CCos_t) / (n*CCos_i + CCos_t) ) )**2
rh = ( ABS( (CCos_i - n*CCos_t) / (CCos_i + n*CCos_t) ) )**2
r = (rv + rh) / TWO
END FUNCTION Fresnel_Reflectance
!------------------------------------------------------------
! Obtain IR refractive index
! Input: Frequency - wavenumber cm-1, valid range 500 - 3500 cm-1
! Return: complex refractive index
!------------------------------------------------------------
FUNCTION Ref_Index(Frequency) RESULT(ref)
REAL(fp), INTENT(IN) :: Frequency
COMPLEX(fp) :: ref
!-------------------------------------------------------------------------
! Refractive index of water from Wieliczka (1989), added
! salinity and chlorinity CORRECTIONS from Friedman (1969). The resolution
! of the Wieliczka data set is reduced to 20 cm-1. The frequancy range
! of the Friedman starts at 666.67 cm-1. For frequency < 666.67 cm-1
! the value at 666.67 is used.
!-------------------------------------------------------------------------
INTEGER, PARAMETER :: nf = 151
! array holding wavenumbers at which the refractive indexes are given
REAL(fp), PARAMETER :: freq(nf) = (/ &
500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0,&
700.0, 720.0, 740.0, 760.0, 780.0, 800.0, 820.0, 840.0, 860.0, 880.0,&
900.0, 920.0, 940.0, 960.0, 980.0, 1000.0, 1020.0, 1040.0, 1060.0, 1080.0,&
1100.0, 1120.0, 1140.0, 1160.0, 1180.0, 1200.0, 1220.0, 1240.0, 1260.0, 1280.0,&
1300.0, 1320.0, 1340.0, 1360.0, 1380.0, 1400.0, 1420.0, 1440.0, 1460.0, 1480.0,&
1500.0, 1520.0, 1540.0, 1560.0, 1580.0, 1600.0, 1620.0, 1640.0, 1660.0, 1680.0,&
1700.0, 1720.0, 1740.0, 1760.0, 1780.0, 1800.0, 1820.0, 1840.0, 1860.0, 1880.0,&
1900.0, 1920.0, 1940.0, 1960.0, 1980.0, 2000.0, 2020.0, 2040.0, 2060.0, 2080.0,&
2100.0, 2120.0, 2140.0, 2160.0, 2180.0, 2200.0, 2220.0, 2240.0, 2260.0, 2280.0,&
2300.0, 2320.0, 2340.0, 2360.0, 2380.0, 2400.0, 2420.0, 2440.0, 2460.0, 2480.0,&
2500.0, 2520.0, 2540.0, 2560.0, 2580.0, 2600.0, 2620.0, 2640.0, 2660.0, 2680.0,&
2700.0, 2720.0, 2740.0, 2760.0, 2780.0, 2800.0, 2820.0, 2840.0, 2860.0, 2880.0,&
2900.0, 2920.0, 2940.0, 2960.0, 2980.0, 3000.0, 3020.0, 3040.0, 3060.0, 3080.0,&
3100.0, 3120.0, 3140.0, 3160.0, 3180.0, 3200.0, 3220.0, 3240.0, 3260.0, 3280.0,&
3300.0, 3320.0, 3340.0, 3360.0, 3380.0, 3400.0, 3420.0, 3440.0, 3460.0, 3480.0,&
3500.0/)
! real part of the refractive index
REAL(fp), PARAMETER :: nr(nf) = (/ &
1.5300, 1.5050, 1.4752, 1.4427, 1.4111, 1.3787, 1.3468, 1.3167, 1.2887, 1.2620,&
1.2348, 1.2077, 1.1823, 1.1596, 1.1417, 1.1268, 1.1152, 1.1143, 1.1229, 1.1327,&
1.1520, 1.1630, 1.1770, 1.1920, 1.2050, 1.2170, 1.2280, 1.2371, 1.2450, 1.2523,&
1.2580, 1.2636, 1.2680, 1.2740, 1.2780, 1.2820, 1.2850, 1.2880, 1.2910, 1.2940,&
1.2970, 1.3000, 1.3020, 1.3040, 1.3070, 1.3080, 1.3100, 1.3120, 1.3150, 1.3170,&
1.3200, 1.3230, 1.3280, 1.3350, 1.3490, 1.3530, 1.3430, 1.3080, 1.2630, 1.2460,&
1.2410, 1.2540, 1.2670, 1.2760, 1.2810, 1.2900, 1.2950, 1.3000, 1.3040, 1.3070,&
1.3100, 1.3130, 1.3160, 1.3180, 1.3200, 1.3210, 1.3220, 1.3240, 1.3240, 1.3250,&
1.3260, 1.3250, 1.3250, 1.3250, 1.3250, 1.3250, 1.3260, 1.3270, 1.3280, 1.3280,&
1.3290, 1.3300, 1.3310, 1.3330, 1.3350, 1.3360, 1.3380, 1.3400, 1.3410, 1.3430,&
1.3450, 1.3480, 1.3500, 1.3520, 1.3550, 1.3580, 1.3600, 1.3610, 1.3640, 1.3660,&
1.3680, 1.3710, 1.3740, 1.3770, 1.3810, 1.3850, 1.3890, 1.3940, 1.3980, 1.4020,&
1.4070, 1.4120, 1.4180, 1.4240, 1.4310, 1.4370, 1.4440, 1.4510, 1.4590, 1.4660,&
1.4720, 1.4780, 1.4840, 1.4830, 1.4820, 1.4750, 1.4690, 1.4580, 1.4390, 1.4230,&
1.4130, 1.3960, 1.3810, 1.3560, 1.3290, 1.2970, 1.2600, 1.2170, 1.1860, 1.1590,&
1.1410/)
! imaginary part of the refractive index
REAL(fp), PARAMETER :: ni(nf) = (/ &
0.3874, 0.4031, 0.4149, 0.4213, 0.4243, 0.4227, 0.4180, 0.4105, 0.3997, 0.3877,&
0.3731, 0.3556, 0.3341, 0.3079, 0.2803, 0.2489, 0.2144, 0.1776, 0.1505, 0.1243,&
0.0970, 0.0822, 0.0693, 0.0593, 0.0517, 0.0469, 0.0432, 0.0404, 0.0384, 0.0369,&
0.0357, 0.0353, 0.0352, 0.0348, 0.0343, 0.0340, 0.0338, 0.0334, 0.0329, 0.0330,&
0.0329, 0.0326, 0.0323, 0.0320, 0.0318, 0.0317, 0.0316, 0.0317, 0.0323, 0.0329,&
0.0337, 0.0352, 0.0394, 0.0454, 0.0538, 0.0711, 0.1040, 0.1239, 0.1144, 0.0844,&
0.0539, 0.0364, 0.0264, 0.0177, 0.0145, 0.0127, 0.0117, 0.0108, 0.0104, 0.0103,&
0.0104, 0.0105, 0.0108, 0.0114, 0.0122, 0.0129, 0.0135, 0.0143, 0.0150, 0.0155,&
0.0156, 0.0155, 0.0152, 0.0148, 0.0140, 0.0132, 0.0123, 0.0115, 0.0106, 0.0098,&
0.0091, 0.0085, 0.0079, 0.0073, 0.0068, 0.0064, 0.0060, 0.0056, 0.0052, 0.0049,&
0.0046, 0.0043, 0.0041, 0.0039, 0.0038, 0.0037, 0.0036, 0.0036, 0.0037, 0.0038,&
0.0040, 0.0042, 0.0046, 0.0050, 0.0055, 0.0061, 0.0070, 0.0080, 0.0092, 0.0106,&
0.0124, 0.0146, 0.0174, 0.0197, 0.0259, 0.0350, 0.0429, 0.0476, 0.0541, 0.0655,&
0.0793, 0.0921, 0.1058, 0.1232, 0.1450, 0.1656, 0.1833, 0.1990, 0.2161, 0.2267,&
0.2369, 0.2487, 0.2606, 0.2740, 0.2819, 0.2854, 0.2826, 0.2765, 0.2637, 0.2451,&
0.2230/)
REAL(fp), PARAMETER :: df = freq(2) - freq(1)
INTEGER :: idx
REAL(fp) :: c
IF(Frequency < freq(1))THEN
ref = CMPLX( nr(1), ni(1), fp )
ELSE IF( Frequency > freq(nf) )THEN
ref = CMPLX( nr(nf), ni(nf), fp )
ELSE
! Linear interpolation
idx = INT((Frequency - freq(1))/df) + 1 ! find the starting index
c = (Frequency - freq(idx))/(freq(idx+1) - freq(idx))
ref = CMPLX( nr(idx) + c*(nr(idx+1) - nr(idx)), &
ni(idx) + c*(ni(idx+1) - ni(idx)), &
fp )
END IF
END FUNCTION Ref_Index
END MODULE CRTM_IR_Water_SfcOptics