diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90
index 6cd379ea9c..f1da664a7a 100644
--- a/config_src/coupled_driver/ocean_model_MOM.F90
+++ b/config_src/coupled_driver/ocean_model_MOM.F90
@@ -157,6 +157,7 @@ module ocean_model_mod
   logical :: use_ice_shelf    ! If true, the ice shelf model is enabled.
   logical :: icebergs_apply_rigid_boundary  ! If true, the icebergs can change ocean bd condition.
   real :: kv_iceberg          ! The viscosity of the icebergs in m2/s (for ice rigidity)
+  real :: Lat_fusion          ! Latent heat of fusio
   real :: density_iceberg     ! A typical density of icebergs in kg/m3 (for ice rigidity) 
   type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL()
   logical :: restore_salinity ! If true, the coupled MOM driver adds a term to
@@ -293,6 +294,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in)
                  "The viscosity of the icebergs",  units="m2 s-1",default=1.0e10)
     call get_param(param_file, mod, "DENSITY_ICEBERGS",  OS%density_iceberg, &
                   "A typical density of icebergs.", units="kg m-3", default=917.0)
+    call get_param(param_file, mod, "Lat_fusion",  OS%Lat_fusion, &
+                 "Latent heat of fusion for ice",  units="J kg-1",default=3.34e5)
   endif
 
   OS%press_to_z = 1.0/(Rho0*G_Earth)
@@ -420,7 +423,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
     endif
     if (OS%icebergs_apply_rigid_boundary)  then
       !This assumes that the iceshelf and ocean are on the same grid. I hope this is true
-      call add_berg_flux_to_shelf(OS%grid, OS%fluxes,OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg)
+      call add_berg_flux_to_shelf(OS%grid, OS%fluxes,OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%Lat_fusion, OS%State, time_step)
     endif
     ! Indicate that there are new unused fluxes.
     OS%fluxes%fluxes_used = .false.
@@ -434,7 +437,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
     endif
     if (OS%icebergs_apply_rigid_boundary)  then
      !This assumes that the iceshelf and ocean are on the same grid. I hope this is true
-     call add_berg_flux_to_shelf(OS%grid, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg)
+     call add_berg_flux_to_shelf(OS%grid, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%Lat_fusion, OS%State, time_step)
     endif
   
     call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid, weight)
@@ -503,15 +506,19 @@ end subroutine update_ocean_model
 ! </DESCRIPTION>
 !
 
-subroutine add_berg_flux_to_shelf(G, fluxes, use_ice_shelf, density_ice, kv_ice)
+subroutine add_berg_flux_to_shelf(G, fluxes, use_ice_shelf, density_ice, kv_ice, Lat_fusion, state, time_step)
   type(ocean_grid_type),              intent(inout)    :: G
   type(forcing),                      intent(inout) :: fluxes
+  type(surface),                      intent(inout) :: state
   logical,                            intent(in) :: use_ice_shelf
   real, intent(in) :: kv_ice       ! The viscosity of ice, in m2 s-1.
   real, intent(in) :: density_ice  ! A typical density of ice, in kg m-3.
+  real, intent(in) :: Lat_fusion   ! The latent heat of fusion, in J kg-1.
+  real, intent(in) :: time_step   ! The latent heat of fusion, in J kg-1.
 ! Arguments:
 !  (in)      fluxes - A structure of surface fluxes that may be used.
 !  (in)      G - The ocean's grid structure.
+  real :: fraz          ! refreezing rate in kg m-2 s-1
   integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed
@@ -559,6 +566,33 @@ subroutine add_berg_flux_to_shelf(G, fluxes, use_ice_shelf, density_ice, kv_ice)
                                    max(fluxes%mass_berg(i,j), fluxes%mass_berg(i,j+1)))
     enddo ; enddo
     call pass_vector(fluxes%frac_shelf_u, fluxes%frac_shelf_v, G%domain, TO_ALL, CGRID_NE)
+    
+    !Zero'ing out other fluxes under the tabular icebergs 
+    do j=jsd,jed ; do i=isd,ied
+        if (fluxes%frac_shelf_h(i,j) > 0.5) then  !Only applying for ice shelf covering most of cell
+            
+          if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
+          if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
+          if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
+          if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
+
+          ! Add frazil formation diagnosed by the ocean model (J m-2) in the
+          ! form of surface layer evaporation (kg m-2 s-1). Update lprec in the
+          ! control structure for diagnostic purposes.
+
+          if (associated(state%frazil)) then
+            fraz = state%frazil(i,j) / time_step / Lat_fusion
+            if (associated(fluxes%evap)) fluxes%evap(i,j) = fluxes%evap(i,j) - fraz
+            !CS%lprec(i,j)=CS%lprec(i,j) - fraz
+            state%frazil(i,j) = 0.0
+          endif
+
+          !Alon: Should these be set to zero too?
+          if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
+          if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
+          if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
+        endif
+    enddo ; enddo
 
 end subroutine add_berg_flux_to_shelf