Separate frazil for ice sheets and sea ice

config_src/drivers/FMS_cap/ocean_model_MOM.F90

@@ -56,7 +56,7 @@ module ocean_model_mod
 use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
 use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces
 use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
-use MOM_ice_shelf, only : ice_sheet_calving_to_ocean_sfc
+use MOM_ice_shelf, only : ice_sheet_calving_to_ocean_sfc, adjust_ice_sheet_frazil
 use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init
 use MOM_wave_interface, only: Update_Surface_Waves
 use iso_fortran_env, only : int64
@@ -413,6 +413,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
 
     call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+    if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+      call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
     call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
 
   endif
@@ -682,6 +685,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
     call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles)
   endif
 
+  !only ,ale ice-shelf frazil adjustments if sfc_state%frazil was updated (do_thermo=True)
+  if (do_thermo .and. OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
 ! Translate state into Ocean.
 !  call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US, &
 !                                   OS%fluxes%p_surf_full, OS%press_to_z)
@@ -990,6 +997,9 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc)
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+  if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
 
 end subroutine ocean_model_init_sfc

config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90

@@ -53,6 +53,7 @@ module MOM_ocean_model_mct
 use MOM_verticalGrid,         only : verticalGrid_type
 use MOM_ice_shelf,            only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
 use MOM_ice_shelf,            only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
+use MOM_ice_shelf,            only : adjust_ice_sheet_frazil
 use MOM_coupler_types,        only : coupler_1d_bc_type, coupler_2d_bc_type
 use MOM_coupler_types,        only : coupler_type_spawn, coupler_type_write_chksums
 use MOM_coupler_types,        only : coupler_type_initialized, coupler_type_copy_data
@@ -404,6 +405,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
     call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+    if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+      call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
     call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
   endif
 
@@ -661,6 +665,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
     call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles)
   endif
 
+  !only make ice-shelf frazil adjustments if sfc_state%frazil was updated (do_thermo=True)
+  if (do_thermo .and. OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
 ! Translate state into Ocean.
 !  call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, &
 !                                   Ice_ocean_boundary%p, OS%press_to_z)
@@ -975,6 +983,9 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc)
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+  if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
 
 end subroutine ocean_model_init_sfc

config_src/drivers/STALE_mct_cap/ocn_comp_mct.F90

@@ -31,7 +31,7 @@ module ocn_comp_mct
 use MOM_variables,        only: surface
 use MOM_domains,          only: MOM_infra_init
 use MOM_restart,          only: save_restart
-use MOM_ice_shelf,        only: ice_shelf_save_restart
+use MOM_ice_shelf,        only: ice_shelf_save_restart, adjust_ice_sheet_frazil
 use MOM_domains,          only: num_pes, root_pe, pe_here
 use MOM_grid,             only: ocean_grid_type, get_global_grid_size
 use MOM_error_handler,    only: MOM_error, FATAL, is_root_pe, WARNING
@@ -781,6 +781,9 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc)
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+  if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
 
 end subroutine ocean_model_init_sfc

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -49,6 +49,7 @@ module MOM_ocean_model_nuopc
 use MOM_verticalGrid,        only : verticalGrid_type
 use MOM_ice_shelf,           only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
 use MOM_ice_shelf,           only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
+use MOM_ice_shelf,           only : adjust_ice_sheet_frazil
 use MOM_coupler_types,       only : coupler_1d_bc_type, coupler_2d_bc_type
 use MOM_coupler_types,       only : coupler_type_spawn, coupler_type_write_chksums
 use MOM_coupler_types,       only : coupler_type_initialized, coupler_type_copy_data
@@ -431,11 +432,18 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
     call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+    if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+      call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
     call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
 
   endif
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
+
+  if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
 ! get number of processors and PE list for stocasthci physics initialization
   call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, &
                  "If true, then stochastically perturb the thermodynamic "//&
@@ -712,6 +720,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
     endif
   endif
 
+  !only make ice-shelf frazil adjustments if sfc_state%frazil was updated (do_thermo=True)
+  if (do_thermo .and. OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
 ! Translate state into Ocean.
 !  call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, &
 !                                   Ice_ocean_boundary%p, OS%press_to_z)
@@ -1031,6 +1043,9 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc)
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
+  if (OS%use_ice_shelf .and. allocated(OS%sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(OS%sfc_state, OS%fluxes, OS%Ice_shelf_CSp)
+
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
 
 end subroutine ocean_model_init_sfc

config_src/drivers/solo_driver/MOM_driver.F90

@@ -51,7 +51,7 @@ program MOM6
   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_ice_shelf,       only : ice_shelf_query, adjust_ice_sheet_frazil
   use MOM_ice_shelf_initialize, only : initialize_ice_SMB
   use MOM_interpolate,     only : time_interp_external_init
   use MOM_io,              only : file_exists, open_ASCII_file, close_file
@@ -316,6 +316,9 @@ program MOM6
 
   call extract_surface_state(MOM_CSp, sfc_state)
 
+  if (use_ice_shelf .and. allocated(sfc_state%frazil)) &
+    call adjust_ice_sheet_frazil(sfc_state, fluxes, Ice_shelf_CSp)
+
   call surface_forcing_init(Time, grid, US, param_file, diag, &
                             surface_forcing_CSp, tracer_flow_CSp)
   call callTree_waypoint("done surface_forcing_init")

src/ice_shelf/MOM_ice_shelf.F90

@@ -82,6 +82,7 @@ module MOM_ice_shelf
 public ice_shelf_save_restart, solo_step_ice_shelf, add_shelf_forces
 public initialize_ice_shelf_fluxes, initialize_ice_shelf_forces
 public ice_sheet_calving_to_ocean_sfc
+public adjust_ice_sheet_frazil
 
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
@@ -210,7 +211,7 @@ module MOM_ice_shelf
              id_tfreeze = -1, id_tfl_shelf = -1, &
              id_thermal_driving = -1, id_haline_driving = -1, &
              id_u_ml = -1, id_v_ml = -1, id_sbdry = -1, &
-             id_h_shelf = -1, id_dhdt_shelf, id_h_mask = -1, &
+             id_h_shelf = -1, id_dhdt_shelf, id_h_mask = -1, id_frazil = -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, &
@@ -881,7 +882,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
 
     !Add frazil formation
     if (add_frazil .and. (ISS%hmask(i,j) == 1 .or. ISS%hmask(i,j) == 2)) &
-      ISS%water_flux(i,j) = ISS%water_flux(i,j) - sfc_state%frazil(i,j) * I_dt_LHF
+      ISS%water_flux(i,j) = ISS%water_flux(i,j) - ISS%frazil(i,j) * I_dt_LHF
     fluxes%iceshelf_melt(i,j) = ISS%water_flux(i,j)
   enddo ; enddo ! i- and j-loops
 
@@ -974,10 +975,14 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
   if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf, ISS%h_shelf, CS%diag)
   if (CS%id_dhdt_shelf > 0) call post_data(CS%id_dhdt_shelf, ISS%dhdt_shelf, CS%diag)
   if (CS%id_h_mask > 0) call post_data(CS%id_h_mask,ISS%hmask,CS%diag)
+  if (CS%id_frazil > 0) call post_data(CS%id_frazil,ISS%frazil,CS%diag)
   if (CS%active_shelf_dynamics) &
       call process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott)
   call disable_averaging(CS%diag)
 
+  !reset used frazil
+  if (add_frazil) ISS%frazil(:,:) = 0.0
+
   call cpu_clock_end(id_clock_shelf)
 
   if (CS%debug) call MOM_forcing_chksum("End of shelf calc flux", fluxes, G, CS%US, haloshift=0)
@@ -993,6 +998,59 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
 
 end subroutine shelf_calc_flux
 
+!> Copies frazil from the ocean surface state to the ice sheet state. Removes frazil that will
+!! be used by the ice sheet from the ocean surface state
+subroutine adjust_ice_sheet_frazil(sfc_state_in, fluxes_in, CS)
+  type(surface), target,  intent(inout) :: sfc_state_in !< A structure containing fields that
+                                                 !! describe the surface state of the ocean.  The
+                                                 !! intent is only inout to allow for halo updates.
+  type(forcing),  target, intent(in)    :: fluxes_in !< structure containing pointers to any
+                                                 !! possible thermodynamic or mass-flux forcing fields.
+  type(ice_shelf_CS),     pointer       :: CS    !< A pointer to the control structure returned
+                                                 !! by a previous call to initialize_ice_shelf.
+  ! Local variables
+  type(ocean_grid_type), pointer :: G => NULL()  !< The grid structure used by the ice shelf.
+  type(ice_shelf_state), pointer :: ISS => NULL() !< A structure with elements that describe
+                                                 !! the ice-shelf state
+  type(surface), pointer :: sfc_state => NULL()
+  type(forcing), pointer :: fluxes => NULL()
+  integer :: i,j,is,ie,js,je
+
+  G => CS%grid ; ISS => CS%ISS
+
+  if (CS%rotate_index) then
+    allocate(sfc_state)
+    call rotate_surface_state(sfc_state_in, sfc_state, G, CS%turns)
+    allocate(fluxes)
+    call allocate_forcing_type(fluxes_in, G, fluxes, turns=CS%turns)
+    call rotate_forcing(fluxes_in, fluxes, CS%turns)
+  else
+    sfc_state => sfc_state_in
+    fluxes => fluxes_in
+  endif
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+
+  do j=js,je ; do i=is,ie
+    !Copy frazil to the ice sheet module where ice sheet is present.
+    !No scaling to account for partial ice-sheet cells is necessary here, as
+    !this is taken care of when applied to the ice sheet.
+    if (fluxes%frac_shelf_h(i,j)>0.0) ISS%frazil(i,j) = sfc_state%frazil(i,j)
+    !Remove the frazil that is used by the ice sheet from sfc_state%frazil
+    !The sfc_state%frazil is sent to the sea-ice module
+    sfc_state%frazil(i,j) = sfc_state%frazil(i,j) * (1.0-fluxes%frac_shelf_h(i,j))
+  enddo; enddo
+
+  if (CS%rotate_index) then
+    call rotate_surface_state(sfc_state, sfc_state_in, G, -CS%turns)
+    ! call rotate_forcing(fluxes, fluxes_in, -CS%turns)
+    call deallocate_surface_state(sfc_state)
+    deallocate(sfc_state)
+    call deallocate_forcing_type(fluxes)
+    deallocate(fluxes)
+  endif
+end subroutine adjust_ice_sheet_frazil
+
 function integrate_over_ice_sheet_area(G, ISS, var, unscale, hemisphere) result(var_out)
   type(ocean_grid_type), intent(in) :: G  !< The grid structure used by the ice shelf.
   type(ice_shelf_state), intent(in) :: ISS  !< A structure with elements that describe the ice-shelf state
@@ -1964,7 +2022,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,
                               "ice sheet/shelf thickness", "m", conversion=US%Z_to_m)
   call register_restart_field(ISS%melt_mask, "melt_mask", .false., CS%restart_CSp, &
                               "Mask that is >0 where ice-shelf melting is allowed", "none")
-
   if (CS%calve_ice_shelf_bergs) then
     call register_restart_field(ISS%calving, "shelf_calving", .true., CS%restart_CSp, &
                                 "Calving flux from ice shelf into icebergs", "kg m-2", conversion=US%RZ_to_kg_m2)
@@ -2126,6 +2183,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,
       'Heat conduction into ice shelf', 'W m-2', conversion=-US%QRZ_T_to_W_m2)
   CS%id_ustar_shelf = register_diag_field('ice_shelf_model', 'ustar_shelf', CS%diag%axesT1, CS%Time, &
       'Fric vel under shelf', 'm/s', conversion=US%Z_to_m*US%s_to_T)
+  CS%id_frazil = register_diag_field('ice_shelf_model', 'frazil', CS%diag%axesT1, CS%Time, &
+     'Frazil heat rejected by the ocean', 'J m-2', conversion=US%Q_to_J_kg*US%RZ_to_kg_m2)
   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', conversion=1.0)

src/ice_shelf/MOM_ice_shelf_state.F90

@@ -48,10 +48,10 @@ module MOM_ice_shelf_state
                                !! ocean-ice interface [R Z T-1 ~> kg m-2 s-1].
     tflux_shelf => NULL(), &   !< The downward diffusive heat flux in the ice
                                !! shelf at the ice-ocean interface [Q R Z T-1 ~> W m-2].
-
     tfreeze => NULL(), &       !< The freezing point potential temperature
                                !! at the ice-ocean interface [C ~> degC].
-
+    frazil => NULL(), &        !< Accumulated heating [J m-2] from frazil formation in the ocean
+                               !! under ice-shelf cells
     !only active when calve_ice_shelf_bergs=true:
     calving => NULL(), &       !< The mass flux per unit area of the ice shelf to convert to
                                !! bergs [R Z T-1 ~> kg m-2 s-1].
@@ -88,6 +88,7 @@ subroutine ice_shelf_state_init(ISS, G)
   allocate(ISS%tflux_shelf(isd:ied,jsd:jed), source=0.0 )
   allocate(ISS%tfreeze(isd:ied,jsd:jed), source=0.0 )
 
+  allocate(ISS%frazil(isd:ied,jsd:jed), source=0.0 )
   allocate(ISS%calving(isd:ied,jsd:jed), source=0.0 )
   allocate(ISS%calving_hflx(isd:ied,jsd:jed), source=0.0 )
 end subroutine ice_shelf_state_init
@@ -102,7 +103,7 @@ subroutine ice_shelf_state_end(ISS)
   deallocate(ISS%mass_shelf, ISS%area_shelf_h, ISS%h_shelf, ISS%dhdt_shelf, ISS%hmask)
 
   deallocate(ISS%tflux_ocn, ISS%water_flux, ISS%salt_flux, ISS%tflux_shelf)
-  deallocate(ISS%tfreeze)
+  deallocate(ISS%tfreeze, ISS%frazil)
 
   deallocate(ISS%calving, ISS%calving_hflx)