[MITgcm-devel] upcoming changes in seaice_growth.F

[Gael Forget]

Hi Ian,

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.

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? 

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

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