[MITgcm-devel] upcoming changes in seaice_growth.F
Gael Forget
gforget at MIT.EDU
Fri May 27 07:16:15 EDT 2011
Hi Ian,
I believe I understand how you want those things handled. It never appeared to
me that it was worth having two different treatments of limit cases (one for
solve4temp and one for the area increments), since I have never seen those
type of things have a significant influence. I assume you have, so let's have two.
We don't need to have more "croppings" versions than necessary though. So I still
want to use the same in FENTY_AREA_... and in the rest of the evolution branch,
which I intend to implement as follows.
I had implemented the change we had agreed upon when you were here as
tmpscal1 = MAX(areaMin,AREApreTH(I,J))
hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
tmpscal2 = HEFFpreTH(I,J)/tmpscal1
#ifdef SEAICE_GROWTH_LEGACY
heffActual(I,J) = MAX(tmpscal2,hiceMin)
#else
heffActual(I,J) = sqrt(tmpscal2*tmpscal2 + hiceMin*hiceMin)
#endif
and use heffActual throughout as it was done until FENTY_AREA_...
Now we want
tmpscal1 = MAX(areaMin,AREApreTH(I,J))
hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
tmpscal2 = HEFFpreTH(I,J)/tmpscal1
#ifdef SEAICE_GROWTH_LEGACY
heffActual(I,J) = MAX(tmpscal2,hiceMin)
#else
heffActual(I,J) = sqrt(tmpscal2*tmpscal2 + hiceMin*hiceMin)
tmpscal1 = sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hiceMin*hiceMin)
recip_heffActual(I,J) = AREApreTH(I,J) / tmpscal1
#endif
and use " * recip_heffActual " in the area increments, throughout the evolution branch.
And as I said in my email last night, I will also CPP OUT the other bloc
that used areaMin ("treat the case of very small area").
all good?
Cheers,
Gael
On May 26, 2011, at 11:03 PM, Ian Fenty wrote:
> 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
>>> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20110527/54b69abc/attachment-0001.htm>
More information about the MITgcm-devel
mailing list