flwout set to zero at ice free points

[dabail10]

In ice_flux.F90 (and icedrv_flux.F90 in icepack) we set flwout to zero so that it can be accumulated in merge_fluxes. However, there is a timestepping bug here such that if you have a point that goes from ice free to ice covered in a single timestep, then flwout will not be computed in merge_fluxes (because aicen_init is all zero) and stays zero. However, the coupler fluxes are passed based on aice, so for a single timestep flwout is passed to the coupler with a zero value. This causes a crash in CAM and likely is a problem in other atmospheric models. I'm not sure why this was never a problem before.

The solution is to remove the line:

flwout = c0

in ice_flux.F90 and icedrv_flux.F90 and only do this just before the call to merge_fluxes based on aice_init. Working on it now.

[dabail10]

I have a solution. Let me know if you think this is ok. Basically, after merge_fluxes where lwout is computed, I make sure it never goes above the open ocean value (negative).

      ! Make sure flwout is not zero.

      TsfK = sst + Tffresh
      flwout = min(flwout,-stefan_boltzmann * TsfK**4)

[anton-seaice]

@kieranricardo - this might be relevant to you ?

[dabail10]

@peverwhee

[dabail10]

Actually I have this backwards. The reasonable range of temperatures in sea ice covered regions would be:

-80C to 0C or 193.15 to 273.15K.

So, this means a range of 79 W/m^2 to 316 W/m2. So, this should be:

flwout = min(flwout, -79)

since flwout is negative.

[dabail10]

Actually a better option is to add a check in scale_fluxes something like:

if (aice > puny .and. abs(flwout) < puny) then
   flwout = flwout_ocn
endif

This makes sure to catch the one special case of a point going from ice-free to ice-covered and the flwout is zero.

Dave

[dabail10]

I have a PR in for this now, but it occurs to me that there may be similar problems with fswabs, fsens, flat ... while a zero value is valid in these fields, we should technically be setting points that go from ice free to ice covered to have an open ocean value for these. For example, you can look in init_coupler_flux for some of these initial settings for different seasons.

[anton-seaice]

This doesn't solve you problem Dave, but I am slightly confused about the current behavior here.

In scale_fluxes, if there is some ice (aice>0) then the exported longwave is the radiation from the ice only and doesn't include the ocean. However if there is no ice, then the exported longwave is the radiation from the ocean.

These seems inconsistent. If the ice area is small, then the exported longwave is small (as it doesn't contain ocean longwave). If the ice area is zero, the longwave history output is large!

This is in scale_fluxes, so I think it only impacts the history output ? Why isn't

https://github.com/CICE-Consortium/CICE/blob/6aa677ec49b415298b6368a19e5122f9c9d49f8e/cicecore/cicedyn/general/ice_flux.F90#L1263-L1267

flwout just c0 ?

[dabail10]

This was super confusing! Took me a long time to figure this out. Let me try to explain.

1. In init_coupler_flux we have:

      flwout  (:,:,:) = -stefan_boltzmann*Tffresh**4
                        ! in case atm model diagnoses Tsfc from flwout

So, for all points this is set to the lwout equivalent to 0C (273.15K)

2. However in init_flux_atm:

flwout(:,:,:) = c0

This is everywhere so that flwout is computed as flwout = flwout + aicen(n)*flwoutn(n) for n=1,ncat in merge_fluxes. Since we are accumulating over subgridscale categories, this has to be initialized to zero.

3. The call to merge_fluxes is only done for aicen_init(n) > puny. So, flwout is only computed when there is ice at the beginning of the step.

4. However, it is possible that ice grows during the step and in scale_fluxes we do:

flwout = flwout / aice

where aice is the ice area at the end of the step. So, if flwout was not computed at the beginning of the step, then flwout is zero which is not physically possible.

Does this make sense?

[apcraig]

It seems like there might be two different issues.

1. What's on the history file. We can define this however we like. Ice free areas can have flwout whatever. This assumes the history file output is not being used to validate conservation. (see point 2).

2. What's being coupled. So, this seems important for conservation. Shouldn't the implementation be based entirely on conservation. The lw flux passed out at the ice surface should be consistent with the lw flux "used" in the ice model. If the ice fraction is zero, the lw flux value plays no role. But ultimately, the lw flux should be conserved. When ice goes from/to zero fraction, what lw flux is used in the ice model and are we conserving when we pass to the coupler? If lw flux is zero when aice > 0, and we fix that by setting lw flux to the ocean temperature lw flux, that seems wrong. If lw flux is non-zero when aice=0 (ice fraction -> 0), that seems wrong too. I wondered whether maybe the lags in the system were taking into account missing lw flux, but I don't think they can be. It could be that we have a conservation error that is reduced with this fix, but it still seems like we're not conserving. Is there an implementation that eliminates all conservation errors?

[dabail10]

This is not really a conservation issue. What the atmosphere is receiving is:

flwout_atm = aice flwout_ice + (1-aice) flwout_ocn

Because aice is non-zero, this means the atmosphere can receive a zero flwout from the ice meaning a temperature of zero degree Kelvin! This is completely unphysical. You might think that aice is very small as it is the growth over one timestep, but in CESM we are seeing the ice area go from zero to 0.8 in a single timestep over all of Hudson Bay!

[apcraig]

It is a conservation error as you describe it. If aice is 0.0001 when ice first grows, one could argue the error is small. If aice is 0.8, that's interesting and significant, but in either case, it seems like the system is not conserving. Changing the flw from zero to something associated with the ocean temperature may be reducing the error, but it's not conserving the heat. The flw*aice that the coupler receives should match exactly the long wave flux "used" in sea ice in order to conserve. Going to and from aice=0 may be a problem that we cannot overcome give the current implementation constraints, so minimizing the error may be the best option. But I'm pretty sure there is a conservation error.

I thought we were coupling the ice model each ice timestep to avoid these kinds of problems, in part by avoiding having to average (or similar) an ice fraction when some timesteps have zero ice fraction, but for other reasons related to ice state switching to/from zero ice fraction. I'm a little surprised this error has been in the implementation for so long. It must be quite small most of the time and especially when assessing conservation over longer periods.

[dabail10]

I agree this has been there a long time! They only recently added a check to the atmospheric radiation for this case. You do raise a good point in that the flwout should be computed based on the first temperature of the ice that forms in that timestep. Here we could use the Tsfc that is computed in the thermodynamics.