Ice dynamics (#80)

* In MOM_ice_shelf_dynamics.F90 changes are  made to calc_shelf_visc() and calc_shelf_driving_stress() to account for an irregular quadrilateral grid. In MOM_ice_shelf_initialize.F90 changes are made to initialize_ice_thickness_from_file() to correct masks at initialization.

* Changed indentation

* Changes are made to 'calc_shelf_visc()` to make computations of the velocity derivatives rotation-invariant. Changes in `update_ice_shelf()` utilize G%IareaT(:,:) instead of 1/G%areaT(:,:).

* Allocate shelf_sfc_mass_flux for dynamic ice shelf configurations

* Conditional allocation of ice shelf fluxes and data override call for sfc mass flux

* get rid of excessive line length

* call data override init for ice shelves

* A new subroutine change_thickness_using_precip() changes the ice thickness in response to the atmospheric mass flux

* Conditional calls to data_override

  - Introduces a new flag: DATA_OVERRIDE_SHELF_FLUXES, which
    if set to True, and ACTIVE_SHELF_DYNAMICS is also True, will
    enable data_override capability for the surface mass deposition
    (only avaialble for MOSAIC grid types)

* Changes are made to extrapolate_metric() in MOM_grid_initialize.F90 to avoid negative values at the southern boundary of the domain. This is done for grids extending over Antarctica.

* Reversed changes in MOM_grid_initialize.F90 to debug a regression test.

* Uncommented lines reading the ice velocity from file

* Changes to scaling factors and conversions of the 'sfc_mass_flux' parameter

* Remove unnecessary diagnostics for ice shelves from MOM_forcing_type

* Correct rescaling for data override of ice shelf sfc mass fluxes

* Add diagnostic rescaling for ice shelf mass fluxes

* Uncommented assignement of 'area_shelf_h(i,j)=G%areaT(i,j)' in MOM_ice_shelf_initialize.F90 and modified the four halo updates in MOM_ice_shelf.F90.

Co-authored-by: Matthew Harrison <Matthew.Harrison@noaa.gov>

config_src/drivers/solo_driver/MOM_driver.F90

@@ -26,6 +26,7 @@ program MOM_main
 
   use MOM_cpu_clock,       only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
   use MOM_cpu_clock,       only : CLOCK_COMPONENT
+  use MOM_data_override,   only : data_override_init
   use MOM_diag_mediator,   only : enable_averaging, disable_averaging, diag_mediator_end
   use MOM_diag_mediator,   only : diag_ctrl, diag_mediator_close_registration
   use MOM,                 only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
@@ -47,6 +48,7 @@ program MOM_main
   use MOM_ice_shelf,       only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS
   use MOM_ice_shelf,       only : shelf_calc_flux, add_shelf_forces, ice_shelf_save_restart
   use MOM_ice_shelf,       only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces
+  use MOM_ice_shelf,       only : ice_shelf_query
   use MOM_interpolate,     only : time_interp_external_init
   use MOM_io,              only : file_exists, open_ASCII_file, close_file
   use MOM_io,              only : check_nml_error, io_infra_init, io_infra_end
@@ -176,6 +178,8 @@ program MOM_main
   type(surface_forcing_CS),  pointer :: surface_forcing_CSp => NULL()
   type(write_cputime_CS),    pointer :: write_CPU_CSp => NULL()
   type(ice_shelf_CS),        pointer :: ice_shelf_CSp => NULL()
+  logical                            :: override_shelf_fluxes !< If true, and shelf dynamics are active,
+                                        !! the data_override feature is enabled (only for MOSAIC grid types)
   type(wave_parameters_cs),  pointer :: waves_CSp => NULL()
   type(MOM_restart_CS),      pointer :: &
     restart_CSp => NULL()     !< A pointer to the restart control structure
@@ -302,6 +306,8 @@ program MOM_main
     ! when using an ice shelf
     call initialize_ice_shelf_fluxes(ice_shelf_CSp, grid, US, fluxes)
     call initialize_ice_shelf_forces(ice_shelf_CSp, grid, US, forces)
+    call ice_shelf_query(ice_shelf_CSp, grid, data_override_shelf_fluxes=override_shelf_fluxes)
+    if (override_shelf_fluxes) call data_override_init(Ocean_Domain_in=grid%domain%mpp_domain)
   endif
 
 

src/core/MOM_forcing_type.F90

@@ -170,6 +170,8 @@ module MOM_forcing_type
                                  !! exactly 0 away from shelves or on land.
   real, pointer, dimension(:,:) :: iceshelf_melt => NULL() !< Ice shelf melt rate (positive)
                                  !! or freezing (negative) [R Z T-1 ~> kg m-2 s-1]
+  real, pointer, dimension(:,:) :: shelf_sfc_mass_flux => NULL() !< Ice shelf surface mass flux
+                                 !! deposition from the atmosphere. [R Z T-1 ~> kg m-2 s-1]
 
   ! Scalars set by surface forcing modules
   real :: vPrecGlobalAdj = 0.     !< adjustment to restoring vprec to zero out global net [kg m-2 s-1]
@@ -377,7 +379,6 @@ module MOM_forcing_type
   ! Iceberg + Ice shelf diagnostic handles
   integer :: id_ustar_ice_cover = -1
   integer :: id_frac_ice_cover = -1
-
   ! wave forcing diagnostics handles.
   integer :: id_lamult = -1
   !>@}
@@ -2099,6 +2100,12 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
       fluxes%iceshelf_melt(i,j)  = flux_tmp%iceshelf_melt(i,j)
     enddo ; enddo
   endif
+  if (associated(fluxes%shelf_sfc_mass_flux) &
+                 .and. associated(flux_tmp%shelf_sfc_mass_flux)) then
+    do i=isd,ied ; do j=jsd,jed
+      fluxes%shelf_sfc_mass_flux(i,j)  = flux_tmp%shelf_sfc_mass_flux(i,j)
+    enddo ; enddo
+  endif
   if (associated(fluxes%frac_shelf_h) .and. associated(flux_tmp%frac_shelf_h)) then
     do i=isd,ied ; do j=jsd,jed
       fluxes%frac_shelf_h(i,j)  = flux_tmp%frac_shelf_h(i,j)
@@ -2928,7 +2935,8 @@ end subroutine forcing_diagnostics
 
 !> Conditionally allocate fields within the forcing type
 subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
-                                     shelf, iceberg, salt, fix_accum_bug, cfc, waves)
+                                  shelf, iceberg, salt, fix_accum_bug, cfc, waves, &
+                                  shelf_sfc_accumulation)
   type(ocean_grid_type), intent(in) :: G       !< Ocean grid structure
   type(forcing),      intent(inout) :: fluxes  !< A structure containing thermodynamic forcing fields
   logical, optional,     intent(in) :: water   !< If present and true, allocate water fluxes
@@ -2942,14 +2950,21 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
                                                !! accumulation of ustar_gustless
   logical, optional,     intent(in) :: cfc     !< If present and true, allocate cfc fluxes
   logical, optional,     intent(in) :: waves   !< If present and true, allocate wave fields
+  logical, optional,     intent(in) :: shelf_sfc_accumulation !< If present and true, and shelf is true,
+                                               !! then allocate surface flux deposition from the atmosphere
+                                               !! over ice shelves and ice sheets.
 
   ! Local variables
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
   logical :: heat_water
+  logical :: shelf_sfc_acc
 
   isd  = G%isd   ; ied  = G%ied    ; jsd  = G%jsd   ; jed  = G%jed
   IsdB = G%IsdB  ; IedB = G%IedB   ; JsdB = G%JsdB  ; JedB = G%JedB
 
+  shelf_sfc_acc=.false.
+  if (present(shelf_sfc_accumulation)) shelf_sfc_acc=shelf_sfc_accumulation
+
   call myAlloc(fluxes%ustar,isd,ied,jsd,jed, ustar)
   call myAlloc(fluxes%ustar_gustless,isd,ied,jsd,jed, ustar)
 
@@ -2987,9 +3002,13 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
 
   call myAlloc(fluxes%p_surf,isd,ied,jsd,jed, press)
 
-  call myAlloc(fluxes%frac_shelf_h,isd,ied,jsd,jed, shelf)
-  call myAlloc(fluxes%ustar_shelf,isd,ied,jsd,jed, shelf)
-  call myAlloc(fluxes%iceshelf_melt,isd,ied,jsd,jed, shelf)
+  ! These fields should only be allocated if ice shelf is enabled.
+  if (present(shelf)) then; if (shelf) then
+    call myAlloc(fluxes%frac_shelf_h,isd,ied,jsd,jed, shelf)
+    call myAlloc(fluxes%ustar_shelf,isd,ied,jsd,jed, shelf)
+    call myAlloc(fluxes%iceshelf_melt,isd,ied,jsd,jed, shelf)
+    if (shelf_sfc_acc) call myAlloc(fluxes%shelf_sfc_mass_flux,isd,ied,jsd,jed, shelf_sfc_acc)
+  endif; endif
 
   !These fields should only on allocated when iceberg area is being passed through the coupler.
   call myAlloc(fluxes%ustar_berg,isd,ied,jsd,jed, iceberg)
@@ -3257,6 +3276,8 @@ subroutine deallocate_forcing_type(fluxes)
   if (associated(fluxes%ustar_tidal))          deallocate(fluxes%ustar_tidal)
   if (associated(fluxes%ustar_shelf))          deallocate(fluxes%ustar_shelf)
   if (associated(fluxes%iceshelf_melt))        deallocate(fluxes%iceshelf_melt)
+  if (associated(fluxes%shelf_sfc_mass_flux)) &
+                                               deallocate(fluxes%shelf_sfc_mass_flux)
   if (associated(fluxes%frac_shelf_h))         deallocate(fluxes%frac_shelf_h)
   if (associated(fluxes%ustar_berg))           deallocate(fluxes%ustar_berg)
   if (associated(fluxes%area_berg))            deallocate(fluxes%area_berg)
@@ -3355,6 +3376,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns)
     call rotate_array(fluxes_in%frac_shelf_h, turns, fluxes%frac_shelf_h)
     call rotate_array(fluxes_in%ustar_shelf, turns, fluxes%ustar_shelf)
     call rotate_array(fluxes_in%iceshelf_melt, turns, fluxes%iceshelf_melt)
+    call rotate_array(fluxes_in%shelf_sfc_mass_flux, turns, fluxes%shelf_sfc_mass_flux)
   endif
 
   if (do_iceberg) then
@@ -3603,6 +3625,7 @@ subroutine homogenize_forcing(fluxes, G, GV, US)
     call homogenize_field_t(fluxes%frac_shelf_h, G)
     call homogenize_field_t(fluxes%ustar_shelf, G, tmp_scale=US%Z_to_m*US%s_to_T)
     call homogenize_field_t(fluxes%iceshelf_melt, G, tmp_scale=US%RZ_T_to_kg_m2s)
+    call homogenize_field_t(fluxes%shelf_sfc_mass_flux, G, tmp_scale=US%RZ_T_to_kg_m2s)
   endif
 
   if (do_iceberg) then

src/ice_shelf/MOM_ice_shelf.F90

@@ -9,6 +9,7 @@ module MOM_ice_shelf
 use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
 use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_ROUTINE
 use MOM_coms,                 only : num_PEs
+use MOM_data_override,       only : data_override
 use MOM_diag_mediator, only    : MOM_diag_ctrl=>diag_ctrl
 use MOM_IS_diag_mediator, only : post_data=>post_IS_data
 use MOM_IS_diag_mediator, only : register_diag_field=>register_MOM_IS_diag_field, safe_alloc_ptr
@@ -160,6 +161,8 @@ module MOM_ice_shelf
                                          !! equation of state to use.
   logical :: active_shelf_dynamics       !< True if the ice shelf mass changes as a result
                                          !! the dynamic ice-shelf model.
+  logical :: data_override_shelf_fluxes  !< True if the ice shelf surface mass fluxes can be
+                                         !! written using the data_override feature (only for MOSAIC grids)
   logical :: override_shelf_movement     !< If true, user code specifies the shelf movement
                                          !! instead of using the dynamic ice-shelf mode.
   logical :: isthermo                    !< True if the ice shelf can exchange heat and
@@ -188,7 +191,8 @@ module MOM_ice_shelf
              id_h_shelf = -1, id_h_mask = -1, &
              id_surf_elev = -1, id_bathym = -1, &
              id_area_shelf_h = -1, &
-             id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1
+             id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1, &
+             id_shelf_sfc_mass_flux = -1
   !>@}
 
   integer :: id_read_mass !< An integer handle used in time interpolation of
@@ -321,6 +325,11 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
   G => CS%grid ; US => CS%US
   ISS => CS%ISS
 
+  if (CS%data_override_shelf_fluxes .and. CS%active_shelf_dynamics) then
+       call data_override(G%Domain, 'shelf_sfc_mass_flux', fluxes_in%shelf_sfc_mass_flux, CS%Time, &
+                          scale=US%kg_m2s_to_RZ_T)
+  endif
+
   if (CS%rotate_index) then
     allocate(sfc_state)
     call rotate_surface_state(sfc_state_in, CS%Grid_in, sfc_state, CS%Grid, CS%turns)
@@ -730,6 +739,15 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
       call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, &
                    scale=US%RZ_to_kg_m2)
     endif
+
+   call change_thickness_using_precip(CS, ISS, G, US, fluxes,US%s_to_T*time_step, Time)
+
+   if (CS%debug) then
+     call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, scale=US%Z_to_m)
+     call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using surf acc", G%HI, haloshift=0, &
+                  scale=US%RZ_to_kg_m2)
+   endif
+
   endif
 
   if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0)
@@ -753,6 +771,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
   if (CS%id_shelf_mass > 0) call post_data(CS%id_shelf_mass, ISS%mass_shelf, CS%diag)
   if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, ISS%area_shelf_h, CS%diag)
   if (CS%id_ustar_shelf > 0) call post_data(CS%id_ustar_shelf, fluxes%ustar_shelf, CS%diag)
+  if (CS%id_shelf_sfc_mass_flux > 0) call post_data(CS%id_shelf_sfc_mass_flux, fluxes%shelf_sfc_mass_flux, CS%diag)
+
   if (CS%id_melt > 0) call post_data(CS%id_melt, fluxes%iceshelf_melt, CS%diag)
   if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (sfc_state%sst-ISS%tfreeze), CS%diag)
   if (CS%id_Sbdry > 0) call post_data(CS%id_Sbdry, Sbdry, CS%diag)
@@ -1265,7 +1285,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
   allocate(CS%Grid)
   call MOM_domains_init(CS%Grid%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,&
        domain_name='MOM_Ice_Shelf_in')
-!  allocate(CS%Grid_in%HI)
+  !allocate(CS%Grid_in%HI)
   !call hor_index_init(CS%Grid%Domain, CS%Grid%HI, param_file, &
   !     local_indexing=.not.global_indexing)
   call MOM_grid_init(CS%Grid, param_file, CS%US)
@@ -1354,6 +1374,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
                  "If true, user provided code specifies the ice-shelf "//&
                  "movement instead of the dynamic ice model.", default=.false.)
     CS%active_shelf_dynamics = .not.CS%override_shelf_movement
+    call get_param(param_file, mdl, "DATA_OVERRIDE_SHELF_FLUXES", &
+                  CS%data_override_shelf_fluxes, &
+                 "If true, the data override feature is used to write "//&
+                 "the surface mass flux deposition. This option is only "//&
+                 "available for MOSAIC grid types.", default=.false.)
     call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", CS%GL_regularize, &
                  "If true, regularize the floatation condition at the "//&
                  "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
@@ -1815,6 +1840,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
   if (CS%active_shelf_dynamics) then
     CS%id_h_mask = register_diag_field('ice_shelf_model', 'h_mask', CS%diag%axesT1, CS%Time, &
        'ice shelf thickness mask', 'none')
+    CS%id_shelf_sfc_mass_flux = register_diag_field('ice_shelf_model', 'sfc_mass_flux', CS%diag%axesT1, CS%Time, &
+       'ice shelf surface mass flux deposition from atmosphere', 'none', conversion=US%RZ_T_to_kg_m2s)
   endif
   call MOM_IS_diag_mediator_close_registration(CS%diag)
 
@@ -1845,7 +1872,7 @@ subroutine initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in)
    ! when SHELF_THERMO = True. These fluxes are necessary if one wants to
    ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode).
     call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., &
-         press=.true., water=CS%isthermo, heat=CS%isthermo)
+         press=.true., water=CS%isthermo, heat=CS%isthermo, shelf_sfc_accumulation = CS%active_shelf_dynamics)
   else
     call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
     call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., press=.true.)
@@ -1975,6 +2002,57 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim)
   end select
 
 end subroutine initialize_shelf_mass
+!> This subroutine applies net accumulation/ablation at the top surface to the dynamic ice shelf.
+!>>acc_rate[m-s]=surf_mass_flux/density_ice is ablation/accumulation rate
+!>>positive for accumulation negative for ablation
+subroutine change_thickness_using_precip(CS, ISS, G, US, fluxes,time_step, Time)
+  type(ice_shelf_CS),    intent(in)    :: CS  !< A pointer to the ice shelf control structure
+  type(ocean_grid_type), intent(inout) :: G  !< The ocean's grid structure.
+  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
+                                              !! the ice-shelf state
+  type(forcing),         intent(in)    :: fluxes  !< A structure of surface fluxes that
+                                                  !! includes surface mass flux
+  type(time_type),       intent(in)    :: Time !< The current model time
+  type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
+  real,                  intent(in)    :: time_step !< The time step for this update [T ~> s].
+
+  ! locals
+  integer :: i, j
+  real ::I_rho_ice
+
+  I_rho_ice = 1.0 / CS%density_ice
+
+  !update time
+!  CS%Time = Time
+
+
+!    CS%time_step = time_step
+    ! update surface mass flux  rate
+!    if (CS%surf_mass_flux_from_file) call update_surf_mass_flux(G, US, CS, ISS, Time)
+
+  do j=G%jsc,G%jec ; do i=G%isc,G%iec
+    if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 2)) then
+
+      if (-fluxes%shelf_sfc_mass_flux(i,j) * time_step  < ISS%h_shelf(i,j)) then
+        ISS%h_shelf(i,j) = ISS%h_shelf(i,j) + fluxes%shelf_sfc_mass_flux(i,j) * time_step * I_rho_ice
+      else
+        ! the ice is about to ablate, so set thickness, area, and mask to zero
+        ! NOTE: this is not mass conservative should maybe scale salt & heat flux for this cell
+        ISS%h_shelf(i,j) = 0.0
+        ISS%hmask(i,j) = 0.0
+        ISS%area_shelf_h(i,j) = 0.0
+      endif
+      ISS%mass_shelf(i,j) = ISS%h_shelf(i,j) * CS%density_ice
+    endif
+  enddo ; enddo
+
+  call pass_var(ISS%area_shelf_h, G%domain, complete=.false.)
+  call pass_var(ISS%h_shelf, G%domain, complete=.false.)
+  call pass_var(ISS%hmask, G%domain, complete=.false.)
+  call pass_var(ISS%mass_shelf, G%domain, complete=.true.)
+
+end subroutine change_thickness_using_precip
+
 
 !> Updates the ice shelf mass using data from a file.
 subroutine update_shelf_mass(G, US, CS, ISS, Time)
@@ -2032,12 +2110,13 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time)
 end subroutine update_shelf_mass
 
 !> Save the ice shelf restart file
-subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf)
+subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes)
   type(ice_shelf_CS),         pointer    :: CS !< ice shelf control structure
   type(ocean_grid_type), intent(in)      :: G  !< A pointer to an ocean grid control structure.
   real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: frac_shelf_h !< Ice shelf area fraction [nodim].
   real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: mass_shelf !<Ice shelf mass [R Z -> kg m-2]
-
+  logical, optional                      :: data_override_shelf_fluxes !< If true, shelf fluxes can be written using
+                                               !! the data_override capability (only for MOSAIC grids)
 
   integer :: i, j
 
@@ -2055,6 +2134,11 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf)
     enddo ; enddo
   endif
 
+  if (present(data_override_shelf_fluxes)) then
+     data_override_shelf_fluxes=.false.
+     if (CS%active_shelf_dynamics) data_override_shelf_fluxes = CS%data_override_shelf_fluxes
+  endif
+
 end subroutine ice_shelf_query
 
 !> Save the ice shelf restart file