Skip to content

Add analytic LandEM LAI/canopy Jacobians + land-surface TL/AD/K unit tests #282

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
BenjaminTJohnson wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Add analytic LandEM LAI/canopy Jacobians + land-surface TL/AD/K unit tests #282

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
017eacf
initial commit of canopy_water_content TL/AD
BenjaminTJohnson Jan 10, 2026
b479f58
updated canopy water Jacobian in MW land SfcOptics - add TL and K-mat…
BenjaminTJohnson Jan 10, 2026
ed8d422
would be useful if I actually added the tests, here they are.
BenjaminTJohnson Jan 12, 2026
acd161c
Add analytic LandEM LAI jacobians and land surface tests
BenjaminTJohnson Jan 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
201 changes: 187 additions & 14 deletions src/SfcOptics/CRTM_MW_Land_SfcOptics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ MODULE CRTM_MW_Land_SfcOptics
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 NESDIS_LandEM_Module, ONLY: NESDIS_LandEM
USE NESDIS_LandEM_Module, ONLY: NESDIS_LandEM, &
NESDIS_LandEM_LAI_Derivative
! Disable implicit typing
IMPLICIT NONE

Expand Down Expand Up @@ -78,6 +79,8 @@ MODULE CRTM_MW_Land_SfcOptics
INTEGER, PARAMETER :: DWARF_TREES_SHRUBS_GROUNDCOVER = 10
INTEGER, PARAMETER :: BARE_SOIL = 11
INTEGER, PARAMETER :: CULTIVATIONS = 12
REAL(fp), PARAMETER :: FREQUENCY_CUTOFF = 80.0_fp ! GHz
REAL(fp), PARAMETER :: DEFAULT_EMISSIVITY = 0.95_fp


! --------------------------------------
Expand Down Expand Up @@ -188,11 +191,10 @@ FUNCTION Compute_MW_Land_SfcOptics( &
INTEGER :: err_stat
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_MW_Land_SfcOptics'
REAL(fp), PARAMETER :: FREQUENCY_CUTOFF = 80.0_fp ! GHz
REAL(fp), PARAMETER :: DEFAULT_EMISSIVITY = 0.95_fp
! Local variables
CHARACTER(ML) :: msg
INTEGER :: i
REAL(fp) :: lai_eff


! Set up
Expand All @@ -219,6 +221,7 @@ FUNCTION Compute_MW_Land_SfcOptics( &

! Compute the surface optical parameters
IF ( SC(SensorIndex)%Frequency(ChannelIndex) < FREQUENCY_CUTOFF ) THEN
lai_eff = MAX(Surface%Lai + Surface%Canopy_Water_Content, ZERO)
! Frequency is low enough for the model
DO i = 1, SfcOptics%n_Angles
CALL NESDIS_LandEM(SfcOptics%Angle(i), & ! Input, Degree
Expand All @@ -227,7 +230,7 @@ FUNCTION Compute_MW_Land_SfcOptics( &
Surface%Vegetation_Fraction, & ! Input
Surface%Soil_Temperature, & ! Input, K
Surface%Land_Temperature, & ! Input, K
Surface%Lai, & ! Input, Leaf Area Index
lai_eff, & ! Input, Leaf Area Index + canopy water
Surface%Soil_Type, & ! Input, Soil Type (1 - 9)
Surface%Vegetation_Type, & ! Input, Vegetation Type (1 - 13)
ZERO, & ! Input, Snow depth, mm
Expand Down Expand Up @@ -260,11 +263,50 @@ END FUNCTION Compute_MW_Land_SfcOptics
!
! This function is a wrapper for third party code.
!
! NB: CURRENTLY THIS IS A STUB FUNCTION AS THERE ARE NO TL
! COMPONENTS IN THE MW LAND SFCOPTICS COMPUTATIONS.
! This implementation includes analytic derivatives with respect to the
! effective LAI (LAI + canopy water content).
!
! CALLING SEQUENCE:
! Error_Status = Compute_MW_Land_SfcOptics_TL( SfcOptics_TL )
! Error_Status = Compute_MW_Land_SfcOptics_TL( Surface , &
! SfcOptics , &
! Surface_TL , &
! SensorIndex , &
! ChannelIndex, &
! SfcOptics_TL )
!
! 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 forward surface
! optical properties.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_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)
!
! SensorIndex: Sensor index id.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ChannelIndex: Channel index id.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
! SfcOptics_TL: Structure containing the tangent-linear surface
Expand Down Expand Up @@ -293,25 +335,69 @@ END FUNCTION Compute_MW_Land_SfcOptics
!----------------------------------------------------------------------------------

FUNCTION Compute_MW_Land_SfcOptics_TL( &
Surface , & ! Input
SfcOptics , & ! Input
Surface_TL , & ! Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
SfcOptics_TL) & ! TL Output
RESULT ( err_stat )
! Arguments
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface
TYPE(CRTM_SfcOptics_type), INTENT(IN) :: SfcOptics
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface_TL
INTEGER, INTENT(IN) :: SensorIndex
INTEGER, INTENT(IN) :: ChannelIndex
TYPE(CRTM_SfcOptics_type), INTENT(IN OUT) :: SfcOptics_TL
! Function result
INTEGER :: err_stat
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_MW_Land_SfcOptics_TL'
! Local variables
INTEGER :: i
REAL(fp) :: lai_sum
REAL(fp) :: lai_eff
REAL(fp) :: lai_eff_tl
REAL(fp) :: d_emiss_h
REAL(fp) :: d_emiss_v
REAL(fp) :: frequency


! Set up
err_stat = SUCCESS


! Compute the tangent-linear surface optical parameters
! ***No TL models yet, so default TL output is zero***
SfcOptics_TL%Reflectivity = ZERO
SfcOptics_TL%Emissivity = ZERO
frequency = SC(SensorIndex)%Frequency(ChannelIndex)
IF ( frequency >= FREQUENCY_CUTOFF ) RETURN

lai_sum = Surface%Lai + Surface%Canopy_Water_Content
IF ( lai_sum <= ZERO ) RETURN
lai_eff = MAX(lai_sum, ZERO)
lai_eff_tl = Surface_TL%Lai + Surface_TL%Canopy_Water_Content
IF ( lai_eff_tl == ZERO ) RETURN

DO i = 1, SfcOptics%n_Angles
CALL NESDIS_LandEM_LAI_Derivative(SfcOptics%Angle(i), & ! Input, Degree
frequency, & ! Input, GHz
Surface%Soil_Moisture_Content, & ! Input, g.cm^-3
Surface%Vegetation_Fraction, & ! Input
Surface%Soil_Temperature, & ! Input, K
Surface%Land_Temperature, & ! Input, K
lai_eff, & ! Input, Effective LAI
Surface%Soil_Type, & ! Input, Soil Type (1 - 9)
Surface%Vegetation_Type, & ! Input, Vegetation Type (1 - 13)
ZERO, & ! Input, Snow depth, mm
d_emiss_h, & ! Output, H component
d_emiss_v ) ! Output, V component

SfcOptics_TL%Emissivity(i,2) = d_emiss_h * lai_eff_tl
SfcOptics_TL%Emissivity(i,1) = d_emiss_v * lai_eff_tl
SfcOptics_TL%Reflectivity(i,2,i,2) = -SfcOptics_TL%Emissivity(i,2)
SfcOptics_TL%Reflectivity(i,1,i,1) = -SfcOptics_TL%Emissivity(i,1)
END DO

END FUNCTION Compute_MW_Land_SfcOptics_TL

Expand All @@ -329,13 +415,32 @@ END FUNCTION Compute_MW_Land_SfcOptics_TL
!
! This function is a wrapper for third party code.
!
! NB: CURRENTLY THIS IS A STUB FUNCTION AS THERE ARE NO AD
! COMPONENTS IN THE MW LAND SFCOPTICS COMPUTATIONS.
! This implementation includes analytic derivatives with respect to the
! effective LAI (LAI + canopy water content).
!
! CALLING SEQUENCE:
! Error_Status = Compute_MW_Land_SfcOptics_AD( SfcOptics_AD )
! Error_Status = Compute_MW_Land_SfcOptics_AD( Surface , &
! SfcOptics , &
! SfcOptics_AD, &
! SensorIndex , &
! ChannelIndex, &
! Surface_AD )
!
! 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 forward surface
! optical properties.
! UNITS: N/A
! TYPE: CRTM_SfcOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SfcOptics_AD: Structure containing the adjoint surface optical
! properties required for the adjoint radiative
! transfer calculation.
Expand All @@ -345,6 +450,18 @@ END FUNCTION Compute_MW_Land_SfcOptics_TL
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN OUT)
!
! SensorIndex: Sensor index id.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ChannelIndex: Channel index id.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
! Error_Status: The return value is an integer defining the error status.
! The error codes are defined in the Message_Handler module.
Expand All @@ -364,25 +481,81 @@ END FUNCTION Compute_MW_Land_SfcOptics_TL
!----------------------------------------------------------------------------------

FUNCTION Compute_MW_Land_SfcOptics_AD( &
SfcOptics_AD) & ! AD Input
Surface , & ! Input
SfcOptics , & ! Input
SfcOptics_AD, & ! AD Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
Surface_AD ) & ! AD Output
RESULT( err_stat )
! Arguments
TYPE(CRTM_Surface_type), INTENT(IN) :: Surface
TYPE(CRTM_SfcOptics_type), INTENT(IN) :: SfcOptics
TYPE(CRTM_SfcOptics_type), INTENT(IN OUT) :: SfcOptics_AD
INTEGER, INTENT(IN) :: SensorIndex
INTEGER, INTENT(IN) :: ChannelIndex
TYPE(CRTM_Surface_type), INTENT(IN OUT) :: Surface_AD
! Function result
INTEGER :: err_stat
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_MW_Land_SfcOptics_AD'
! Local variables
INTEGER :: i
REAL(fp) :: lai_sum
REAL(fp) :: lai_eff
REAL(fp) :: d_emiss_h
REAL(fp) :: d_emiss_v
REAL(fp) :: emiss_h_ad
REAL(fp) :: emiss_v_ad
REAL(fp) :: lai_ad
REAL(fp) :: frequency


! Set up
err_stat = SUCCESS


! Compute the adjoint surface optical parameters
! ***No AD models yet, so there is no impact on AD result***
frequency = SC(SensorIndex)%Frequency(ChannelIndex)
IF ( frequency >= FREQUENCY_CUTOFF ) THEN
SfcOptics_AD%Reflectivity = ZERO
SfcOptics_AD%Emissivity = ZERO
RETURN
END IF

lai_sum = Surface%Lai + Surface%Canopy_Water_Content
IF ( lai_sum <= ZERO ) THEN
SfcOptics_AD%Reflectivity = ZERO
SfcOptics_AD%Emissivity = ZERO
RETURN
END IF
lai_eff = MAX(lai_sum, ZERO)
DO i = 1, SfcOptics%n_Angles
CALL NESDIS_LandEM_LAI_Derivative(SfcOptics%Angle(i), & ! Input, Degree
frequency, & ! Input, GHz
Surface%Soil_Moisture_Content, & ! Input, g.cm^-3
Surface%Vegetation_Fraction, & ! Input
Surface%Soil_Temperature, & ! Input, K
Surface%Land_Temperature, & ! Input, K
lai_eff, & ! Input, Effective LAI
Surface%Soil_Type, & ! Input, Soil Type (1 - 9)
Surface%Vegetation_Type, & ! Input, Vegetation Type (1 - 13)
ZERO, & ! Input, Snow depth, mm
d_emiss_h, & ! Output, H component
d_emiss_v ) ! Output, V component

emiss_h_ad = SfcOptics_AD%Emissivity(i,2) - SfcOptics_AD%Reflectivity(i,2,i,2)
emiss_v_ad = SfcOptics_AD%Emissivity(i,1) - SfcOptics_AD%Reflectivity(i,1,i,1)
lai_ad = (emiss_h_ad * d_emiss_h) + (emiss_v_ad * d_emiss_v)
Surface_AD%Lai = Surface_AD%Lai + lai_ad
Surface_AD%Canopy_Water_Content = Surface_AD%Canopy_Water_Content + lai_ad
SfcOptics_AD%Emissivity(i,2) = ZERO
SfcOptics_AD%Emissivity(i,1) = ZERO
SfcOptics_AD%Reflectivity(i,2,i,2) = ZERO
SfcOptics_AD%Reflectivity(i,1,i,1) = ZERO
END DO

SfcOptics_AD%Reflectivity = ZERO
SfcOptics_AD%Emissivity = ZERO

END FUNCTION Compute_MW_Land_SfcOptics_AD

Expand Down
14 changes: 12 additions & 2 deletions src/SfcOptics/CRTM_SfcOptics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1295,7 +1295,12 @@ FUNCTION CRTM_Compute_SfcOptics_TL( &
Microwave_Land: IF( Surface%Land_Coverage > ZERO) THEN

! Compute the surface optics
Error_Status = Compute_MW_Land_SfcOptics_TL( SfcOptics_TL )
Error_Status = Compute_MW_Land_SfcOptics_TL( Surface , &
SfcOptics , &
Surface_TL , &
SensorIndex , &
ChannelIndex , &
SfcOptics_TL )
IF ( Error_Status /= SUCCESS ) THEN
WRITE( Message,'("Error computing MW land SfcOptics_TL at ",&
&"channel index ",i0)' ) ChannelIndex
Expand Down Expand Up @@ -2222,7 +2227,12 @@ FUNCTION CRTM_Compute_SfcOptics_AD( &
(Reflectivity_AD(1:nZ,1:2,1:nZ,1:2)*Surface%Land_Coverage)

! Compute the surface optics adjoints
Error_Status = Compute_MW_Land_SfcOptics_AD( SfcOptics_AD )
Error_Status = Compute_MW_Land_SfcOptics_AD( Surface , &
SfcOptics , &
SfcOptics_AD , &
SensorIndex , &
ChannelIndex , &
Surface_AD )
IF ( Error_Status /= SUCCESS ) THEN
WRITE( Message,'("Error computing MW land SfcOptics_AD at ",&
&"channel index ",i0)' ) ChannelIndex
Expand Down
Loading