Why is Surface_Reflectivity always zero in RTSolution

[chengdang]

Why is Surface_Reflectivity always zero in RTSolution

Note that for the data structure RTSolution, Surface Reflectivity is always zero (in most cases) regardless of channels.

This is because in function Assign_Common_Input

CRTMv3/src/RTSolution/Common_RTSolution.f90

Line 189 in 9901de1

FUNCTION Assign_Common_Input( &

1. The weight of angle "Cosine_Sensor_Zenith" is set as zero:

      SfcOptics%Index_Sat_Ang        = RTV%n_Angles 
      RTV%COS_Angle( RTV%n_Angles )  = GeometryInfo%Cosine_Sensor_Zenith
      RTV%COS_Weight( RTV%n_Angles ) = ZERO

2. Then the surface reflectivity is weighted by Gaussian weight, with the last angle being Cosine_Sensor_Zenith if Cosine_Sensor_Zenith is not equal to any Gaussian angles:

      DO j = 1, SfcOptics%n_Angles
        SfcOptics%Reflectivity(1:SfcOptics%n_Angles,1,j,1) = SfcOptics%Direct_Reflectivity(j,1)*SfcOptics%Weight(j)
      END DO

3. The computed reflectivity is later assigned to RT solution as the final output:

...
    RTSolution%Surface_Reflectivity = SfcOptics%Reflectivity( SfcOptics%Index_Sat_Ang, 1, SfcOptics%Index_Sat_Ang, 1 )
...

Because SfcOptics%Weight(SfcOptics%n_Angles) is zero, SfcOptics%Reflectivity( SfcOptics%Index_Sat_Ang, 1, SfcOptics%Index_Sat_Ang, 1 ) is always zero regardless of solvers, channels, and surface type

Exception

The only exception would be if Cosine_Sensor_Zenith is almost equal to COS_Angle from function Double_Gauss_Quadrature. In this case, RTV%COS_Weight( RTV%Index_Sat_Ang ) != ZERO, but in this case, the surface reflectivity value is wrong.

      ! Search for a matching angle
      DO i = 1, RTV%n_Streams
        IF( ABS( RTV%COS_Angle(i) - GeometryInfo%Cosine_Sensor_Zenith ) < ANGLE_THRESHOLD ) THEN
          SfcOptics%Index_Sat_Ang = i         ! Identify the sensor zenith angle
          RTV%n_Angles = RTV%n_Streams
          ! *****FLAW*****
          ! CAN WE REPLACE THIS CONSTRUCT WITH ANOTHER
          GO TO 100
          ! *****FLAW*****
        END IF
      END DO

[BenjaminTJohnson]

if this is for visible only, is the solar zenith angle appropriately set?

[chengdang]

The solar angle is set correctly for the TOA boundary condition, so it is not a problem

[chengdang]

[will revisit this with more tests]

Is this a problem for TOA results, for solar channels, with ADA solver?

In ADA module, in particular subroutine CRTM_AMOM_layer. The layer results are almost always weighted by RTV%COS_Weight:

         RTV%s_Layer_Refl(i,j,KL) = c * bb(i,j) * COS_Weight(j)
         RTV%s_Layer_Trans(i,j,KL) = c * ff(i,j) * COS_Weight(j)

Technically, we want the the Gaussian weight of this extra angle (Cosine_Sensor_Zenith) to be zero, because function Double_Gauss_Quadrature (in function Assign_Common_Input) has already populated a set of angles (number of streams) with total weights equal to 1, for ADA calculation.

 CALL Double_Gauss_Quadrature( RTV%n_Streams, RTV%COS_Angle, RTV%COS_Weight )

Which begs the question, what's the purpose of adding this extra angle when GeometryInfo%Cosine_Sensor_Zenith does not equal to any Gaussian angles? It looks like we are just adding a bunch of zeros to ADA calculation?

If we search Index_Sat_Ang, this index is also used for the calculation of direct reflectivity, which is computed directly from various LUTs and not weighted by COS_Weight. In subroutine CRTM_ADA:

...
    nZ = RTV%n_Angles * RTV%n_Stokes
...
    IF( RTV%Solar_Flag_true ) THEN
      RTV%s_Level_Rad_UP(1:nZ,n_Layers ) = RTV%s_Level_Rad_UP(1:nZ,n_Layers )+direct_reflectivity(1:nZ)* &
        RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_opt(n_Layers)/RTV%COS_SUN)
    END IF

Consider RTV%n_Stokes = 1 for now, nZ = RTV%n_Angles = SfcOptics%Index_Sat_Ang. At the surface, we have s_Level_Rad_UP for direct light at Cosine_Sensor_Zenith.

At TOA, eventually CRTM assign upward radiance based on Index_Sat_Ang as well (function Assign_Common_Output):

    n1 = (SfcOptics%Index_Sat_Ang-1)*RTV%n_Stokes + 1
    ! ADA and SOI specific assignments
    IF( RTV%Scattering_RT ) THEN

      ! Assign radiance output from ADA or SOI
      IF ( RTV%aircraft%rt ) THEN
        Radiance(:) = RTV%s_Level_Rad_UP(n1:n1-1+RTV%n_Stokes, RTV%aircraft%idx)
      ELSE
        Radiance(:) = RTV%s_Level_Rad_UP(n1:n1-1+RTV%n_Stokes, 0)
      END IF

Again for RTV%n_Stokes = 1, n1 = SfcOptics%Index_Sat_Ang, so the radiance at TOA would be the attenuated direct reflection from surface at angle Index_Sat_Ang, and average hemispheric diffuse radiance (??? need to check ADA solver for this), which is CRTM best approximation to space observations.

[chengdang]

So, to conclude, I think we could update the calculation of RTSolution%Surface_Reflectivity.

RTSolution%Surface_Reflectivity should be a weighted sum of diffuse reflectance and direct reflectance, where the weights are determined by the fractions of diffuse and direct radiance reaching the surface. However, given the design of CRTM and the ADA solver, downward diffuse radiance at the surface is not computed in most applications. Therefore, the simplest approach would be to assign the direct_reflectivity at Cosine_Sensor_Zenith to RTSolution%Surface_Reflectivity.

Alternatively, we could add a note to the CRTM User Guide indicating that users should use the direct reflectivity as the reference quantity, which can be obtained from the CRTM output as: direct_reflectivity = 1 - RTSolution%Surface_Emissivity