Noah MP Glacier: should snowalb_bats_glacier be called even if cosz<0?

[climbfuji]

Description

Tracking a bug in NEPTUNE/CCPP physics where I am getting a FPE invalid from Noah MP Glacier line 974 because:

DH DEBUG: sl1, sl2, cosz:   0.5000000E+00   0.4000000E+01  -0.2500000E+00

just before

cf1=((1.+sl1)/(1.+sl2*cosz)-sl1)

According to the note in

ccpp-physics/physics/SFC_Models/Land/Noahmp/module_sf_noahmp_glacier.F90

Line 969 in 5a36313

! when cosz > 0

this calculation is only valid for cosz>0. Looking through the entire file, there are notes all over the place that cosz should be between 0 and 1, but there are also if statements, for example in glacier_atm, that set swdown to zero if cosz < 0.

I am hoping for some feedback/suggestions from the NoahMP developers before I spend hours debugging this further (and fixing it the wrong way). Thanks!

[dustinswales]

That sure looks like a bug to me.

[rhaesung]

@HelinWei-NOAA @barlage Could you please provide any suggestions or feedback on this issue?

[HelinWei-NOAA]

I agree it is a bug. My suggestion is to add if (cosz > 0) block before calling both snowalb_bats_glacier and snowalb_class_glacier at line 845

[climbfuji]

Thanks @HelinWei-NOAA and everyone else. This is resolved by NCAR#1091 in NCAR ccpp-physics main.

[barlage]

Though I don't have an issue with this if statement proposed in NCAR#1091, this should probably be dealt with in the same way that it is in the non-glacier code (similar to what @HelinWei-NOAA said above). Unfortunately, that is done with a goto statement which I assume we want to avoid and should be replaced with an if block around the albedo calculation section, i.e., replace the goto 100 with then and 100 continue with end if.

However, I see another potential issue with this as a solution in the glacier code since the ice albedo is fixed. This means the albedo at night will be >0 and will then get replaced at the driver level. And will then get used on the time step when the sun comes up (if the land model is not called before the radiation time step). Unless the initialization of sag,fsa,fsr are moved up above the albedo calculation.

So, apologies for thinking through this here, but probably the solution proposed by @HelinWei-NOAA is good for now and we can deal with the timing issue later.