Merge pull request #4 from OCO/ticket/2957-l_rad-u_sign_fix

Delivery is done, so we can go ahead and merge this in. Fix l_rad coxmunk U-sign conventions
nasa · Jun 4, 2021 · 0f85a04 · 0f85a04
2 parents 62a2c3d + d3f8d29
commit 0f85a04
Show file tree

Hide file tree

Showing 3 changed files with 130 additions and 46 deletions.
diff --git a/lib/Implementation/l_rad_fortran/l_surface.F90 b/lib/Implementation/l_rad_fortran/l_surface.F90
@@ -1,6 +1,9 @@
 module l_surface_m
 implicit none
 
+!  10/26/20. Revision for BRDF consistency, R. Spurr
+!    -- All changes marked by "10/26/20. BRDF Upgrade"
+
 PUBLIC
 
 contains
@@ -9,6 +12,9 @@ module l_surface_m
 ! set nspars to 3, spars(1) to ws, spars(2) to ri and spars(3) to shadow
 ! factor (set to 1.d0 for including shadowing).
 
+! 10/26/20. BRDF Upgrade. 
+!   -- This NOTE is not correct any more !!!!
+
       subroutine L_R1_glint_exact &
        (nstokes,nspars,& !I
         xj,xi,phi,ws,ri,alb,sfac,& !I ! V. Natraj, 8/17/2010
@@ -314,14 +320,21 @@ subroutine gisscoxmunk_vfunction &
       AF21 = AF21*AF21
       AF22 = AF22*AF22
 
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1
+!       Formerly, R1(3) = + ( CTTPT+CTPPP ) * DCOEFF , Now R1(3) = - ( CTTPT+CTPPP ) * DCOEFF
+
       R1(1)=(AF11+AF12+AF21+AF22)*AF
       R1(2)=(AF11-AF22+AF12-AF21)*AF
       IF (NSTOKES .EQ. 3) THEN
         CTTPT = CF11*CF21
         CTPPP = CF12*CF22
+
+        R1(NSTOKES) = - (CTTPT+CTPPP)*DCOEFF 
+
 !        R1(NSTOKES) = (-CTTPT-CTPPP)*DCOEFF
-        R1(NSTOKES) = (CTTPT+CTPPP)*DCOEFF ! Change sign of U since this code uses the opposite sign convention as we do 
-                                           ! in 2OS and LIDORT
+!        R1(NSTOKES) = (CTTPT+CTPPP)*DCOEFF ! Change sign of U since this code uses the opposite sign convention as we do 
+!                                           ! in 2OS and LIDORT
       ENDIF
 
 !  No Shadow code if not flagged
@@ -590,24 +603,34 @@ subroutine gisscoxmunk_vfunction_plus &
       L_AF21 = 2.0d0 * CF21 * L_CF21
       L_AF22 = 2.0d0 * CF22 * L_CF22
 
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1
+!       Formerly, R1(3) = + ( CTTPT+CTPPP ) * DCOEFF , Now R1(3) = - ( CTTPT+CTPPP ) * DCOEFF
+
       R1(1)=(AF11+AF12+AF21+AF22)*AF   ! CN2
       R1(2)=(AF11-AF22+AF12-AF21)*AF   ! CN2
       IF (NSTOKES .EQ. 3) THEN
         CTTPT = CF11*CF21             ! CN2
         CTPPP = CF12*CF22             ! CN2
-!        R1(NSTOKES) = (-CTTPT-CTPPP)*DCOEFF   ! CN2
-        R1(NSTOKES) = (CTTPT+CTPPP)*DCOEFF ! Change sign of U since this code uses the opposite sign convention as we do 
-                                           ! in 2OS and LIDORT
+
+        R1(NSTOKES) = - (CTTPT+CTPPP)*DCOEFF 
+!        R1(NSTOKES) = (-CTTPT-CTPPP)*DCOEFF
+
+
+
       ENDIF
 
+! 10/26/20. BRDF Upgrade. TWO BUGS HERE, Connected to the Stokes-U component
+!   1. Former Code for LS_R1(3,2) was signed correctly, even though R1(3) was incorrectly signed
+!   2. Ls_R1(3,2) should have DCOEFF multiplier, not AF. Was a factor of 2 too small
+
       !derivs wrt ri
       Ls_R1(1,2) = (L_AF11+L_AF12+L_AF21+L_AF22)*AF
       Ls_R1(2,2) = (L_AF11-L_AF22+L_AF12-L_AF21)*AF
       IF (NSTOKES .EQ. 3) &
-!      Ls_R1(3,2) = -(L_CF11*CF21+CF11*L_CF21+&
-!                     L_CF12*CF22+CF12*L_CF22)*AF
-      Ls_R1(3,2) = (L_CF11*CF21+CF11*L_CF21+&   ! Change sign of U since this code uses the opposite sign convention as we 
-                    L_CF12*CF22+CF12*L_CF22)*AF ! do in 2OS and LIDORT
+
+!         Ls_R1(NSTOKES,2) = -(L_CF11*CF21+CF11*L_CF21+L_CF12*CF22+CF12*L_CF22)*AF
+         Ls_R1(NSTOKES,2) = -(L_CF11*CF21+CF11*L_CF21+L_CF12*CF22+CF12*L_CF22)*DCOEFF
 
 !  Derivative before shadow effect
 
@@ -879,6 +902,11 @@ subroutine landbrdf_vfunction &
       AF22 = AF22*AF22
 
       FACTOR = 0.5d0/DMOD
+
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1
+!       Formerly, R1(3) = + ( CTTPT+CTPPP ) * FACTOR , Now R1(3) = - ( CTTPT+CTPPP ) * FACTOR
+
       R1(1) = (AF11+AF12+AF21+AF22) * FACTOR
       R1(2) = (AF11-AF12+AF21-AF22) * FACTOR
 
@@ -889,8 +917,11 @@ subroutine landbrdf_vfunction &
         CTTPT=CF11*CF21
         CTPPP=CF12*CF22
         FACTOR = 1.d0/DMOD
+
+        R1(3) = - ( CTTPT+CTPPP ) * FACTOR
+
 !        R1(3)  = (-CTTPT-CTPPP) * FACTOR
-        R1(3)  = (CTTPT+CTPPP) * FACTOR ! Change sign for U just as in Cox-Munk
+!        R1(3)  = (CTTPT+CTPPP) * FACTOR ! Change sign for U just as in Cox-Munk
       ENDIF
 
 !  Set the H-function
@@ -948,9 +979,14 @@ subroutine landbrdf_vfunction &
 
 !  This is just the Rahman Kernel.........different name !!
 
+! 10/26/20. BRDF Upgrade. 
+!   incident and reflected zenith angles swapped here, so swap them in the call
+!   Does not matter for this kernel, as scalar only.
+!   XJ, SXJ, XI, SXI  ==> XI, SXI, XJ, SXJ 
+
       CALL rahman_function_2os &
             ( 3, PARS(2:4),& !I
-              XJ, SXJ, XI, SXI,& !I
+              XI, SXI, XJ, SXJ,& !I
               CKPHI_REF, SKPHI_REF,& !I
               RAHMAN_KERNEL )  !O
 
@@ -968,6 +1004,10 @@ subroutine rahman_function_2os &
               XJ, SXJ, XI, SXI, CPHI, SKPHI,& !I
               RAHMAN_KERNEL ) !O
 
+! 10/26/20. BRDF Upgrade. 
+!   incident and reflected zenith angles swapped before this kernel is called
+!   so do not need to swap them here. Does not matter for this kernel, as scalar only.
+
 !  Revision. 24 October 2007.
 !  --------------------------
 
@@ -1252,6 +1292,10 @@ subroutine landbrdf_vfunction_plus &
       AF21 = AF21*AF21
       AF22 = AF22*AF22
 
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1
+!       Formerly, R1(3) = + ( CTTPT+CTPPP ) * FACTOR , Now R1(3) = - ( CTTPT+CTPPP ) * FACTOR
+
       FACTOR = 0.5d0/DMOD
       R1(1) = (AF11+AF12+AF21+AF22) * FACTOR
       R1(2) = (AF11-AF12+AF21-AF22) * FACTOR
@@ -1263,8 +1307,11 @@ subroutine landbrdf_vfunction_plus &
         CTTPT=CF11*CF21
         CTPPP=CF12*CF22
         FACTOR = 1.d0/DMOD
+        R1(3)  =  - ( CTTPT+CTPPP) * FACTOR
+
 !        R1(3)  = (-CTTPT-CTPPP) * FACTOR
-        R1(3)  = (CTTPT+CTPPP) * FACTOR ! Change sign for U just as in Cox-Munk
+!       R1(3)  = (CTTPT+CTPPP) * FACTOR ! Change sign for U just as in Cox-Munk
+
       ENDIF
 
 !  Set the H-function
@@ -1326,9 +1373,14 @@ subroutine landbrdf_vfunction_plus &
 
 !  This is just the Rahman Kernel.........different name !!
 
+! 10/26/20. BRDF Upgrade. 
+!   incident and reflected zenith angles swapped here, so swap them in the call
+!   Does not matter for this kernel, as scalar only.
+!   XJ, SXJ, XI, SXI  ==> XI, SXI, XJ, SXJ 
+
       CALL rahman_function_2os_plus &
             ( 3, PARS(2:4), DO_DERIV_PARS,& !I
-              XJ, SXJ, XI, SXI,& !I
+              XI, SXI, XJ, SXJ,& !I
               CKPHI_REF, SKPHI_REF,& !I
               RAHMAN_KERNEL, RAHMAN_DERIVATIVES ) !O
 
@@ -1356,6 +1408,10 @@ subroutine rahman_function_2os_plus &
               XJ, SXJ, XI, SXI, CPHI, SKPHI,& !I
               RAHMAN_KERNEL, RAHMAN_DERIVATIVES ) !O
 
+! 10/26/20. BRDF Upgrade. 
+!   incident and reflected zenith angles swapped before this kernel is called
+!   so do not need to swap them here. Does not matter for this kernel, as scalar only.
+
       IMPLICIT NONE
 
       DOUBLE PRECISION SMALL

diff --git a/lib/Implementation/l_rad_fortran/l_surface_fourier.F90 b/lib/Implementation/l_rad_fortran/l_surface_fourier.F90
@@ -1,12 +1,15 @@
 module l_surface_fourier_m
 
+!  10/26/20. Revision for BRDF consistency, R. Spurr
+!    -- All changes marked by "10/26/20. BRDF Upgrade"
+
 USE l_surface_m
 implicit none
 PUBLIC
 
 contains
 
-! NOT : For Lambertian, set nspars to 1 and spars(1) to Asurf. For glint, 
+! NOTE : For Lambertian, set nspars to 1 and spars(1) to Asurf. For glint, 
 ! set nspars to 3, spars(1) to ws, spars(2) to ri and spars(3) to shadow &
 ! fac or (set to 1.d0 for including shadowing).
 
@@ -203,12 +206,19 @@ SUBROUTINE GISSCOXMUNK_FOURIER &
 !      R1M(3,1)= (-CTTPT-CTPPP)*DCOEFF
 !      R1M(3,2)= (-CTTPT+CTPPP)*DCOEFF
 
-! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1M
 
-      R1M(1,3)= (CTTTP+CPTPP)*DCOEFF
-      R1M(2,3)= (CTTTP-CPTPP)*DCOEFF
-      R1M(3,1)= (CTTPT+CTPPP)*DCOEFF
-      R1M(3,2)= (CTTPT-CTPPP)*DCOEFF
+      R1M(1,3) =  - ( CTTTP+CPTPP ) * DCOEFF   ! Sine
+      R1M(2,3) =  - ( CTTTP-CPTPP ) * DCOEFF   ! Sine
+      R1M(3,1) =  - ( CTTPT+CTPPP ) * DCOEFF   ! Sine
+      R1M(3,2) =  - ( CTTPT-CTPPP ) * DCOEFF   ! Sine
+
+! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+!      R1M(1,3)= (CTTTP+CPTPP)*DCOEFF
+!      R1M(2,3)= (CTTTP-CPTPP)*DCOEFF
+!      R1M(3,1)= (CTTPT+CTPPP)*DCOEFF
+!      R1M(3,2)= (CTTPT-CTPPP)*DCOEFF
 
       R1M(3,3)= (CTTPP+CTPPT)*DCOEFF
       R1M(4,4)= (CTTPP-CTPPT)*DCOEFF
@@ -464,18 +474,30 @@ SUBROUTINE GISSCOXMUNK_FOURIER_PLUS &
 !      R1M(3,1)= (-CTTPT-CTPPP)*DCOEFF
 !      R1M(3,2)= (-CTTPT+CTPPP)*DCOEFF
 
-! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1M
+
+      R1M(1,3) =  - ( CTTTP+CPTPP ) * DCOEFF   ! Sine
+      R1M(2,3) =  - ( CTTTP-CPTPP ) * DCOEFF   ! Sine
+      R1M(3,1) =  - ( CTTPT+CTPPP ) * DCOEFF   ! Sine
+      R1M(3,2) =  - ( CTTPT-CTPPP ) * DCOEFF   ! Sine
 
-      R1M(1,3)= (CTTTP+CPTPP)*DCOEFF
-      R1M(2,3)= (CTTTP-CPTPP)*DCOEFF
-      R1M(3,1)= (CTTPT+CTPPP)*DCOEFF
-      R1M(3,2)= (CTTPT-CTPPP)*DCOEFF
+! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+!      R1M(1,3)= (CTTTP+CPTPP)*DCOEFF
+!      R1M(2,3)= (CTTTP-CPTPP)*DCOEFF
+!      R1M(3,1)= (CTTPT+CTPPP)*DCOEFF
+!      R1M(3,2)= (CTTPT-CTPPP)*DCOEFF
 
       R1M(3,3)= (CTTPP+CTPPT)*DCOEFF
       R1M(4,4)= (CTTPP-CTPPT)*DCOEFF
 
 !  Derivatives wrt ri, V. Natraj, 8/17/2010
 
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1M
+!   Original coding was correct, corresponds to sign-change introduced above
+!     - this was a bug in the code.......!!!
+
       LS_R1M(1,3,2)= (-L_CTTTP-L_CPTPP)*DCOEFF
       LS_R1M(2,3,2)= (-L_CTTTP+L_CPTPP)*DCOEFF
       LS_R1M(3,1,2)= (-L_CTTPT-L_CTPPP)*DCOEFF
@@ -1245,17 +1267,20 @@ SUBROUTINE LANDBRDF_FOURIER &
       CPTPP=CF21*CF22
 
       FACTOR = 1.d0/DMOD
-!      R1M(1,3)= (-CTTTP-CPTPP)*FACTOR
-!      R1M(2,3)= (-CTTTP+CPTPP)*FACTOR
-!      R1M(3,1)= (-CTTPT-CTPPP)*FACTOR
-!      R1M(3,2)= (-CTTPT+CTPPP)*FACTOR
 
-! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1M
 
-      R1M(1,3)= (CTTTP+CPTPP)*FACTOR
-      R1M(2,3)= (CTTTP-CPTPP)*FACTOR
-      R1M(3,1)= (CTTPT+CTPPP)*FACTOR
-      R1M(3,2)= (CTTPT-CTPPP)*FACTOR
+      R1M(1,3) =  - ( CTTTP+CPTPP ) * FACTOR   ! Sine
+      R1M(2,3) =  - ( CTTTP-CPTPP ) * FACTOR   ! Sine
+      R1M(3,1) =  - ( CTTPT+CTPPP ) * FACTOR   ! Sine
+      R1M(3,2) =  - ( CTTPT-CTPPP ) * FACTOR   ! Sine
+
+! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+!      R1M(1,3)= (CTTTP+CPTPP)*FACTOR
+!      R1M(2,3)= (CTTTP-CPTPP)*FACTOR
+!      R1M(3,1)= (CTTPT+CTPPP)*FACTOR
+!      R1M(3,2)= (CTTPT-CTPPP)*FACTOR
 
       R1M(3,3)= (CTTPP+CTPPT)*FACTOR
       R1M(4,4)= (CTTPP-CTPPT)*FACTOR
@@ -1510,17 +1535,20 @@ SUBROUTINE LANDBRDF_FOURIER_PLUS &
       CPTPP=CF21*CF22
 
       FACTOR = 1.d0/DMOD
-!      R1M(1,3)= (-CTTTP-CPTPP)*FACTOR
-!      R1M(2,3)= (-CTTTP+CPTPP)*FACTOR
-!      R1M(3,1)= (-CTTPT-CTPPP)*FACTOR
-!      R1M(3,2)= (-CTTPT+CTPPP)*FACTOR
 
-! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+! 10/26/20. BRDF Upgrade. 
+!   Sign Switch for the stokes-U contributions to R1M
+
+      R1M(1,3) =  - ( CTTTP+CPTPP ) * FACTOR   ! Sine
+      R1M(2,3) =  - ( CTTTP-CPTPP ) * FACTOR   ! Sine
+      R1M(3,1) =  - ( CTTPT+CTPPP ) * FACTOR   ! Sine
+      R1M(3,2) =  - ( CTTPT-CTPPP ) * FACTOR   ! Sine
 
-      R1M(1,3)= (CTTTP+CPTPP)*FACTOR
-      R1M(2,3)= (CTTTP-CPTPP)*FACTOR
-      R1M(3,1)= (CTTPT+CTPPP)*FACTOR
-      R1M(3,2)= (CTTPT-CTPPP)*FACTOR
+! Change the sign of the sine terms to account for the opposite sign convention compared to 2OS and LIDORT
+!      R1M(1,3)= (CTTTP+CPTPP)*FACTOR
+!      R1M(2,3)= (CTTTP-CPTPP)*FACTOR
+!      R1M(3,1)= (CTTPT+CTPPP)*FACTOR
+!      R1M(3,2)= (CTTPT-CTPPP)*FACTOR
 
       R1M(3,3)= (CTTPP+CTPPT)*FACTOR
       R1M(4,4)= (CTTPP-CTPPT)*FACTOR

diff --git a/lib/Implementation/l_rad_rt_test.cc b/lib/Implementation/l_rad_rt_test.cc
@@ -306,8 +306,8 @@ BOOST_AUTO_TEST_CASE(reflectance_first_order)
   wn = 12929.94, 12930.30;
   Array<double, 2> stokes_expect(2, 3);
   stokes_expect = 
-    0.0079613385061960122852, -0.001546138923946610913, -0.00013167264462035855243,
-    0.0079615098335930490486, -0.0015462247310868153412, -0.00013171900922801001578;
+    0.0079613385061960122852, -0.0015461389239466109130, -0.00055660191445895516611,
+    0.0079615098335930490486, -0.0015462247310868153412, -0.00055663280457999967689;
   BOOST_CHECK_MATRIX_CLOSE_TOL(rt_first_order->stokes(wn, 0), 
                                stokes_expect, check_tol);
   BOOST_CHECK_MATRIX_CLOSE_TOL
@@ -328,8 +328,8 @@ BOOST_AUTO_TEST_CASE(reflectance_second_order)
   wn = 12929.94, 12930.30;
   Array<double, 2> stokes_expect(2, 3);
   stokes_expect = 
-      0.0079672428053023970629, -0.0023774729200916268519, -0.00028391678080464213037,
-      0.0079674137933562587388, -0.002377560330641938658, -0.00028396880453497954679;
+    0.0079681282901736590757, -0.0024250036984472403759, -0.00074540308734331189606,
+    0.0079682994187089527943, -0.0024250909594306377728, -0.00074543862544189246414;
   BOOST_CHECK_MATRIX_CLOSE_TOL(rt_second_order->stokes(wn, 0), 
                                stokes_expect, check_tol);
   BOOST_CHECK_MATRIX_CLOSE_TOL