From 28fdbebf4a3b98485c1e977344fd402f29f6a0b6 Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Thu, 19 Jan 2023 14:59:26 -0500 Subject: [PATCH 1/3] Fix for rare instability in (probabilistic) seabed stress (#810) * Modified doc for Dupont et al ref * Minor modif to v_i calculation in seabed_stress_factor_prob to prevent rare instability --- cicecore/cicedyn/dynamics/ice_dyn_shared.F90 | 11 ++++++----- cicecore/cicedyn/general/ice_init.F90 | 2 +- doc/source/master_list.bib | 10 ++++++++++ doc/source/science_guide/sg_dynamics.rst | 2 +- doc/source/user_guide/ug_case_settings.rst | 2 +- 5 files changed, 19 insertions(+), 8 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 index a12e6fddd..50f1aae6e 100644 --- a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 @@ -59,7 +59,7 @@ module ice_dyn_shared yield_curve , & ! 'ellipse' ('teardrop' needs further testing) visc_method , & ! method for viscosity calc at U points (C, CD grids) seabed_stress_method ! method for seabed stress calculation - ! LKD: Lemieux et al. 2015, probabilistic: Dupont et al. in prep. + ! LKD: Lemieux et al. 2015, probabilistic: Dupont et al. 2022 real (kind=dbl_kind), parameter, public :: & u0 = 5e-5_dbl_kind, & ! residual velocity for seabed stress (m/s) @@ -1347,8 +1347,9 @@ end subroutine seabed_stress_factor_LKD ! a normal distribution with sigma_b = 2.5d0. An improvement would ! be to provide the distribution based on high resolution data. ! -! Dupont, F. Dumont, D., Lemieux, J.F., Dumas-Lefebvre, E., Caya, A. -! in prep. +! Dupont, F., D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, A. Caya (2022). +! A probabilistic seabed-ice keel interaction model, The Cryosphere, 16, +! 1963-1977. ! ! authors: D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, F. Dupont ! @@ -1481,13 +1482,13 @@ subroutine seabed_stress_factor_prob (nx_block, ny_block, & do n =1, ncat v_i = v_i + vcat(n)**2 / (max(acat(n), puny)) enddo - v_i = v_i - m_i**2 + v_i = max((v_i - m_i**2), puny) mu_i = log(m_i/sqrt(c1 + v_i/m_i**2)) ! parameters for the log-normal sigma_i = sqrt(log(c1 + v_i/m_i**2)) ! max thickness associated with percentile of log-normal PDF - ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al.) + ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al. 2022) x_kmax = exp(mu_i + sqrt(c2*sigma_i)*1.9430d0) diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index d56ad002e..4c8fb1fee 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -392,7 +392,7 @@ subroutine input_data dyscale = 1.0_dbl_kind ! user defined rectgrid y-grid scale factor (e.g., 1.02) close_boundaries = .false. ! true = set land on edges of grid seabed_stress= .false. ! if true, seabed stress for landfast is on - seabed_stress_method = 'LKD'! LKD = Lemieux et al 2015, probabilistic = Dupont et al. in prep + seabed_stress_method = 'LKD'! LKD = Lemieux et al 2015, probabilistic = Dupont et al. 2022 k1 = 7.5_dbl_kind ! 1st free parameter for landfast parameterization k2 = 15.0_dbl_kind ! 2nd free parameter (N/m^3) for landfast parametrization alphab = 20.0_dbl_kind ! alphab=Cb factor in Lemieux et al 2015 diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index a2da9b9f8..9e387efb9 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -1082,6 +1082,16 @@ @article{Bouchat22 year = {2022} } +@Article{Dupont22, + author = {F. Dupont and D. Dumont and J.F. Lemieux and E. Dumas-Lefebvre and A. Caya}, + title = "{A probabilistic seabed-ice keel interaction model}", + journal = TC, + year = {2022}, + volume = {16}, + pages = {1963-1977}, + url = {https://doi.org/10.5194/tc-16-1963-2022} +} + @Article{Tsujino18, author = "H. Tsujino and S. Urakawa and R.J. Small and W.M. Kim and S.G. Yeager and et al.", title = "{JRA‐55 based surface dataset for driving ocean–sea‐ice models (JRA55‐do)}", diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index e6b918538..6b269d453 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -382,7 +382,7 @@ The value of :math:`k_1` can be changed at runtime using the namelist variable ` This more sophisticated grounding parameterization computes the seabed stress based on the probability of contact between the ice thickness distribution -(ITD) and the seabed. Multi-thickness category models such as CICE typically use a +(ITD) and the seabed :cite:`Dupont22`. 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 events. diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 587adcd56..9906fba87 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -485,7 +485,7 @@ dynamics_nml "``revised_evp``", "logical", "use revised EVP formulation", "``.false.``" "``seabed_stress``", "logical", "use seabed stress parameterization for landfast ice", "``.false.``" "``seabed_stress_method``", "``LKD``", "linear keel draft method :cite:`Lemieux16`", "``LKD``" - "", "``probabilistic``", "probability of contact method (Dupont et al., in prep)", "" + "", "``probabilistic``", "probability of contact method :cite:`Dupont22`", "" "``ssh_stress``", "``coupled``", "computed from coupled sea surface height gradient", "``geostrophic``" "", "``geostropic``", "computed from ocean velocity", "" "``threshold_hw``", "real", "Max water depth for grounding (see :cite:`Amundrud04`)", "30." From b946a95ea46096d739b4b0d725d0eab81a53fbee Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Tue, 31 Jan 2023 20:17:26 -0700 Subject: [PATCH 2/3] Some small CESM updates. (#812) * Fix OMP setup * Update meshgrid * Small updates for CESM * Add change for UFS --- cicecore/cicedyn/infrastructure/ice_grid.F90 | 1 - cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 | 6 ++++++ cicecore/shared/ice_fileunits.F90 | 2 -- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index b775c21f2..0d56d3400 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -507,7 +507,6 @@ subroutine init_grid2 ! Diagnose OpenMP thread schedule, force order in output !----------------------------------------------------------------- -! This code does not work in CESM. Needs to be investigated further. #if defined (_OPENMP) !$OMP PARALLEL DO ORDERED PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 index 7019f7128..cdfbac87a 100644 --- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 @@ -749,6 +749,9 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, & ! if (ndim3 == ncat .and. ncat>1) then if (ndim3 == ncat .and. ndims == 3) then call pio_read_darray(File, vardesc, iodesc3d_ncat, work, status) +#ifdef CESMCOUPLED + where (work == PIO_FILL_DOUBLE) work = c0 +#endif if (present(field_loc)) then do n=1,ndim3 call ice_HaloUpdate (work(:,:,n,:), halo_info, & @@ -758,6 +761,9 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, & ! elseif (ndim3 == 1) then elseif (ndim3 == 1 .and. ndims == 2) then call pio_read_darray(File, vardesc, iodesc2d, work, status) +#ifdef CESMCOUPLED + where (work == PIO_FILL_DOUBLE) work = c0 +#endif if (present(field_loc)) then call ice_HaloUpdate (work(:,:,1,:), halo_info, & field_loc, field_type) diff --git a/cicecore/shared/ice_fileunits.F90 b/cicecore/shared/ice_fileunits.F90 index 7e425e5e7..72a40f513 100644 --- a/cicecore/shared/ice_fileunits.F90 +++ b/cicecore/shared/ice_fileunits.F90 @@ -81,10 +81,8 @@ module ice_fileunits integer (kind=int_kind), public :: & nu_diag = ice_stdout ! diagnostics output file, unit number may be overwritten -#ifdef CESMCOUPLED logical (kind=log_kind), public :: & nu_diag_set = .false. ! flag to indicate whether nu_diag is already set -#endif integer (kind=int_kind), public :: & ice_IOUnitsMinUnit = 11, & ! do not use unit numbers below From d73bb8b064e062217afc54e33faa9a0247be50e3 Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Thu, 2 Mar 2023 15:15:46 -0700 Subject: [PATCH 3/3] Add time_period_freq (#816) * Add time_period_freq to history file metadata --- .../io/io_netcdf/ice_history_write.F90 | 22 ++++++++++++++++++- .../io/io_pio2/ice_history_write.F90 | 18 +++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 index 019ab8ce9..d85ec5e3c 100644 --- a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 @@ -48,7 +48,7 @@ subroutine ice_write_hist (ns) use ice_blocks, only: nx_block, ny_block use ice_broadcast, only: broadcast_scalar use ice_calendar, only: msec, timesecs, idate, idate0, write_ic, & - histfreq, days_per_year, use_leap_years, dayyr, & + histfreq, histfreq_n, days_per_year, use_leap_years, dayyr, & hh_init, mm_init, ss_init use ice_communicate, only: my_task, master_task use ice_domain, only: distrb_info @@ -86,6 +86,7 @@ subroutine ice_write_hist (ns) integer (kind=int_kind), dimension(6) :: dimidex real (kind=dbl_kind) :: ltime2 character (char_len) :: title + character (char_len) :: time_period_freq = 'none' character (char_len_long) :: ncfile(max_nstrm) real (kind=dbl_kind) :: secday, rad_to_deg @@ -682,6 +683,25 @@ subroutine ice_write_hist (ns) if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: global attribute date2') + select case (histfreq(ns)) + case ("y", "Y") + write(time_period_freq,'(a,i0)') 'year_',histfreq_n(ns) + case ("m", "M") + write(time_period_freq,'(a,i0)') 'month_',histfreq_n(ns) + case ("d", "D") + write(time_period_freq,'(a,i0)') 'day_',histfreq_n(ns) + case ("h", "H") + write(time_period_freq,'(a,i0)') 'hour_',histfreq_n(ns) + case ("1") + write(time_period_freq,'(a,i0)') 'step_',histfreq_n(ns) + end select + + if (.not.write_ic .and. trim(time_period_freq) /= 'none') then + status = nf90_put_att(ncid,nf90_global,'time_period_freq',trim(time_period_freq)) + if (status /= nf90_noerr) call abort_ice(subname// & + 'ERROR: global attribute time_period_freq') + endif + title = 'CF-1.0' status = & nf90_put_att(ncid,nf90_global,'conventions',title) diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 index 6407d8c76..a697a98d5 100644 --- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 @@ -76,6 +76,7 @@ subroutine ice_write_hist (ns) integer (kind=int_kind), dimension(6) :: dimidex real (kind= dbl_kind) :: ltime2 character (char_len) :: title + character (char_len) :: time_period_freq = 'none' character (char_len_long) :: ncfile(max_nstrm) integer (kind=int_kind) :: iotype @@ -649,6 +650,23 @@ subroutine ice_write_hist (ns) write(title,'(a,i6)') 'seconds elapsed into model date: ',msec status = pio_put_att(File,pio_global,'comment3',trim(title)) + select case (histfreq(ns)) + case ("y", "Y") + write(time_period_freq,'(a,i0)') 'year_',histfreq_n(ns) + case ("m", "M") + write(time_period_freq,'(a,i0)') 'month_',histfreq_n(ns) + case ("d", "D") + write(time_period_freq,'(a,i0)') 'day_',histfreq_n(ns) + case ("h", "H") + write(time_period_freq,'(a,i0)') 'hour_',histfreq_n(ns) + case ("1") + write(time_period_freq,'(a,i0)') 'step_',histfreq_n(ns) + end select + + if (.not.write_ic .and. trim(time_period_freq) /= 'none') then + status = pio_put_att(File,pio_global,'time_period_freq',trim(time_period_freq)) + endif + title = 'CF-1.0' status = & pio_put_att(File,pio_global,'conventions',trim(title))