[MITgcm-devel] conservation // sublimation

[Ian Fenty]

Ice-people,

After discussing with Martin, I took a crack at fixing the non-conservation bug that Gael found in the sublimation code.  The approach I took, which works and is physically sound, is to cap the latent heat flux over the ice in solve4temp so that sublimation cannot remove more than the total mass of snow and ice present over a time step.  

A maximum latent heat flux (F_lh_max) is calculated in seaice_growth using the true (not regularized) ice and snow thickness and passed as an argument to solve4temp.  In the Newton-Rhapson loop, if F_lh is capped, the derivative of F_lh with respect to tsurfLoc is set to zero so that both numerical convergence and energy conservation are achieved. 

#ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
c     if the latent heat flux implied by tsurfLoc exceeds
c     F_lh_max, cap F_lh and decouple the flux magnitude from TICE 
          if (F_lh(I,J) .GT. F_lh_max(I,J)) THEN 
             F_lh(I,J)  = F_lh_max(I,J)
             dqhice_dTice = ZERO
          endif
#endif

In addition, the order in which thickness tendency terms are applied in seaice_growth must be changed such that loss of ice and snow by sublimation occurs before ice-ocean melting (if ice and snow are removed by ice-ocean fluxes prior to sublimation we can end up with nonzero r_FWbySublim).

Unlike other ice fluxes, the latent heat flux does actually have an upper bound which we simply overlooked in the past (which was a good approximation being that actual latent heat fluxes << maximum latent heat fluxes).  By bounding F_lh and switching the order of the d_HEFF operations, we are ensured no latent heat flux residuals and therefore no need for complicated retroactive accounting of energy and freshwater fluxes in QNET and EmPmR.  

Comments?  

-Ian


On Jun 7, 2011, at 9:29 PM, Gael Forget wrote:

> Hi Martin et al.,
> 
> as you may have notice, I did a slight re-organization of the
> SEAICE_ADD_SUBLIMATION_TO_FWBUDGET code.
> 
> The reason why I did this now is that, once we double check it,
> I think it should become the default for the EVOLUTION branch (i.e. 
> I would replace #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
> with a simple #ifndef SEAICE_GROWTH_LEGACY) and then
> include it in the seaice tracer stuff (seaice_tracer_phys.F).
> 
> I tested the changes with global_ocean.cs32x15, after adding 
> SEAICE_ADD_SUBLIMATION_TO_FWBUDGET, over one year. 
> I checked that the experiment thus tested all of the 'if ...' blocs, and it believe it 
> did. Then I checked that the results remained identical as I edited the code, 
> which they did. So I am pretty confident I did not change a thing. Am I mistaken?
> 
> However, one point got my attention, which I think is a bug that is likely to break 
> conservation. Indeed in the case when the pre-determined sublimation increment goes 
> beyond all of the ice/snow that is present, we take the rest of a_FWbySublim (now called 
> r_FWbySublim in analogy with the rest of the code) form the ocean, consistent with
> evaporation not sublimation. So I think we need to do something about the different latent 
> heats. It took me a while, and I had to ask for Jean-Michel's help, but I think we figured 
> the correct bug fix.  I wrote it down as a comment (line 1499 to 1503).
> 
> What do you think?
> 
> Cheers,
> Gael
>