diff --git a/src/SfcOptics/CRTM_MW_Land_SfcOptics.f90 b/src/SfcOptics/CRTM_MW_Land_SfcOptics.f90 index 47f8248d..6951684e 100644 --- a/src/SfcOptics/CRTM_MW_Land_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_MW_Land_SfcOptics.f90 @@ -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 @@ -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 ! -------------------------------------- @@ -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 @@ -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 @@ -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 @@ -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 @@ -293,15 +335,32 @@ 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 @@ -309,9 +368,36 @@ FUNCTION Compute_MW_Land_SfcOptics_TL( & ! 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 @@ -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. @@ -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. @@ -364,15 +481,34 @@ 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 @@ -380,9 +516,46 @@ FUNCTION Compute_MW_Land_SfcOptics_AD( & ! 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 diff --git a/src/SfcOptics/CRTM_SfcOptics.f90 b/src/SfcOptics/CRTM_SfcOptics.f90 index ea3565c2..c1b7f10c 100644 --- a/src/SfcOptics/CRTM_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_SfcOptics.f90 @@ -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 @@ -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 diff --git a/src/SfcOptics/NESDIS_Emissivity/NESDIS_LandEM_Module.f90 b/src/SfcOptics/NESDIS_Emissivity/NESDIS_LandEM_Module.f90 index 3662d293..1dc74107 100644 --- a/src/SfcOptics/NESDIS_Emissivity/NESDIS_LandEM_Module.f90 +++ b/src/SfcOptics/NESDIS_Emissivity/NESDIS_LandEM_Module.f90 @@ -17,6 +17,7 @@ MODULE NESDIS_LandEM_Module PRIVATE ! Procedures PUBLIC :: NESDIS_LandEM + PUBLIC :: NESDIS_LandEM_LAI_Derivative ! Parameters PUBLIC :: ZERO PUBLIC :: POINT1 @@ -233,6 +234,168 @@ END SUBROUTINE NESDIS_LandEM + SUBROUTINE NESDIS_LandEM_LAI_Derivative(Angle, & ! Input + Frequency, & ! Input + Soil_Moisture_Content, & ! Input + Vegetation_Fraction, & ! Input + Soil_Temperature, & ! Input + t_skin, & ! Input + Lai, & ! Input + Soil_Type, & ! Input + Vegetation_Type, & ! Input + Snow_Depth, & ! Input + dEmissivity_H, & ! Output + dEmissivity_V) ! Output + ! Arguments + REAL(fp), intent(in) :: Angle + REAL(fp), intent(in) :: Frequency + REAL(fp), intent(in) :: Soil_Moisture_Content + REAL(fp), intent(in) :: Vegetation_Fraction + REAL(fp), intent(in) :: Soil_Temperature + REAL(fp), intent(in) :: t_skin + REAL(fp), intent(in) :: Lai + INTEGER, intent(in) :: Soil_Type + INTEGER, intent(in) :: Vegetation_Type + REAL(fp), intent(in) :: Snow_Depth + REAL(fp), intent(out) :: dEmissivity_H + REAL(fp), intent(out) :: dEmissivity_V + ! Local parameters + REAL(fp), PARAMETER :: rhos = 2.65_fp + REAL(fp), PARAMETER, dimension(0:9) :: frac_sand = (/ 0.80_fp, & + 0.92_fp, 0.10_fp, 0.20_fp, 0.51_fp, 0.50_fp, & + 0.35_fp, 0.60_fp, 0.42_fp, 0.92_fp /) + REAL(fp), PARAMETER, dimension(0:9) :: frac_clay = (/ 0.20_fp, & + 0.06_fp, 0.34_fp, 0.63_fp, 0.14_fp, 0.43_fp, & + 0.34_fp, 0.28_fp, 0.085_fp, 0.06_fp /) + REAL(fp), PARAMETER, dimension(0:9) :: rhob_soil = (/ 1.48_fp, & + 1.68_fp, 1.27_fp, 1.21_fp, 1.48_fp, 1.31_fp, & + 1.32_fp, 1.40_fp, 1.54_fp, 1.68_fp /) + REAL(fp), PARAMETER, dimension(0:13) :: veg_rho = (/ 0.33_fp, & + 0.40_fp, 0.40_fp, 0.40_fp, 0.40_fp, 0.40_fp, & + 0.25_fp, 0.25_fp, 0.40_fp, 0.40_fp, 0.40_fp, & + 0.40_fp, 0.33_fp, 0.33_fp /) + REAL(fp), PARAMETER, dimension(0:13) :: veg_mge = (/ 0.50_fp, & + 0.45_fp, 0.45_fp, 0.45_fp, 0.40_fp, 0.40_fp, & + 0.30_fp, 0.35_fp, 0.30_fp, 0.30_fp, 0.40_fp, & + 0.30_fp, 0.50_fp, 0.40_fp /) + REAL(fp), PARAMETER, dimension(0:13) :: leaf_th = (/ 0.07_fp, & + 0.18_fp, 0.18_fp, 0.18_fp, 0.18_fp, 0.18_fp, & + 0.12_fp, 0.12_fp, 0.12_fp, 0.12_fp, 0.12_fp, & + 0.12_fp, 0.15_fp, 0.12_fp /) + ! Local variables + REAL(fp) :: mv, veg_frac, theta, theta_t, mu, sigma, vlai, mge, rhoveg + REAL(fp) :: leaf_thick, t_soil, sand, clay, rhob, local_snow_depth + REAL(fp) :: gv, gh, ssalb_h, ssalb_v, tau_h, tau_v, tau_coeff, tv, th + REAL(fp) :: r21_h, r21_v, r23_h, r23_v, t21_h, t21_v + REAL(fp) :: alfa_v, alfa_h, kk_h, kk_v, gamma_h, gamma_v, beta_v, beta_h + REAL(fp) :: fact1, fact2, gsect0, gsect1_h, gsect1_v, gsect2_h, gsect2_v + REAL(fp) :: d_fact1, d_fact2, d_gsect2_h, d_gsect2_v + REAL(fp) :: a_h, b_h, a_v, b_v + REAL(fp) :: d_a_h, d_b_h, d_a_v, d_b_v + REAL(fp) :: d_esh_dtau, d_esv_dtau + REAL(fp) :: esh, esv + COMPLEX(fp) :: esoil, eveg, eair + + dEmissivity_H = ZERO + dEmissivity_V = ZERO + + eair = CMPLX(ONE,-ZERO,fp) + theta = Angle*PI/180.0_fp + + mv = Soil_Moisture_Content + veg_frac = Vegetation_Fraction + t_soil = Soil_Temperature + sand = frac_sand(Soil_Type) + clay = frac_clay(Soil_Type) + rhob = rhob_soil(Soil_Type) + local_snow_depth = Snow_Depth + + if ( (t_soil <= 100.0_fp .OR. t_soil >= 350.0_fp) .AND. & + (t_skin >= 100.0_fp .AND. t_skin <= 350.0_fp) ) t_soil = t_skin + + mv = MAX(MIN(mv,ONE),ZERO) + + IF (local_snow_depth > POINT1) RETURN + + veg_frac = MAX(MIN(veg_frac,ONE),ZERO) + mu = COS(theta) + sigma = POINT5 + + vlai = Lai*veg_frac + mge = veg_mge(Vegetation_Type) + rhoveg = veg_rho(Vegetation_Type) + leaf_thick = leaf_th(Vegetation_Type) + + r21_h = ZERO + r21_v = ZERO + t21_h = ONE + t21_v = ONE + + CALL Soil_Diel(Frequency, t_soil, mv, rhob, rhos, sand, clay, esoil) + theta_t = ASIN(REAL(SIN(theta)*SQRT(eair)/SQRT(esoil),fp)) + CALL Reflectance(eair, esoil, theta, theta_t, r23_v, r23_h) + CALL Roughness_Reflectance(Frequency, sigma, r23_v, r23_h) + CALL Canopy_Diel(Frequency, mge, eveg, rhoveg) + CALL Canopy_Optic(vlai,Frequency,theta,eveg,leaf_thick,gv,gh,ssalb_v,ssalb_h,tau_v,tau_h, & + tv, th) + + ! Sensitivity to LAI enters through canopy optical depth (tau). + tau_coeff = POINT5*(TWO-tv-th) + + alfa_h = SQRT((ONE - ssalb_h)/(ONE - gh*ssalb_h)) + kk_h = SQRT((ONE - ssalb_h)*(ONE - gh*ssalb_h))/mu + beta_h = (ONE - alfa_h)/(ONE + alfa_h) + gamma_h = (beta_h - r23_h)/(ONE-beta_h*r23_h) + + alfa_v = SQRT((ONE-ssalb_v)/(ONE - gv*ssalb_v)) + kk_v = SQRT((ONE-ssalb_v)*(ONE - gv*ssalb_v))/mu + beta_v = (ONE - alfa_v)/(ONE + alfa_v) + gamma_v = (beta_v -r23_v)/(ONE-beta_v*r23_v) + + fact1=gamma_h*EXP(-TWO*kk_h*tau_h) + fact2=gamma_v*EXP(-TWO*kk_v*tau_v) + + gsect0 =(EXP(C_2*frequency/t_skin) -ONE)/(EXP(C_2*frequency/t_soil) -ONE) + + gsect1_h=(ONE-r23_h)*(gsect0-ONE) + gsect2_h=((ONE-beta_h*beta_h)/(ONE-beta_h*r23_h))*EXP(-kk_h*tau_h) + + gsect1_v=(ONE-r23_v)*(gsect0-ONE) + gsect2_v=((ONE-beta_v*beta_v)/(ONE-beta_v*r23_v))*EXP(-kk_h*tau_v) + + a_h = (ONE - beta_h)*(ONE + fact1)+gsect1_h*gsect2_h + b_h = ONE-beta_h*r21_h-(beta_h-r21_h)*fact1 + + a_v = (ONE - beta_v)*(ONE + fact2)+gsect1_v*gsect2_v + b_v = ONE-beta_v*r21_v-(beta_v-r21_v)*fact2 + + esh = t21_h*a_h/b_h + esv = t21_v*a_v/b_v + + d_fact1 = -TWO*kk_h*fact1 + d_fact2 = -TWO*kk_v*fact2 + d_gsect2_h = -kk_h*gsect2_h + d_gsect2_v = -kk_h*gsect2_v + + d_a_h = (ONE - beta_h)*d_fact1 + gsect1_h*d_gsect2_h + d_b_h = -(beta_h - r21_h)*d_fact1 + d_esh_dtau = t21_h*(d_a_h*b_h - a_h*d_b_h)/(b_h*b_h) + + d_a_v = (ONE - beta_v)*d_fact2 + gsect1_v*d_gsect2_v + d_b_v = -(beta_v - r21_v)*d_fact2 + d_esv_dtau = t21_v*(d_a_v*b_v - a_v*d_b_v)/(b_v*b_v) + + dEmissivity_H = d_esh_dtau * tau_coeff * veg_frac + dEmissivity_V = d_esv_dtau * tau_coeff * veg_frac + + if (esh < EMISSH_DEFAULT .OR. esh > ONE) dEmissivity_H = ZERO + if (esv < EMISSV_DEFAULT .OR. esv > ONE) dEmissivity_V = ZERO + + END SUBROUTINE NESDIS_LandEM_LAI_Derivative + + + + subroutine SnowEM_Default(frequency,ts, depth,Emissivity_V,Emissivity_H) @@ -314,10 +477,11 @@ end subroutine SnowEM_Default subroutine Canopy_Optic(vlai,frequency,theta,esv,d,gv,gh,& - ssalb_v,ssalb_h,tau_v, tau_h) + ssalb_v,ssalb_h,tau_v, tau_h, tv_out, th_out) REAL(fp) :: frequency,theta,d,vlai,ssalb_v,ssalb_h,tau_v,tau_h,gv, gh, mu + REAL(fp), OPTIONAL, INTENT(OUT) :: tv_out, th_out COMPLEX(fp) :: ix,k0,kz0,kz1,rhc,rvc,esv,expval1,factt,factrvc,factrhc REAL(fp) :: rh,rv,th,tv REAL(fp), PARAMETER :: threshold = 0.999_fp @@ -352,6 +516,9 @@ subroutine Canopy_Optic(vlai,frequency,theta,esv,d,gv,gh,& ssalb_v = MIN((rv+rh)/(TWO-tv-th),threshold) ssalb_h = ssalb_v + IF (PRESENT(tv_out)) tv_out = tv + IF (PRESENT(th_out)) th_out = th + end subroutine Canopy_Optic diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 692dad4c..01cc803c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -234,6 +234,24 @@ list( APPEND OMPoverChannels_Sensor_Ids atms_n21 ) +list( APPEND Land_Surface_Sensor_Ids + atms_n21 +) + +list( APPEND Land_Surface_TL_Tests + Canopy_Water + LAI +) + +list( APPEND Land_Surface_K_Tests + Canopy_Water + LAI +) + +list( APPEND Land_Surface_AD_Tests + LAI_Canopy_Water +) + list (APPEND common_tests Simple AOD @@ -294,6 +312,39 @@ add_test(NAME test_Unit_TL_TEST COMMAND $) set_tests_properties(test_Unit_TL_TEST PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") +foreach(testtype IN LISTS Land_Surface_TL_Tests) + string(TOLOWER ${testtype} testtype_lower) + add_executable(Unit_TL_${testtype}_TEST mains/unit/Unit_Test/test_TL_${testtype_lower}.f90) + target_link_libraries(Unit_TL_${testtype}_TEST PRIVATE crtm) + foreach(sensor_id IN LISTS Land_Surface_Sensor_Ids) + add_test(NAME test_Unit_TL_${testtype}_${sensor_id} + COMMAND $ "${sensor_id}") + set_tests_properties(test_Unit_TL_${testtype}_${sensor_id} PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + endforeach() +endforeach() + +foreach(testtype IN LISTS Land_Surface_K_Tests) + string(TOLOWER ${testtype} testtype_lower) + add_executable(Unit_K_Matrix_${testtype}_TEST mains/unit/Unit_Test/test_k_matrix_${testtype_lower}.f90) + target_link_libraries(Unit_K_Matrix_${testtype}_TEST PRIVATE crtm) + foreach(sensor_id IN LISTS Land_Surface_Sensor_Ids) + add_test(NAME test_Unit_K_Matrix_${testtype}_${sensor_id} + COMMAND $ "${sensor_id}") + set_tests_properties(test_Unit_K_Matrix_${testtype}_${sensor_id} PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + endforeach() +endforeach() + +foreach(testtype IN LISTS Land_Surface_AD_Tests) + string(TOLOWER ${testtype} testtype_lower) + add_executable(Unit_AD_${testtype}_TEST mains/unit/Unit_Test/test_AD_${testtype_lower}.f90) + target_link_libraries(Unit_AD_${testtype}_TEST PRIVATE crtm) + foreach(sensor_id IN LISTS Land_Surface_Sensor_Ids) + add_test(NAME test_Unit_AD_${testtype}_${sensor_id} + COMMAND $ "${sensor_id}") + set_tests_properties(test_Unit_AD_${testtype}_${sensor_id} PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + endforeach() +endforeach() + add_executable(Unit_Aerosol_Bypass mains/unit/Unit_Test/test_Aerosol_Bypass.f90) target_link_libraries(Unit_Aerosol_Bypass PRIVATE crtm) add_test(NAME test_Unit_Aerosol_Bypass diff --git a/test/mains/unit/Unit_Test/test_AD_lai_canopy_water.f90 b/test/mains/unit/Unit_Test/test_AD_lai_canopy_water.f90 new file mode 100644 index 00000000..7dbbedcc --- /dev/null +++ b/test/mains/unit/Unit_Test/test_AD_lai_canopy_water.f90 @@ -0,0 +1,224 @@ +! +! test_AD_lai_canopy_water.f90 +! +! Program to test the CRTM adjoint response for LAI and canopy water content. +! +PROGRAM test_AD_lai_canopy_water + + ! ============================================================================ + ! **** ENVIRONMENT SETUP FOR RTM USAGE **** + ! + USE CRTM_Module + ! Disable all implicit typing + IMPLICIT NONE + ! ============================================================================ + + + ! ---------- + ! Parameters + ! ---------- + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_AD_lai_canopy_water' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + INTEGER, PARAMETER :: N_PROFILES = 2 + INTEGER, PARAMETER :: N_LAYERS = 92 + INTEGER, PARAMETER :: N_ABSORBERS = 2 + INTEGER, PARAMETER :: N_CLOUDS = 0 + INTEGER, PARAMETER :: N_AEROSOLS = 0 + INTEGER, PARAMETER :: N_SENSORS = 1 + INTEGER, PARAMETER :: TEST_CHANNEL = 1 + INTEGER, PARAMETER :: TEST_PROFILE = 2 + REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp + REAL(fp), PARAMETER :: SCAN_ANGLE = 26.37293341421_fp + REAL(fp), PARAMETER :: TOLERANCE = 1.0e-5_fp + + + ! --------- + ! Variables + ! --------- + CHARACTER(256) :: Message + CHARACTER(256) :: Version + CHARACTER(256) :: Sensor_Id + INTEGER :: Error_Status + INTEGER :: Allocate_Status + INTEGER :: n_Channels + REAL(fp) :: Perturb_Lai + REAL(fp) :: Perturb_Canopy + REAL(fp) :: Dot_TL + REAL(fp) :: Dot_AD + REAL(fp) :: Dot_Diff + REAL(fp) :: Dot_Norm + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(N_SENSORS) + TYPE(CRTM_Geometry_type) :: Geometry(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES) + TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_TL(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_AD(:,:) + TYPE(CRTM_Atmosphere_type) :: Atmosphere_TL(N_PROFILES) + TYPE(CRTM_Surface_type) :: Surface_TL(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atmosphere_AD(N_PROFILES) + TYPE(CRTM_Surface_type) :: Surface_AD(N_PROFILES) + + + CALL CRTM_Version( Version ) + CALL Program_Message( PROGRAM_NAME, & + 'Program to test LAI and canopy water AD response.', & + 'CRTM Version: '//TRIM(Version) ) + + CALL GET_COMMAND_ARGUMENT(1, Sensor_Id) + IF ( LEN_TRIM(Sensor_Id) == 0 ) Sensor_Id = 'atms_n21' + Sensor_Id = ADJUSTL(Sensor_Id) + WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id) + + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo , & + File_Path=COEFFICIENTS_PATH ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error initializing CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo)) + + ALLOCATE( RTSolution( n_Channels, N_PROFILES ), & + RTSolution_TL( n_Channels, N_PROFILES ), & + RTSolution_AD( n_Channels, N_PROFILES ), & + STAT = Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating structure arrays' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN + Message = 'Error allocating CRTM Atmosphere structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atmosphere_TL, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_TL)) ) THEN + Message = 'Error allocating CRTM Atmosphere_TL structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atmosphere_AD, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_AD)) ) THEN + Message = 'Error allocating CRTM Atmosphere_AD structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL Load_Atm_Data() + CALL Load_Sfc_Data() + + CALL CRTM_Geometry_SetValue( Geometry, & + Sensor_Zenith_Angle = ZENITH_ANGLE, & + Sensor_Scan_Angle = SCAN_ANGLE ) + + Sfc(:)%Land_Coverage = 1.0_fp + Sfc(:)%Water_Coverage = 0.0_fp + Sfc(:)%Snow_Coverage = 0.0_fp + Sfc(:)%Ice_Coverage = 0.0_fp + Sfc(:)%Canopy_Water_Content = 0.15_fp + Sfc(:)%Lai = 2.0_fp + + CALL CRTM_Atmosphere_Zero( Atmosphere_TL ) + CALL CRTM_Surface_Zero( Surface_TL ) + Perturb_Lai = 0.01_fp + Perturb_Canopy = -0.02_fp + Surface_TL(TEST_PROFILE)%Lai = Perturb_Lai + Surface_TL(TEST_PROFILE)%Canopy_Water_Content = Perturb_Canopy + + Error_Status = CRTM_Tangent_Linear( Atm , & + Sfc , & + Atmosphere_TL , & + Surface_TL , & + Geometry , & + ChannelInfo , & + RTSolution , & + RTSolution_TL ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Tangent-linear Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Zero( Atmosphere_AD ) + CALL CRTM_Surface_Zero( Surface_AD ) + RTSolution_AD%Radiance = ZERO + RTSolution_AD%Brightness_Temperature = ZERO + RTSolution_AD(TEST_CHANNEL,TEST_PROFILE)%Radiance = ONE + + Error_Status = CRTM_Adjoint( Atm , & + Sfc , & + RTSolution_AD, & + Geometry, & + ChannelInfo, & + Atmosphere_AD, & + Surface_AD, & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Adjoint Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Dot_TL = RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + Dot_AD = (Surface_AD(TEST_PROFILE)%Lai * Surface_TL(TEST_PROFILE)%Lai) + & + (Surface_AD(TEST_PROFILE)%Canopy_Water_Content * Surface_TL(TEST_PROFILE)%Canopy_Water_Content) + Dot_Diff = ABS(Dot_TL - Dot_AD) + Dot_Norm = MAX(ABS(Dot_TL), ONE) + + WRITE(*,*) 'TL dot:', Dot_TL + WRITE(*,*) 'AD dot:', Dot_AD + WRITE(*,*) 'Diff:', Dot_Diff + + Error_Status = CRTM_Destroy( ChannelInfo ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error destroying CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Destroy(Atmosphere_AD) + CALL CRTM_Atmosphere_Destroy(Atmosphere_TL) + CALL CRTM_Atmosphere_Destroy(Atm) + + DEALLOCATE( RTSolution, RTSolution_TL, RTSolution_AD, & + STAT = Allocate_Status ) + + IF ( ABS(Dot_TL) <= 0.0_fp ) THEN + Message = 'TL dot product is zero for LAI/canopy perturbations' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + IF ( (Dot_Diff / Dot_Norm) < TOLERANCE ) THEN + STOP 0 + END IF + + WRITE(*,*) 'FAIL diff/norm=', Dot_Diff / Dot_Norm, ' TOL=', TOLERANCE + STOP 1 + +CONTAINS + + INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Load_Sfc_Data.inc' + +END PROGRAM test_AD_lai_canopy_water diff --git a/test/mains/unit/Unit_Test/test_TL_canopy_water.f90 b/test/mains/unit/Unit_Test/test_TL_canopy_water.f90 new file mode 100644 index 00000000..2084662f --- /dev/null +++ b/test/mains/unit/Unit_Test/test_TL_canopy_water.f90 @@ -0,0 +1,199 @@ +! +! test_TL_canopy_water.f90 +! +! Program to test the CRTM tangent-linear response to canopy water content. +! +PROGRAM test_TL_canopy_water + + ! ============================================================================ + ! **** ENVIRONMENT SETUP FOR RTM USAGE **** + ! + USE CRTM_Module + ! Disable all implicit typing + IMPLICIT NONE + ! ============================================================================ + + + ! ---------- + ! Parameters + ! ---------- + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_TL_canopy_water' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + INTEGER, PARAMETER :: N_PROFILES = 2 + INTEGER, PARAMETER :: N_LAYERS = 92 + INTEGER, PARAMETER :: N_ABSORBERS = 2 + INTEGER, PARAMETER :: N_CLOUDS = 0 + INTEGER, PARAMETER :: N_AEROSOLS = 0 + INTEGER, PARAMETER :: N_SENSORS = 1 + INTEGER, PARAMETER :: TEST_CHANNEL = 1 + INTEGER, PARAMETER :: TEST_PROFILE = 2 + REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp + REAL(fp), PARAMETER :: SCAN_ANGLE = 26.37293341421_fp + REAL(fp), PARAMETER :: TOLERANCE = 0.2_fp + + + ! --------- + ! Variables + ! --------- + CHARACTER(256) :: Message + CHARACTER(256) :: Version + CHARACTER(256) :: Sensor_Id + INTEGER :: Error_Status + INTEGER :: Allocate_Status + INTEGER :: n_Channels + REAL(fp) :: Perturbation + REAL(fp) :: Ratio + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(N_SENSORS) + TYPE(CRTM_Geometry_type) :: Geometry(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES) + TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_Perturb(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_TL(:,:) + TYPE(CRTM_Atmosphere_type) :: Atmosphere_TL(N_PROFILES) + TYPE(CRTM_Surface_type) :: Surface_TL(N_PROFILES) + + + CALL CRTM_Version( Version ) + CALL Program_Message( PROGRAM_NAME, & + 'Program to test canopy water content TL response.', & + 'CRTM Version: '//TRIM(Version) ) + + CALL GET_COMMAND_ARGUMENT(1, Sensor_Id) + IF ( LEN_TRIM(Sensor_Id) == 0 ) Sensor_Id = 'atms_n21' + Sensor_Id = ADJUSTL(Sensor_Id) + WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id) + + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo , & + File_Path=COEFFICIENTS_PATH ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error initializing CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo)) + + ALLOCATE( RTSolution( n_Channels, N_PROFILES ), & + RTSolution_Perturb( n_Channels, N_PROFILES ), & + RTSolution_TL( n_Channels, N_PROFILES ), & + STAT = Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating structure arrays' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN + Message = 'Error allocating CRTM Atmosphere structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atmosphere_TL, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_TL)) ) THEN + Message = 'Error allocating CRTM Atmosphere_TL structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL Load_Atm_Data() + CALL Load_Sfc_Data() + + CALL CRTM_Geometry_SetValue( Geometry, & + Sensor_Zenith_Angle = ZENITH_ANGLE, & + Sensor_Scan_Angle = SCAN_ANGLE ) + + Sfc(:)%Land_Coverage = 1.0_fp + Sfc(:)%Water_Coverage = 0.0_fp + Sfc(:)%Snow_Coverage = 0.0_fp + Sfc(:)%Ice_Coverage = 0.0_fp + Sfc(:)%Canopy_Water_Content = 0.15_fp + Sfc(:)%Lai = 2.0_fp + + CALL CRTM_Atmosphere_Zero( Atmosphere_TL ) + CALL CRTM_Surface_Zero( Surface_TL ) + Perturbation = 0.01_fp + Surface_TL(TEST_PROFILE)%Canopy_Water_Content = Perturbation + + Error_Status = CRTM_Tangent_Linear( Atm , & + Sfc , & + Atmosphere_TL , & + Surface_TL , & + Geometry , & + ChannelInfo , & + RTSolution , & + RTSolution_TL ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Tangent-linear Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Sfc(TEST_PROFILE)%Canopy_Water_Content = Sfc(TEST_PROFILE)%Canopy_Water_Content + Perturbation + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution_Perturb ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in perturbed CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + IF ( ABS(RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance) <= 0.0_fp ) THEN + WRITE(*,*) 'TL radiance:', RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + Message = 'TL radiance is zero for canopy water perturbation' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Ratio = ( RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Radiance - & + RTSolution(TEST_CHANNEL,TEST_PROFILE)%Radiance ) / & + RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'Base radiance:', RTSolution(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'Perturb radiance:', RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'TL radiance:', RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'Ratio: ', Ratio + + Error_Status = CRTM_Destroy( ChannelInfo ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error destroying CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Destroy(Atmosphere_TL) + CALL CRTM_Atmosphere_Destroy(Atm) + + DEALLOCATE( RTSolution, RTSolution_TL, RTSolution_Perturb, & + STAT = Allocate_Status ) + + IF ( ABS(1.0_fp - Ratio) < TOLERANCE ) THEN + STOP 0 + END IF + + WRITE(*,*) 'FAIL abs(1 - Ratio)=', ABS(1.0_fp - Ratio), ' TOL=', TOLERANCE + STOP 1 + +CONTAINS + + INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Load_Sfc_Data.inc' + +END PROGRAM test_TL_canopy_water diff --git a/test/mains/unit/Unit_Test/test_TL_lai.f90 b/test/mains/unit/Unit_Test/test_TL_lai.f90 new file mode 100644 index 00000000..5b4ff722 --- /dev/null +++ b/test/mains/unit/Unit_Test/test_TL_lai.f90 @@ -0,0 +1,199 @@ +! +! test_TL_lai.f90 +! +! Program to test the CRTM tangent-linear response to LAI. +! +PROGRAM test_TL_lai + + ! ============================================================================ + ! **** ENVIRONMENT SETUP FOR RTM USAGE **** + ! + USE CRTM_Module + ! Disable all implicit typing + IMPLICIT NONE + ! ============================================================================ + + + ! ---------- + ! Parameters + ! ---------- + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_TL_lai' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + INTEGER, PARAMETER :: N_PROFILES = 2 + INTEGER, PARAMETER :: N_LAYERS = 92 + INTEGER, PARAMETER :: N_ABSORBERS = 2 + INTEGER, PARAMETER :: N_CLOUDS = 0 + INTEGER, PARAMETER :: N_AEROSOLS = 0 + INTEGER, PARAMETER :: N_SENSORS = 1 + INTEGER, PARAMETER :: TEST_CHANNEL = 1 + INTEGER, PARAMETER :: TEST_PROFILE = 2 + REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp + REAL(fp), PARAMETER :: SCAN_ANGLE = 26.37293341421_fp + REAL(fp), PARAMETER :: TOLERANCE = 0.2_fp + + + ! --------- + ! Variables + ! --------- + CHARACTER(256) :: Message + CHARACTER(256) :: Version + CHARACTER(256) :: Sensor_Id + INTEGER :: Error_Status + INTEGER :: Allocate_Status + INTEGER :: n_Channels + REAL(fp) :: Perturbation + REAL(fp) :: Ratio + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(N_SENSORS) + TYPE(CRTM_Geometry_type) :: Geometry(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES) + TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_Perturb(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_TL(:,:) + TYPE(CRTM_Atmosphere_type) :: Atmosphere_TL(N_PROFILES) + TYPE(CRTM_Surface_type) :: Surface_TL(N_PROFILES) + + + CALL CRTM_Version( Version ) + CALL Program_Message( PROGRAM_NAME, & + 'Program to test LAI TL response.', & + 'CRTM Version: '//TRIM(Version) ) + + CALL GET_COMMAND_ARGUMENT(1, Sensor_Id) + IF ( LEN_TRIM(Sensor_Id) == 0 ) Sensor_Id = 'atms_n21' + Sensor_Id = ADJUSTL(Sensor_Id) + WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id) + + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo , & + File_Path=COEFFICIENTS_PATH ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error initializing CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo)) + + ALLOCATE( RTSolution( n_Channels, N_PROFILES ), & + RTSolution_Perturb( n_Channels, N_PROFILES ), & + RTSolution_TL( n_Channels, N_PROFILES ), & + STAT = Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating structure arrays' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN + Message = 'Error allocating CRTM Atmosphere structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atmosphere_TL, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_TL)) ) THEN + Message = 'Error allocating CRTM Atmosphere_TL structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL Load_Atm_Data() + CALL Load_Sfc_Data() + + CALL CRTM_Geometry_SetValue( Geometry, & + Sensor_Zenith_Angle = ZENITH_ANGLE, & + Sensor_Scan_Angle = SCAN_ANGLE ) + + Sfc(:)%Land_Coverage = 1.0_fp + Sfc(:)%Water_Coverage = 0.0_fp + Sfc(:)%Snow_Coverage = 0.0_fp + Sfc(:)%Ice_Coverage = 0.0_fp + Sfc(:)%Canopy_Water_Content = 0.0_fp + Sfc(:)%Lai = 2.0_fp + + CALL CRTM_Atmosphere_Zero( Atmosphere_TL ) + CALL CRTM_Surface_Zero( Surface_TL ) + Perturbation = 0.01_fp + Surface_TL(TEST_PROFILE)%Lai = Perturbation + + Error_Status = CRTM_Tangent_Linear( Atm , & + Sfc , & + Atmosphere_TL , & + Surface_TL , & + Geometry , & + ChannelInfo , & + RTSolution , & + RTSolution_TL ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Tangent-linear Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Sfc(TEST_PROFILE)%Lai = Sfc(TEST_PROFILE)%Lai + Perturbation + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution_Perturb ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in perturbed CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + IF ( ABS(RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance) <= 0.0_fp ) THEN + WRITE(*,*) 'TL radiance:', RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + Message = 'TL radiance is zero for LAI perturbation' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Ratio = ( RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Radiance - & + RTSolution(TEST_CHANNEL,TEST_PROFILE)%Radiance ) / & + RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'Base radiance:', RTSolution(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'Perturb radiance:', RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'TL radiance:', RTSolution_TL(TEST_CHANNEL,TEST_PROFILE)%Radiance + WRITE(*,*) 'Ratio: ', Ratio + + Error_Status = CRTM_Destroy( ChannelInfo ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error destroying CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Destroy(Atmosphere_TL) + CALL CRTM_Atmosphere_Destroy(Atm) + + DEALLOCATE( RTSolution, RTSolution_TL, RTSolution_Perturb, & + STAT = Allocate_Status ) + + IF ( ABS(1.0_fp - Ratio) < TOLERANCE ) THEN + STOP 0 + END IF + + WRITE(*,*) 'FAIL abs(1 - Ratio)=', ABS(1.0_fp - Ratio), ' TOL=', TOLERANCE + STOP 1 + +CONTAINS + + INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Load_Sfc_Data.inc' + +END PROGRAM test_TL_lai diff --git a/test/mains/unit/Unit_Test/test_k_matrix_canopy_water.f90 b/test/mains/unit/Unit_Test/test_k_matrix_canopy_water.f90 new file mode 100644 index 00000000..f190c701 --- /dev/null +++ b/test/mains/unit/Unit_Test/test_k_matrix_canopy_water.f90 @@ -0,0 +1,210 @@ +! +! test_k_matrix_canopy_water.f90 +! +! Program to test the CRTM K-matrix response to canopy water content. +! +PROGRAM test_k_matrix_canopy_water + + ! ============================================================================ + ! **** ENVIRONMENT SETUP FOR RTM USAGE **** + ! + USE CRTM_Module + ! Disable all implicit typing + IMPLICIT NONE + ! ============================================================================ + + + ! ---------- + ! Parameters + ! ---------- + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_k_matrix_canopy_water' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + INTEGER, PARAMETER :: N_PROFILES = 2 + INTEGER, PARAMETER :: N_LAYERS = 92 + INTEGER, PARAMETER :: N_ABSORBERS = 2 + INTEGER, PARAMETER :: N_CLOUDS = 0 + INTEGER, PARAMETER :: N_AEROSOLS = 0 + INTEGER, PARAMETER :: N_SENSORS = 1 + INTEGER, PARAMETER :: TEST_CHANNEL = 1 + INTEGER, PARAMETER :: TEST_PROFILE = 2 + REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp + REAL(fp), PARAMETER :: SCAN_ANGLE = 26.37293341421_fp + REAL(fp), PARAMETER :: TOLERANCE = 0.2_fp + + + ! --------- + ! Variables + ! --------- + CHARACTER(256) :: Message + CHARACTER(256) :: Version + CHARACTER(256) :: Sensor_Id + INTEGER :: Error_Status + INTEGER :: Allocate_Status + INTEGER :: n_Channels + REAL(fp) :: Perturbation + REAL(fp) :: Ratio + REAL(fp) :: FD_Deriv + REAL(fp) :: K_Deriv + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(N_SENSORS) + TYPE(CRTM_Geometry_type) :: Geometry(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES) + TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_Perturb(:,:) + + TYPE(CRTM_Atmosphere_type), ALLOCATABLE :: Atmosphere_K(:,:) + TYPE(CRTM_Surface_type) , ALLOCATABLE :: Surface_K(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_K(:,:) + + + CALL CRTM_Version( Version ) + CALL Program_Message( PROGRAM_NAME, & + 'Program to test canopy water content K-matrix response.', & + 'CRTM Version: '//TRIM(Version) ) + + CALL GET_COMMAND_ARGUMENT(1, Sensor_Id) + IF ( LEN_TRIM(Sensor_Id) == 0 ) Sensor_Id = 'atms_n21' + Sensor_Id = ADJUSTL(Sensor_Id) + WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id) + + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo , & + File_Path=COEFFICIENTS_PATH ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error initializing CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo)) + + ALLOCATE( RTSolution( n_Channels, N_PROFILES ), & + RTSolution_Perturb( n_Channels, N_PROFILES ), & + RTSolution_K( n_Channels, N_PROFILES ), & + Atmosphere_K( n_Channels, N_PROFILES ), & + Surface_K( n_Channels, N_PROFILES ), & + STAT = Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating structure arrays' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN + Message = 'Error allocating CRTM Atmosphere structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atmosphere_K, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_K)) ) THEN + Message = 'Error allocating CRTM Atmosphere_K structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL Load_Atm_Data() + CALL Load_Sfc_Data() + + CALL CRTM_Geometry_SetValue( Geometry, & + Sensor_Zenith_Angle = ZENITH_ANGLE, & + Sensor_Scan_Angle = SCAN_ANGLE ) + + Sfc(:)%Land_Coverage = 1.0_fp + Sfc(:)%Water_Coverage = 0.0_fp + Sfc(:)%Snow_Coverage = 0.0_fp + Sfc(:)%Ice_Coverage = 0.0_fp + Sfc(:)%Canopy_Water_Content = 0.15_fp + Sfc(:)%Lai = 2.0_fp + + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Perturbation = 0.01_fp + Sfc(TEST_PROFILE)%Canopy_Water_Content = Sfc(TEST_PROFILE)%Canopy_Water_Content + Perturbation + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution_Perturb ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in perturbed CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + Sfc(TEST_PROFILE)%Canopy_Water_Content = Sfc(TEST_PROFILE)%Canopy_Water_Content - Perturbation + + CALL CRTM_Atmosphere_Zero( Atmosphere_K ) + CALL CRTM_Surface_Zero( Surface_K ) + + RTSolution_K%Radiance = ZERO + RTSolution_K%Brightness_Temperature = ONE + + Error_Status = CRTM_K_Matrix( Atm , & + Sfc , & + RTSolution_K, & + Geometry , & + ChannelInfo , & + Atmosphere_K, & + Surface_K , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM K_Matrix Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + FD_Deriv = ( RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature - & + RTSolution(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature ) / Perturbation + K_Deriv = Surface_K(TEST_CHANNEL,TEST_PROFILE)%Canopy_Water_Content + + IF ( ABS(K_Deriv) <= 0.0_fp ) THEN + WRITE(*,*) 'K-matrix canopy water derivative:', K_Deriv + Message = 'K-matrix canopy water derivative is zero' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Ratio = FD_Deriv / K_Deriv + WRITE(*,*) 'Base BT:', RTSolution(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature + WRITE(*,*) 'Perturb BT:', RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature + WRITE(*,*) 'FD dTB/dx:', FD_Deriv + WRITE(*,*) 'K dTB/dx:', K_Deriv + WRITE(*,*) 'Ratio: ', Ratio + + Error_Status = CRTM_Destroy( ChannelInfo ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error destroying CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Destroy(Atmosphere_K) + CALL CRTM_Atmosphere_Destroy(Atm) + + DEALLOCATE( RTSolution, RTSolution_Perturb, RTSolution_K, & + Atmosphere_K, Surface_K, STAT = Allocate_Status ) + + IF ( ABS(1.0_fp - Ratio) < TOLERANCE ) THEN + STOP 0 + END IF + + WRITE(*,*) 'FAIL abs(1 - Ratio)=', ABS(1.0_fp - Ratio), ' TOL=', TOLERANCE + STOP 1 + +CONTAINS + + INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Load_Sfc_Data.inc' + +END PROGRAM test_k_matrix_canopy_water diff --git a/test/mains/unit/Unit_Test/test_k_matrix_lai.f90 b/test/mains/unit/Unit_Test/test_k_matrix_lai.f90 new file mode 100644 index 00000000..d1d41f95 --- /dev/null +++ b/test/mains/unit/Unit_Test/test_k_matrix_lai.f90 @@ -0,0 +1,210 @@ +! +! test_k_matrix_lai.f90 +! +! Program to test the CRTM K-matrix response to LAI. +! +PROGRAM test_k_matrix_lai + + ! ============================================================================ + ! **** ENVIRONMENT SETUP FOR RTM USAGE **** + ! + USE CRTM_Module + ! Disable all implicit typing + IMPLICIT NONE + ! ============================================================================ + + + ! ---------- + ! Parameters + ! ---------- + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_k_matrix_lai' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + INTEGER, PARAMETER :: N_PROFILES = 2 + INTEGER, PARAMETER :: N_LAYERS = 92 + INTEGER, PARAMETER :: N_ABSORBERS = 2 + INTEGER, PARAMETER :: N_CLOUDS = 0 + INTEGER, PARAMETER :: N_AEROSOLS = 0 + INTEGER, PARAMETER :: N_SENSORS = 1 + INTEGER, PARAMETER :: TEST_CHANNEL = 1 + INTEGER, PARAMETER :: TEST_PROFILE = 2 + REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp + REAL(fp), PARAMETER :: SCAN_ANGLE = 26.37293341421_fp + REAL(fp), PARAMETER :: TOLERANCE = 0.2_fp + + + ! --------- + ! Variables + ! --------- + CHARACTER(256) :: Message + CHARACTER(256) :: Version + CHARACTER(256) :: Sensor_Id + INTEGER :: Error_Status + INTEGER :: Allocate_Status + INTEGER :: n_Channels + REAL(fp) :: Perturbation + REAL(fp) :: Ratio + REAL(fp) :: FD_Deriv + REAL(fp) :: K_Deriv + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(N_SENSORS) + TYPE(CRTM_Geometry_type) :: Geometry(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES) + TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_Perturb(:,:) + + TYPE(CRTM_Atmosphere_type), ALLOCATABLE :: Atmosphere_K(:,:) + TYPE(CRTM_Surface_type) , ALLOCATABLE :: Surface_K(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_K(:,:) + + + CALL CRTM_Version( Version ) + CALL Program_Message( PROGRAM_NAME, & + 'Program to test LAI K-matrix response.', & + 'CRTM Version: '//TRIM(Version) ) + + CALL GET_COMMAND_ARGUMENT(1, Sensor_Id) + IF ( LEN_TRIM(Sensor_Id) == 0 ) Sensor_Id = 'atms_n21' + Sensor_Id = ADJUSTL(Sensor_Id) + WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id) + + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo , & + File_Path=COEFFICIENTS_PATH ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error initializing CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo)) + + ALLOCATE( RTSolution( n_Channels, N_PROFILES ), & + RTSolution_Perturb( n_Channels, N_PROFILES ), & + RTSolution_K( n_Channels, N_PROFILES ), & + Atmosphere_K( n_Channels, N_PROFILES ), & + Surface_K( n_Channels, N_PROFILES ), & + STAT = Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating structure arrays' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN + Message = 'Error allocating CRTM Atmosphere structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atmosphere_K, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_K)) ) THEN + Message = 'Error allocating CRTM Atmosphere_K structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL Load_Atm_Data() + CALL Load_Sfc_Data() + + CALL CRTM_Geometry_SetValue( Geometry, & + Sensor_Zenith_Angle = ZENITH_ANGLE, & + Sensor_Scan_Angle = SCAN_ANGLE ) + + Sfc(:)%Land_Coverage = 1.0_fp + Sfc(:)%Water_Coverage = 0.0_fp + Sfc(:)%Snow_Coverage = 0.0_fp + Sfc(:)%Ice_Coverage = 0.0_fp + Sfc(:)%Canopy_Water_Content = 0.0_fp + Sfc(:)%Lai = 2.0_fp + + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Perturbation = 0.01_fp + Sfc(TEST_PROFILE)%Lai = Sfc(TEST_PROFILE)%Lai + Perturbation + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo , & + RTSolution_Perturb ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in perturbed CRTM Forward Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + Sfc(TEST_PROFILE)%Lai = Sfc(TEST_PROFILE)%Lai - Perturbation + + CALL CRTM_Atmosphere_Zero( Atmosphere_K ) + CALL CRTM_Surface_Zero( Surface_K ) + + RTSolution_K%Radiance = ZERO + RTSolution_K%Brightness_Temperature = ONE + + Error_Status = CRTM_K_Matrix( Atm , & + Sfc , & + RTSolution_K, & + Geometry , & + ChannelInfo , & + Atmosphere_K, & + Surface_K , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM K_Matrix Model' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + FD_Deriv = ( RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature - & + RTSolution(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature ) / Perturbation + K_Deriv = Surface_K(TEST_CHANNEL,TEST_PROFILE)%Lai + + IF ( ABS(K_Deriv) <= 0.0_fp ) THEN + WRITE(*,*) 'K-matrix LAI derivative:', K_Deriv + Message = 'K-matrix LAI derivative is zero' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + Ratio = FD_Deriv / K_Deriv + WRITE(*,*) 'Base BT:', RTSolution(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature + WRITE(*,*) 'Perturb BT:', RTSolution_Perturb(TEST_CHANNEL,TEST_PROFILE)%Brightness_Temperature + WRITE(*,*) 'FD dTB/dx:', FD_Deriv + WRITE(*,*) 'K dTB/dx:', K_Deriv + WRITE(*,*) 'Ratio: ', Ratio + + Error_Status = CRTM_Destroy( ChannelInfo ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error destroying CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Destroy(Atmosphere_K) + CALL CRTM_Atmosphere_Destroy(Atm) + + DEALLOCATE( RTSolution, RTSolution_Perturb, RTSolution_K, & + Atmosphere_K, Surface_K, STAT = Allocate_Status ) + + IF ( ABS(1.0_fp - Ratio) < TOLERANCE ) THEN + STOP 0 + END IF + + WRITE(*,*) 'FAIL abs(1 - Ratio)=', ABS(1.0_fp - Ratio), ' TOL=', TOLERANCE + STOP 1 + +CONTAINS + + INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Load_Sfc_Data.inc' + +END PROGRAM test_k_matrix_lai