CAN be a negative values for ice_area_new (I don’t know what that would physically signify.
If that is the case then bot_band, which gets set to -1 in https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.glacier/src/glacier.c#L106, never gets updated and it results in a negative array index in https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.glacier/src/glacier.c#L119-L130 in either line 121 (if GLACIER_OVERFLOW is TRUE) or line 129 (if GLACIER_OVERFLOW is FALSE).
But negative array indices are not allowed I hear you protest - agreed, but since we do a lot of pointer arithmetic in C, I don’t think it is formally disallowed. The negative array index ends up clobbering something that should not be clobbered. Note that when I set GLACIER_OUTFLOW to FALSE and then replace line 129 with
snow[iveg][0].gl_overflow = 0.; # instead of snow[iveg][bot_band].gl_overflow = 0.;
then the model runs to completion on h2o (UW machine)
Note that this is not the solution. @jhamman: Should ice_area_new be allowed to be negative or how should that be handled to prevent mass balance errors? If we do not allow it to go negative, then what should the corresponding value of bot_band be?
From @jhamman:
I don’t think so. if (ice_vol > 0.0) (line 81) there should be a positive ice area.
GLACIER_OVERFLOW should only be non-zero (positive) when there are ice volume (area) exceeds the grid cell area (normalized). (edited)
From @bartnijssen
Agreed - note that ice_area_temp = pow((ice_vol / BAHR_C), (1 / BAHR_LAMBDA)); (l.87) is not negative.
The negative value happens in the time scaling a few lines later:
The exponent evaluates to a value very close to 1 (1.0000137) when you run at a three hour timestep, so the equation reduces to:
ice_area_new = -0.0000137*ice_area_old + 1.0000137*ice_area_temp
which ends up being negative if ice_area_temp is very small compared to ice_area_old. This would imply a big relative decrease in a glacier in a single time step, so perhaps the edge case when the glacier is very small??
From @jhamman
yes, okay. I think we’ll just need work around this edge case. Easiest is to skip time scaling when negative values would occur. Alternatively we can just move the rest of the glacier to the snow pack in this case but that will be harder to deal with.
From @anewman89
If the glacier is really small, it seems reasonable to assume no time scaling, as the scaling constant could be considered infinite. Little volume/mass would imply little change in area, just reduction in depth?
From @jhamman
I think this would be a fine way to move forward.
Via slack discussion from @bartnijssen
The problem happens in: https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.glacier/src/glacier.c in
int gl_volume_area(...). It turns out that the result of https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.glacier/src/glacier.c#L89-L91CAN be a negative values for
ice_area_new(I don’t know what that would physically signify. If that is the case thenbot_band, which gets set to-1in https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.glacier/src/glacier.c#L106, never gets updated and it results in a negative array index in https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.glacier/src/glacier.c#L119-L130 in either line 121 (ifGLACIER_OVERFLOWisTRUE) or line 129 (ifGLACIER_OVERFLOWisFALSE). But negative array indices are not allowed I hear you protest - agreed, but since we do a lot of pointer arithmetic in C, I don’t think it is formally disallowed. The negative array index ends up clobbering something that should not be clobbered. Note that when I setGLACIER_OUTFLOWtoFALSEand then replace line 129 withsnow[iveg][0].gl_overflow = 0.; # instead of snow[iveg][bot_band].gl_overflow = 0.;then the model runs to completion on h2o (UW machine) Note that this is not the solution. @jhamman: Shouldice_area_newbe allowed to be negative or how should that be handled to prevent mass balance errors? If we do not allow it to go negative, then what should the corresponding value ofbot_bandbe?From @jhamman: I don’t think so.
if (ice_vol > 0.0)(line 81) there should be a positive ice area.GLACIER_OVERFLOWshould only be non-zero (positive) when there are ice volume (area) exceeds the grid cell area (normalized). (edited)From @bartnijssen Agreed - note that
ice_area_temp = pow((ice_vol / BAHR_C), (1 / BAHR_LAMBDA));(l.87) is not negative. The negative value happens in the time scaling a few lines later:The exponent evaluates to a value very close to 1 (1.0000137) when you run at a three hour timestep, so the equation reduces to:
ice_area_new = -0.0000137*ice_area_old + 1.0000137*ice_area_tempwhich ends up being negative ifice_area_tempis very small compared toice_area_old. This would imply a big relative decrease in a glacier in a single time step, so perhaps the edge case when the glacier is very small??From @jhamman yes, okay. I think we’ll just need work around this edge case. Easiest is to skip time scaling when negative values would occur. Alternatively we can just move the rest of the glacier to the snow pack in this case but that will be harder to deal with.
From @anewman89 If the glacier is really small, it seems reasonable to assume no time scaling, as the scaling constant could be considered infinite. Little volume/mass would imply little change in area, just reduction in depth?
From @jhamman I think this would be a fine way to move forward.