Add probabilistic grounding scheme for landfast ice #565

JFLemieux73 · 2021-02-16T16:15:25Z

For detailed information about submitting Pull Requests (PRs) to the CICE-Consortium,
please refer to: https://github.com/CICE-Consortium/About-Us/wiki/Resource-Index#information-for-developers

PR checklist

We think that the expression "seabed stress" is more appropriate than "basal stress" We have therefore changed the code and the documentation.

…iles)

… prep

… by seabedstress

eclare108213

This looks okay to me, just reading through the code changes. There are a lot of potential divides-by-zero or logarithm issues. How extensively has this been tested?

Shouldn't the new flag kseabed be added to the three tests for which seabedstress=T? And at least one of them be set for the new parameterization?

JFLemieux73 · 2021-02-17T13:35:36Z

Hi @eclare108213 ,

I don't think we should have problems with division by zero...as long as there is some sea ice. The calculations are only done if the total concentration is greater than 0.05 [ if (atot .gt. 0.05_dbl_kind) ], We also have (max(acat(n), puny)) where there could be a div by zero.

What do you mean by logarithm issues?

In terms of tests, we ran 1 year of gx3 and 1 year of gx1. Using our in-house code (cice5 with modifications including the same seabed2 algorithm) we ran many 10 year hindcast simulations (resolution of 0.25 degree). No crash.

I could add a new test in another PR.

apcraig

Should kseabed be a character string rather than "1" and "2"? I think we'd like to work towards that implementation when we can.

There should be an abort somewhere if kseabed is not set to a valid option.

The value of kseabed should be written out in ice_init.F90 as a readable character string. (probably only when seabedstress is true)

kseabed needs to be added to the namelist documention.

apcraig · 2021-02-17T18:06:22Z

I think we should also add a test with kseabed=2 assuming it's ready to use. I can provide some guidance on what and how to do this if needed.

JFLemieux73 · 2021-02-17T18:58:57Z

Thanks @apcraig.

kseabed is inspired from kstrength, kdyn, ktherm...I will think about good character strings to replace these.

I will do the other changes requested.

Just to clarify...If I just push the same branch, will this update the PR?

apcraig · 2021-02-17T19:01:23Z

@JFLemieux73, You can push to the same branch. The PR will update and all the checks will automatically rerun including documentation generation. It's very typical that a PR is updated as part of the draft and review process in this way. Thanks.

phil-blain

I noted a few things below. One general comment about naming: in the spirit of writing self-documenting code, I would prefer mean_thickness_seabed_stress_factor and probabilistic_seabed_stress_factor over seabed1_stress_factor and seabed2_stress_factor.

phil-blain · 2021-02-17T18:20:41Z

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

+          case (1)
+             call seabed1_stress_factor (nx_block,         ny_block,       &
+                                        icellu  (iblk),                   &
+                                        indxui(:,iblk),   indxuj(:,iblk), &
+                                        vice(:,:,iblk),   aice(:,:,iblk), &
+                                        hwater(:,:,iblk), Tbu(:,:,iblk))
+          case (2)
+             call seabed2_stress_factor (nx_block,         ny_block,                   &
+                                        icellt(iblk), indxti(:,iblk), indxtj(:,iblk), &
+                                        icellu(iblk), indxui(:,iblk), indxuj(:,iblk), &
+                                        aicen(:,:,:,iblk), vicen(:,:,:,iblk),         &
+                                        hwater(:,:,iblk), Tbu(:,:,iblk))


Here and for seabed2 below, the first ampersand is not aligned with those below:

Suggested change

case (1)

call seabed1_stress_factor (nx_block, ny_block, &

icellu (iblk), &

indxui(:,iblk), indxuj(:,iblk), &

vice(:,:,iblk), aice(:,:,iblk), &

hwater(:,:,iblk), Tbu(:,:,iblk))

case (2)

call seabed2_stress_factor (nx_block, ny_block, &

icellt(iblk), indxti(:,iblk), indxtj(:,iblk), &

icellu(iblk), indxui(:,iblk), indxuj(:,iblk), &

aicen(:,:,:,iblk), vicen(:,:,:,iblk), &

hwater(:,:,iblk), Tbu(:,:,iblk))

case (1)

call seabed1_stress_factor (nx_block, ny_block, &

icellu (iblk), &

indxui(:,iblk), indxuj(:,iblk), &

vice(:,:,iblk), aice(:,:,iblk), &

hwater(:,:,iblk), Tbu(:,:,iblk))

case (2)

call seabed2_stress_factor (nx_block, ny_block, &

icellt(iblk), indxti(:,iblk), indxtj(:,iblk), &

icellu(iblk), indxui(:,iblk), indxuj(:,iblk), &

aicen(:,:,:,iblk), vicen(:,:,:,iblk), &

hwater(:,:,iblk), Tbu(:,:,iblk))

phil-blain · 2021-02-17T18:28:00Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+      call icepack_query_parameters(puny_out=puny)
+
+      Tbt=c0
+      sigma_b = 2.5d0 ! Standard deviation of bathymetry. 


this does not change in the subroutine so could be a parameter

phil-blain · 2021-02-17T18:32:32Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+
+         atot = sum(aicen(i,j,1:ncat))
+
+         if (atot .gt. 0.05_dbl_kind .and. hwater(i,j) .lt. max_depth) then


I think we prefer to use < and > in new code instead of .lt. and .gt.

phil-blain · 2021-02-17T19:08:08Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+            x_k = (/(wid_i*(i*c1-0.5d0), i=1, ncat_i)/) 
+            y_n = (/( (mu_b-c3*sigma_b)+(i*c1-0.5d0)*(c2*c3*sigma_b/ncat_b), i=1, ncat_b )/)


I think these formula would be more legible with a little more whitespace

phil-blain · 2021-02-17T19:10:07Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+!            x95  = exp(mu_i + sqrt(c2*sigma_i)*1.1631d0)
+!            x99  = exp(mu_i + sqrt(c2*sigma_i)*1.6450d0)
+!            x995 = exp(mu_i + sqrt(c2*sigma_i)*1.8214d0)
+!            x996 = exp(mu_i + sqrt(c2*sigma_i)*1.8753d0)
+            x997 = exp(mu_i + sqrt(c2*sigma_i)*1.9430d0)
+!            x998 = exp(mu_i + sqrt(c2*sigma_i)*2.0352d0)
+!            x999 = exp(mu_i + sqrt(c2*sigma_i)*2.1851d0)
+!            xinf = exp(mu_i + sqrt(c2*sigma_i)*4.4981d0) ! 0.9999999999


I don't think we want to keep those commented lines. If we want to keep the ability to tweak the factor in the square root more easily, why not declare a local parameter variable ?

I agree it's better not to have commented lines like this if they aren't ever expected to be used again. Presumably Dupont et al documents the procedure and the values that were tried, so this information shouldn't need to be kept in the code.

phil-blain · 2021-02-17T19:11:41Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+            do n =1, ncat_i
+               if (x_k(n) .gt. x_kmax) P_x(n)=c0
+            enddo


minor: I think maybe this loop could be rewritten with where statement

and > instead of .gt.

and below, < for .le., == for .eq.

phil-blain · 2021-02-17T19:16:58Z

cicecore/drivers/direct/hadgem3/CICE_InitMod.F90

@@ -71,7 +71,7 @@ subroutine cice_init
      use ice_domain, only: init_domain_blocks
      use ice_domain_size, only: ncat, nfsd
      use ice_dyn_eap, only: init_eap, alloc_dyn_eap
-      use ice_dyn_shared, only: kdyn, init_dyn, basalstress, alloc_dyn_shared
+      use ice_dyn_shared, only: kdyn, init_dyn, seabedstress, alloc_dyn_shared


This variable is not used in that driver, so we might as well remove it from that line:

Suggested change

use ice_dyn_shared, only: kdyn, init_dyn, seabedstress, alloc_dyn_shared

use ice_dyn_shared, only: kdyn, init_dyn, alloc_dyn_shared

phil-blain · 2021-02-17T19:18:41Z

doc/source/science_guide/sg_dynamics.rst

@@ -236,8 +236,31 @@ pending further testing.
 Seabed stress
 ***************

-The parameterization for the seabed stress is described in :cite:`Lemieux16`. The components of the basal seabed stress are 
-:math:`\tau_{bx}=C_bu` and :math:`\tau_{by}=C_bv`, where :math:`C_b` is a coefficient expressed as
+CICE now includes two options for calculating the seabed stress,


minor: I would not use words such as "now" (or "new", a little bit below), since it is not clear to the reader when is "now". The documentation reflects the actual state of the code, so just "CICE includes two options..." would be better in my opinion.

phil-blain · 2021-02-17T19:21:04Z

doc/source/science_guide/sg_dynamics.rst

+on the probability of contact between the ice thickness distribution
+(ITD) and the seabed. Multi-thickness category models such as CICE typically use a
+few thickness categories (5-10). This crude representation of the ITD
+does not resolve the tail of the ITD which is crucial for grounding


Suggested change

does not resolve the tail of the ITD which is crucial for grounding

does not resolve the tail of the ITD, which is crucial for grounding

@phil-blain you seem to have good ideas...Right now we choose the method with the integer kseabed. @apcraig would prefer character strings...any suggestions (should be short I guess)?

I think mean_thickness and probabilistic would be fit well.

Note: you should be able to click the "Accept suggestion" button on the suggestions I made above, and then they will be committed to your branch (on GitHub). You will then be able to pull these changes locally if you want to tweak the PR further.

doc/source/user_guide/ug_testing.rst

Co-authored-by: Philippe Blain <levraiphilippeblain@gmail.com>

JFLemieux73 · 2021-02-18T19:19:06Z

I reran the test suite with the latest changes:

269 measured results of 269 total results
269 of 269 tests PASSED
0 of 269 tests PENDING
0 of 269 tests MISSING data
0 of 269 tests FAILED

I am almost done...

JFLemieux73 · 2021-02-18T21:46:35Z

@apcraig I have made the changes requested. The only thing missing is an additional test for the new seabed stress method. Is it ok if we do it in another PR?

eclare108213

I have a few suggestions for improvements, nothing major. A test for the new parameterization can be added in a separate PR, although it would be nice to have it in this one.

eclare108213 · 2021-02-18T22:05:03Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+!            x95  = exp(mu_i + sqrt(c2*sigma_i)*1.1631d0)
+!            x99  = exp(mu_i + sqrt(c2*sigma_i)*1.6450d0)
+!            x995 = exp(mu_i + sqrt(c2*sigma_i)*1.8214d0)
+!            x996 = exp(mu_i + sqrt(c2*sigma_i)*1.8753d0)
+            x997 = exp(mu_i + sqrt(c2*sigma_i)*1.9430d0)
+!            x998 = exp(mu_i + sqrt(c2*sigma_i)*2.0352d0)
+!            x999 = exp(mu_i + sqrt(c2*sigma_i)*2.1851d0)
+!            xinf = exp(mu_i + sqrt(c2*sigma_i)*4.4981d0) ! 0.9999999999


I agree it's better not to have commented lines like this if they aren't ever expected to be used again. Presumably Dupont et al documents the procedure and the values that were tried, so this information shouldn't need to be kept in the code.

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

doc/source/science_guide/sg_dynamics.rst

eclare108213 · 2021-02-18T22:26:38Z

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

+            do n =1, ncat_i
+               if (x_k(n) .gt. x_kmax) P_x(n)=c0
+            enddo


and > instead of .gt.

and below, < for .le., == for .eq.

…4 months (gx1)

JFLemieux73 · 2021-02-19T21:37:29Z

Ok I have modified the routine following the comments of @eclare108213.

Just to be safe I reran the test suite:

269 measured results of 269 total results
269 of 269 tests PASSED
0 of 269 tests PENDING
0 of 269 tests MISSING data
0 of 269 tests FAILED

Is it ok?

eclare108213 · 2021-02-19T21:47:18Z

Looks good to me. Let's let @apcraig update his review. Thanks JF!

apcraig

I can add a test separately or we can do it now. Is the "probabilistic" option ready to test? Right now alt01, alt03, and boxrestore seem to have seabed_stress set to true. I assume that means they are all using the LKD method? How about if we add "seabed_stress_method='probabilistic'" to alt01 or alt03? That will change the answers there, but that's OK. Is there any reason why that might be a problem or an unreasonable way to proceed?

apcraig · 2021-02-22T18:17:13Z

One other thing, seabed_stress_method still needs to be documented in the namelist section (doc/source/case_settings/ug_case_settings.rst). Also, the seabedstress documentation in the same section needs to be moved down in the appropriate location for alphabetization.

JFLemieux73 · 2021-02-22T18:28:28Z

Hi @apcraig , Yes "probabilistic" is ready to test. Replacing LKD by "seabed_stress_method='probabilistic'" for alt01 or alt03 is a good plan. I could do it. I guess I just have to modify the set_nml.alt file?

apcraig · 2021-02-22T18:51:32Z

Exactly, just modify the appropriate set_nml file. You should be able to run a quick test with the updated file leveraging -s to specify the option. If you have questions, let me know. thanks!

JFLemieux73 · 2021-02-22T19:36:54Z

Ok the new test seems to work.

apcraig

Looks good, thanks JF. I will merge this afternoon unless there are other comments.

apcraig · 2021-02-24T23:01:20Z

I just ran a full test suite. The changes to the seabed_stress_method from LKD to probabilistic in the alt01 case did not change any answers.

I think this suggests we are not getting any grounding in the alt01 short test runs. Does that sound plausible? It also suggests we may not be getting any grounding in any short runs, so we're not actually testing the grounding scheme in our tests even when seabed_stress = .true.

Any thoughts?

JFLemieux73 · 2021-02-25T13:46:21Z

I think the ice is too thin to get grounding. As the runs are short it does not have time to get thicker.

Just to be clear, these simulations start on 1 January and run just for a few days, right?

If you look here:
#559
you can see that there is grounding in April (after 3 months). The easiest way to know if there is grounding is to look at the taubx and tauby fields.

I could investigate what is going on in the alt01 test.

eclare108213 · 2021-02-25T16:09:13Z

This is a good reason to start from new restart files that have included the landfast ice parameterization.

JFLemieux73 · 2021-02-26T14:17:53Z

I setup the test:
./cice.setup -c test_alt01 -g gx1 -m daley -e intel -p 8x1 -s alt01

taubx and tauby are indeed zero.

JFLemieux73 · 2021-02-26T14:19:30Z

I think the best way to test both LKD and probabilistic schemes would be to restart in April. The ice would then be thick enough to get grounding.

JFLemieux73 · 2021-03-09T18:48:53Z

@apcraig I would like to submit a PR with minor modifications to the documentation and ice_in. Should I work at the same time on improving our tests for grounding (i.e. restart in April) or you had planned to do that?

apcraig · 2021-03-09T19:31:06Z

@JFLemieux73, please feel free to go ahead and submit a PR now. I think we should wait to add the new "April" tests until we have the new ics from JRA55 as well as the new time manager. Let me create an issue for this so we don't forget.

JFLemieux73 added 18 commits February 15, 2021 15:08

In the process of adding seabed2 parameterization

8f49435

Added calls foe icepack_query_parameters

73cdca4

Modified modules used in subroutine to match latest code (it now comp…

d7a9eec

…iles)

Changed expression basal stress to seabed stress

3a4c2b4

Replaced basal by better expression seabed at many places

e0aaf74

seabed everywhere...basalstress logical replaced by seabedstress

c9222d0

Added kseabed in namelist for choosing seabed stress method

1f6f09a

Changed gacc to gravit (from icepack) and use of alphab from ice_in

41ef64e

Cleaned up the seabed2 code, added comments

c99a3df

Changing the doc for new grounding scheme

001914f

Almost done with the doc

257e7ac

Minor changes to the doc

774d2b1

Modifs to the doc

8f320e9

Value of sigmab was not the right one with respect to Dupont et el in…

a8bcbed

… prep

For some rebase did not work completely...replced basalstress logical…

c6426b3

… by seabedstress

Other issue with the rebase for ice_init...fixed

4bfe465

Minor mistake corrected...it compiles

16e9b1d

changed basalstress to seabedstress in options for tests

cdc82d1

eclare108213 requested review from eclare108213, phil-blain and apcraig February 16, 2021 21:40

eclare108213 added Dynamics Type: Feature labels Feb 16, 2021

eclare108213 reviewed Feb 16, 2021

View reviewed changes

apcraig requested changes Feb 17, 2021

View reviewed changes

phil-blain reviewed Feb 17, 2021

View reviewed changes

phil-blain reviewed Feb 18, 2021

View reviewed changes

doc/source/user_guide/ug_testing.rst Outdated Show resolved Hide resolved

Update doc/source/user_guide/ug_testing.rst

9f689ee

Co-authored-by: Philippe Blain <levraiphilippeblain@gmail.com>

Cosmetic changes to ice_init to align prints

852658e

eclare108213 approved these changes Feb 18, 2021

View reviewed changes

JFLemieux73 added 2 commits February 19, 2021 14:22

Some modifs requested for PR: same taub after 4 months (gx1)

b6bb0a4

Done with requested changes to new grounding method: same taub after …

0a930db

…4 months (gx1)

apcraig reviewed Feb 22, 2021

View reviewed changes

Modifs to doc and new test of prob method

173d99e

apcraig approved these changes Feb 22, 2021

View reviewed changes

apcraig merged commit 3c516c8 into CICE-Consortium:master Feb 22, 2021

JFLemieux73 deleted the Seabed_stress2 branch February 25, 2021 13:47

apcraig mentioned this pull request Feb 27, 2021

Implement new grounding parameterization for landfast ice #559

Closed

apcraig mentioned this pull request Mar 9, 2021

Add new landfast tests #573

Closed


		atot = sum(aicen(i,j,1:ncat))

		if (atot .gt. 0.05_dbl_kind .and. hwater(i,j) .lt. max_depth) then

		x_k = (/(wid_i(ic1-0.5d0), i=1, ncat_i)/)
		y_n = (/( (mu_b-c3sigma_b)+(ic1-0.5d0)(c2c3*sigma_b/ncat_b), i=1, ncat_b )/)

	use ice_dyn_shared, only: kdyn, init_dyn, seabedstress, alloc_dyn_shared
	use ice_dyn_shared, only: kdyn, init_dyn, alloc_dyn_shared

	does not resolve the tail of the ITD which is crucial for grounding
	does not resolve the tail of the ITD, which is crucial for grounding

Add probabilistic grounding scheme for landfast ice #565

Add probabilistic grounding scheme for landfast ice #565

Conversation

JFLemieux73 commented Feb 16, 2021

PR checklist

eclare108213 left a comment

Choose a reason for hiding this comment

JFLemieux73 commented Feb 17, 2021 • edited Loading

apcraig left a comment

Choose a reason for hiding this comment

apcraig commented Feb 17, 2021

JFLemieux73 commented Feb 17, 2021

apcraig commented Feb 17, 2021

phil-blain left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

phil-blain Feb 17, 2021 • edited Loading

Choose a reason for hiding this comment

JFLemieux73 commented Feb 18, 2021

JFLemieux73 commented Feb 18, 2021

eclare108213 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

JFLemieux73 commented Feb 19, 2021

eclare108213 commented Feb 19, 2021

apcraig left a comment

Choose a reason for hiding this comment

apcraig commented Feb 22, 2021

JFLemieux73 commented Feb 22, 2021

apcraig commented Feb 22, 2021

JFLemieux73 commented Feb 22, 2021

apcraig left a comment

Choose a reason for hiding this comment

apcraig commented Feb 24, 2021

JFLemieux73 commented Feb 25, 2021

eclare108213 commented Feb 25, 2021

JFLemieux73 commented Feb 26, 2021

JFLemieux73 commented Feb 26, 2021

JFLemieux73 commented Mar 9, 2021

apcraig commented Mar 9, 2021

JFLemieux73 commented Feb 17, 2021 •

edited

Loading

phil-blain Feb 17, 2021 •

edited

Loading