Initial implementation of the scaled coxmunk BRDF.

Author: gmcgarragh

Makefile.in

@@ -2297,6 +2297,59 @@ function ConfigCommon.lambertian_retrieval:initial_guess()
    return ig 
 end
 
+------------------------------------------------------------
+--- Lambertian ground state vector component and initial guess
+------------------------------------------------------------
+
+ConfigCommon.brdf_scale_retrieval = CreatorMultiSpec:new {}
+
+function ConfigCommon.brdf_scale_retrieval:create()
+   local num_coeff = self:apriori_v(0):rows()
+   local num_spec = self.config.number_pixel:rows()
+
+   local ap = Blitz_double_array_2d(num_spec, num_coeff)
+   local flag = Blitz_bool_array_2d(num_spec, num_coeff)
+
+   for i = 1, num_spec do
+       ap:set(i-1, Range.all(), self:apriori_v(i - 1))
+       flag:set(i-1, Range.all(), self:retrieval_flag(i))
+   end
+
+   local lambertian = GroundBrdfWeight(ap, flag, 
+                                       self.config.common.band_reference,
+                                       self.config.common.desc_band_name,
+                                       self.scaled_brdf_name)
+   return lambertian
+end
+
+function ConfigCommon.brdf_scale_retrieval:initial_guess()
+   local ig = CompositeInitialGuess()
+   for i=1,self.config.number_pixel:rows() do
+       local flag = self:retrieval_flag(i)
+
+       local ap = self:apriori_v(i - 1) 
+       local apcov_in = self:covariance_v(i - 1)
+
+       -- Copy input covariance to modified one which repeats the last covariance value if
+       -- the number of degrees is defined
+       local apcov_mod = Blitz_double_array_2d(ap:rows(), ap:rows())
+       apcov_mod:set(Range.all(), Range.all(), 0.0)
+       src_idx = 0
+       for dst_idx = 0, ap:rows()-1 do
+          apcov_mod:set(dst_idx, dst_idx, apcov_in(src_idx, src_idx))
+          if dst_idx < apcov_in:rows()-1 then
+             src_idx = src_idx + 1
+          end
+       end
+
+       local band_ig = InitialGuessValue()
+       band_ig:apriori_subset(flag, ap)
+       band_ig:apriori_covariance_subset(flag, apcov_mod)
+       ig:add_builder(band_ig)
+   end
+   return ig 
+end
+
 ------------------------------------------------------------
 --- Creates a GroundOco only using a Lambertian retrieval
 ------------------------------------------------------------
@@ -2345,7 +2398,7 @@ function ConfigCommon.ground_coxmunk:register_output(ro)
 end
 
 ------------------------------------------------------------
---- Creates a GroundOco using both Lambertioan and Coxmunk 
+--- Creates a GroundOco using both Lambertian and Coxmunk 
 --- retrieval
 ------------------------------------------------------------
 
@@ -2364,6 +2417,26 @@ function ConfigCommon.ground_coxmunk_plus_lamb:register_output(ro)
    ro:push_back(GroundCoxmunkPlusLambertianOutput.create(self.config.coxmunk, self.config.coxmunk_lambertian, self.config.common.hdf_band_name))
 end
 
+------------------------------------------------------------
+--- Creates a GroundOco using both Lambertian and Coxmunk 
+--- retrieval
+------------------------------------------------------------
+
+ConfigCommon.ground_coxmunk_scaled = CompositeCreator:new {}
+
+function ConfigCommon.ground_coxmunk_scaled:sub_object_key()
+   -- Ordering is important here or else the statevector is setup incorrectly
+   return { "coxmunk", "coxmunk_scaled" }
+end
+
+function ConfigCommon.ground_coxmunk_scaled:create_parent_object(sub_object)
+   return GroundCoxmunkScaled.create(self.config.coxmunk, self.config.coxmunk_scaled)
+end
+
+function ConfigCommon.ground_coxmunk_scaled:register_output(ro)
+   ro:push_back(GroundCoxmunkScaledOutput.create(self.config.l1b, self.config.coxmunk, self.config.coxmunk_scaled, self.config.common.hdf_band_name))
+end
+
 ------------------------------------------------------------
 --- Modify the GroundBrdf a_priori, scaling by the 
 --- kernel value

input/common/config/config_common.lua

@@ -471,6 +471,15 @@ GeocarbBaseConfig = GeocarbConfig:new {
                creator = ConfigCommon.lambertian_retrieval,
             },
 
+            -- Lambertian component of coxmunk + lambertian
+            coxmunk_scaled = {
+               apriori = ConfigCommon.hdf_apriori_i("Ground/Coxmunk_Scaled"),
+               covariance = ConfigCommon.hdf_covariance_i("Ground/Coxmunk_Scaled"),
+               retrieve_bands = { true, true, true, true },
+               scaled_brdf_name = "CoxMunk",
+               creator = ConfigCommon.brdf_scale_retrieval,
+            },
+
             -- Brdf vegetative kernel with Rahman retrieved parameters
             brdf_veg = {
                apriori = ConfigCommon.brdf_veg_apriori("Ground/Brdf"),

input/geocarb/config/geocarb_base_config.lua

@@ -357,6 +357,28 @@ function GeocarbConfig.ground_from_ground_type:get_creator()
    end
 end
 
+GeocarbConfig.ground_from_ground_type_scaled = DispatchCreator:new()
+
+function GeocarbConfig.ground_from_ground_type_scaled:get_creator()
+   local ground_type = self.config:ground_type_name()
+
+   if (ground_type == "lambertian") then
+      return ConfigCommon.ground_lambertian
+   elseif (ground_type == "brdf_soil") then
+      return ConfigCommon.ground_brdf_soil
+   elseif (ground_type == "brdf_veg") then
+      return ConfigCommon.ground_brdf_veg
+   elseif (ground_type == "coxmunk") then
+      if(self.config.using_radiance_scaling ~= nil) then
+         return ConfigCommon.ground_coxmunk
+      else
+         return ConfigCommon.ground_coxmunk_scaled
+      end
+   else
+      error("Invalid ground type value: " .. ground_type)
+   end
+end
+
 ------------------------------------------------------------
 --- Create lambertian ground initial guess from radiance
 --- but the other parts from the static HDF file

input/geocarb/config/geocarb_config.lua

@@ -466,6 +466,15 @@ OcoBaseConfig = OcoConfig:new {
                creator = ConfigCommon.lambertian_retrieval,
             },
 
+            -- Lambertian component of coxmunk + lambertian
+            coxmunk_scaled = {
+               apriori = ConfigCommon.hdf_apriori_i("Ground/Coxmunk_Scaled"),
+               covariance = ConfigCommon.hdf_covariance_i("Ground/Coxmunk_Scaled"),
+               retrieve_bands = { true, true, true },
+               scaled_brdf_name = "CoxMunk",
+               creator = ConfigCommon.brdf_scale_retrieval,
+            },
+
             -- Brdf vegetative kernel with Rahman retrieved parameters
             brdf_veg = {
                apriori = ConfigCommon.brdf_veg_apriori("Ground/BrdfQuadratic"),
@@ -526,7 +535,8 @@ OcoBaseConfig = OcoConfig:new {
             creator = ConfigCommon.absorber_creator,
             gases = {"CO2", "H2O", "O2"},
             CO2 = {
-               apriori = ConfigCommon.co2_profile_file_apriori,
+               apriori = ConfigCommon.reference_co2_apriori_met_apriori,
+--             apriori = ConfigCommon.co2_profile_file_apriori,
                covariance = ConfigCommon.hdf_covariance("Gas/CO2"),
                absco = "v5.0_alpha/raw/co2_devi2015_wco2scale=nist_sco2scale=unity.hdf",
                table_scale = {1.0, 1.0, 1.004},

input/geocarb/input/l2_geocarb_static_input.h5

@@ -374,6 +374,28 @@ function OcoConfig.ground_from_ground_type:get_creator()
    end
 end
 
+OcoConfig.ground_from_ground_type_scaled = DispatchCreator:new()
+
+function OcoConfig.ground_from_ground_type_scaled:get_creator()
+   local ground_type = self.config:ground_type_name()
+
+   if (ground_type == "lambertian") then
+      return ConfigCommon.ground_lambertian
+   elseif (ground_type == "brdf_soil") then
+      return ConfigCommon.ground_brdf_soil
+   elseif (ground_type == "brdf_veg") then
+      return ConfigCommon.ground_brdf_veg
+   elseif (ground_type == "coxmunk") then
+      if(self.config.using_radiance_scaling ~= nil) then
+         return ConfigCommon.ground_coxmunk
+      else
+         return ConfigCommon.ground_coxmunk_scaled
+      end
+   else
+      error("Invalid ground type value: " .. ground_type)
+   end
+end
+
 ------------------------------------------------------------
 --- Create lambertian ground initial guess from radiance
 --- but the other parts from the static HDF file

input/oco/config/oco_base_config.lua

@@ -4,7 +4,7 @@ module brdf_functions_m
 
     USE LIDORT_pars, only : fpk, DEG_TO_RAD, MAXBEAMS, MAXSTREAMS_BRDF, MAX_BRDF_KERNELS, &
                             MAXSTHALF_BRDF, MAX_BRDF_PARAMETERS, ZERO, ONE, TWO, PIE, &
-                            MAX_USER_STREAMS, MAX_USER_RELAZMS, &
+                            MAX_USER_STREAMS, MAX_USER_RELAZMS, COXMUNK_IDX, &
                             BREONVEG_IDX, BREONSOIL_IDX, RAHMAN_IDX
 
     implicit none
@@ -28,7 +28,18 @@ real(kind=c_double) function exact_brdf_value_soil_f(params, sza, vza, azm) bind
 
         exact_brdf_value_soil_f = exact_brdf_value_f(BREONSOIL_IDX, params, sza, vza, azm)
     end function exact_brdf_value_soil_f
- 
+
+    real(kind=c_double) function exact_brdf_value_coxmunk(params, sza, vza, azm) bind(c)
+        real(kind=c_double), intent(in) :: params(5)
+        real(kind=c_double), intent(in) :: sza
+        real(kind=c_double), intent(in) :: vza
+        real(kind=c_double), intent(in) :: azm
+
+        exact_brdf_value_coxmunk = exact_brdf_value_cm_f(params, sza, vza, azm)
+    end function exact_brdf_value_coxmunk
+
+ !------------------------------------------------------------------------------
+
     real(kind=c_double) function exact_brdf_value_f(breon_type, params, sza, vza, azm) bind(c)
 
       USE brdf_sup_routines_m, only : BRDF_FUNCTION
@@ -200,6 +211,168 @@ real(kind=c_double) function exact_brdf_value_f(breon_type, params, sza, vza, az
 
     END FUNCTION exact_brdf_value_f
 
+ !------------------------------------------------------------------------------
+
+    real(kind=c_double) function exact_brdf_value_cm_f(params, sza, vza, azm) bind(c)
+
+      USE brdf_sup_routines_m, only : BRDF_FUNCTION
+
+!  Parameters for lrad
+
+      INTEGER, PARAMETER :: NSTOKES = 3
+      INTEGER, PARAMETER :: NSPARS = 5
+
+      real(kind=c_double), intent(in) :: params(NSPARS)
+      real(kind=c_double), intent(in) :: sza
+      real(kind=c_double), intent(in) :: vza
+      real(kind=c_double), intent(in) :: azm
+
+!  Number and index-list of bidirectional functions
+
+      INTEGER   :: N_BRDF_KERNELS
+      INTEGER   :: WHICH_BRDF ( MAX_BRDF_KERNELS )
+
+!  Parameters required for Kernel families
+
+      INTEGER   :: N_BRDF_PARAMETERS ( MAX_BRDF_KERNELS )
+      REAL(fpk) :: BRDF_PARAMETERS   ( MAX_BRDF_KERNELS, MAX_BRDF_PARAMETERS )
+
+!  Input kernel amplitude factors
+
+      REAL(fpk) :: BRDF_FACTORS ( MAX_BRDF_KERNELS )
+
+!  Cos and since of SZAs
+
+      REAL(fpk)   :: SZASURCOS(MAXBEAMS)
+      REAL(fpk)   :: SZASURSIN(MAXBEAMS)
+
+!  Local angle control
+
+      INTEGER    :: NBEAMS
+      INTEGER    :: N_USER_STREAMS
+      INTEGER    :: N_USER_RELAZMS
+
+!  Angles
+
+      REAL(fpk) :: BEAM_SZAS         (MAXBEAMS)
+      REAL(fpk) :: USER_ANGLES_INPUT (MAX_USER_STREAMS)
+      REAL(fpk) :: USER_RELAZMS      (MAX_USER_RELAZMS)
+      REAL(fpk) :: USER_STREAMS      (MAX_USER_STREAMS)
+      REAL(fpk) :: USER_SINES        (MAX_USER_STREAMS)
+      REAL(fpk) :: PHIANG            (MAX_USER_RELAZMS)
+      REAL(fpk) :: COSPHI            (MAX_USER_RELAZMS)
+      REAL(fpk) :: SINPHI            (MAX_USER_RELAZMS)
+
+!  Exact DB values from LIDORT
+
+      REAL(fpk)  :: EXACTDB_BRDFUNC ( MAX_USER_STREAMS, MAX_USER_RELAZMS, MAXBEAMS )
+
+!  Lrad variables
+
+      INTEGER   :: HFUNCTION_INDEX
+      real(fpk) :: SPARS(NSPARS)
+      real(fpk) :: XI, SXI, XJ, SXJ
+      real(fpk) :: CKPHI_REF, SKPHI_REF
+      real(fpk) :: R1(NSTOKES)
+
+!  Calculated BRDF
+
+      REAL(fpk)  :: EXACT_CALC_BRDF
+
+!  Help
+
+      INTEGER    :: IB, K, UI, IA
+      REAL(fpk)  :: MUX
+
+!  Set values for the basic LIDORT variables
+
+      N_BRDF_KERNELS = 1
+      WHICH_BRDF(1) = COXMUNK_IDX
+      N_BRDF_PARAMETERS(1) = 2
+      BRDF_PARAMETERS = ZERO
+      ! Breon refactive index squared, same as value hardcoded inside lrad
+      BRDF_PARAMETERS(1,1) = params(1)
+      BRDF_PARAMETERS(1,2) = params(2)
+      BRDF_FACTORS(1) = 1.
+      NBEAMS = 1
+      BEAM_SZAS(1) = sza
+      N_USER_STREAMS = 1
+      USER_ANGLES_INPUT(1) = vza
+      N_USER_RELAZMS = 1
+      USER_RELAZMS(1) = azm
+
+!  Set values for the basic l_rad variables
+
+      spars(1) = 1.
+      spars(2) = BRDF_PARAMETERS(1,1)
+      spars(3) = BRDF_PARAMETERS(1,2)
+      spars(4) = 0.
+      spars(5) = 0.
+
+!  Cosine and since of SZA
+
+      DO IB = 1, NBEAMS
+         MUX =  COS(BEAM_SZAS(IB)*DEG_TO_RAD)
+         SZASURCOS(IB) = MUX
+         SZASURSIN(IB) = SQRT(1.0_fpk-MUX*MUX)
+      ENDDO
+
+      XI = SZASURCOS(1)
+      SXI = SZASURSIN(1)
+
+!  Cosine and since of VZA
+
+      DO UI = 1, N_USER_STREAMS
+         MUX =  COS(USER_ANGLES_INPUT(UI)*DEG_TO_RAD)
+         USER_STREAMS(UI) = MUX
+         USER_SINES(UI) = SQRT(1.0_fpk-MUX*MUX)
+      ENDDO
+
+      XJ = USER_STREAMS(1)
+      SXJ = USER_SINES(1)
+
+!  Cosine and since of AZM
+
+      DO IA = 1, N_USER_RELAZMS
+         PHIANG(IA) = USER_RELAZMS(IA)*DEG_TO_RAD
+         COSPHI(IA) = COS(PHIANG(IA))
+         SINPHI(IA) = SIN(PHIANG(IA))
+      ENDDO
+
+      CKPHI_REF = COSPHI(1)
+      SKPHI_REF = SINPHI(1)
+
+!  Exact BRDF (LIDORT)
+!  -------------------
+
+      EXACT_CALC_BRDF = ZERO
+
+!  Kernel loop
+
+      DO K = 1, N_BRDF_KERNELS
+
+         DO IA = 1, N_USER_RELAZMS
+            DO IB = 1, NBEAMS
+               DO UI = 1, N_USER_STREAMS
+                  CALL BRDF_FUNCTION &
+                     ( MAX_BRDF_PARAMETERS, WHICH_BRDF(K), N_BRDF_PARAMETERS(K), BRDF_PARAMETERS(K,:), & ! Inputs
+                       SZASURCOS(IB), SZASURSIN(IB), USER_STREAMS(UI),         & ! Inputs
+                       USER_SINES(UI), PHIANG(IA), COSPHI(IA), SINPHI(IA),     & ! Inputs
+                       EXACTDB_BRDFUNC(UI,IA,IB) )                               ! Output
+               ENDDO
+            ENDDO
+         ENDDO
+
+         EXACT_CALC_BRDF = EXACT_CALC_BRDF + BRDF_FACTORS(K) * EXACTDB_BRDFUNC(1,1,1)
+
+      ENDDO
+
+      exact_brdf_value_cm_f = EXACT_CALC_BRDF
+
+    END FUNCTION exact_brdf_value_cm_f
+
+ !------------------------------------------------------------------------------
+
     real(kind=c_double) function black_sky_albedo_veg_f(params, sza) bind(c)
         real(kind=c_double), intent(in) :: params(5)
         real(kind=c_double), intent(in) :: sza