diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90
index 8d04ba19c7..c86a225b5b 100644
--- a/config_src/drivers/nuopc_cap/mom_cap.F90
+++ b/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
!! |
!!
!!
+!! | Foxx_rofs |
+!! kg m-2 s-1 |
+!! runoff |
+!! mass flux of submarine melt. Only allocated in Ice_ocean_boundary if
+!! advertised by the coupler. |
+!! |
+!!
+!!
!! | Forr_rofl_glc |
!! kg m-2 s-1 |
!! runoff |
diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90
index 1accf7472a..47bf52818c 100644
--- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90
+++ b/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
diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
index b8bb87bff3..65b506fc3d 100644
--- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
+++ b/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
diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
index c83a47d498..3e62abfc44 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/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(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
@@ -670,6 +694,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
G%areaT(i,j)
net_FW2(i,j) = net_FW(i,j) / G%areaT(i,j)
enddo ; enddo
+ ! fluxes%srunoff is only associated when submarine-melt coupling is actually configured,
+ ! unlike the other water fluxes above, so it needs its own guard here.
+ if (associated(fluxes%srunoff)) then
+ do j=js,je ; do i=is,ie
+ net_FW(i,j) = net_FW(i,j) + fluxes%srunoff(i,j) * G%areaT(i,j)
+ net_FW2(i,j) = net_FW(i,j) / G%areaT(i,j)
+ enddo ; enddo
+ endif
if (CS%adjust_net_fresh_water_by_scaling) then
call adjust_area_mean_to_zero(net_FW2, G, fluxes%netFWGlobalScl, unscale=US%RZ_T_to_kg_m2s)
@@ -1149,9 +1181,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 +1249,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 +1623,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
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index 90a483ec30..9827f3dc05 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -30,7 +30,7 @@ module MOM_forcing_type
#include
-public extractFluxes1d, extractFluxes2d, optics_type
+public extractFluxes1d, extractFluxes2d, optics_type, distribute_srunoff
public MOM_forcing_chksum, MOM_mech_forcing_chksum
public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, find_ustar
public forcing_accumulate, fluxes_accumulate
@@ -114,7 +114,9 @@ module MOM_forcing_type
latent => NULL(), & !< latent [Q R Z T-1 ~> W m-2] (typically < 0)
sens => NULL(), & !< sensible [Q R Z T-1 ~> W m-2] (typically negative)
seaice_melt_heat => NULL(), & !< sea ice and snow melt or formation [Q R Z T-1 ~> W m-2] (typically negative)
- heat_added => NULL() !< additional heat flux from SST restoring or flux adjustments [Q R Z T-1 ~> W m-2]
+ heat_added => NULL(), & !< additional heat flux from SST restoring or flux adjustments [Q R Z T-1 ~> W m-2]
+ latent_srunoff => NULL() !< latent from submarine melt (srunoff), which can be distributed over depth
+ !! and so is not included in latent (typically < 0) [Q R Z T-1 ~> W m-2]
! components of latent heat fluxes used for diagnostic purposes
real, pointer, dimension(:,:) :: &
@@ -131,6 +133,7 @@ module MOM_forcing_type
vprec => NULL(), & !< virtual liquid precip associated w/ SSS restoring [R Z T-1 ~> kg m-2 s-1]
lrunoff => NULL(), & !< liquid river runoff entering ocean [R Z T-1 ~> kg m-2 s-1]
frunoff => NULL(), & !< frozen river runoff (calving) entering ocean [R Z T-1 ~> kg m-2 s-1]
+ srunoff => NULL(), & !< submarine melt entering ocean [R Z T-1 ~> kg m-2 s-1]
lrunoff_glc => NULL(), & !< liquid river glacier runoff entering ocean [R Z T-1 ~> kg m-2 s-1]
frunoff_glc => NULL(), & !< frozen river glacier runoff entering ocean [R Z T-1 ~> kg m-2 s-1]
seaice_melt => NULL() !< snow/seaice melt (positive) or formation (negative) [R Z T-1 ~> kg m-2 s-1]
@@ -151,6 +154,9 @@ module MOM_forcing_type
heat_content_vprec => NULL(), & !< heat content associated with virtual >0 precip [Q R Z T-1 ~> W m-2]
heat_content_lrunoff => NULL(), & !< heat content associated with liquid runoff [Q R Z T-1 ~> W m-2]
heat_content_frunoff => NULL(), & !< heat content associated with frozen runoff [Q R Z T-1 ~> W m-2]
+ heat_content_srunoff => NULL(), & !< heat content associated with submarine melt [Q R Z T-1 ~> W m-2]
+ heat_content_eff_srunoff => NULL(), & !< effective heat content associated with submarine melt,
+ !! i.e. heat content plus e.g. latent heat [Q R Z T-1 ~> W m-2]
heat_content_lrunoff_glc => NULL(), & !< heat content associated with liquid runoff [Q R Z T-1 ~> W m-2]
heat_content_frunoff_glc => NULL(), & !< heat content associated with frozen runoff [Q R Z T-1 ~> W m-2]
heat_content_massout => NULL(), & !< heat content associated with mass leaving ocean [Q R Z T-1 ~> W m-2]
@@ -225,6 +231,17 @@ module MOM_forcing_type
!! average of salt_left_behind
real :: C_p !< heat capacity of seawater [Q C-1 ~> J kg-1 degC-1].
!! C_p is is the same value as in thermovar_ptrs_type.
+ real :: latent_heat_fusion !< The latent heat of fusion, the same value as
+ !! CS%latent_heat_fusion in surface_forcing_CS [J kg-1].
+ logical :: srunoff_latent_heat = .false. !< 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 = 0.0 !< 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. The same value as
+ !! CS%srunoff_depth_max in surface_forcing_CS [H ~> m or kg m-2].
! arrays needed in the some tracer modules, e.g., MOM_CFC_cap
real, pointer, dimension(:,:) :: &
@@ -336,6 +353,7 @@ module MOM_forcing_type
integer :: id_precip = -1, id_vprec = -1
integer :: id_lprec = -1, id_fprec = -1
integer :: id_lrunoff = -1, id_frunoff = -1
+ integer :: id_srunoff = -1
integer :: id_lrunoff_glc = -1, id_frunoff_glc = -1
integer :: id_net_massout = -1, id_net_massin = -1
integer :: id_massout_flux = -1, id_massin_flux = -1
@@ -346,6 +364,7 @@ module MOM_forcing_type
integer :: id_total_precip = -1, id_total_vprec = -1
integer :: id_total_lprec = -1, id_total_fprec = -1
integer :: id_total_lrunoff = -1, id_total_frunoff = -1
+ integer :: id_total_srunoff = -1
integer :: id_total_lrunoff_glc = -1, id_total_frunoff_glc = -1
integer :: id_total_net_massout = -1, id_total_net_massin = -1
integer :: id_total_seaice_melt = -1
@@ -361,9 +380,10 @@ module MOM_forcing_type
integer :: id_sw = -1, id_lw = -1
integer :: id_sw_vis = -1, id_sw_nir = -1
integer :: id_lat_evap = -1, id_lat_frunoff = -1
- integer :: id_lat_frunoff_glc = -1
+ integer :: id_lat_frunoff_glc = -1, id_lat_srunoff = -1
integer :: id_lat = -1, id_lat_fprec = -1
integer :: id_heat_content_lrunoff = -1, id_heat_content_frunoff = -1
+ integer :: id_heat_content_srunoff = -1, id_heat_content_eff_srunoff = -1
integer :: id_heat_content_lrunoff_glc= -1, id_heat_content_frunoff_glc= -1
integer :: id_heat_content_lprec = -1, id_heat_content_fprec = -1
integer :: id_heat_content_cond = -1, id_heat_content_surfwater = -1
@@ -378,9 +398,11 @@ module MOM_forcing_type
integer :: id_total_sens = -1, id_total_LwLatSens = -1
integer :: id_total_sw = -1, id_total_lw = -1
integer :: id_total_lat_evap = -1, id_total_lat_frunoff = -1
+ integer :: id_total_lat_srunoff = -1
integer :: id_total_lat_frunoff_glc = -1
integer :: id_total_lat = -1, id_total_lat_fprec = -1
integer :: id_total_heat_content_lrunoff = -1, id_total_heat_content_frunoff = -1
+ integer :: id_total_heat_content_srunoff = -1, id_total_heat_content_eff_srunoff = -1
integer :: id_total_heat_content_lrunoff_glc= -1, id_total_heat_content_frunoff_glc=-1
integer :: id_total_heat_content_lprec = -1, id_total_heat_content_fprec = -1
integer :: id_total_heat_content_cond = -1, id_total_heat_content_surfwater = -1
@@ -1030,6 +1052,147 @@ subroutine extractFluxes2d(G, GV, US, fluxes, optics, nsw, dt, FluxRescaleDepth,
end subroutine extractFluxes2d
+!> Distribute the mass and (optionally) heat content of submarine melt (srunoff) uniformly by
+!! thickness over the upper fluxes%srunoff_depth_max of a column. If fluxes%srunoff_depth_max is 0,
+!! the distribution depth for a column instead falls back to that column's own top layer thickness,
+!! so that srunoff is deposited entirely within the top layer.
+!! This routine is called from both:
+!! - applyBoundaryFluxesInOut, to update thickness and temperature and dilute salinity.
+!! - applyTracerBoundaryFluxesInOut, to dilute a tracer concentration.
+subroutine distribute_srunoff(G, GV, dt, fluxes, j, h2d, C2d, T2d, C_p, US, &
+ g_Hconv2, cTKE, dSV_dT, dSV_dS, dHeat_total, dTempxPmE_total)
+ type(ocean_grid_type), intent(in) :: G !< Grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ real, intent(in) :: dt !< Time step [T ~> s]
+ type(forcing), intent(in) :: fluxes !< Surface fluxes, including srunoff
+ integer, intent(in) :: j !< j-index to work on
+ real, dimension(SZI_(G),SZK_(GV)), &
+ intent(inout) :: h2d !< Layer thicknesses [H ~> m or kg m-2]
+ real, dimension(SZI_(G),SZK_(GV)), &
+ intent(inout) :: C2d !< A concentration diluted by srunoff's
+ !! fresh water, e.g. salinity or a tracer's own
+ !! concentration [CU]
+ real, dimension(SZI_(G),SZK_(GV)), &
+ optional, intent(inout) :: T2d !< Layer temperatures [C ~> degC]
+ real, optional, intent(in) :: C_p !< The heat capacity of seawater
+ !! [Q C-1 ~> J kg-1 degC-1]
+ type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
+ real, optional, intent(in) :: g_Hconv2 !< A conversion factor for use in the
+ !! TKE calculation.
+ !! [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2]
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+ optional, intent(inout) :: cTKE !< Turbulent kinetic energy requirement
+ !! to mix forcing through each layer.
+ !! [R Z3 T-2 ~> J m-2]
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+ optional, intent(in) :: dSV_dT !< Partial derivative of specific
+ !! volume with potential temperature.
+ !! [R-1 C-1 ~> m3 kg-1 degC-1]
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+ optional, intent(in) :: dSV_dS !< Partial derivative of specific
+ !! volume with salinity.
+ !! [R-1 S-1 ~> m3 kg-1 ppt-1]
+ real, dimension(SZI_(G)), &
+ optional, intent(out) :: dHeat_total !< The total heat added to each
+ !! column by srunoff.
+ !! [Q H ~> J m-2]
+ real, dimension(SZI_(G)), &
+ optional, intent(out) :: dTempxPmE_total !< The total srunoff mass
+ !! added to each column,
+ !! weighted by ambient
+ !! temperature.
+ !! [C H ~> degC m or degC kg m-2]
+
+ ! Local variables
+ real :: dM_tot ! Total mass of srunoff to be distributed this timestep [H ~> m or kg m-2]
+ real :: dM_k ! The portion of dM_tot applied to a layer [H ~> m or kg m-2]
+ real :: dHeat ! The heat added to the column so far, accumulated for diagnostics [Q H ~> J m-2]
+ real :: dTempxPmE ! The ambient-temperature-weighted mass added to the column so far,
+ ! accumulated for diagnostics [C H ~> degC m or degC kg m-2]
+ real :: depth ! The requested depth over which srunoff is distributed [H ~> m or kg m-2]
+ real :: total_h ! The full thickness of the water column [H ~> m or kg m-2]
+ real :: eff_depth ! The depth actually spanned in this column, capped at total_h so that the
+ ! full amount of srunoff is always distributed even in a column shallower
+ ! than depth [H ~> m or kg m-2]
+ real :: h_above ! Cumulative thickness above the current layer, before adding srunoff [H ~> m or kg m-2]
+ real :: h_in_range ! The portion of a layer's thickness within eff_depth [H ~> m or kg m-2]
+ real :: T_eff ! The effective temperature carried by srunoff mass added to a layer [C ~> degC]
+ real :: hOld ! The layer thickness before adding srunoff mass [H ~> m or kg m-2]
+ real :: I_Habs ! The inverse of eff_depth [H-1 ~> m-1 or m2 kg-1]
+ logical :: do_heat ! If true, T2d and C_p are present and get updated
+ logical :: do_energetics ! If true, cTKE and its supporting arguments are also present
+ integer :: i, k, is, ie, nz
+
+ is = G%isc ; ie = G%iec ; nz = GV%ke
+ do_heat = present(T2d) .and. present(C_p)
+ do_energetics = do_heat .and. present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS) &
+ .and. present(g_Hconv2)
+ if (do_heat .and. fluxes%srunoff_latent_heat .and. .not.present(US)) call MOM_error(FATAL, &
+ "distribute_srunoff: US must be present when T2d is present and "//&
+ "fluxes%srunoff_latent_heat is true.")
+
+ do i=is,ie
+ if (present(dHeat_total)) dHeat_total(i) = 0.0
+ if (present(dTempxPmE_total)) dTempxPmE_total(i) = 0.0
+
+ if ((G%mask2dT(i,j) > 0.0) .and. (fluxes%srunoff(i,j) > 0.0)) then
+
+ dM_tot = GV%RZ_to_H * dt * fluxes%srunoff(i,j)
+
+ ! When fluxes%srunoff_depth_max is 0 (the default), fall back to this column's own top layer
+ ! thickness as the distribution depth, so srunoff is deposited within the top layer.
+ depth = fluxes%srunoff_depth_max
+ if (depth <= 0.0) depth = max(h2d(i,1), GV%H_subroundoff)
+
+ ! Cap the depth actually used at the full column thickness, so that a column shallower
+ ! than the requested depth still receives all of dM_tot.
+ ! Note, with the currently implementation, if the column is zero thickness any srunoff
+ ! is just dropped.
+ total_h = 0.0
+ do k=1,nz ; total_h = total_h + h2d(i,k) ; enddo
+ eff_depth = min(depth, total_h)
+ I_Habs = 0.0 ; if (eff_depth > 0.0) I_Habs = 1.0 / eff_depth
+
+ h_above = 0.0
+ dHeat = 0.0
+ dTempxPmE = 0.0
+
+ ! Spread dM_tot over the upper eff_depth of the column, in proportion to how much of each
+ ! layer's thickness falls within that range.
+ do k=1,nz
+ h_in_range = max(0.0, min(h_above + h2d(i,k), eff_depth) - h_above)
+ h_above = h_above + h2d(i,k)
+ if (h_in_range > 0.0) then
+ dM_k = dM_tot * (h_in_range * I_Habs)
+ hOld = h2d(i,k)
+
+ if (do_heat) then
+ T_eff = T2d(i,k)
+ if (fluxes%srunoff_latent_heat) &
+ T_eff = T_eff - fluxes%latent_heat_fusion*US%J_kg_to_Q/C_p
+ dHeat = dHeat + C_p*dM_k*T_eff
+ dTempxPmE = dTempxPmE + dM_k*T2d(i,k)
+
+ if (do_energetics) then
+ ! Calculate the energy required to mix the srunoff added to this layer.
+ cTKE(i,j,k) = cTKE(i,j,k) + 0.5*g_Hconv2*(hOld*dM_k) * &
+ ((T2d(i,k) - T_eff) * dSV_dT(i,j,k) + C2d(i,k) * dSV_dS(i,j,k))
+ endif
+ endif
+
+ h2d(i,k) = hOld + dM_k
+ if (do_heat) T2d(i,k) = (hOld*T2d(i,k) + dM_k*T_eff) / h2d(i,k)
+ C2d(i,k) = (hOld*C2d(i,k)) / h2d(i,k)
+ endif
+ enddo ! k-loop for layer distribution
+
+ if (present(dHeat_total)) dHeat_total(i) = dHeat
+ if (present(dTempxPmE_total)) dTempxPmE_total(i) = dTempxPmE
+
+ endif
+ enddo
+
+end subroutine distribute_srunoff
!> This routine calculates surface buoyancy flux by adding up the heat, FW & salt fluxes.
!! These are actual fluxes, with units of stuff per time. Setting dt=1 in the call to
@@ -1363,6 +1526,9 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
if (associated(fluxes%latent_frunoff_diag)) &
call hchksum(fluxes%latent_frunoff_diag, mesg//" fluxes%latent_frunoff_diag", G%HI, &
haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
+ if (associated(fluxes%latent_srunoff)) &
+ call hchksum(fluxes%latent_srunoff, mesg//" fluxes%latent_srunoff", G%HI, &
+ haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
if (associated(fluxes%latent_frunoff_glc_diag)) &
call hchksum(fluxes%latent_frunoff_glc_diag, mesg//" fluxes%latent_frunoff_glc_diag", G%HI, &
haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
@@ -1402,6 +1568,8 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s)
if (associated(fluxes%frunoff_glc)) &
call hchksum(fluxes%frunoff_glc, mesg//" fluxes%frunoff_glc", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s)
+ if (associated(fluxes%srunoff)) &
+ call hchksum(fluxes%srunoff, mesg//" fluxes%srunoff", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s)
if (associated(fluxes%heat_content_lrunoff)) &
call hchksum(fluxes%heat_content_lrunoff, mesg//" fluxes%heat_content_lrunoff", G%HI, &
haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
@@ -1414,6 +1582,12 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
if (associated(fluxes%heat_content_frunoff_glc)) &
call hchksum(fluxes%heat_content_frunoff_glc, mesg//" fluxes%heat_content_frunoff_glc", G%HI, &
haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
+ if (associated(fluxes%heat_content_srunoff)) &
+ call hchksum(fluxes%heat_content_srunoff, mesg//" fluxes%heat_content_srunoff", G%HI, &
+ haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
+ if (associated(fluxes%heat_content_eff_srunoff)) &
+ call hchksum(fluxes%heat_content_eff_srunoff, mesg//" fluxes%heat_content_eff_srunoff", G%HI, &
+ haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
if (associated(fluxes%heat_content_lprec)) &
call hchksum(fluxes%heat_content_lprec, mesg//" fluxes%heat_content_lprec", G%HI, &
haloshift=hshift, unscale=US%QRZ_T_to_W_m2)
@@ -1517,6 +1691,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg)
call locMsg(fluxes%latent_fprec_diag,'latent_fprec_diag')
call locMsg(fluxes%latent_frunoff_diag,'latent_frunoff_diag')
call locMsg(fluxes%latent_frunoff_glc_diag,'latent_frunoff_glc_diag')
+ call locMsg(fluxes%latent_srunoff,'latent_srunoff')
call locMsg(fluxes%sens,'sens')
call locMsg(fluxes%evap,'evap')
call locMsg(fluxes%lprec,'lprec')
@@ -1532,10 +1707,13 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg)
call locMsg(fluxes%lrunoff_glc,'lrunoff_glc')
call locMsg(fluxes%frunoff,'frunoff')
call locMsg(fluxes%frunoff_glc,'frunoff_glc')
+ call locMsg(fluxes%srunoff,'srunoff')
call locMsg(fluxes%heat_content_lrunoff,'heat_content_lrunoff')
call locMsg(fluxes%heat_content_lrunoff_glc,'heat_content_lrunoff_glc')
call locMsg(fluxes%heat_content_frunoff,'heat_content_frunoff')
call locMsg(fluxes%heat_content_frunoff_glc,'heat_content_frunoff_glc')
+ call locMsg(fluxes%heat_content_srunoff,'heat_content_srunoff')
+ call locMsg(fluxes%heat_content_eff_srunoff,'heat_content_eff_srunoff')
call locMsg(fluxes%heat_content_lprec,'heat_content_lprec')
call locMsg(fluxes%heat_content_fprec,'heat_content_fprec')
call locMsg(fluxes%heat_content_vprec,'heat_content_vprec')
@@ -1653,7 +1831,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
! surface mass flux maps
handles%id_prcme = register_diag_field('ocean_model', 'PRCmE', diag%axesT1, Time, &
- 'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', &
+ 'Net surface water flux (precip+melt+lrunoff+srunoff+ice calving-evap)', &
'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, &
standard_name='water_flux_into_sea_water', cmor_field_name='wfo', &
cmor_standard_name='water_flux_into_sea_water',cmor_long_name='Water Flux Into Sea Water')
@@ -1709,6 +1887,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
cmor_standard_name='water_flux_into_sea_water_from_rivers', &
cmor_long_name='Water Flux into Sea Water From Rivers')
+ handles%id_srunoff = register_diag_field('ocean_model', 'srunoff', diag%axesT1, Time, &
+ 'Submarine melt into ocean', &
+ units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s)
+
if (present(use_glc_runoff)) then
handles%id_frunoff_glc = register_diag_field('ocean_model', 'frunoff_glc', diag%axesT1, Time, &
'Frozen glacier runoff (calving) and iceberg melt into ocean', &
@@ -1730,11 +1912,13 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s)
handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, &
- 'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', &
+ 'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation). '//&
+ 'Excludes coupled submarine melt (srunoff).', &
'kg m-2', conversion=diag%GV%H_to_kg_m2)
handles%id_massin_flux = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, &
- 'Net mass flux of freshwater into the ocean (used in boundary flux calculation)', &
+ 'Net mass flux of freshwater into the ocean (used in boundary flux calculation). '//&
+ 'Excludes coupled submarine melt (srunoff).', &
'kg m-2', conversion=diag%GV%H_to_kg_m2)
!=========================================================================
@@ -1803,6 +1987,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
cmor_standard_name='water_flux_into_sea_water_from_rivers_area_integrated', &
cmor_long_name='Water Flux into Sea Water From Rivers Area Integrated')
+ handles%id_total_srunoff = register_scalar_field('ocean_model', 'total_srunoff', Time, diag, &
+ long_name='Area integrated submarine melt into ocean', &
+ units='kg s-1', conversion=US%RZL2_to_kg*US%s_to_T)
+
if (present(use_glc_runoff)) then
handles%id_total_frunoff_glc = register_scalar_field('ocean_model', 'total_frunoff_glc', Time, diag, &
long_name='Area integrated frozen glacier runoff (calving) & iceberg melt into ocean', &
@@ -1873,6 +2061,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
'W m-2', conversion=US%QRZ_T_to_W_m2, &
standard_name='temperature_flux_due_to_runoff_expressed_as_heat_flux_into_sea_water')
+ handles%id_heat_content_srunoff = register_diag_field('ocean_model', 'heat_content_srunoff', &
+ diag%axesT1, Time, 'Heat content (relative to 0C) of submarine melt into ocean', &
+ 'W m-2', conversion=US%QRZ_T_to_W_m2)
+
+ handles%id_heat_content_eff_srunoff = register_diag_field('ocean_model', 'heat_content_eff_srunoff', &
+ diag%axesT1, Time, 'Effective heat content (relative to 0C) of submarine melt into ocean, '//&
+ 'i.e. heat content plus e.g. latent heat', 'W m-2', conversion=US%QRZ_T_to_W_m2)
+
if (present(use_glc_runoff)) then
handles%id_heat_content_frunoff_glc = register_diag_field('ocean_model', 'heat_content_frunoff_glc', &
diag%axesT1, Time, 'Heat content (relative to 0C) of solid glacier runoff into ocean', &
@@ -1927,7 +2123,8 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
cmor_long_name='Heat Content (relative to 0degC) of Water Leaving Ocean via Evaporation and Ice Formation')
handles%id_heat_content_massin = register_diag_field('ocean_model', 'heat_content_massin', &
- diag%axesT1, Time,'Heat content (relative to 0degC) of net mass entering ocean ocean',&
+ diag%axesT1, Time,'Heat content (relative to 0degC) of net mass entering ocean ocean. '//&
+ 'Excludes coupled submarine melt (srunoff).',&
'W m-2', conversion=US%QRZ_T_to_W_m2)
handles%id_net_heat_coupler = register_diag_field('ocean_model', 'net_heat_coupler', &
@@ -1992,6 +2189,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
'Latent heat flux into ocean due to melting of frozen glacier runoff', 'W m-2', conversion=US%QRZ_T_to_W_m2)
endif
+ handles%id_lat_srunoff = register_diag_field('ocean_model', 'latent_srunoff', diag%axesT1, Time, &
+ 'Latent heat flux into ocean due to melting of coupled submarine melt (srunoff)', &
+ 'W m-2', conversion=US%QRZ_T_to_W_m2)
+
handles%id_sens = register_diag_field('ocean_model', 'sensible', diag%axesT1, Time, &
'Sensible heat flux into ocean', 'W m-2', conversion=US%QRZ_T_to_W_m2, &
standard_name='surface_downward_sensible_heat_flux', &
@@ -2031,6 +2232,16 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
cmor_long_name= &
'Temperature Flux due to Runoff Expressed as Heat Flux into Sea Water Area Integrated')
+ handles%id_total_heat_content_srunoff = register_scalar_field('ocean_model', &
+ 'total_heat_content_srunoff', Time, diag, &
+ long_name='Area integrated heat content (relative to 0C) of submarine melt', &
+ units='W', conversion=US%QRZ_T_to_W_m2*US%L_to_m**2)
+
+ handles%id_total_heat_content_eff_srunoff = register_scalar_field('ocean_model', &
+ 'total_heat_content_eff_srunoff', Time, diag, &
+ long_name='Area integrated effective heat content (relative to 0C) of submarine melt', &
+ units='W', conversion=US%QRZ_T_to_W_m2*US%L_to_m**2)
+
if (present(use_glc_runoff)) then
handles%id_total_heat_content_frunoff_glc = register_scalar_field('ocean_model', &
'total_heat_content_frunoff_glc', Time, diag, &
@@ -2167,6 +2378,11 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
units='W', conversion=US%QRZ_T_to_W_m2*US%L_to_m**2) ! todo: update cmor names
endif
+ handles%id_total_lat_srunoff = register_scalar_field('ocean_model', &
+ 'total_lat_srunoff', Time, diag, &
+ long_name='Area integrated latent heat flux due to melting coupled submarine melt (srunoff)', &
+ units='W', conversion=US%QRZ_T_to_W_m2*US%L_to_m**2)
+
handles%id_total_sens = register_scalar_field('ocean_model', &
'total_sens', Time, diag, &
long_name='Area integrated downward sensible heat flux', &
@@ -2446,6 +2662,19 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
fluxes%salt_flux(i,j) = wt1*fluxes%salt_flux(i,j) + wt2*flux_tmp%salt_flux(i,j)
enddo ; enddo
+ ! Unlike the other water-group and heat fields above, fluxes%srunoff and fluxes%latent_srunoff
+ ! are only associated when submarine-melt coupling is actually configured, so they need their
+ ! own guards here.
+ if (associated(fluxes%srunoff) .and. associated(flux_tmp%srunoff)) then
+ do j=js,je ; do i=is,ie
+ fluxes%srunoff(i,j) = wt1*fluxes%srunoff(i,j) + wt2*flux_tmp%srunoff(i,j)
+ enddo ; enddo
+ endif
+ if (associated(fluxes%latent_srunoff) .and. associated(flux_tmp%latent_srunoff)) then
+ do j=js,je ; do i=is,ie
+ fluxes%latent_srunoff(i,j) = wt1*fluxes%latent_srunoff(i,j) + wt2*flux_tmp%latent_srunoff(i,j)
+ enddo ; enddo
+ endif
if (associated(fluxes%salt_flux_in) .and. associated(flux_tmp%salt_flux_in)) then
do j=js,je ; do i=is,ie
fluxes%salt_flux_in(i,j) = wt1*fluxes%salt_flux_in(i,j) + wt2*flux_tmp%salt_flux_in(i,j)
@@ -2504,6 +2733,17 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
fluxes%heat_content_frunoff(i,j) = wt1*fluxes%heat_content_frunoff(i,j) + wt2*flux_tmp%heat_content_frunoff(i,j)
enddo ; enddo
endif
+ if (associated(fluxes%heat_content_srunoff) .and. associated(flux_tmp%heat_content_srunoff)) then
+ do j=js,je ; do i=is,ie
+ fluxes%heat_content_srunoff(i,j) = wt1*fluxes%heat_content_srunoff(i,j) + wt2*flux_tmp%heat_content_srunoff(i,j)
+ enddo ; enddo
+ endif
+ if (associated(fluxes%heat_content_eff_srunoff) .and. associated(flux_tmp%heat_content_eff_srunoff)) then
+ do j=js,je ; do i=is,ie
+ fluxes%heat_content_eff_srunoff(i,j) = wt1*fluxes%heat_content_eff_srunoff(i,j) + &
+ wt2*flux_tmp%heat_content_eff_srunoff(i,j)
+ enddo ; enddo
+ endif
if (associated(fluxes%heat_content_lrunoff_glc) .and. associated(flux_tmp%heat_content_lrunoff_glc)) then
do j=js,je ; do i=is,ie
fluxes%heat_content_lrunoff_glc(i,j) = wt1*fluxes%heat_content_lrunoff_glc(i,j) + &
@@ -2687,6 +2927,9 @@ subroutine get_net_mass_forcing(fluxes, G, US, net_mass_src)
if (associated(fluxes%frunoff)) then ; do j=js,je ; do i=is,ie
net_mass_src(i,j) = net_mass_src(i,j) + fluxes%frunoff(i,j)
enddo ; enddo ; endif
+ if (associated(fluxes%srunoff)) then ; do j=js,je ; do i=is,ie
+ net_mass_src(i,j) = net_mass_src(i,j) + fluxes%srunoff(i,j)
+ enddo ; enddo ; endif
if (associated(fluxes%lrunoff_glc)) then ; do j=js,je ; do i=is,ie
net_mass_src(i,j) = net_mass_src(i,j) + fluxes%lrunoff_glc(i,j)
enddo ; enddo ; endif
@@ -2854,6 +3097,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
if (associated(fluxes%evap)) res(i,j) = res(i,j) + fluxes%evap(i,j)
if (associated(fluxes%lrunoff)) res(i,j) = res(i,j) + fluxes%lrunoff(i,j)
if (associated(fluxes%frunoff)) res(i,j) = res(i,j) + fluxes%frunoff(i,j)
+ if (associated(fluxes%srunoff)) res(i,j) = res(i,j) + fluxes%srunoff(i,j)
if (associated(fluxes%lrunoff_glc)) res(i,j) = res(i,j) + fluxes%lrunoff_glc(i,j)
if (associated(fluxes%frunoff_glc)) res(i,j) = res(i,j) + fluxes%frunoff_glc(i,j)
if (associated(fluxes%vprec)) res(i,j) = res(i,j) + fluxes%vprec(i,j)
@@ -2902,6 +3146,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
if (associated(fluxes%fprec)) res(i,j) = res(i,j) + fluxes%fprec(i,j)
if (associated(fluxes%lrunoff)) res(i,j) = res(i,j) + fluxes%lrunoff(i,j)
if (associated(fluxes%frunoff)) res(i,j) = res(i,j) + fluxes%frunoff(i,j)
+ if (associated(fluxes%srunoff)) res(i,j) = res(i,j) + fluxes%srunoff(i,j)
if (associated(fluxes%lrunoff_glc)) res(i,j) = res(i,j) + fluxes%lrunoff_glc(i,j)
if (associated(fluxes%frunoff_glc)) res(i,j) = res(i,j) + fluxes%frunoff_glc(i,j)
@@ -3023,6 +3268,14 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
endif
endif
+ if (associated(fluxes%srunoff)) then
+ if (handles%id_srunoff > 0) call post_data(handles%id_srunoff, fluxes%srunoff, diag)
+ if (handles%id_total_srunoff > 0) then
+ total_mass_flux = global_area_integral(fluxes%srunoff, G, tmp_scale=US%RZ_T_to_kg_m2s)
+ call post_data(handles%id_total_srunoff, total_mass_flux, diag)
+ endif
+ endif
+
if (associated(fluxes%seaice_melt)) then
if (handles%id_seaice_melt > 0) call post_data(handles%id_seaice_melt, fluxes%seaice_melt, diag)
if (handles%id_total_seaice_melt > 0) then
@@ -3061,6 +3314,20 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
call post_data(handles%id_total_heat_content_frunoff_glc, total_heat_flux, diag)
endif
+ if ((handles%id_heat_content_srunoff > 0) .and. associated(fluxes%heat_content_srunoff)) &
+ call post_data(handles%id_heat_content_srunoff, fluxes%heat_content_srunoff, diag)
+ if ((handles%id_total_heat_content_srunoff > 0) .and. associated(fluxes%heat_content_srunoff)) then
+ total_heat_flux = global_area_integral(fluxes%heat_content_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_total_heat_content_srunoff, total_heat_flux, diag)
+ endif
+
+ if ((handles%id_heat_content_eff_srunoff > 0) .and. associated(fluxes%heat_content_eff_srunoff)) &
+ call post_data(handles%id_heat_content_eff_srunoff, fluxes%heat_content_eff_srunoff, diag)
+ if ((handles%id_total_heat_content_eff_srunoff > 0) .and. associated(fluxes%heat_content_eff_srunoff)) then
+ total_heat_flux = global_area_integral(fluxes%heat_content_eff_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_total_heat_content_eff_srunoff, total_heat_flux, diag)
+ endif
+
if ((handles%id_heat_content_lprec > 0) .and. associated(fluxes%heat_content_lprec)) &
call post_data(handles%id_heat_content_lprec, fluxes%heat_content_lprec, diag)
if ((handles%id_total_heat_content_lprec > 0) .and. associated(fluxes%heat_content_lprec)) then
@@ -3116,6 +3383,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
res(i,j) = 0.0
if (associated(fluxes%LW)) res(i,j) = res(i,j) + fluxes%lw(i,j)
if (associated(fluxes%latent)) res(i,j) = res(i,j) + fluxes%latent(i,j)
+ if (associated(fluxes%latent_srunoff)) res(i,j) = res(i,j) + fluxes%latent_srunoff(i,j)
if (associated(fluxes%sens)) res(i,j) = res(i,j) + fluxes%sens(i,j)
if (associated(fluxes%SW)) res(i,j) = res(i,j) + fluxes%sw(i,j)
if (associated(fluxes%seaice_melt_heat)) res(i,j) = res(i,j) + fluxes%seaice_melt_heat(i,j)
@@ -3145,6 +3413,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j)
if (associated(fluxes%heat_content_frunoff)) &
res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j)
+ ! The effective heat content of submarine melt includes its latent contribution.
+ if (associated(fluxes%heat_content_eff_srunoff)) &
+ res(i,j) = res(i,j) + fluxes%heat_content_eff_srunoff(i,j)
if (associated(fluxes%heat_content_lrunoff_glc)) &
res(i,j) = res(i,j) + fluxes%heat_content_lrunoff_glc(i,j)
if (associated(fluxes%heat_content_frunoff_glc)) &
@@ -3183,6 +3454,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
res(i,j) = 0.0
if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j)
if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j)
+ if (associated(fluxes%heat_content_srunoff)) res(i,j) = res(i,j) + fluxes%heat_content_srunoff(i,j)
if (associated(fluxes%heat_content_lrunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff_glc(i,j)
if (associated(fluxes%heat_content_frunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff_glc(i,j)
if (associated(fluxes%heat_content_lprec)) res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j)
@@ -3208,6 +3480,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
res(i,j) = 0.0
if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j)
if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j)
+ if (associated(fluxes%heat_content_srunoff)) res(i,j) = res(i,j) + fluxes%heat_content_srunoff(i,j)
if (associated(fluxes%heat_content_lrunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff_glc(i,j)
if (associated(fluxes%heat_content_frunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff_glc(i,j)
enddo ; enddo
@@ -3230,6 +3503,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
do j=js,je ; do i=is,ie
res(i,j) = (fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)
enddo ; enddo
+ if (associated(fluxes%latent_srunoff)) then ; do j=js,je ; do i=is,ie
+ res(i,j) = res(i,j) + fluxes%latent_srunoff(i,j)
+ enddo ; enddo ; endif
call post_data(handles%id_LwLatSens, res, diag)
endif
@@ -3238,6 +3514,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
do j=js,je ; do i=is,ie
res(i,j) = (fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)
enddo ; enddo
+ if (associated(fluxes%latent_srunoff)) then ; do j=js,je ; do i=is,ie
+ res(i,j) = res(i,j) + fluxes%latent_srunoff(i,j)
+ enddo ; enddo ; endif
total_heat_flux = global_area_integral(res, G, tmp_scale=US%QRZ_T_to_W_m2)
call post_data(handles%id_total_LwLatSens, total_heat_flux, diag)
endif
@@ -3247,6 +3526,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
do j=js,je ; do i=is,ie
res(i,j) = ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j))
enddo ; enddo
+ if (associated(fluxes%latent_srunoff)) then ; do j=js,je ; do i=is,ie
+ res(i,j) = res(i,j) + fluxes%latent_srunoff(i,j)
+ enddo ; enddo ; endif
ave_heat_flux = global_area_mean(res, G, tmp_scale=US%QRZ_T_to_W_m2)
call post_data(handles%id_LwLatSens_ga, ave_heat_flux, diag)
endif
@@ -3283,16 +3565,33 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
call post_data(handles%id_lw_ga, ave_heat_flux, diag)
endif
- if ((handles%id_lat > 0) .and. associated(fluxes%latent)) then
- call post_data(handles%id_lat, fluxes%latent, diag)
- endif
- if ((handles%id_total_lat > 0) .and. associated(fluxes%latent)) then
- total_heat_flux = global_area_integral(fluxes%latent, G, tmp_scale=US%QRZ_T_to_W_m2)
- call post_data(handles%id_total_lat, total_heat_flux, diag)
- endif
- if ((handles%id_lat_ga > 0) .and. associated(fluxes%latent)) then
- ave_heat_flux = global_area_mean(fluxes%latent, G, tmp_scale=US%QRZ_T_to_W_m2)
- call post_data(handles%id_lat_ga, ave_heat_flux, diag)
+ if (associated(fluxes%latent_srunoff)) then
+ if ((handles%id_lat > 0) .or. (handles%id_total_lat > 0) .or. (handles%id_lat_ga > 0)) then
+ ! fluxes%latent is guaranteed to be associated here too, since latent_srunoff
+ ! is only allocated when heat is true, which is also required for fluxes%latent.
+ do j=js,je ; do i=is,ie
+ res(i,j) = fluxes%latent(i,j) + fluxes%latent_srunoff(i,j)
+ enddo ; enddo
+ if (handles%id_lat > 0) call post_data(handles%id_lat, res, diag)
+ if (handles%id_total_lat > 0) then
+ total_heat_flux = global_area_integral(res, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_total_lat, total_heat_flux, diag)
+ endif
+ if (handles%id_lat_ga > 0) then
+ ave_heat_flux = global_area_mean(res, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_lat_ga, ave_heat_flux, diag)
+ endif
+ endif
+ elseif (associated(fluxes%latent)) then
+ if (handles%id_lat > 0) call post_data(handles%id_lat, fluxes%latent, diag)
+ if (handles%id_total_lat > 0) then
+ total_heat_flux = global_area_integral(fluxes%latent, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_total_lat, total_heat_flux, diag)
+ endif
+ if (handles%id_lat_ga > 0) then
+ ave_heat_flux = global_area_mean(fluxes%latent, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_lat_ga, ave_heat_flux, diag)
+ endif
endif
if ((handles%id_lat_evap > 0) .and. associated(fluxes%latent_evap_diag)) then
@@ -3327,6 +3626,14 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
call post_data(handles%id_total_lat_frunoff_glc, total_heat_flux, diag)
endif
+ if ((handles%id_lat_srunoff > 0) .and. associated(fluxes%latent_srunoff)) then
+ call post_data(handles%id_lat_srunoff, fluxes%latent_srunoff, diag)
+ endif
+ if (handles%id_total_lat_srunoff > 0 .and. associated(fluxes%latent_srunoff)) then
+ total_heat_flux = global_area_integral(fluxes%latent_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
+ call post_data(handles%id_total_lat_srunoff, total_heat_flux, diag)
+ endif
+
if ((handles%id_sens > 0) .and. associated(fluxes%sens)) then
call post_data(handles%id_sens, fluxes%sens, diag)
endif
@@ -3455,7 +3762,7 @@ end subroutine forcing_diagnostics
subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
shelf, iceberg, salt, fix_accum_bug, fix_excess_salt_accum_bug, &
cfc, marbl, waves, shelf_sfc_accumulation, lamult, hevap, &
- ice_ncat, tau_mag)
+ ice_ncat, tau_mag, srunoff)
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
@@ -3483,6 +3790,8 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
!! via coupler.
integer, optional, intent(in) :: ice_ncat !< number of ice categories
logical, optional, intent(in) :: tau_mag !< If present and true, allocate tau_mag and related fields
+ logical, optional, intent(in) :: srunoff !< If present and true, allocate fields for
+ !! submarine melt
! Local variables
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
@@ -3512,6 +3821,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
call myAlloc(fluxes%vprec,isd,ied,jsd,jed, water)
call myAlloc(fluxes%lrunoff,isd,ied,jsd,jed, water)
call myAlloc(fluxes%frunoff,isd,ied,jsd,jed, water)
+ call myAlloc(fluxes%srunoff,isd,ied,jsd,jed, srunoff)
call myAlloc(fluxes%lrunoff_glc,isd,ied,jsd,jed, water)
call myAlloc(fluxes%frunoff_glc,isd,ied,jsd,jed, water)
call myAlloc(fluxes%seaice_melt,isd,ied,jsd,jed, water)
@@ -3543,6 +3853,12 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
call myAlloc(fluxes%heat_content_massin,isd,ied,jsd,jed, enthalpy_mom)
endif ; endif
+ if (present(heat) .and. present(srunoff)) then ; if (heat .and. srunoff) then
+ call myAlloc(fluxes%latent_srunoff,isd,ied,jsd,jed, .true.)
+ call myAlloc(fluxes%heat_content_srunoff,isd,ied,jsd,jed, .true.)
+ call myAlloc(fluxes%heat_content_eff_srunoff,isd,ied,jsd,jed, .true.)
+ endif ; endif
+
call myAlloc(fluxes%p_surf,isd,ied,jsd,jed, press)
! These fields should only be allocated if ice shelf is enabled.
@@ -3629,6 +3945,15 @@ subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes, turns)
call myAlloc(fluxes%buoy, G%isd, G%ied, G%jsd, G%jed, &
associated(fluxes_ref%buoy))
+ call myAlloc(fluxes%srunoff, G%isd, G%ied, G%jsd, G%jed, &
+ associated(fluxes_ref%srunoff))
+ call myAlloc(fluxes%heat_content_srunoff, G%isd, G%ied, G%jsd, G%jed, &
+ associated(fluxes_ref%heat_content_srunoff))
+ call myAlloc(fluxes%heat_content_eff_srunoff, G%isd, G%ied, G%jsd, G%jed, &
+ associated(fluxes_ref%heat_content_eff_srunoff))
+ call myAlloc(fluxes%latent_srunoff, G%isd, G%ied, G%jsd, G%jed, &
+ associated(fluxes_ref%latent_srunoff))
+
call myAlloc(fluxes%BBL_tidal_dis, G%isd, G%ied, G%jsd, G%jed, &
associated(fluxes_ref%BBL_tidal_dis))
call myAlloc(fluxes%ustar_tidal, G%isd, G%ied, G%jsd, G%jed, &
@@ -3842,10 +4167,13 @@ subroutine deallocate_forcing_type(fluxes)
if (associated(fluxes%latent_fprec_diag)) deallocate(fluxes%latent_fprec_diag)
if (associated(fluxes%latent_frunoff_diag)) deallocate(fluxes%latent_frunoff_diag)
if (associated(fluxes%latent_frunoff_glc_diag)) deallocate(fluxes%latent_frunoff_glc_diag)
+ if (associated(fluxes%latent_srunoff)) deallocate(fluxes%latent_srunoff)
if (associated(fluxes%sens)) deallocate(fluxes%sens)
if (associated(fluxes%heat_added)) deallocate(fluxes%heat_added)
if (associated(fluxes%heat_content_lrunoff)) deallocate(fluxes%heat_content_lrunoff)
if (associated(fluxes%heat_content_frunoff)) deallocate(fluxes%heat_content_frunoff)
+ if (associated(fluxes%heat_content_srunoff)) deallocate(fluxes%heat_content_srunoff)
+ if (associated(fluxes%heat_content_eff_srunoff)) deallocate(fluxes%heat_content_eff_srunoff)
if (associated(fluxes%heat_content_lrunoff_glc)) deallocate(fluxes%heat_content_lrunoff_glc)
if (associated(fluxes%heat_content_frunoff_glc)) deallocate(fluxes%heat_content_frunoff_glc)
if (associated(fluxes%heat_content_lprec)) deallocate(fluxes%heat_content_lprec)
@@ -3860,6 +4188,7 @@ subroutine deallocate_forcing_type(fluxes)
if (associated(fluxes%vprec)) deallocate(fluxes%vprec)
if (associated(fluxes%lrunoff)) deallocate(fluxes%lrunoff)
if (associated(fluxes%frunoff)) deallocate(fluxes%frunoff)
+ if (associated(fluxes%srunoff)) deallocate(fluxes%srunoff)
if (associated(fluxes%lrunoff_glc)) deallocate(fluxes%lrunoff_glc)
if (associated(fluxes%frunoff_glc)) deallocate(fluxes%frunoff_glc)
if (associated(fluxes%seaice_melt)) deallocate(fluxes%seaice_melt)
@@ -4008,6 +4337,15 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns)
endif
! The following fields are handled by drivers rather than control flags.
+ if (associated(fluxes_in%srunoff)) &
+ call rotate_array(fluxes_in%srunoff, turns, fluxes%srunoff)
+ if (associated(fluxes_in%heat_content_srunoff)) &
+ call rotate_array(fluxes_in%heat_content_srunoff, turns, fluxes%heat_content_srunoff)
+ if (associated(fluxes_in%heat_content_eff_srunoff)) &
+ call rotate_array(fluxes_in%heat_content_eff_srunoff, turns, fluxes%heat_content_eff_srunoff)
+ if (associated(fluxes_in%latent_srunoff)) &
+ call rotate_array(fluxes_in%latent_srunoff, turns, fluxes%latent_srunoff)
+
if (associated(fluxes_in%sw_vis_dir)) &
call rotate_array(fluxes_in%sw_vis_dir, turns, fluxes%sw_vis_dir)
if (associated(fluxes_in%sw_vis_dif)) &
@@ -4052,6 +4390,9 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns)
fluxes%fluxes_used = fluxes_in%fluxes_used
fluxes%dt_buoy_accum = fluxes_in%dt_buoy_accum
fluxes%C_p = fluxes_in%C_p
+ fluxes%latent_heat_fusion = fluxes_in%latent_heat_fusion
+ fluxes%srunoff_latent_heat = fluxes_in%srunoff_latent_heat
+ fluxes%srunoff_depth_max = fluxes_in%srunoff_depth_max
! NOTE: gustless_accum_bug is set during allocation
fluxes%num_msg = fluxes_in%num_msg
@@ -4224,6 +4565,10 @@ subroutine homogenize_forcing(fluxes, G, GV, US)
call homogenize_field_t(fluxes%vprec, G, tmp_scale=US%RZ_T_to_kg_m2s)
call homogenize_field_t(fluxes%lrunoff, G, tmp_scale=US%RZ_T_to_kg_m2s)
call homogenize_field_t(fluxes%frunoff, G, tmp_scale=US%RZ_T_to_kg_m2s)
+ ! Unlike the other water-group fields above, fluxes%srunoff is only associated when
+ ! submarine-melt coupling is actually configured, so it needs its own guard here.
+ if (associated(fluxes%srunoff)) &
+ call homogenize_field_t(fluxes%srunoff, G, tmp_scale=US%RZ_T_to_kg_m2s)
call homogenize_field_t(fluxes%lrunoff_glc, G, tmp_scale=US%RZ_T_to_kg_m2s)
call homogenize_field_t(fluxes%frunoff_glc, G, tmp_scale=US%RZ_T_to_kg_m2s)
call homogenize_field_t(fluxes%seaice_melt, G, tmp_scale=US%RZ_T_to_kg_m2s)
@@ -4247,6 +4592,11 @@ subroutine homogenize_forcing(fluxes, G, GV, US)
call homogenize_field_t(fluxes%latent_frunoff_glc_diag, G, tmp_scale=US%QRZ_T_to_W_m2)
endif
+ ! Unlike the other heat fluxes above, fluxes%latent_srunoff is only associated when
+ ! submarine-melt coupling is actually configured, so it needs its own guard here.
+ if (associated(fluxes%latent_srunoff)) &
+ call homogenize_field_t(fluxes%latent_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
+
if (do_salt) call homogenize_field_t(fluxes%salt_flux, G, tmp_scale=US%RZ_T_to_kg_m2s)
if (do_heat .and. do_water) then
@@ -4256,6 +4606,13 @@ subroutine homogenize_forcing(fluxes, G, GV, US)
call homogenize_field_t(fluxes%heat_content_vprec, G, tmp_scale=US%QRZ_T_to_W_m2)
call homogenize_field_t(fluxes%heat_content_lrunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
call homogenize_field_t(fluxes%heat_content_frunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
+ ! Unlike the other heat_content fields above, fluxes%heat_content_srunoff and
+ ! fluxes%heat_content_eff_srunoff are only associated when submarine-melt coupling is
+ ! actually configured, so they need their own guards here.
+ if (associated(fluxes%heat_content_srunoff)) &
+ call homogenize_field_t(fluxes%heat_content_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
+ if (associated(fluxes%heat_content_eff_srunoff)) &
+ call homogenize_field_t(fluxes%heat_content_eff_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2)
call homogenize_field_t(fluxes%heat_content_lrunoff_glc, G, tmp_scale=US%QRZ_T_to_W_m2)
call homogenize_field_t(fluxes%heat_content_frunoff_glc, G, tmp_scale=US%QRZ_T_to_W_m2)
call homogenize_field_t(fluxes%heat_content_massout, G, tmp_scale=US%QRZ_T_to_W_m2)
diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90
index fba33fe0ee..c241aff51c 100644
--- a/src/diagnostics/MOM_sum_output.F90
+++ b/src/diagnostics/MOM_sum_output.F90
@@ -1061,6 +1061,13 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
(((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + (fluxes%lrunoff(i,j) + fluxes%lrunoff_glc(i,j))) + &
(fluxes%fprec(i,j) + (fluxes%frunoff(i,j) + fluxes%frunoff_glc(i,j)))))
enddo ; enddo
+ ! fluxes%srunoff is only associated when submarine-melt coupling is actually configured,
+ ! unlike the other water-group fields above, so it needs its own guard here.
+ if (associated(fluxes%srunoff)) then
+ do j=js,je ; do i=is,ie
+ FW_in(i,j) = FW_in(i,j) + dt*G%areaT(i,j)*fluxes%srunoff(i,j)
+ enddo ; enddo
+ endif
else
call MOM_error(WARNING, &
"accumulate_net_input called with associated evap field, but no precip field.")
@@ -1102,6 +1109,14 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
fluxes%heat_content_cond(i,j) + fluxes%heat_content_fprec(i,j) + &
fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j))
enddo ; enddo
+ ! fluxes%heat_content_eff_srunoff is only associated when submarine-melt coupling is actually
+ ! configured, unlike the other heat_content fields above, so it needs its own guard here.
+ ! The effective heat content of submarine melt includes its latent contribution.
+ if (associated(fluxes%heat_content_eff_srunoff)) then
+ do j=js,je ; do i=is,ie
+ heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * fluxes%heat_content_eff_srunoff(i,j)
+ enddo ; enddo
+ endif
elseif (associated(tv%TempxPmE)) then
do j=js,je ; do i=is,ie
heat_in(i,j) = heat_in(i,j) + (tv%C_p * G%areaT(i,j)) * tv%TempxPmE(i,j)
@@ -1112,6 +1127,15 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
enddo ; enddo
endif
+ ! Without coupler enthalpies, only the sensible heat content of submarine melt is included above,
+ ! so add the rest of the heat it applies (e.g. the latent heat of melting).
+ if (.not. associated(fluxes%heat_content_evap)) then
+ if (associated(fluxes%heat_content_eff_srunoff)) then ; do j=js,je ; do i=is,ie
+ heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * &
+ (fluxes%heat_content_eff_srunoff(i,j) - fluxes%heat_content_srunoff(i,j))
+ enddo ; enddo ; endif
+ endif
+
! The following heat sources may or may not be used.
if (associated(tv%internal_heat)) then
do j=js,je ; do i=is,ie
diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90
index 25c5fdef14..1f1eaa664c 100644
--- a/src/parameterizations/vertical/MOM_diabatic_aux.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90
@@ -15,7 +15,7 @@ module MOM_diabatic_aux
use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
-use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint
+use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint, distribute_srunoff
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : thickness_to_dz
use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
@@ -762,7 +762,12 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
mixing_depth, & ! The mixing depth for brine plumes [H ~> m or kg m-2]
- total_h ! Total thickness of the water column [H ~> m or kg m-2]
+ total_h, & ! Total thickness of the water column [H ~> m or kg m-2]
+ dHeat_srunoff, & ! The heat added to each column by srunoff, from distribute_srunoff, for
+ ! the heat_content_eff_srunoff diagnostic [Q H ~> J m-2]
+ dTempxPmE_srunoff ! The ambient-temperature-weighted mass added to each column by srunoff,
+ ! from distribute_srunoff, for the tv%TempxPmE and heat_content_srunoff
+ ! diagnostics [C H ~> degC m or degC kg m-2]
real, dimension(SZI_(G), SZK_(GV)) :: &
h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
! dz, & ! Layer thicknesses in depth units [Z ~> m]
@@ -869,7 +874,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
!$OMP drhodt,drhods,pen_sw_bnd_rate, &
!$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, &
!$OMP mixing_depth,A_brine,fraction_left_brine, &
- !$OMP plume_fraction,dK,total_h) &
+ !$OMP plume_fraction,dK,total_h,dHeat_srunoff, &
+ !$OMP dTempxPmE_srunoff) &
!$OMP firstprivate(SurfPressure,plume_flux)
do j=js,je
! Work in vertical slices for efficiency
@@ -990,6 +996,42 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
endif
enddo
+ if (associated(fluxes%srunoff)) then
+ ! Add srunoff's contribution to netmassinout_rate/net_heat_rate, used for the surface
+ ! buoyancy flux. This means that for the purposes of the surface buoyancy flux, the mass
+ ! and heat associated with srunoff are treated as if they were applied at the surface,
+ ! analogous to how penetrating SW is handled. srunoff enters at the ambient temperature, so
+ ! only its latent heat of melting changes the temperature and is included in net_heat_rate.
+ if (calculate_buoyancy) then
+ do i=is,ie
+ if ((G%mask2dT(i,j) > 0.0) .and. (fluxes%srunoff(i,j) > 0.0)) then
+ netmassinout_rate(i) = netmassinout_rate(i) + GV%RZ_to_H * fluxes%srunoff(i,j)
+ if (associated(fluxes%latent_srunoff)) &
+ netheat_rate(i) = netheat_rate(i) + fluxes%latent_srunoff(i,j) * (GV%RZ_to_H / tv%C_p)
+ endif
+ enddo ! i
+ endif
+
+ ! Distribute srunoff's mass and heat over a range of depths, if requested,
+ ! updating temperature and salinity together. This is done before loops A and B
+ ! below so loop B sees the srunoff-fattened column, making grounding less likely.
+ call distribute_srunoff(G, GV, dt, fluxes, j, h2d, tv%S(:,j,:), T2d=T2d, C_p=tv%C_p, US=US, &
+ g_Hconv2=g_Hconv2, cTKE=cTKE, dSV_dT=dSV_dT, dSV_dS=dSV_dS, &
+ dHeat_total=dHeat_srunoff, dTempxPmE_total=dTempxPmE_srunoff)
+
+ ! Add srunoff's contributions to relevant diagnostics. heat_content_eff_srunoff is all of
+ ! the heat applied (i.e. heat content plus e.g. latent heat), while heat_content_srunoff is
+ ! only the heat content of the melt water at the ambient temperatures.
+ do i=is,ie
+ if (associated(fluxes%heat_content_eff_srunoff)) &
+ fluxes%heat_content_eff_srunoff(i,j) = dHeat_srunoff(i) * GV%H_to_RZ / dt
+ if (associated(fluxes%heat_content_srunoff)) &
+ fluxes%heat_content_srunoff(i,j) = tv%C_p * dTempxPmE_srunoff(i) * GV%H_to_RZ / dt
+ if (associated(tv%TempxPmE)) &
+ tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + dTempxPmE_srunoff(i) * GV%H_to_RZ
+ enddo
+ endif
+
! Apply the surface boundary fluxes in three steps:
! A/ update mass, temp, and salinity due to all terms except mass leaving
! ocean (and corresponding outward heat content), and ignoring penetrative SW.
@@ -1240,7 +1282,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
.false., .true., T2d, Pen_SW_bnd)
endif
-
! Step D/ copy updated thickness and temperature
! 2d slice now back into model state.
do k=1,nz ; do i=is,ie
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 5389de3e1d..287751c18c 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -333,6 +333,16 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
"diabatic was called with a negative timestep.")
+ if (associated(fluxes%srunoff)) then
+ ! KPP code does not account for submarine-melt contributions to the surface buoyancy fluxes.
+ if (CS%useKPP) call MOM_error(FATAL, &
+ "MOM_diabatic_driver: submarine melt (srunoff) coupling is not supported with USE_KPP=True.")
+
+ ! Distribution of submarine melt does not handle negative values, which is not physical anyway.
+ if (any(fluxes%srunoff(G%isc:G%iec,G%jsc:G%jec) < 0.0)) call MOM_error(FATAL, &
+ "MOM_diabatic_driver: coupled submarine melt (srunoff) contains negative values.")
+ endif
+
showCallTree = callTree_showQuery()
! Offer diagnostics of various state variables at the start of diabatic
@@ -1999,6 +2009,10 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
showCallTree = callTree_showQuery()
if (showCallTree) call callTree_enter("layered_diabatic(), MOM_diabatic_driver.F90")
+ if (associated(fluxes%srunoff)) call MOM_error(FATAL, &
+ "layered_diabatic: submarine melt (srunoff) coupling requires the ALE algorithm; "//&
+ "it is not supported with the layered (non-ALE) diabatic driver.")
+
! set equivalence between the same bits of memory for these arrays
eaml => eatr ; ebml => ebtr
@@ -3334,7 +3348,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
"and applied as either incoming or outgoing depending on the sign of the net. "//&
"If false, the net incoming fresh water flux is added to the model and "//&
"thereafter the net outgoing is removed from the topmost non-vanished "//&
- "layers of the updated state.", default=.true.)
+ "layers of the updated state. This does not apply to submarine melt (srunoff) "//&
+ "which is handled separately.", default=.true.)
call get_param(param_file, mdl, "DEBUG", CS%debug, &
"If true, write out verbose debugging data.", &
diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90
index ae6d98e3a7..32b30784d4 100644
--- a/src/tracer/MOM_tracer_diabatic.F90
+++ b/src/tracer/MOM_tracer_diabatic.F90
@@ -9,7 +9,7 @@ module MOM_tracer_diabatic
use MOM_grid, only : ocean_grid_type
use MOM_verticalGrid, only : verticalGrid_type
-use MOM_forcing_type, only : forcing
+use MOM_forcing_type, only : forcing, distribute_srunoff
use MOM_error_handler, only : MOM_error, FATAL, WARNING
implicit none ; private
@@ -530,6 +530,11 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, evap_CFL_lim
netMassIn(i) = fluxes%netMassIn(i,j)
enddo
+ ! Distribute submarine melt (srunoff) mass over a range of depths, if requested, diluting
+ ! this tracer accordingly. This is done before loops A and B below so loop B sees the
+ ! srunoff-fattened column, making grounding less likely.
+ if (associated(fluxes%srunoff)) call distribute_srunoff(G, GV, dt, fluxes, j, h2d, Tr2d)
+
! Apply the surface boundary fluxes in three steps:
! A/ update concentration from mass entering the ocean
! B/ update concentration from mass leaving ocean.