Add handling of coupled submarine melt

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -940,6 +940,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Sa_pslv"        , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl"      , "will provide") !-> liquid runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"      , "will provide") !-> ice runoff
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofs"      , "will provide") !-> submarine melt
   if (cesm_coupled) then
     call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofl_glc"  , "will provide") !-> liquid glc runoff
     call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofi_glc"  , "will provide") !-> frozen glc runoff
@@ -1391,6 +1392,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
     call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", &
                            ice_ocean_boundary=Ice_ocean_boundary, mesh=Emesh, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call allocate_optional_iob_fields(importState, Ice_ocean_boundary, isc, iec, jsc, jec)
 
     call MOM_RealizeFields(exportState, fldsFrOcn_num, fldsFrOcn, "Ocn export", mesh=Emesh, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -1682,6 +1684,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
     call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", &
          ice_ocean_boundary=Ice_ocean_boundary, grid=gridIn, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call allocate_optional_iob_fields(importState, Ice_ocean_boundary, isc, iec, jsc, jec)
 
     call MOM_RealizeFields(exportState, fldsFrOcn_num, fldsFrOcn, "Ocn export", grid=gridOut, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -2523,6 +2526,24 @@ subroutine State_SetScalar(value, scalar_id, State, mytask, scalar_name, scalar_
 
 end subroutine State_SetScalar
 
+!> Allocate ice_ocean_boundary fields whose presence depends on whether the coupler
+!! actually connects the corresponding field, rather than on a compile-time flag.
+subroutine allocate_optional_iob_fields(importState, ice_ocean_boundary, isc, iec, jsc, jec)
+  type(ESMF_State)             , intent(in)    :: importState !< Import ESMF_State object.
+  type(ice_ocean_boundary_type), intent(inout) :: ice_ocean_boundary !< Fields from the mediator.
+  integer                      , intent(in)    :: isc !< Start i-index of the computational domain.
+  integer                      , intent(in)    :: iec !< End i-index of the computational domain.
+  integer                      , intent(in)    :: jsc !< Start j-index of the computational domain.
+  integer                      , intent(in)    :: jec !< End j-index of the computational domain.
+  character(len=*), parameter :: subname = '(MOM_cap:allocate_optional_iob_fields)'
+
+  if (NUOPC_IsConnected(importState, fieldName='Foxx_rofs')) then
+    if (is_root_pe()) write(stdout,*) subname // " Allocating ice_ocean_boundary%srunoff"
+    allocate(ice_ocean_boundary%srunoff(isc:iec,jsc:jec), source=0.0)
+  endif
+
+end subroutine allocate_optional_iob_fields
+
 !> Realize the import and export fields using either a grid or a mesh.
 subroutine MOM_RealizeFields(state, nfields, field_defs, tag, ice_ocean_boundary, grid, mesh, rc)
   type(ESMF_State)             , intent(inout)           :: state !< ESMF_State object for
@@ -3082,6 +3103,14 @@ end subroutine shr_log_setLogUnit
 !!     <td></td>
 !! </tr>
 !! <tr>
+!!     <td>Foxx_rofs</td>
+!!     <td>kg m-2 s-1</td>
+!!     <td>runoff</td>
+!!     <td>mass flux of submarine melt. Only allocated in Ice_ocean_boundary if
+!!     advertised by the coupler.</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
 !!     <td>Forr_rofl_glc</td>
 !!     <td>kg m-2 s-1</td>
 !!     <td>runoff</td>

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -239,6 +239,13 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
        isc, iec, jsc, jec, ice_ocean_boundary%frunoff, areacor=med2mod_areacor, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
+  ! submarine melt. This is not preset to 0.0 since it is only allocated if connected.
+  if ( associated(ice_ocean_boundary%srunoff) ) then
+    call state_getimport(importState, 'Foxx_rofs',  &
+         isc, iec, jsc, jec, ice_ocean_boundary%srunoff, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  endif
+
   ! liquid glc runoff
   if ( associated(ice_ocean_boundary%lrunoff_glc) ) then
     ice_ocean_boundary%lrunoff_glc (:,:) = 0._ESMF_KIND_R8

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -390,7 +390,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
                               do_integrals=.true., gas_fields_ocn=gas_fields_ocn, &
                               use_meltpot=use_melt_pot, use_marbl_tracers=OS%use_MARBL)
 
-  call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, &
+  call surface_forcing_init(Time_in, OS%grid, OS%GV, OS%US, param_file, OS%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves)
 
   if (OS%use_ice_shelf)  then

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -37,6 +37,7 @@ module MOM_surface_forcing_nuopc
 use MOM_spatial_means,    only : adjust_area_mean_to_zero
 use MOM_unit_scaling,     only : unit_scale_type
 use MOM_variables,        only : surface
+use MOM_verticalGrid,     only : verticalGrid_type
 use user_revise_forcing,  only : user_alter_forcing, user_revise_forcing_init
 use user_revise_forcing,  only : user_revise_forcing_CS
 use iso_fortran_env,      only : int64
@@ -86,6 +87,14 @@ module MOM_surface_forcing_nuopc
   logical :: use_marbl_tracers  !< enables the MARBL tracer package.
   logical :: enthalpy_cpl       !< Controls if enthalpy terms are provided by the coupler or computed
                                 !! internally.
+  logical :: srunoff_latent_heat !< If true, give submarine melt (srunoff) an effective temperature
+                                !! offset by latent_heat_fusion/C_p, following the Gade (1979)
+                                !! meltwater mixing line, so that melting the ice consumes latent
+                                !! heat. If false, srunoff enters the ocean at the ambient
+                                !! temperature.
+  real :: srunoff_depth_max     !< The maximum depth over which the mass and heat from coupled
+                                !! submarine melt are spread, distributed uniformly by thickness. If 0,
+                                !! srunoff is applied entirely within the top layer [H ~> m or kg m-2].
   real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied
                                 !! from an input file.
@@ -166,6 +175,7 @@ module MOM_surface_forcing_nuopc
 type, public :: ice_ocean_boundary_type
   real, pointer, dimension(:,:) :: lrunoff           =>NULL() !< liquid runoff [km m-2 s-1]
   real, pointer, dimension(:,:) :: frunoff           =>NULL() !< ice runoff [km m-2 s-1]
+  real, pointer, dimension(:,:) :: srunoff           =>NULL() !< submarine melt [km m-2 s-1]
   real, pointer, dimension(:,:) :: lrunoff_glc       =>NULL() !< liquid glc runoff via rof [km m-2 s-1]
   real, pointer, dimension(:,:) :: frunoff_glc       =>NULL() !< frozen glc runoff via rof [km m-2 s-1]
   real, pointer, dimension(:,:) :: u_flux            =>NULL() !< i-direction wind stress [Pa]
@@ -329,7 +339,10 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., &
                                press=.true., fix_accum_bug=.not.CS%ustar_gustless_bug, &
                                cfc=CS%use_CFC, marbl=CS%use_marbl_tracers, hevap=CS%enthalpy_cpl, &
-                               tau_mag=.true., ice_ncat=IOB%ice_ncat)
+                               tau_mag=.true., ice_ncat=IOB%ice_ncat, srunoff=associated(IOB%srunoff))
+    fluxes%latent_heat_fusion = CS%latent_heat_fusion
+    fluxes%srunoff_latent_heat = CS%srunoff_latent_heat
+    fluxes%srunoff_depth_max = CS%srunoff_depth_max
     call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
@@ -503,6 +516,15 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff(i-i0,j-j0) * G%mask2dT(i,j)
     endif
 
+    ! submarine melt flux
+    if (associated(IOB%srunoff)) then
+      fluxes%srunoff(i,j) = kg_m2_s_conversion * IOB%srunoff(i-i0,j-j0) * G%mask2dT(i,j)
+      ! Total latent heat extracted by submarine melt (not depth-dependent)
+      if (CS%srunoff_latent_heat) &
+        fluxes%latent_srunoff_diag(i,j) = - G%mask2dT(i,j) * &
+            IOB%srunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
+    endif
+
     ! add liquid glc runoff flux via rof
     if (associated(IOB%lrunoff_glc)) then
       fluxes%lrunoff_glc(i,j) = kg_m2_s_conversion * IOB%lrunoff_glc(i-i0,j-j0) * G%mask2dT(i,j)
@@ -557,6 +579,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%latent_frunoff_glc_diag(i,j) = fluxes%latent_frunoff_glc_diag(i,j) - G%mask2dT(i,j) * &
           IOB%frunoff_glc(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
     endif
+    ! Submarine melt (srunoff) is not included in fluxes%latent since it is depth-dependent not a
+    ! surface flux. The latent heat associated with srunoff is included in latent heat diagnostics
     if (associated(IOB%q_flux)) then
       fluxes%latent(i,j)           = fluxes%latent(i,j) + &
           IOB%q_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_vapor
@@ -1149,9 +1173,10 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
 end subroutine forcing_save_restart
 
 !> Initialize the surface forcing, including setting parameters and allocating permanent memory.
-subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp, use_waves)
+subroutine surface_forcing_init(Time, G, GV, US, param_file, diag, CS, restore_salt, restore_temp, use_waves)
   type(time_type),          intent(in)    :: Time !< The current model time
   type(ocean_grid_type),    intent(in)    :: G    !< The ocean's grid structure
+  type(verticalGrid_type),  intent(in)    :: GV   !< The ocean's vertical grid structure
   type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
   type(param_file_type),    intent(in)    :: param_file !< A structure to parse for run-time parameters
   type(diag_ctrl), target,  intent(inout) :: diag !< A structure that is used to regulate
@@ -1216,6 +1241,15 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
                  "The latent heat of fusion.", units="J/kg", default=hlf)
   call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", CS%latent_heat_vapor, &
                  "The latent heat of fusion.", units="J/kg", default=hlv)
+  call get_param(param_file, mdl, "SRUNOFF_LATENT_HEAT", CS%srunoff_latent_heat, &
+                 "If true, give submarine melt (srunoff) an effective temperature offset by "//&
+                 "LATENT_HEAT_FUSION/C_p, following the Gade (1979) meltwater mixing line, "//&
+                 "so that melting the ice consumes latent heat. If false, srunoff enters "// &
+                 "the ocean at the ambient temperature.", default=.false.)
+  call get_param(param_file, mdl, "SRUNOFF_DEPTH_MAX", CS%srunoff_depth_max, &
+                 "The maximum depth over which the mass and heat from coupled submarine melt are spread, "//&
+                 "distributed uniformly by thickness. If 0, srunoff is applied entirely "//&
+                 "within the top layer.", units="m", default=0.0, scale=GV%m_to_H)
   call get_param(param_file, mdl, "MAX_P_SURF", CS%max_p_surf, &
                  "The maximum surface pressure that can be exerted by the "//&
                  "atmosphere and floating sea-ice or ice shelves. This is "//&
@@ -1581,6 +1615,9 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   chks = field_chksum( iobt%fprec          ) ; if (root) write(outunit,100) 'iobt%fprec          ', chks
   chks = field_chksum( iobt%lrunoff        ) ; if (root) write(outunit,100) 'iobt%lrunoff        ', chks
   chks = field_chksum( iobt%frunoff        ) ; if (root) write(outunit,100) 'iobt%frunoff        ', chks
+  if (associated(iobt%srunoff)) then
+    chks = field_chksum( iobt%srunoff      ) ; if (root) write(outunit,100) 'iobt%srunoff        ', chks
+  endif
   chks = field_chksum( iobt%lrunoff_glc    ) ; if (root) write(outunit,100) 'iobt%lrunoff_glc    ', chks
   chks = field_chksum( iobt%frunoff_glc    ) ; if (root) write(outunit,100) 'iobt%frunoff_glc    ', chks
   chks = field_chksum( iobt%p              ) ; if (root) write(outunit,100) 'iobt%p              ', chks