@ifenty great!
Could you turn this into a "draft" PR? That way it's clear we are still working on this. I think you need to be author to do so, I don't see the button at my end.
Have you checked in a FWD/AD simulation that this change does not affect the result any more that any other numerical change? I would expect differences in the simulations of order <1% in variables like sea ice thickness, concentration, surface temperature, etc.

There are a few things to think about:

• TICES is still used in seaice_growth_adx.F to decide if precipitation will be snow or rain. This happens before calling solve4temp (which would compute new surface temperature). We'll have to move this part of the code to after calling solve4temp if we really want to get rid of the dependency on TICES of the previous timestep.
• There is a place in seaice_reg_ridge.F where TICES is reset if the ice thickness is below a certain threshold. I am not sure how to handle this cleanly. One option would be to have a CPP flag (like SEAICE_RESET_TSURF or TICES) and apply the change of seaice_growth_adx also to seaice_growth, just within this CPP flag?
• we could think of removing the pickup of TICES, if a flag like SEAICE_RESET_TSURF is defined.

I added the label "work in progress" which has similar meaning as "draft" as far as I know.

Miraculously, resetting ticeInMult = celsius2k does not affect the AD gradient, but only the FD gradients (from 16 to 9 digits of agreement):

Y Y Y Y 16>16< 9 16 16 16 16 16 22 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 16 16 16 pass  offline_exf_seaice  (e=0, w=33, lfd=1, dop=1)
Y Y Y Y 16>16<16 10 10  6  6 16 16  6  6 16 22  5  5 22 22 22 22 pass  offline_exf_seaice.thsice

* There is a place in seaice_reg_ridge.F where TICES is reset if the ice thickness is below a certain threshold. I am not sure how to handle this cleanly. One option would be to have a CPP flag (like SEAICE_RESET_TSURF or TICES) and apply the change of seaice_growth_adx also to seaice_growth, just within this CPP flag?

I guess the larger question is do we need to carry around the ice/snow surface temperature outside of seaice_growth_adx? it doesn't play any role in seaice_reg_ridge (e.g., sea-ice viscosity is not a function of surface temperature) or anywhere else besides being a diagnostic variable.

Also, can you @mjlosch explain what is the difference between kind=isbyte and byte =isbyte ?

it's the same, but "byte" is deprecated.

if ticeMultIn is set to celsius2k as first-guess and is updated to the EBM temperature after solve4temp, do we even need both ticeMultIn and ticeMultOut? maybe just a single ticeMult that is updated in solve4temp?

I agree, but that's even independent of TICES.

I guess the larger question is do we need to carry around the ice/snow surface temperature outside of seaice_growth_adx? it doesn't play any role in seaice_reg_ridge (e.g., sea-ice viscosity is not a function of surface temperature) or anywhere else besides being a diagnostic variable.

There's still "seaice_growth.F" which requires (so far) that this field is available, and in "cost_test_ice.F", we also need this as a global field. That's why I suggested a CPP flag: if we use this new method, TICES doesn't need to be in the pickup. The it's only needed in the pickup, as a diagnostics (which could be done locally in seaice_growth_adx, although this implies a staggering in time with respect to HEFF, etc.) and in "cost_ice_test.F"

I guess the larger question is do we need to carry around the ice/snow surface temperature outside of seaice_growth_adx? it doesn't play any role in seaice_reg_ridge (e.g., sea-ice viscosity is not a function of surface temperature) or anywhere else besides being a diagnostic variable.

There's still "seaice_growth.F" which requires (so far) that this field is available, and in "cost_test_ice.F", we also need this as a global field. That's why I suggested a CPP flag: if we use this new method, TICES doesn't need to be in the pickup. The it's only needed in the pickup, as a diagnostics (which could be done locally in seaice_growth_adx, although this implies a staggering in time with respect to HEFF, etc.) and in "cost_ice_test.F"

@mjlosch Oh, yes. another reason to carry TICES is so that ice/snow surface temperature can be used as a cost term. We should just set TICES equal to ticeMultOut after solve4temp.

@ifenty I turned this into a draft to make it clear that we are still working on this. There are few things to sort out, some of the points are repeated from a previous post.

• TICES is still used in seaice_growth_adx.F to decide if precipitation will be snow or rain. This happens before calling solve4temp (which would compute new surface temperature). We'll have to move this part of the code to after calling solve4temp if we really want to get rid of the dependency on TICES of the previous timestep.
• There is a place in seaice_reg_ridge.F where TICES is reset if the ice thickness is below a certain threshold. I am not sure how to handle this cleanly. One option would be to have a CPP flag (like SEAICE_RESET_TSURF or SEAICE_RESET_TICES) and apply the change of seaice_growth_adx also to seaice_growth, just within this CPP flag?
• we could think of removing the pickup of TICES, if a flag like SEAICE_RESET_TSURF is defined.
• what are the effects/differences in a long run (order 1-10 years of simulation)? I assume that they are very small, but we should be able to quantify them. It would be great if @ifenty could explore this a little so that we can give some numbers or even a plot here.
• Do we want to restrict these changes to seaice_growth_adx.F, or also add them to seaice_growth.F for consistency? This would have the advantage that we can cleanly remove resetting TICES in seaice_reg_ridge.F

I propose this:

1. we get replace ticeInMult and tIceOutMult by a single ticeMult = TICES in seaice_growth.F and seaice_growth_adx.F, then pass it to solve4temp, copy the result to TICES, and then diagnose TICES after seaice_growth (or at the end), but no longer in seaice_diagnostics_state.F (this will change the diagnostics of TICES, but not the results)
2. We add a flag SEAICE_RESET_TICES. If set, TICES is reset in seaice_growth_adx.F, including using TICES to decide over snow or rain, i.e. use the surface temperature of ice AFTER solve4temp (will change only offline_exf_seaice)
3. We also add this flag to seaice_growth.F (will change many results at some low level)
4. decide if SEAICE_RESET_TICES should be the default or if it should even be removed again. This will affect all setups out there that use seaice, but only very slightly, as we can probably show that the differences between resetting TICES or not are small. Alternatively, we could skip the SEAICE_RESET_TICES flag and immediately go for resetting TICES after make sure that the effects are small.
5. cleanup: remove resetting TICES according to HEFF in seaice_reg_ridge.F. and remove a lot of store directives.

The intended outcome is that we break the dependency of TICES across the time loop cycles, with minimal impact on the results; this will help move PR #468 forward. TICES is still there as a global diagnostics to be used in any cost function that we want to have. We keep it in the pickups for now (although not necessary, alternatively we could have a CPP flag SEAICE_TICES_IN_PICKUPS or so).

@jm-c what do you think about this plan?

@ifenty I just realised in a different context that TICES (copied to ticeMultIn -> TSURFin) is used in S/R seaice_solve4temp to set the albedo.

This means that setting ticeMultIn = 0degC will not work as it will corrupt the albedo computations, and we still need TICES from the previous time step. If we want to proceed here we need to have some new ideas. What do you think?

@mjlosch I looked at the code and see the issue that you are raising. What if we allowed the albedo (and thereby the absorbedSW) to change within the iterations? We could start the iterations with the surface temperature set to 0 Celsius (melting condition) with its corresponding melting condition albedo. Then if the surface temperature in a later iteration changed to be <0 C, the albedo would change to its freezing condition value. from that iteration on, we would be in the same situation as if we had inherited the 'melting or freezing' condition information from the last iteration of TICES.

Perhaps another way of doing solving this problem is to have a single, trial iteration starting from T=0. If the new T is < 0 then we can say that the surface albedo has to be from a freezing condition and then proceed with the iterations without allowing albedos to change. Doing that would only make sense if the iterative temperature updates never switched directions (i.e., never went from 0, to <0, then back up to >=0).

What do you think?

@ifenty I also thought about including the albedos in the iterative loop. I foresee that this may jeopardise convergence as it adds (non-smooth) nonlinearity. Are you willing to test this?

This second option appears to much of a hack to me, but again, one can give it a try.

@mjlosch Sure, I can test it out. My intuition is that if we always start from a melting condition (Tsurf = 0C) then we'll be OK.

@ifenty should we pursue this PR any further, especially in the light of the dependence of the albedos on TICES?

@ifenty This is just a measure of the effect of incorrect albedo computations in this branch. I compared a pair of 3-year runs using this branch and the master branch. This branch has a slightly larger sea-ice concentration cost (by 4%) than the master branch.

@owang01 thanks for checking. I added the albedo computation into the iteration loop. If @ifenty is right, this should bring down your cost again and the differences between this branch and the default will be minimal. But I think we also need to check the differences in long forward integrations (O(100y)).

@mjlosch Thank you for implementing this. As I mentioned in Monday's meeting, @ifenty and I are also implementing something to solve the albedo computation. See https://github.com/owang01/MITgcm/tree/tices_reset_in_seaice_growth_adx. A three-year forward run with our implementation has almost the same sea-ice concentration cost (increase by 0.04%) as the master branch.

The way how it is implemented is as follows. We start using dry albedo for the iterations. During the iterations, if the sea ice surface temperature (tsurfLoc) is ever warmer than melting temperature, we set a switch (tsurfabovemelt as a function of i,j) to TRUE. After the iterations, if the switch tsurfabovemelt == FALSE, we will continue using dry albedo. Otherwise, we use wet albedo and compute absorbedSW, IcePenetSW, and F_ia based on wet albedo. So the fluxes would be consistent with the albedo we are using.

We are also doing a pair of 28-year forward runs again the master branch to see how large sea-ice would differ.

For 28-year runs, the sea-ice concentration cost difference between this branch and master is < 1% . (As pointed out by Martin, sea-ice concentration cost may not be a good matrix for testing the implementation).

@mjlosch https://github.com/owang01/MITgcm/tree/tices_reset_in_seaice_growth_adx has the latest modification to seaice_solve4temp.F we implemented.