I rebased on ocean/develop, compiled, and this passes global ocean test cases.

@maltrud please look at the files changed tab and see if these look like the right code changes.

An extra line needs to be added here:
https://github.com/MPAS-Dev/MPAS-Model/blob/ocean/develop/src/core_ocean/shared/mpas_ocn_diagnostics.F#L1832
to modify surfaceFluxAttenuationCoefficientPrecEvapIce in the same way as surfaceFluxAttenuationCoefficient under land ice. Or we would need to add yet another surfaceFluxAttenuationCoefficientLandIce and another set of loops...

@xylar , for the first question, the effect of land ice is already included in activeTracersSurfaceFlux(indexTempFlux,iCell) so doesn't need to be handled explicitly.

for the second question, evap and precip are forced to have water entering/leaving at the local SST (instead of at 0C as it was before). we account for that added heat with an additional surface heat flux in activeTracersSurfaceFlux(indexTempFlux,iCell). but this flux doesn't change buoyancy (only heat content) so we have to remove it from the buoyancy flux.

@maltrud, your comment refers to KPP? Okay, I'll defer to you and desist on that.

The issue remains that surfaceFluxAttenuationCoefficientPrecEvapIce needs to be modified under land ice in the same way that we modify surfaceFluxAttenuationCoefficient. Otherwise, land ice freshwater fluxes will no longer be distributed as expected.

@xylar , for the third question, i chose not to touch land ice fluxes. i don't think that what i've done would change the land ice flux treatment at all. if you believe that it would be beneficial to treat land ice fluxes differently than they are now following a similar procedure, that's up to you.

@xylar ok i just saw your message about needing to modify stuff under land ice--i'll check that.

i chose not to touch land ice fluxes.

Because landIceFreshwaterFluxes are added to surfaceThicknessFluxes, which were prerviously distributed with surfaceFluxAttenuationCoefficient, and now are distributed with surfaceFluxAttenuationCoefficientPrecEvapIce, I very much do expect land-ice fluxes to be affected.

yes, i was incorrect before. instead of saying i chose not to modify land ice treatment, i should have said i chose not to exempt land ice from being treated the same as other (non-runoff) thickness fluxes. if you'd like to add another attenuation scale for land ice fluxes, go ahead and do that. this PR is for water cycle, so ice cavities won't be enabled. for testing with ice cavities enabled, one can set the attenuation scale back to 0.001 to get the previous behaviour.

@maltrud, if you look carefully at the line I posted before, which I will share again:
https://github.com/MPAS-Dev/MPAS-Model/blob/ocean/develop/src/core_ocean/shared/mpas_ocn_diagnostics.F#L1832
we currently modify surfaceFluxAttenuationCoefficient to have a different value under land ice than elsewhere. The current code (before this PR) thus distributes land ice fluxes (including thickness fluxes) over 10 m even if they are distributed over 0.001 m elsewhere.

With this PR, that functionality is broken unless the same treatment is added for the new surfaceFluxAttenuationCoefficientPrecEvapIce. I really must insist that some careful testing be done to make sure this doesn't break ice-shelf cases before this can go in, regardless of whether it is intended for water cycle runs. It isn't acceptable to me to break the cryosphere configuration in haste to get new features in for water cycle.

ok, there are too many issues to be solved before the freeze. maybe this will get in later, maybe not.

@maltrud, I'm sorry for my tone last night. As @mark-petersen kindly pointed out to me, I probably shouldn't be responding to pull requests at 1 am.

It seems like @vanroekel was fine with the KPP changes made here (the remaining concerns about fluxes from land ice and sea ice seem unrelated to this PR).

The only remaining change I needed was the addition of this one line:

            surfaceFluxAttenuationCoefficientPrecEvapIce(iCell) = config_land_ice_flux_attenuation_coefficient

below the line:

            surfaceFluxAttenuationCoefficient(iCell) = config_land_ice_flux_attenuation_coefficient

at
https://github.com/MPAS-Dev/MPAS-Model/blob/ocean/develop/src/core_ocean/shared/mpas_ocn_diagnostics.F#L1832

I got frustrated because it didn't seem like I was making myself very clear on this and that the runs with ice shelves might get broken in the haste to get the PR in before the code freeze. But the fix is simple enough that a quick nightly test could confirm that no harm is done once the above line is added.

If you wan to push for this to happen before the code freeze, I'd support that. I'll do what I can to make sure ice-shelf runs aren't negatively affected.

It was still worth having a PR with the code and issues raised for later reference. No harm done. Mat, don't worry about it now, but after you return we should discuss what is worth merging versus adding to a development list for later.

@maltrud, could you please comment here on what simulations you ran with the thickness flux spreading, what parameter values you used, and what differences you saw? Were there improvements that would indicate that thickness flux spreading should be on by default in V2 simulations?

@mark-petersen i did simulations with an attenuation scale 0.001, 10., and 50. meters. there were significant improvements in the 50m case, but i'm no longer fully confident in the results due to possible energy conservation issues. they're likely quite small errors, but i need to fix that before definitively saying it should be implemented in v2, especially since it won't be ready by july 15.