Shortwave issues in the CICE6 log

[cecilehannay]

The run cam_dev/MOM6/CICE6 (#115) shows some large sw errors or the order of E-02 in the cice6 log (see: #115 (comment))

arwt incoming sw (W) = 6.85618158455401250E+13 1.43538638248901375E+14
arwt absorbed sw (W) = 6.94156657060417812E+13 1.43187650989964625E+14
arwt swdn error = 1.23005355032898781E-02 -2.45123972989368718E-03

We run two short experiments:
005.cam6: cam6/MOM6/CICE6
005.cice5: cam_dev/MOM6/CICE5

that points to the cam parameterization reordering causing these errors as the error in 005.cam6 is of the order of E-16

arwt incoming sw (W)   =    5.21296342785727903E+07   1.63098032321698850E+15 
arwt absorbed sw (W)   =    5.21296342785727829E+07   1.63098032321698850E+15 
arwt swdn error        =   -1.42924093808776020E-16   0.00000000000000000E+00

---

Gokhan and Gustavo looked at the heat budget tables from cases cam6.005 and cice5.005. There is an issue in cice5.005 (see below). This is a showstopper they think we should solve that before proceeding with fully coupled runs.

hnetsw looks ok in the cam6.005:
/glade/scratch/hannay/archive/b.cesm3_cam041_mom.B1850WcMOM.ne30_L58_t061.cam6.005/logs/diags.log.4204015.chadmin1.ib0.cheyenne.ucar.edu.220518-144444.gz

On the other hand, the same term does not look good in cice5.005 (hnetsw ~ 0.27):
/glade/scratch/hannay/archive/b.cesm3_cam041_mom.B1850WcMOM.ne30_L58_t061.cice5.005/logs/diags.log.4210639.chadmin1.ib0.cheyenne.ucar.edu.220518-161801.gz

---

Let's see if I can describe this using equations. For simplicity, let's assume the coupling is hourly.

The old CAM ordering:

1. We receive Fswdn at step 0 and we have an intial guess at albedos. So, our net shortwave is (1-alpha(0))*Fswdn(0)

2. The sea ice state is updated and at the end of the step we compute alpha(1) and (1-alpha(1))*Fswdn(0)

3. We send the albedos to CAM and it then computes Fsdwn(1) and sends it to the ice.

4. So we compute the mismatch between (1-alpha(1))*Fswdn(0) and (1-alpha(1))*Fswdn(1). Then we rescale the absorbed shortwave to account for this mismatch. The error reported in the log is the mismatch.

5. We recompute (1-alpha(2))*Fswdn(1) and send the albedos back to CAM ... Note the the computation of absorbed shortwave and albedos is done together and cannot be separated.

WIth the new CAM ordering:

1. We receive Fswdn at step 0 as before and using the initial guess at albedos we have (1-alpha(0))*Fswdn(0).

2. The sea ice state is updated and we compute (1-alpha(1))*Fswdn(0) and send the albedos back to CAM.

3. Now CAM resends Fswdn(0) ...

4. So, then we compute (1-alpha(1))*Fswdn(0). This should identically match what we had at the end of last timestep.

5. Then we update the sea ice state and recompute (1-alpha(2))*Fswdn(0) and we send albedos.

6. CAM then sends Fswdn(1) (which used the alpha(1) albedos).

7. We then compute (1-alpha(2))*Fswdn(1) and compute the mismatch with (1-alpha(2))*Fswdn(0). So, the difference is that we are using the incoming shortwave from n and n-1 here. Where in the old way we were using n+1 and n whereas the albedo is at n+1 always.

So, with the new ordering, we are an extra hour off. Now, if we move the albedo/shortwave calculation to the beginning of the sea ice timestep, then internally the sea ice albedo and shortwave would be in sync. However, CAM would be receiving albedos from one hour earlier. We had the rescaling in CICE to deal with this one-hour lag, which is now a two-hour lag and hence the rescaling error is much larger.

Does this sort of make sense? The system is still conserving energy with the rescaling, it is just a matter of how much rescaling we have to do.

---

A couple questions:

After CICE estimates the absorbed solar at n+1 w/ (1-alpha(n+1))Fswdn(n), you say that only the albedo's are sent to CAM. Then what's the point of estimating absorbed solar if what you need to provide is the albedo's to CAM? Is it solely to provide the correction at the next time-step?

You say or imply a couple times that the albedo's are sent to CAM and then it is used to provide Fsdwn to CICE. It's counterintuitive to me that Fsdwn would be dependent on albedo ... but I'll take your word for it if you say so.

One correction is that for the reordered physics, we actually call the radiation twice in nstep=0. Once right before the CPL step so as to provide something to initialize the surface models ... and then again in the new radiation location near the end of the time loop (note we don't apply that first radiation call to the CAM state --it's solely to provide fluxes for the surface models). So, maybe consider subscripting nstep=0 like this: Fswdn(0a) and Fswdn(0b). At nstep=1, you will receive Fswdn(0b), and so you will need to do the correction since this is different than the absorbed solar provided in nstep=0, Fsdwn(0a).

You are correctly indexing Fsdwn with time-subscripts that are offset by 1 compared w/ the default new ordering. One could argue you should be subscripting alpha similarly; CICE should be using the nextsw_cday to compute a zenith angle from computing the albedo. nextsw_cday is consistent with the Fswdn you receive. It is the calendar day pertaining to when the Fsdwn was computed.

Your saying that the correction is to adjust for the difference is Fsdwn between n+1 and n. And in the new ordering, the correction is the difference between n and n-1. How would that lead to larger corrections if It's still a difference of one time-step you are correcting for. I'm afraid I'm missing something fundamental about how you see this leading to larger corrections.

---

So, it’s not that we are estimating net shortwave. We cannot decouple the albedos from the net shortwave. This is the end result of the call to our radiation scheme. We have to correct once we receive new incoming shortwave. CAM uses the surface albedos from all components to get the upward shortwave in the calculation. Not sure if it ever was the case where CAM received the upward shortwave. I realized later that step 0 is special. I should have started with step 1. I believe that it is a difference of two hours instead of one. I think the error is in the mismatch of the albedos and the incoming shortwave.

---

Within each gridcell, the ice model includes 5 ice thickness categories that have vertically resolved temperature. We perform a radiative transfer calculation for each ice category and possibly 3 different surface types within each category (bare ice, snow-covered ice, ponded ice). This then provides the absorbed SW within each vertical layer of each ice category (which we need for our ice thermodynamics) and a consistent albedo of each surface type for each category. So this is not just zenith angle dependent but also dependent on the ice state (ponding, snow, thickness) that is changing at each timestep. A single category area-weighted albedo is then provided to the coupler for the ice-covered portion of the gridcell. As mentioned below, I believe that a surface albedo is needed for CAM to perform its radiative transfer. Sorry if this is more information than needed but I wanted to clarify why we do our own SW radiative transfer within the ice model. I believe that the land model may do something similar.

I expect that the CAM radiative transfer does require a surface albedo for its simulation (and so for the calculation of Fswdn). For example, multiple reflections between the surface and clouds, which depend on the surface albedo, are important for the resulting surface Fswdn. I'm guessing that these impacts are included in the CAM RT. Someone please correct me if I'm wrong.

---

Thanks for the additional information. It is helpful. I can see why the albedo's are indexed with the current time-step, since they depend on the state at the current time-step.

Using Dave's terminology and indexing, in the default ordering, the magnitude of the correction at time step n is:
diff = (1-alpha(n))*Fswdn(n-1) - (1-alpha(n))*Fswdn(n)
For the new ordering, it is:
diff = (1-alpha(n))*Fswdn(n-2) - (1-alpha(n))*Fswdn(n-1)
Is this accurate? If so, I don't see how the magnitude of "diff" is different in these two scenarios.

---

Thanks for the additional information. It is helpful. I can see why the albedo's are indexed with the current time-step, since they depend on the state at the current time-step.

Using Dave's terminology and indexing, in the default ordering, the magnitude of the correction at time step n is:
diff = (1-alpha(n))*Fswdn(n-1) - (1-alpha(n))*Fswdn(n)
For the new ordering, it is:
diff = (1-alpha(n))*Fswdn(n-2) - (1-alpha(n))*Fswdn(n-1)
Is this accurate? If so, I don't see how the magnitude of "diff" is different in these two scenarios.

---

The observation is that the SW corrections in CICE are larger with new ordering, regardless if I understand it.

Dave, you had mentioned that there is a limitation to the correction, in particular if Fswdn at "n" is 0 (night-time), and "n+1" is non-zero (daytime), then the correction at n+1 won't work. If that is true, then if the corrections are larger with the new ordering, then this "leak" of energy due to day/night transitions is also larger. Does that seem reasonable?

Because in the grand scheme of things, it is this 0.27 W/m2 imbalance in the "hnetsw" that is the actual problem. If the CICE corrections were perfect, than we wouldn't have this problem. The problem is that the corrections aren't perfect and there's always an energy leak, and the leak becomes intolerable with the new physics ordering.

Dave and Marika, are we going to try and do a run with your first guess as to the fix -- moving the albedo calculation to the top of CICE time-loop?

Returning to Gokhan's point. In the worst case scenario, we can go back to the old physics ordering ... but w/ a dirtier fix to handle the issues that the reordering was intended to solve. It's dirty because it will apply a correction to the surface stresses when they are applied to the atmosphere. That is, the atmosphere would see a different surface stress than the rest of the CESM components. Not ideal IMO, but this is what NorESM does and it does seem to work for them.

---

@adamrher, @dabail10,

[cecilehannay]

While we are trying to figure out the sw issue, should be do a coupled simulation for the week-end with:
cam6-zm2 + MOM6 + CICE5
it would allow to test the salinity bug fix.

We could run 15-20 years during the week-end.

[JulioTBacmeister]

Sounds like a good idea.

[cecilehannay]

From @adamrher:

I can summarize how the new ordering increases time truncation errors in radiation. But time-truncation errors are separate from conservation errors, which is the current problem. Should I add a slide showing that the surface models are now getting radiation offset by a time-step in the new ordering?

How sure are we that the 0.27 W/m2 errors are coming from CICE? Do the SW correction in CICE conserve energy? I had thought there are scenarios where this does not conserve energy, when a grid point goes from night to day in one time-step, then the correction cannot be applied.

---

From Marika:
@gustavo-marques: Do we know that the coupler budget error is related to CICE? Or to the CAM SW re-ordering?

[gustavo-marques]

The coupler budget error only occurs when we use CAM SW re-ordering (cam_dev).

[cecilehannay]

-- From Dave:

Just a quick update. Keith Oleson and I have figured out that the land does indeed scale. They save fractions of the total absorbed flux for each layer. Then when new incoming shortwave comes in, the total absorbed is reset and then each layer absorbed is set as a fraction of the total again.

So, in Fortran pseudo code:

fswabs_new = (1-alvdr)*swvdr + ... (for all four components of the shortwave)

fswsfc_new = fswabs_new * (fswsfc / fswabs)

So, this says to set the fraction of the total absorbed for the sfc times the new absorbed. You can actually rewrite this as:

fswsfc_new = fswsfc * (fswabs_new / fswabs)

where fswabs_new / fswabs is the scale factor we use in the sea ice. Now, the trick is if fswabs (the old) is zero. In the sea ice we just set the new components to zero. However, I think it would be more correct to just set 100% of the total absorbed to be absorbed right at the surface. This might be the root of the error. That is, we would be two timesteps from the day-night boundary instead of one and there is more area. I will implement this and see if that changes the 'scaling error'.

---

We are not sure if these budgets can be obtained over shorter intervals than a month. I believe Gustavo will be checking that.

Should we be worried about a cloud-free radiation initialization in general? Large (SW) fluxes can result in significant melting of snow and sea-ice during the first coupling interval. Can there be perhaps a more reasonable approach for this?

---

Does this cloud free radiation in nstep=0 with cam_dev occur in the first time step of RUN_TYPE=hybrid runs?
If so, that sounds problematic. A use-case where this might really matter is short-term prediction runs, which have RUN_TYPE=hybrid.

---

I don't think we need to worry about the cloud-free radiation initialization in general. Nonetheless, I've done several comparisons w/ and w/o reordering, (1) the climate is basically indistinguishable, (2) the solutions follow a to the naked eye an identical trajectory, and then at about 48 hours the solutions begin to noticeably diverge. So the differences due to reordering are similar to perturbation to the system.

In the grand scheme of things, how much harm can we do in a single time-step? And in the context of the clear sky radiation at nstep=0, apparently not much harm. But I now see Keith's email and so I expect some push back!

I believe the clock is reset in hybrid runs, so nstep=0 would trigger the cloud free radiation ... assuming its not reading the radiation off the atm restart file, which it very well could be if rpointer.atm is set to a atm restart file.