Merge remote-tracking branch 'Consortium/main' into feature/updcice

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)
 

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

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

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)

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))

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)

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)}",

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.

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."