Check with Alex

[BaptisteVandecrux]

pySICE/fortran/version 6.1 2022/s.f

Line 569 in 53488af

if (powe.le.0.9) powe=0.

Seems quite drastic

[BaptisteVandecrux]

pySICE/fortran/version 6.1 2022/s.f

Lines 743 to 748 in a6b016d

	if (ctoa(1).lt.THV0)bbasw1=bbasw1*factor
	if (ctoa(1).lt.THV0)bbavis1=bbavis1*factor
	if (ctoa(1).lt.THV0)bbanir1=bbanir1*factor
	if (ctoa(1).lt.THV0)bbasw2=bbasw2*factor
	if (ctoa(1).lt.THV0)bbavis2=bbavis2*factor
	if (ctoa(1).lt.THV0)bbanir2=bbanir2*factor

here ctoa is used to find the fractional pixels, but ctoa was corrected (=brightened up) using that same snow fraction earlier in the code. Additionally, thv0 is a threshold in surface spherical albedo, why is it not compared to the TOA reflectance?

[BaptisteVandecrux]

For the reference here is some info from Alex about the difference between solved spectral albedo:

pySICE/sice.py

Lines 366 to 438 in 17ea753

	def snow_albedo_solved(OLCI_scene, angles, aerosol, atmosphere, snow):
	# solving iteratively the transcendental equation
	# Update 2022: for all pixels!
	snow["alb_sph"] = OLCI_scene.toa * np.nan
	ind_solved = snow.r0.notnull()
	iind_solved = dict(xy=np.arange(len(ind_solved))[ind_solved])
	snow.alb_sph[iind_solved] = 1
	def solver_wrapper(toa, tau, t1t2, r0, u1, u2, albatm, r):
	# it is assumed that albedo is in the range 0.1-1.0
	return zbrent(
	0.1,
	1,
	args=(t1t2, r0, u1, u2, albatm, r, toa),
	max_iter=100,
	tolerance=1e-10,
	)
	solver_wrapper_v = np.vectorize(solver_wrapper)
	# loop over all bands
	for i_channel in range(21):
	snow.alb_sph.sel(band=i_channel)[iind_solved] = solver_wrapper_v(
	OLCI_scene.toa.sel(band=i_channel)[iind_solved],
	aerosol.tau.sel(band=i_channel)[iind_solved],
	atmosphere.t1t2.sel(band=i_channel)[iind_solved],
	snow.r0[iind_solved],
	angles.u1[iind_solved],
	angles.u2[iind_solved],
	atmosphere.albatm.sel(band=i_channel)[iind_solved],
	atmosphere.r.sel(band=i_channel)[iind_solved],
	)
	# ind_bad = snow.alb_sph.sel(band=i_channel) == -999
	# snow["isnow"] = xr.where(ind_bad, -i_channel, snow.isnow)
	ind_bad = snow.alb_sph.sel(band=i_channel) < 0
	snow.alb_sph.loc[dict(band=i_channel)] = xr.where(ind_bad,
	1,
	snow.alb_sph.sel(band=i_channel))
	# some filtering
	snow["alb_sph"] = snow.alb_sph.where(snow.isnow >= 0)
	# ind_neg_alb = (
	# (snow.alb_sph.sel(band=0) < 0)
	# | (snow.alb_sph.sel(band=1) < 0)
	# | (snow.alb_sph.sel(band=2) < 0)
	# )
	# snow["alb_sph"] = xr.where(ind_neg_alb, np.nan, snow["alb_sph"])
	# snow["isnow"] = xr.where(ind_neg_alb, 105, snow.isnow)
	# correcting the retrived spherical albedo for fractional snow cover
	snow["rp"] = snow.factor*(snow.alb_sph ** angles.u1)
	snow["refl"] = (
	snow.factor * snow.r0 * snow.alb_sph ** (angles.u1 * angles.u2 / snow.r0)
	)
	snow["alb_sph"] = snow["factor"] * snow["alb_sph"]
	ind_no_nan = snow["isnow"].notnull()
	snow["isnow"] = xr.where(snow.alb_sph.sel(band=0) > 0.98, 1, snow.isnow)
	snow["isnow"] = xr.where(
	(snow.alb_sph.sel(band=0) <= 0.98) & (snow.factor > 0.99), 2, snow.isnow
	)
	snow["isnow"] = xr.where(
	(snow.alb_sph.sel(band=0) <= 0.98) & (snow.factor <= 0.99), 3, snow.isnow
	)
	snow["isnow"] = snow["isnow"].where(ind_no_nan)
	return OLCI_scene, snow
	@numba.jit(nopython=True, cache=True)
	def f(albedo, t1t2, r0, u1, u2, albatm, r, toa):
	surf = t1t2 * r0 * albedo ** (u1 * u2 / r0) / (1 - albedo * albatm)
	rs = r + surf
	return toa - rs

This part solves the following equation:

Where Rc is the TOA reflectance measured by OLCI, Ra, Ta, ra are the reflecance, transmittance and spherical albedo of the atmosphere, R0 is the reflectance of the non-absorbing underlying surface reflectance, and rs is the surface spherical albedo.

and the directly calculated spectral albedo:

pySICE/sice.py

Lines 522 to 533 in 17ea753

	def snow_albedo_direct(OLCI_scene, angles, aerosol, atmosphere, snow, impurities):
	# direct caluclation including impurities
	# not sure where it is used
	alpha = 1000.0 * 4.0 * np.pi * (bai / wls)
	sdu = impurities.polut * wls ** (-impurities.bm)
	snow["alb_sph_direct"] = snow.factor * np.exp(-np.sqrt((alpha + sdu)*snow.al))
	snow["rp_direct"] = snow.alb_sph_direct ** angles.u1
	snow["refl_direct"] = (
	snow.factor * snow.r0 * snow.alb_sph_direct ** (angles.u1 * angles.u2 / snow.r0)
	)
	return OLCI_scene, snow

It is a new formulation presented in https://doi.org/10.3389/fenvs.2021.644551
and reads as:

Here is the (rephrased) comment from Alex:

The solved spectral albedo is calculated directly from reflectance avoiding any assumptions on impurities, which can be wrong. I trust the solved spectral albedo more...as compared to the direct calculation of spectral albedo.
As a consequence solved spectral albedo is being used for the broadband albedo calculation.

The next point- why I provide directly calculated spectral albedo as spectral output?
This is because I can make it for each OLCI wavelength. This is not possible with the solved albedo.

of course we can give 2 possibilities for the users - so they have 2 versions of spectral labedo (for selected wavelengths or for all wavelengths ( but with assumption on the type of impurities))

However, BBA must be always calculated from the (solved) reflectance-spherical albedo chain (not involving impurity assumption).