diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 02b0b622a3..32a1f75c0c 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -38,6 +38,8 @@ module MOM_variables
     v, &        !< The mixed layer meridional velocity in m s-1.
     sea_lev, &  !< The sea level in m.  If a reduced surface gravity is
                 !! used, that is compensated for in sea_lev.
+    melt_potential, & !< Instantaneous amount of heat that can be used to melt sea ice,
+                      !! in J m-2. This is computed w.r.t. the freezing temperature.
     ocean_mass, &  !< The total mass of the ocean in kg m-2.
     ocean_heat, &  !< The total heat content of the ocean in C kg m-2.
     ocean_salt, &  !< The total salt content of the ocean in kgSalt m-2.
@@ -280,7 +282,7 @@ module MOM_variables
 !> This subroutine allocates the fields for the surface (return) properties of
 !! the ocean model.  Unused fields are unallocated.
 subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
-                                  gas_fields_ocn)
+                                  gas_fields_ocn, use_meltpot)
   type(ocean_grid_type), intent(in)    :: G                !< ocean grid structure
   type(surface),         intent(inout) :: sfc_state        !< ocean surface state type to be allocated.
   logical,     optional, intent(in)    :: use_temperature  !< If true, allocate the space for thermodynamic variables.
@@ -292,8 +294,10 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
                                               !! in the calculation of additional gas or other
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
+  logical,     optional, intent(in)    :: use_meltpot      !< If true, allocate the space for melt potential
 
-  logical :: use_temp, alloc_integ
+  ! local variables
+  logical :: use_temp, alloc_integ, use_melt_potential
   integer :: is, ie, js, je, isd, ied, jsd, jed
   integer :: isdB, iedB, jsdB, jedB
 
@@ -303,6 +307,7 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
 
   use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature
   alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals
+  use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot
 
   if (sfc_state%arrays_allocated) return
 
@@ -317,6 +322,10 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
   allocate(sfc_state%u(IsdB:IedB,jsd:jed)) ; sfc_state%u(:,:) = 0.0
   allocate(sfc_state%v(isd:ied,JsdB:JedB)) ; sfc_state%v(:,:) = 0.0
 
+  if (use_melt_potential) then
+    allocate(sfc_state%melt_potential(isd:ied,jsd:jed)) ; sfc_state%melt_potential(:,:) = 0.0
+  endif
+
   if (alloc_integ) then
     ! Allocate structures for the vertically integrated ocean_mass, ocean_heat,
     ! and ocean_salt.
@@ -342,6 +351,7 @@ subroutine deallocate_surface_state(sfc_state)
 
   if (.not.sfc_state%arrays_allocated) return
 
+  if (allocated(sfc_state%melt_potential)) deallocate(sfc_state%melt_potential)
   if (allocated(sfc_state%SST)) deallocate(sfc_state%SST)
   if (allocated(sfc_state%SSS)) deallocate(sfc_state%SSS)
   if (allocated(sfc_state%sfc_density)) deallocate(sfc_state%sfc_density)