snow density fix #1283
snow density fix #1283
Conversation
| do c = bounds%begc,bounds%endc | ||
| snowdp(c) = snow_depth(c) * frac_sno_eff(c) | ||
| frac_sno_update = frac_sno_eff(c) & | ||
| + tanh(0.1 * (qflx_h2osfc_to_ice(c)+qflx_soliddew_to_top_layer(c))*dtime) & |
There was a problem hiding this comment.
the 0.1 should be the accum_factor from SnowCoverFractionSwensonLawrence2012Mod, but currently it is a private parameter. What's the best way to use that parameter here?
| if (h2osoi_liq(c,lev_top(c)) < 0._r8) then | ||
| ! if (h2osoi_liq(c,lev_top(c)) < 0._r8) then | ||
| ! in some case roundoff level negative values trigger this | ||
| if (h2osoi_liq(c,lev_top(c)) < -1.e-12_r8) then |
There was a problem hiding this comment.
the 0 was changed to -1e-12 b/c I found some cases where <0 was catching values of h2osoi_liq that were negative roundoff values of zero due to the removal of the constrained evaporation from the available moisture being roundoff level different values.
There was a problem hiding this comment.
- I think the appropriate fix for this is to change the
custom_rel_epsilonvalue in the call totruncate_small_values_one_levforh2osoi_liq. You could increase this by one or two orders of magnitude. The issue with the change you made is that in leaves some negative state sitting around, which we don't want.
| else | ||
| rho_avg=200._r8 | ||
| rho_avg = denice |
There was a problem hiding this comment.
I don't know why the old values weren't denice; they should have been.
There was a problem hiding this comment.
- Can you explain this change further? I think I see the rationale for the original code: the 800 value was presumably giving a reasonable bulk density limit for a very dense snow pack, and the 200 value was presumably giving a reasonable bulk density limit when there is a very thin snow pack. I guess I can see the argument for getting rid of the arbitrary 800 value, and instead allowing the snow pack to be as dense as ice. But it doesn't seem like denice is the appropriate value in the
elseblock (though I'm not sure if that value ends up mattering too much).
| endif | ||
|
|
||
| !===================== xm < h2osfc ==================================== | ||
| if(temp1 >= 0._r8) then ! add some frozen water to snow column | ||
|
|
||
| ! update snow depth | ||
| frac_sno_update = frac_sno(c) + tanh(0.1 * (-xm(c)))*(1._r8 - frac_sno(c)) |
There was a problem hiding this comment.
again, 0.1 should be accum_factor
There was a problem hiding this comment.
I have made a few comments inline. However, my biggest comment relates to the question you raised, about how to reference a constant from SnowCoverFractionSwensonLawrence2012Mod. As is often the case, the need to access a private module variable actually indicates a bigger issue, I think. In this case, the problem I see is that it looks like your changes might assume that you're using that snow cover fraction formulation, rather than the one in SnowCoverFractionNiuYang2007Mod or anything that might come along later.
- I think this needs to be reworked somehow so that it works consistently with any snow cover fraction method. This might require introducing one or more new subroutines in the snow cover fraction modules that do what you need.
If that's difficult, I wonder if there is a completely different way to solve the underlying problem (of inconsistency between two history variables) – such as making the H2OSNO history variable based on some new variable that is saved at the same time that snow_depth and snowdp are calculated – i.e., in the middle of the time step, rather than using the H2OSNO value from the end of the time step. I recognize, though, that that might lead to other issues, like H2OSNO now being inconsistent with some other values. Maybe a solution could be to compute the snow density diagnostic in the code at the appropriate time, ensuring that it is using self-consistent values? I'm not sure... just thinking out loud about some ideas if it's too hard to make this work generally across any snow cover fraction method.
|
|
||
| qflx_h2osfc_to_ice => b_waterflux_inst%qflx_h2osfc_to_ice_col , & ! Input: [real(r8) (:) ] conversion of h2osfc to ice (mm H2O /s) | ||
| qflx_soliddew_to_top_layer => b_waterflux_inst%qflx_soliddew_to_top_layer_col, & ! Input: [real(r8) (:) ] rate of solid water deposited on top soil or snow layer (frost) (mm H2O /s) [+] | ||
|
|
There was a problem hiding this comment.
- We generally try to keep the variables from a given instance lumped together in the associate block. So can you please move these added lines to either above or below the waterdiagnostic lines?
There was a problem hiding this comment.
I see two issues: 1) is the diagnostic SNOWDP and snow density correct and 2) is the actual snow depth being correctly calculated. I think that I found a couple of issues with the latter that this pr corrects (i.e the change in snow depth due to phase change and incorporation of frozen surface water into the snowpack). If it was just 1), then I think changing the location of the snowdp assignment could work, but I think 2) needs to include a recalculation of the scf.
| ! add h2osfc to ice conversion flux to snowmelt to trigger | ||
| ! re-calculation of scf when h2osfc freezes (snowmelt | ||
| ! amount is not used in sfc_method, just presence/absence) | ||
| snowmelt(c) = (qflx_h2osfc_to_ice(c) + qflx_snow_drain(c)) * dtime |
There was a problem hiding this comment.
- I'm uncomfortable with this rationale: it seems like this is just begging for someone to later add code that misuses this variable, assuming that it really is snowmelt (since that's what it's called). If you need to make this change, but all you need is a trigger in UpdateSnowDepthAndFrac, then the input variable should be renamed accordingly to not give the false impression that this is actually snow melt. I'm also confused about this, though, because doesn't qflx_h2osfc_to_ice indicate a gain by the snow pack, whereas qflx_snow_drain indicates a loss? Or am I misunderstanding qflx_h2osfc_to_ice?
There was a problem hiding this comment.
This was just used as a trigger to get the SCF to be updated. Normally, it is updated when newsnow > 0 or snowmelt > 0. If neither occurs, then SCF is not updated. I'm not sure, but I think the latter behavior was intentional, and due to undesirable behavior caused by frost deposition or surface water freezing. But now it seems like it would be good to update SCF for changes to due to surface water freezing.
Another way to do this would be to have a second call to scf_method%UpdateSnowDepthAndFrac where the newsnow argument is set to qflx_h2osfc_to_ice and bifall is set to denice. I think that is conceptually what is happening. I can test that.
| if (h2osoi_liq(c,lev_top(c)) < 0._r8) then | ||
| ! if (h2osoi_liq(c,lev_top(c)) < 0._r8) then | ||
| ! in some case roundoff level negative values trigger this | ||
| if (h2osoi_liq(c,lev_top(c)) < -1.e-12_r8) then |
There was a problem hiding this comment.
- I think the appropriate fix for this is to change the
custom_rel_epsilonvalue in the call totruncate_small_values_one_levforh2osoi_liq. You could increase this by one or two orders of magnitude. The issue with the change you made is that in leaves some negative state sitting around, which we don't want.
| else | ||
| rho_avg=200._r8 | ||
| rho_avg = denice |
There was a problem hiding this comment.
- Can you explain this change further? I think I see the rationale for the original code: the 800 value was presumably giving a reasonable bulk density limit for a very dense snow pack, and the 200 value was presumably giving a reasonable bulk density limit when there is a very thin snow pack. I guess I can see the argument for getting rid of the arbitrary 800 value, and instead allowing the snow pack to be as dense as ice. But it doesn't seem like denice is the appropriate value in the
elseblock (though I'm not sure if that value ends up mattering too much).
|
Needs more work... |
Description of changes
Modify updates to snow_depth and snowdp to provide better values for snow density.
Specific notes
Simulations previously created for SnowMIP showed unrealistically large density values when the snowpack was very thin.
This was caused by updates to snow depth that used the value of frac_sno from the beginning of the time step, rather than a value that represented the updates due to changes in snowpack due to incorporation of frozen surface water and frost.
Contributors other than yourself, if any:
CTSM Issues Fixed (include github issue #):
Resolves #1280
Are answers expected to change (and if so in what way)?
yes. snow_depth is calculated differently in some situations, which would impact surface energy budget and albedo.
Any User Interface Changes (namelist or namelist defaults changes)?
Testing performed, if any:
(List what testing you did to show your changes worked as expected)
(This can be manual testing or running of the different test suites)
(Documentation on system testing is here: https://github.com/ESCOMP/ctsm/wiki/System-Testing-Guide)
(aux_clm on cheyenne for intel/gnu and izumi for intel/gnu/nag/pgi is the standard for tags on master)
NOTE: Be sure to check your coding style against the standard
(https://github.com/ESCOMP/ctsm/wiki/CTSM-coding-guidelines) and review
the list of common problems to watch out for
(https://github.com/ESCOMP/CTSM/wiki/List-of-common-problems).