[MITgcm-devel] upcoming changes in seaice_growth.F

[Ian Fenty]

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

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