[MITgcm-devel] upcoming changes in seaice_growth.F

Thu May 26 23:03:36 EDT 2011

Gael,

> thanks for confirming that hiceMin=0.05 is fine for the adjoint.
> I will wait until tomorrow to proceed because Jean Michel
> has other changes underway.

Yes, when implemented as:
heff_star = sqrt(HEFF  * HEFF + hiceMin*hiceMin)
heffActual = heff_star / area_star

then the minimum possible value of the "actual ice thickness" for the 
purposes of the effective conductivty calculation is then 
hiceMin/area_star, i.e. hiceMin when AREA = 1.

I'm emphasizing that we should still use the construct:

d_AREAbyATM(i,J) = HALF*tmpscal1*AREApreTH(I,J) / heff_star(I,J)

and *not*

d_AREAbyATM(i,J) = HALF*tmpscal1*heffActual

because there is no reason to put a minimum bound on AREA in the above 
calculation; AREA can be arbitrarily small.  It's a small HEFF that we 
need to avoid there.

>
> I will add a CPP option around the bloc of code that you dislike,
> and will use your advised sqrt(AREA * AREA + areaMin * areaMin)
> in computing actual thicknesses.
>
> I have one question though. I thought you wanted something like
> 10^-5 for the areaMin default because 0.15 was way too large.
> Do you now suggest 0.1 is best?

What I tried to convey was that for only the purpose of calculating 
heffActual, limiting AREA (the denominator) to 0.1 is fine and that I 
saw no reason to impose a minimum nonzero value on it otherwise.

I will just send you a copy of seaice_growth.F with the changes I 
suggested in my last message so that you can see what I have in mind.

-Ian

Cheers,
> Gael
>
> On May 26, 2011, at 6:00 PM, Ian Fenty wrote:
>
>>
>> On May 25, 2011, at 12:40 PM, Gael Forget wrote:
>>
>>> Hi Ian and co.,
>>>
>>> below I report what Ian and I discussed at the ECCO2 meeting,
>>> which I just got started with seaice_growth.F rev.121 (see below).
>>>
>>> So, to complete the 'merging' process, we planned the following:
>>> 1) a few changes in the EVOLUTION branch
>>>        1.1) reduce the areaMin default
>>>        1.2) always include a_QSWbyATM_cover in QNET
>>>        1.3) use Ian's heff_star sqrt formula in heffActual
>>>        1.4) add the Winton (?) formula for flooding
>>> 1.5) add a check for the areaMax range
>>
>>> 2) remove code duplicates and restore the integrity of the two 
>>> branches structure
>>>   (EVOLUTION: with Ian's codes // LEGACY: for backward compatibilty)
>>
>> The default non-legacy behavior w.r.t. the thickening of existing ice 
>> is now the same as what I put in with:
>> #define FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
>>
>> and w.r.t. the behavior of open water air-sea flux 
>> divergence/convergence, the following are now equivalent:
>> #SEAICE_DO_OPEN_WATER_GROWTH with #FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
>> #SEAICE_DO_OPEN_WATER_MELT with #FENTY_OPEN_WATER_FLUXES_MELT ICE
>>
>> so it seems that the following can be safely removed:
>> #FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
>> #FENTY_OPEN_WATER_FLUXES_MELT_ICE
>>
>> which you actually did it in 1.121.
>>
>>
>>> 3) add somewhere in the code (SEAICE_OPTIONS.h or seaice_check.F may be)
>>>   the combination of CPP options that is advised for adjoint runs.
>>>
>>> With regard to 1.2) I felt it was best not to do it in the LEGACY 
>>> branch,
>>> since it would affect backward reproductibility. Makes sense, right?
>>
>> That makes sense to me.
>>
>>> With regard to 1.4) I argued at the time that we may as well replace the
>>> flooding formula in the EVOLUTION branch (only). But may be we should
>>> keep both and add a CPP option. Opinions?
>>
>> I think that since the non-evolution branch is kept around so that 
>> some old results can be reproduced it is best to leave what is there 
>> alone.
>>
>>>
>>> With regard to 1.1) Ian argued at the time that we don't need
>>> areaMin (old A22) to be a run time parameter. Opinions?
>>
>>> In any case we would leave the 0.15 default for the LEGACY branch.
>>> For the EVOLUTION branch we talked about 10^-5 if I remember right.
>>
>> I still don't understand why we would ever use areaMin > 0 and not 
>> areaMin = 0.  Why do we need to prevent the model from generating 
>> arbitrarily small amounts of AREA?  Maybe I'm missing something 
>> really obvious.
>>
>> The only place where a minimum non-zero value of AREA might be useful 
>> is when we are calculating the "actual" ice thickness (heffActual = 
>> HEFF/AREA) for the heat conduction calculation.  I suggest that we 
>> simply maintain a lower bound of AREA > 0 to keep things physical 
>> when HEFF > 0 and for the calculation of heffActual (and hsnowActual) 
>> to do something like:
>>
>> area_star = sqrt(AREA * AREA + 0.1 * 0.1)
>> heff_star = sqrt(HEFFpreTH  * HEFFpreTH + 0.05 * 0.05)
>>
>> heffActual = heff_star  / area_star
>> hsnowActual = HSNWpreTH / area_star
>>
>> The above caps the max value of heffActual (for the purposes of the 
>> heat conduction only) to 10*HEFF and the minimum value of 0.05 m and 
>> maintains smooth gradients with respect to AREA and HEFF as AREA, 
>> HEFF --> 0.  For actual ice thicknesses of >= 10 m, conductive fluxes 
>> are effectively zero anyway so we don't need to use a number smaller 
>> than 0.01 _d 0 in the area_star calculation.
>>
>> And in case it isn't obvious, a lower bound on heffActual is required 
>> because it is used in a denominator in seaice_solve4temp:
>>
>>     effConduct(I,J) = XKI * XKS /
>> &        (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
>>
>> I experimented with the heff_star and area_star formulations written 
>> above already in my local version and everything worked very well. 
>>  So I propose we do implement that.
>>
>>>
>>> With regard to 2) I left the FENTY_AREA_EXPANSION_CONTRACTION part 
>>> for later,
>>> since 1.3) is a pre-requisite. Ian, you were going to check that 
>>> using 0.05^2 rather
>>> than 0.1^2 in the heff_star sqrt formula is fine in adjoint mode. Is it?
>>
>> I checked that 0.05 m works just as well in the adjoint as 0.1 m in 
>> the contraction bits.  Because it is the denominator which we have to 
>> prevent from getting arbitrarily small to keep gradients bounded, we 
>> can keep the calculation in the following form:
>>
>> d_AREAbyATM(i,J) = HALF*tmpscal1*AREApreTH(I,J) / heff_star(I,J)
>>
>> Reactions?
>>
>> -Ian
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org <mailto: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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20110526/23476cf5/attachment.htm>