Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

icepack_parameters: optionally compute 'dragio' from under-ice roughness #366

Merged
merged 2 commits into from
Jun 22, 2021

Conversation

phil-blain
Copy link
Member

@phil-blain phil-blain commented May 27, 2021

PR checklist

152 measured results of 152 total results
151 of 152 tests PASSED
0 of 152 tests PENDING
1 of 152 tests MISSING data
0 of 152 tests FAILED

The MISS is the test I'm adding which is new so was not in the baseline I'm comparing to (43a909a)

  • How much do the PR code changes differ from the unmodified code?
    • bit for bit
    • different at roundoff level
    • more substantial
  • Does this PR create or have dependencies on CICE or any other models?
    • Yes: I'll then do a CICE PR to add calc_dragio and iceruf_ocn to the CICE namelist
    • No
  • Does this PR add any new test cases?
    • Yes
    • No
  • Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/.)
    • Yes
    • No, does the documentation need to be updated at a later time?
      • Yes
      • No
  • Please provide any additional information or relevant details below:

In some applications, we want the ice-ocean drag coefficient 'dragio' to be computed using an under-ice roughness length and the thickness of the first ocean level instead of imposing its value directly.

Add a new boolean flag 'calc_dragio' and real parameters 'iceruf_ocn', 'thickness_ocn_layer1', where 'iceruf_ocn' is the under-ice roughness length, and 'thickness_ocn_layer1' is the thickness of the first ocean level, to module 'icepack_parameters'.


EDIT 09-06-2021 I've updated the description with the test results and an up-to-date description.

@eclare108213
Copy link
Contributor

@phil-blain
Just a quick thought - please take a look at the form drag parameterization in the code and/or documentation. I think it has something similar to this, e.g. to include the effects of keels. If so, maybe we could generalize the implementation so that it works for your purposes? (It might not be in there in a usable form - my memory is not so good.)
e

@apcraig
Copy link
Contributor

apcraig commented May 27, 2021

I'm curious why this is in icepack at all. It seems you are just storing a value in icepack, but never using it in icepack.

@phil-blain
Copy link
Member Author

phil-blain commented May 28, 2021

@eclare108213, @apcraig

I'm not familiar with the form drag parameterization, but I know that one of our colleague looked at that an made some changes to neutral_drag_coeffs. I will eventually port these changes here also, and so will need iceruf_ocn in icepack_atmo, and that is why I'm implementing it in Icepack and not in CICE.

@eclare108213
Copy link
Contributor

I'm not opposed to this particular modification to the code, but since it's just one step in a larger set of changes, it would be helpful to see the big picture. Do you have something like a design document for the overall approach, or at least a basic outline/description of it?

@phil-blain
Copy link
Member Author

phil-blain commented May 28, 2021

Ok I can share a little bit of the game plan. Right now I have a working NEMO3.6/CICE6 codebase, but with some custom patches to CICE/Icepack to make it work with our version of NEMO:

1- the staggered atmoshepric level PRs
2- the ability to set iceruf_ocn in the CICE namelist, so that NEMO can access the value, and then compute dragio using the depth of the first oceanic level
3- modifications to the hadgem driver. We were historically using the hadgem driver so that's why I've been using it, but we've decided to have our own driver instead, so that will come in a PR probably next week.

After that we will be running NEMO3.6 with a "pristine" version of CICE6, i.e. no custom patches. But not all of our CICE modifications have been re-implemented in CICE6 yet, so as we upgrade our forecasting systems with the new NEMO-CICE executable we will progessively re-implement what we still need (the neutral_drag_coeff work is such an exemple). That will be over the next several months and the plan is that I not do all of this work myself; ideally internal ECCC contributors would re-implement their changes in CICE6 themselves and contribute them to the Consortium (but I might be dreaming here and end up doing most of that myself).

Maybe @JFLemieux73 and @dupontf want to add some more details but that's the gist of it.


So for this here PR it's really just making iceruf_ocn available to NEMO, which will require a subsequent CICE PR if we want to avoid having use icepack_parameters in NEMO. As I said I'm implementing it in Icepack because we might need it in the future, but I could be convinced otherwise.

@eclare108213
Copy link
Contributor

I agree that Icepack is the right place to put this parameter, as opposed to CICE. However, I question whether it needs to be in the ice model at all, since it's essentially just a parameter for an ocean model parameterization. If it were variable, e.g. depending on the ice morphology, then it would make sense to have it in the ice model. (And if that's the plan, then probably some coordination with the form drag implementation is advisable.) But as it is, it doesn't make a lot of sense to house it in the ice model. So you use CICE's namelist as the method for ensuring that the ice and ocean components are using the same parameters?

@phil-blain
Copy link
Member Author

it is used to compute dragio, which is an ice-related parameter, so I would argue that it does have its place in an ice model. Also, when doing sensitivity / tuning studies it is really useful to be able to tweak this value easily, so that's why we want to have it in the CICE namelist.

@eclare108213
Copy link
Contributor

eclare108213 commented Jun 2, 2021

I'm still thinking about this one, and there are several considerations. For example, dragio is not a namelist parameter in either Icepack or CICE -- it's just a parameter set in Icepack that is queried when needed. However, it looks like CESM and the NUOPC couplers change the value (somehow), without it being in namelist. If ECCC also wants to change the value, then we should fix the code to allow that. My take on this is that the current code is poorly implemented.

The follow-on question, then, is whether the value should be constant or variable over the grid and time. The code has an array option (Cdn_ocn) but its use might be limited to just the formdrag parameterization. Maybe that needs to be generalized?

Finally, I still don't think that iceruf_ocn should be a parameter in the sea ice component unless the sea ice component actually uses it somewhere. But this question is really one of coupled model architecture -- which component owns what pieces of information. We made the decision that Icepack would "own" its parameters and be queried by the driver to share a common value -- in that sense, the parameter belongs in Icepack. BUT iceruf_ocn isn't a parameter needed by Icepack or CICE (yet)! We've had several rounds of cleaning up unused variables/parameters, and this essentially would be one of those.

What could be interesting and useful for the broader community, though, is to share the parameterization in which you use iceruf_ocn to compute dragio. The natural home for that, I think, would be in the mixed-layer ocean parameterization. You wouldn't use that parameterization when running coupled, but it would provide context for having iceruf_ocn in Icepack. Would that make sense?

@phil-blain
Copy link
Member Author

Thanks for the feedback Elizabeth. I'll try to add some details.

I'm still thinking about this one, and there are several considerations. For example, dragio is not a namelist parameter in either Icepack or CICE -- it's just a parameter set in Icepack that is queried when needed. However, it looks like CESM and the NUOPC couplers change the value (somehow), without it being in namelist. If ECCC also wants to change the value, then we should fix the code to allow that. My take on this is that the current code is poorly implemented.

I think the approach is simply to call icepack_init_parameters with the specific parameter that you want to tweak, as is done in the cesm1 and cmeps drivers, like here:
https://github.com/CICE-Consortium/CICE/blob/23787406e0951a78f831042177922abe83e327d7/cicecore%2Fdrivers%2Fnuopc%2Fcmeps%2Fice_comp_nuopc.F90#L374

This would be the approach I would take. I don't need to change the way he current CICE/Icepack code is implemented for that.

The follow-on question, then, is whether the value should be constant or variable over the grid and time. The code has an array option (Cdn_ocn) but its use might be limited to just the formdrag parameterization. Maybe that needs to be generalized?

I think this is a separate question.

Finally, I still don't think that iceruf_ocn should be a parameter in the sea ice component unless the sea ice component actually uses it somewhere. But this question is really one of coupled model architecture -- which component owns what pieces of information. We made the decision that Icepack would "own" its parameters and be queried by the driver to share a common value -- in that sense, the parameter belongs in Icepack. BUT iceruf_ocn isn't a parameter needed by Icepack or CICE (yet)! We've had several rounds of cleaning up unused variables/parameters, and this essentially would be one of those.

Yes, I understand that. I also agree it's a little awkward if it's not used anywhere.

What could be interesting and useful for the broader community, though, is to share the parameterization in which you use iceruf_ocn to compute dragio. The natural home for that, I think, would be in the mixed-layer ocean parameterization. You wouldn't use that parameterization when running coupled, but it would provide context for having iceruf_ocn in Icepack. Would that make sense?

We could, but that's not really what I'm working on right now... As you write it would not be used in coupled mode... the parameterization is not very complicated, it's basically the same as what we used in atmo_boundary_layer for the ice/atm interface, although with an additional factor (\lambda^2):

      !compute dragio as a function of z0io and first ocean layer thickness
      zm=gdepw_1d(2)
      zmon2=0.5_wp*zm
      dragio=(vkarmn/log(zmon2/z0io))**2  !dragio at half first layer
      lambda=(zm-z0io)/(zm*(sqrt(dragio)/vkarmn*(log(2._wp)-1._wp+z0io/zm)+1._wp))
      dragio=dragio*lambda**2

Here zm is the thickness of the first ocean layer. This is described in this paper:
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JC010677

However this has only been tested in a coupled context, so I don't know if we want to implement in in the ocean mixed-layer param...

@phil-blain
Copy link
Member Author

Sorry if I sounded harsh above, it was really not my intention. Here is a (slightly edited) copy of an email thread we had yesterday:


@eclare108213 writes:

Hi everyone,
I'm trying to understand the reasons for, and implications of, this PR from Philippe:
#366.
I've looked at the paper [above paper] and I think I understand the reason for the code change -- the stress should be calculated differently when the numerical ocean levels or layers become small enough to lie well within the Ekman layer (or some other turbulence based distance below the ice). Is that correct?

What I'm struggling with is why that needs to be done in the ocean model instead of the ice model. Maybe I'm missing something, but as far as I can tell from the paper, all of the ocean parameters that are used in the new calculation are constants, not variable over the grid or time. We already do similar kinds of calculations for the atmospheric boundary layer, and this appears to be just a variation on the ocean drag coefficient... wouldn't it make more sense to have the ocean calculation also in the ice model? If not, why not?

The numerical instability issue might be a problem, limiting the coupling time step, but if all the parameters are constant, I don't think that would affect the choice of where the calculation is performed, since the ice model is already computing the drag term (or is that not true in your configuration). Is there something else that I'm missing?

I'm thinking that this alternative ice-ocean drag calculation would be a nice option to provide the rest of the community, even if you continue to compute it in the ocean model. If we provide it as an option in Icepack, then the namelist setup would make sense.

[...] I'd really like to understand the issues better, and to include this option in the code if possible.


F. Roy answers:

I've looked at the paper [above paper] and I think I understand the reason for the code change -- the stress should be calculated differently when the numerical ocean levels or layers become small enough to lie well within the Ekman layer (or some other turbulence based distance below the ice). Is that correct?

Correct! In the log layer formulation near the bottom surface of the ice, the stress or square of friction velocity (tau/rho) is constant, but the drag coefficient associated varies depending on the depth of the current you use. In other words, the current beneath the ice varies in the log layer profile but the associated stress is constant. Therefore the drag coefficient must adjust to the depth of the current you use.

What I'm struggling with is why that needs to be done in the ocean model instead of the ice model. Maybe I'm missing something, but as far as I can tell from the paper, all of the ocean parameters that are used in the new calculation are constants, not variable over the grid or time. We already do similar kinds of calculations for the atmospheric boundary layer, and this appears to be just a variation on the ocean drag coefficient... wouldn't it make more sense to have the ocean calculation also in the ice model? If not, why not?

The reason I had put the computation in the ocean model is just because I had the surface ocean model layer available in this part of the code. But the computation could be in the ice model as long as you give the ice model the ocean surface layer thickness prior.

The numerical instability issue might be a problem, limiting the coupling time step, but if all the parameters are constant, I don't think that would affect the choice of where the calculation is performed, since the ice model is already computing the drag term (or is that not true in your configuration). Is there something else that I'm missing?

 No

I'm thinking that this alternative ice-ocean drag calculation would be a nice option to provide the rest of the community, even if you continue to compute it in the ocean model. If we provide it as an option in Icepack, then the namelist setup would make sense.

For group to decide


@eclare108213 answers:

Thanks François, this is helpful. So it sounds like the choice is between the ice model providing its ice roughness to the ocean component, and the ocean model providing its top layer thickness to the ice component? Ice roughness could be ice-state dependent, so from my biased perspective, it makes more sense to have this parameterization in the ice model...
:)

Philippe has been clear that implementing this in the ice code is not part of his workplan. What are the time constraints on getting these changes made on your end? I'm wondering whether we have any flexibility to move your code into the ice component and test it, or if this is one of those things that needs to be done yesterday and fixed later (with the risk that the fix doesn't happen).


I answer:

Thanks Elizabeth for reaching out and François for providing some
context.

I separately came to the same conclusion this morning, that there is no
reason for 'dragio' not to be computed in the ice model, I think it
would make more sense. I did not think outside the box enough last week
and was too much focused on "I need this change so that it works
eaxctly like it does now in our setup".

[...] In the meantime I'll re-work by PR so that 'dragio' can be computed
within Icepack (as an option).


@phil-blain phil-blain force-pushed the under-ice-roughness branch 3 times, most recently from de62844 to fddf286 Compare June 8, 2021 18:35
@phil-blain phil-blain changed the title icepack_parameters: add under-ice roughness length 'iceruf_ocn' icepack_parameters: optionally compute 'dragio' from under-ice roughness Jun 8, 2021
@phil-blain
Copy link
Member Author

I've implemented the computation of dragio from two new parameters, iceruf_ocn and thickness_ocn, the thickness of the first ocean level (names subject to change we don't like them!). I've also changed the PR title to reflect that.

I did not add any tests for this new feature, and did not add the two new parameters to the Icepack namelist, because I'm not 100% sure they would be required / desired. For our needs we only need iceruf_ocn in the namelist, but only in CICE, not Icepack. So this PR has kind of the "minimally required" changes, but I could add more if we decide so.

So I'd like feedback on the testing front: do we want to add a test in Icepack, or it makes more sense to only add a test in CICE? Do we want to do a longer run ? etc.

@eclare108213
Copy link
Contributor

I've implemented the computation of dragio from two new parameters, iceruf_ocn and thickness_ocn, the thickness of the first ocean level (names subject to change we don't like them!). I've also changed the PR title to reflect that.

How about thickness_ocn_sfc or thickness_ocn_layer1? just to be a little more specific.

I did not add any tests for this new feature, and did not add the two new parameters to the Icepack namelist, because I'm not 100% sure they would be required / desired. For our needs we only need iceruf_ocn in the namelist, but only in CICE, not Icepack. So this PR has kind of the "minimally required" changes, but I could add more if we decide so.

I think it's helpful to have parameters like this in the namelist. It's not clear to me how other modeling centers are using namelist versus directly inserting parameter values. I didn't realize that some of you do both. That's a discussion we should have, but not for this PR.

So I'd like feedback on the testing front: do we want to add a test in Icepack, or it makes more sense to only add a test in CICE? Do we want to do a longer run ? etc.

The answers to all of those questions are "yes". It's more important to have a test in Icepack than in CICE, in my opinion, but would be good in both places. A longer run would also be helpful to make sure that it's behaving as expected compared with Francois's paper.

@phil-blain phil-blain force-pushed the under-ice-roughness branch from fddf286 to 2e761e5 Compare June 9, 2021 19:46
@phil-blain
Copy link
Member Author

Thanks for the quick turnaround on feedback, very much appreciated. I've changed thickness_ocn to thickness_ocn_layer1, thanks for the suggestion. I've also added calc_dragio to the Icepack namelist and added a test with that options turned on. I've run the base suite and will update the PR description with the results (all pass, as expected).

@phil-blain phil-blain force-pushed the under-ice-roughness branch 2 times, most recently from 279f5b5 to 8c0bcc0 Compare June 9, 2021 22:31
Copy link
Contributor

@eclare108213 eclare108213 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm much happier with this approach. Have you tried a longer run with calc_dragio=T (either coupled or not) to see if it behaves as expected?

Is iceruf_ocn=0.03 the value that makes dragio equal (or close) to the default value? If not, what value does? (I'm curious)

@phil-blain
Copy link
Member Author

I did not yet do a long run since I was waiting on this PR and the upcoming CICE one to be merged before I rebase my stuff. But I can certainly do that.

The 0.03 value for iceruf_ocn is the one we were using in our systems. I'm not sure exactly how we came to choose it... I'll check. So no, it is not chosen to make dragio close to the default value.

In some applications, we want the ice-ocean drag coefficient 'dragio'
to be computed using an under-ice roughness length and the thickness of
the first ocean level instead of imposing its value directly.

Add a new boolean flag 'calc_dragio' and real parameters 'iceruf_ocn',
'thickness_ocn_layer1', where 'iceruf_ocn' is the under-ice roughness
length, and 'thickness_ocn_layer1' is the thickness of the first ocean
level, to module 'icepack_parameters'.

With 'calc_dragio = .true.', compute the ice-ocean drag coefficient
following [1]. Update the documentation accordingly.

[1] https://doi.org/10.1002/2014JC010677
Add the 'calc_dragio' flag introduced in the previous commit to the
Icepack namelist.

Add a test for the new feature in the base_suite.
@phil-blain phil-blain force-pushed the under-ice-roughness branch from 8c0bcc0 to 7839c16 Compare June 17, 2021 15:46
@phil-blain
Copy link
Member Author

phil-blain commented Jun 17, 2021

The 0.03 value for iceruf_ocn is referenced in Roy et al 2015 and comes from

Maykut, G. A., and M. G. McPhee (1995), Solar heating of the Arctic mixed layer,J. Geophys. Res.,100, 24,691–24,703

I've rebased my stuff and I'm starting a one year coupled run with the changes from this branch just to check everything behaves as it should.

@phil-blain
Copy link
Member Author

I've ran a 2 months coupled run with `calc_dragio=.true.̀ and the new staggered atmospheric level, a configuration which is very close to what we are running for CICE4, and compared with a CICE4 run. Everything seems to behave as expected. So we are good to merge I would say.

Copy link
Contributor

@eclare108213 eclare108213 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Excellent, @phil-blain, thank you for all the extra work on this one. This is great!

@@ -1609,6 +1633,8 @@ subroutine icepack_recompute_constants()

!autodocument_end

real (kind=dbl_kind) :: lambda
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is lambda just a temporary variable, or does it have a physical meaning? units?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is unit-less. I followed our in-house implementation, which in turn follows Roy et al 2015. Lambda is defined in equation (6) (top of p. 7).

@eclare108213 eclare108213 merged commit a80472b into CICE-Consortium:master Jun 22, 2021
@phil-blain
Copy link
Member Author

Thanks @eclare108213 🎊

phil-blain added a commit to phil-blain/CICE that referenced this pull request Jun 22, 2021
In CICE-Consortium/Icepack@a80472b (icepack_parameters: optionally
compute 'dragio' from under-ice roughness (CICE-Consortium/Icepack#366),
2021-06-22), Icepack was updated to optionally compute the ice-ocean
drag coefficicent 'dragio' using an under-ice roughness length and the
thickness of the first ocean level.

Leverage this new feature in CICE by adding 'calc_dragio' and
'iceruf_ocn' to the CICE namelist. Add the new variables to the index in
the documentation and add a test with the new feature (using default
values for 'iceruf_ocn' and 'thickness_ocn_layer1').

As this new feature will mostly be useful in a coupled context, we do
not add 'thickness_ocn_layer1' to the namelist as it is expected that
the ocean model will pass this information to CICE.
apcraig pushed a commit to CICE-Consortium/CICE that referenced this pull request Jun 24, 2021
…#612)

* icepack: optionally compute 'dragio' using under-ice roughness length

In CICE-Consortium/Icepack@a80472b (icepack_parameters: optionally
compute 'dragio' from under-ice roughness (CICE-Consortium/Icepack#366),
2021-06-22), Icepack was updated to optionally compute the ice-ocean
drag coefficicent 'dragio' using an under-ice roughness length and the
thickness of the first ocean level.

Leverage this new feature in CICE by adding 'calc_dragio' and
'iceruf_ocn' to the CICE namelist. Add the new variables to the index in
the documentation and add a test with the new feature (using default
values for 'iceruf_ocn' and 'thickness_ocn_layer1').

As this new feature will mostly be useful in a coupled context, we do
not add 'thickness_ocn_layer1' to the namelist as it is expected that
the ocean model will pass this information to CICE.

* ice_grid: set 'thickness_ocn_layer1' if using 'calc_dragio'

In the previous commit we updated Icepack to allow computing the
ice-ocean drag coefficient 'dragio' using an under-ice roughness length
and the thickness of the first ocean layer, 'thickness_ocean_layer1', a
new Icepack parameter.

In some situations, we have access in CICE to the thicknesses of the
ocean levels, either hard-coded (use_bathymetry = false,
bathymetry_format = default), read from a file (use_bathymetry = true,
bathymetry_format = pop), or generated from the kmt_file
(use_bathymetry = false, bathymetry_format = pop). In these situations,
for consistency set 'thickness_ocean_layer1' in Icepack to the thickness
of the first ocean level, 'thick(1)' if 'calc_dragio' is active.
@phil-blain phil-blain deleted the under-ice-roughness branch February 13, 2024 16:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants