search

ufs-community / ccpp-physics

UFS fork for CCPP
Other
4 stars 34 forks source link

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

Open climbfuji opened 3 weeks ago

climbfuji commented 3 weeks ago

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 https://github.com/ufs-community/ccpp-physics/blob/5a363134a77535f35594e56b58ba1e6141a23d2b/physics/SFC_Models/Land/Noahmp/module_sf_noahmp_glacier.F90#L969

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 commented 3 weeks ago

That sure looks like a bug to me.

rhaesung commented 3 weeks ago

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

HelinWei-NOAA commented 3 weeks ago

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 commented 3 weeks ago

Thanks @HelinWei-NOAA and everyone else. This is resolved by https://github.com/NCAR/ccpp-physics/pull/1091 in NCAR ccpp-physics main.

barlage commented 3 weeks ago

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.