[MITgcm-devel] conservation // sublimation
Ian Fenty
ifenty at gmail.com
Wed Jun 29 17:57:17 EDT 2011
Martin, Gael,
I just checked in my 'a priori' sublimation fix for the EVOLUTON branch. Since Martin's test in the EVOLUTION setup failed even with the 'a posteriori' code as of 6/28, I would ask that the test now be re-run with my code.
If conservation is indeed achieved with my fix (as I'm certain it will be) then I think the 'a posteriori' fix should be wrapped in LEGACY ifdef blocks.
Martin: please check that the diagnostic SIrsSubl ( the "residual" freshwater flux from sublimation ) is aways < epsilon and SIactLHF <= SImaxLHF (the actual latent heat flux is always less than or equal to the maximum latent heat flux)
Gael: I think that defining SEAICE_ADD_SUBLIMATION_TO_FWBUDGET with the 'a priori' fix (provided Martin can confirm that what I checked in works) should be the default in the EVOLUTION branch.
-Ian
On Jun 28, 2011, at 10:27 AM, Gael Forget wrote:
> Martin,
>
>> I have a hard time following the rapid changes to seaice_growth. I just do not have the time to understand new additions the moment they are committed, that's why I seem so silent (and wait for the "call of duty").
> makes sense. no worries.
>
>> Now I have (without thinking too much) rerun my tests as configured with the help of Ian two weeks ago, for legacy and new code after updating to todays (Jun28,10am CET) code with
>> #define SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
>> for both cases and I get this:
> Regarding the diagnostics update vs conservation test, there is one quick thing I would like you to try.
> In your conservation test script I would assume you are using both SIsnPrcp and SIfwSubl.
> This would be in the ice heat budget, when comparing SIatmQnt -(+) oceQnet with a
> delta thickness X fusion heat, an equation to which snow precip and sublimation need to
> be added (or removed depending on which side of the equation you are). Am I mistaken?
> If I am not, I would like you to replace SIfwSubl with SIacSubl there. Any luck?
>
> I am puzzled that the SEAICE_GROWTH_LEGACY switch changes the sublimation term correctness.
> Are you getting much better than residual = 6.511920e-02 kJ/m^2 for the ice heat budget, with
> define SEAICE_GROWTH_LEGACY & undef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET,
> leaving everything else unchanged?
>
> I will try to further look into it by the end of the week.
>
> Ian,
>
> will you be at the office, where I would give you a phone call, say thursday?
>
> thanks for the help,
> Gael
>
>> 1. For SEAICE_GROWTH_LEGACY defined
>>>>> check_conserve
>>> bdir = ../run
>>> global volume budget: ocean
>>> rhosw*dEta = -4.578598e+02, oceFWflx*dt = -4.578598e+02
>>> residual = 1.875833e-12 mm
>>> global volume budget: sea ice
>>> rhoi*dHi+rhos*dHs = 4.138988e+02, net fwFlx *dt = 4.138988e+02
>>> residual = 6.821210e-13 mm
>>> -New- salt budget: ocean
>>> dSalt = -2.330044e-08, oceSflux = 0, SFLUX = -4.075021e-13
>>> residual = -2.330003e-08 g/m^2
>>> -New- heat budget: ocean
>>> rcp*dT = -1.631824e+08, oceQnet= -1.631824e+08, TFLUX = -1.631824e+08
>>> residual = -3.039837e-06 J/m^2
>>> - global heat budget: sea ice
>>> rho*dH*flamb = 1.382422e+02, atmQnt= 3.014246e+02, oceQnet = 1.631824e+02
>>> residual = 2.682209e-13 kJ/m^2
>>
>> 2. for SEAICE_GROWTH_LEGACY undefined (and all sorts of other flags turned on, so this is the ever-evolving "evolution" extasy)
>>>>> check_conserve
>>> bdir = ../run_new
>>> global volume budget: ocean
>>> rhosw*dEta = -4.912088e+00, oceFWflx*dt = -4.912088e+00
>>> residual = -1.072920e-12 mm
>>> global volume budget: sea ice
>>> rhoi*dHi+rhos*dHs = 3.813073e+00, net fwFlx *dt = 3.813073e+00
>>> residual = -4.218847e-14 mm
>>> -New- salt budget: ocean
>>> dSalt = -4.660087e-08, oceSflux = 0, SFLUX = -1.731896e-14
>>> residual = -4.660086e-08 g/m^2
>>> -New- heat budget: ocean
>>> rcp*dT = 1.484117e+07, oceQnet= 1.484117e+07, TFLUX = 1.484117e+07
>>> residual = 1.053512e-05 J/m^2
>>> - global heat budget: sea ice
>>> rho*dH*flamb = -5.943748e+00, atmQnt= -2.085003e+01, oceQnet = -1.484117e+01
>>> residual = 6.511920e-02 kJ/m^2
>>
>> So you see, that the heat budget is broken for the EVOLUTION. I can't tell you why, but this is what Ian and I were getting all the time so that Ian came up with his more robust way of dealing with excess latent heat.
>>
>> Hope that helps.
>>
>> Martin
>>
>> On Jun 28, 2011, at 3:45 AM, Gael Forget wrote:
>>
>>> Ian,
>>>
>>> I am glad you sympathize with the implementation of the 'aposteriori fix'. The implemented fix (as opposed to the
>>> tentative fix I added as a comment in r125) passed my own tests of conservations (based on a modified version of Martin's
>>> earlier script). Yet I would certainly feel more confident if Martin or you could confirm that I got this right, or tell me otherwise.
>>>
>>> I guess I should have figured earlier from your emails last week that you and Martin had tested the tentative commented fix,
>>> and had found it to fail his conservations tests. (which also happened with my own tests yesterday). I was under the
>>> erroneous impression that you 'didn't look' at it much, since your first email was about approaching the problem differently.
>>> Me adding the tentative fix as a comment was just a way to open this thread. I was hoping for feedback from Martin in
>>> regard to conservation tests, for which he is the most expert. (the type of feedback I found in your last email typically).
>>> I understand I was not explicit enough about that when I opened this thread. Sorry for that.
>>>
>>> Anyway, over the week end I decided to give it another look, and I did conservations tests of my own.
>>> My conclusion essentially was that the fix was correct, but misplaced relative to diagnostics. After displacing
>>> the diagnostics and the fix, I felt confident enough that everything was correct and conservative to check it in cvs.
>>> I may very well have missed something though, so I would appreciate if Martin or you could double check.
>>>
>>> Martin,
>>>
>>> could you put the implemented fix to the test? if it was incorrect, could you point me why?
>>>
>>> Ian,
>>>
>>>> Anyway, I'm all for keeping an a posterior treatment (even though it won't do anything) because after spending too much time with the LEGACY solve4temp I decided that it was too much of a pain to try to cap the latent heat flux from within all of those A1, A2, and A3 intermediate variables. So, I'm applying my fix to the EVOLUTION seaice_solve4temp where I cleanly treat F_lh.
>>>
>>> in this context, when you check the 'apriori fix' in cvs, may be it would be best
>>> to bracket it with CPPs, at least for the time being. Would you agree with that?
>>>
>>> Cheers,
>>> Gael
>>>
>>>
>>>
>>>
>>>> A cap on the latent heat flux wasn't a part of the LEGACY experience anyhow so no one's results will change.
>>>>
>>>> -Ian
>>>>
>>>>
>>>> On 6/27/2011 6:22 AM, Gael Forget wrote:
>>>>> Ian,
>>>>>
>>>>> since we agreed that, regardless of your upcoming apriori fix, we
>>>>> need to do something about the r_FWbySublim pathological case
>>>>> aposteriori, I proceeded and included the fix I had proposed at
>>>>> the start of this thread.
>>>>>
>>>>> In second thought, the type of approach Martin had taken (applying the
>>>>> residual fresh water flux as evap -- as an aposteiori fix) makes more sense
>>>>> than an 'if r_FWbySublim > 0 then stop' for several reasons (e.g. why stop
>>>>> since we dont need to?). So I followed that road, and completed it properly.
>>>>>
>>>>> With regard to your upcoming apriori fix via seaice_solve4temp,
>>>>> this does not change anything, except that you dont
>>>>> have to worry about adding an aposteriori test anymore.
>>>>>
>>>>> Cheers,
>>>>> Gael
>>>>>
>>>>> On Jun 21, 2011, at 12:17 PM, Ian Fenty wrote:
>>>>>
>>>>>> Gael,
>>>>>> Our thinking on the subject is almost identical. In fact, the trial version I gave Martin had multicategory support and an almost identical posteriori test. I did overlook legacy solve4temp support, however. I'll commit following JM's checkpoint.
>>>>>> -Ian
>>>>>>
>>>>>> On 6/21/2011 8:12 AM, Gael Forget wrote:
>>>>>>> Hi Ian,
>>>>>>>
>>>>>>> thanks for taking charge of this. Your approach sounds fine.
>>>>>>> You probably thought it through already, but I would like
>>>>>>> to raise two issues before you proceed with the check-in.
>>>>>>>
>>>>>>> The a posteriori fixes (the one in line with Martin's code that
>>>>>>> I had annotated, and the other one that I would have argued for) would have
>>>>>>> been adequate with&without SEAICE_MULTICATEGORY, with&without
>>>>>>> SEAICE_SOLVE4TEMP_LEGACY. Please make sure your fix does too.
>>>>>>>
>>>>>>> Also, in the a priori fix approach (via solve4temp, as you implemented) I see how
>>>>>>> things could go wrong in the future. To remove the risk of the problem re-appearing
>>>>>>> and going undetected again, we now need an a posteriori test in seaice_growth.F.
>>>>>>> Something like "if r_FWbySublim<>0 then print_error and stop" right after the
>>>>>>> sublimation term has been applied. It is probably not great for vectorization,
>>>>>>> but safety should be the highest priority. So please add one such test.
>>>>>>>
>>>>>>> Cheers,
>>>>>>> Gael
>>>>>>>
>>>>>>> On Jun 20, 2011, at 7:40 PM, Ian Fenty wrote:
>>>>>>>
>>>>>>>>
>>>>>>>> Martin,
>>>>>>>>
>>>>>>>> On 6/17/2011 3:56 PM, Martin Losch wrote:
>>>>>>>>> Hi there,
>>>>>>>>>
>>>>>>>>> it looks like, Ian's "crack" is really fixing the heat conservation problem.
>>>>>>>>>
>>>>>>>>> Further I like the conceptual idea: you cannot not sublimate more mass than there is.
>>>>>>>>>
>>>>>>>>> I completely vote for this version.
>>>>>>>>>
>>>>>>>>> Now, what does the adjoint think of this?
>>>>>>>>>
>>>>>>>>
>>>>>>>> I tested today. The adjoint is stable when sublimation is treated with my fix (it may have been stable with the other method but I didn't look). Moreover, including sublimation in the adjoint modifies the dJ/dAQH in a perfectly interpretable way. Basically, the greater the AQH, the lower the latent heat flux/sublimation/mass loss, the greater the ice concentrations at the end of the simulation.
>>>>>>>>
>>>>>>>> So. That's two votes yes and an unknown number of votes abstaining. I put it my fix tomorrow unless there is objection.
>>>>>>>>
>>>>>>>> -Ian
>>>>>>>>
>>>>>>>>> Martin
>>>>>>>>>
>>>>>>>>> On Jun 17, 2011, at 11:06 AM, Ian Fenty wrote:
>>>>>>>>>
>>>>>>>>>> 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
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> _______________________________________________
>>>>>>>>>> MITgcm-devel mailing list
>>>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> MITgcm-devel mailing list
>>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> MITgcm-devel mailing list
>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> MITgcm-devel mailing list
>>>>>>>
>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>
>>>>>> _______________________________________________
>>>>>> MITgcm-devel mailing list
>>>>>> MITgcm-devel at mitgcm.org
>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> MITgcm-devel mailing list
>>>>>
>>>>> MITgcm-devel at mitgcm.org
>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>
>>>> _______________________________________________
>>>> MITgcm-devel mailing list
>>>> MITgcm-devel at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list