From 37a92fd6e5a8264d33430c8682fc689a2e2d1782 Mon Sep 17 00:00:00 2001 From: Dougie Squire Date: Wed, 5 Aug 2026 09:43:55 +1000 Subject: [PATCH 01/11] Add optional submarine melt coupling --- config_src/drivers/nuopc_cap/mom_cap.F90 | 29 +++ .../drivers/nuopc_cap/mom_cap_methods.F90 | 7 + .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 29 +++ src/core/MOM_forcing_type.F90 | 169 ++++++++++++++++-- src/diagnostics/MOM_sum_output.F90 | 7 +- 5 files changed, 223 insertions(+), 18 deletions(-) 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_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index c83a47d498..803c499761 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -86,6 +86,11 @@ 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 :: 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 +171,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] @@ -330,6 +336,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, 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) + fluxes%latent_heat_fusion = CS%latent_heat_fusion + fluxes%srunoff_latent_heat = CS%srunoff_latent_heat 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 +511,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 +574,10 @@ 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. Instead, the latent heat + ! extraction from the ocean due to submarine melt is accounted for as a Gade (1979)-line effective + ! temperature offset via fluxes%heat_content_srunoff in MOM_forcing_type.F90, so that the same + ! code path can be extended to include depth-distribution of submarine melt in the future. 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 @@ -1216,6 +1237,11 @@ 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, "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 +1607,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..74a17d82e7 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -121,7 +121,8 @@ module MOM_forcing_type latent_evap_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from evaporating liquid water (typically < 0) latent_fprec_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting fprec (typically < 0) latent_frunoff_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting frunoff (calving) (typically < 0) - latent_frunoff_glc_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from melting glacier frunoff (typically < 0) + latent_frunoff_glc_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting glacier frunoff (typically < 0) + latent_srunoff_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from submarine melt (srunoff) (typically < 0) ! water mass fluxes into the ocean [R Z T-1 ~> kg m-2 s-1]; these fluxes impact the ocean mass real, pointer, dimension(:,:) :: & @@ -131,6 +132,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 +153,7 @@ 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_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 +228,13 @@ 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. ! arrays needed in the some tracer modules, e.g., MOM_CFC_cap real, pointer, dimension(:,:) :: & @@ -336,6 +346,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 +357,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 +373,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 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 +391,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 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 @@ -636,10 +651,11 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & ! net volume/mass of liquid and solid passing through surface boundary fluxes netMassInOut(i) = dt * (scale * & - (((((((( fluxes%lprec(i,j) & + ((((((((( fluxes%lprec(i,j) & + fluxes%fprec(i,j) ) & + fluxes%evap(i,j) ) & + fluxes%lrunoff(i,j) ) & + + fluxes%srunoff(i,j) ) & + fluxes%lrunoff_glc(i,j)) & + fluxes%vprec(i,j) ) & + fluxes%seaice_melt(i,j)) & @@ -648,10 +664,11 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (do_NMIOr) then ! Repeat the above code without multiplying by a timestep for legacy reasons netMassInOut_rate(i) = (scale * & - (((((((( fluxes%lprec(i,j) & + ((((((((( fluxes%lprec(i,j) & + fluxes%fprec(i,j) ) & + fluxes%evap(i,j) ) & + fluxes%lrunoff(i,j) ) & + + fluxes%srunoff(i,j) ) & + fluxes%lrunoff_glc(i,j)) & + fluxes%vprec(i,j) ) & + fluxes%seaice_melt(i,j)) & @@ -747,6 +764,32 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & endif endif + ! Add explicit heat flux for submarine melt (srunoff). Submarine melt is added at the ambient + ! temperature which is not available at depth in the coupler, so there is no option for a + ! coupler-supplied heat content. The heat content computed here optionally includes the latent + ! heat to melt the ice, accounted for by an effective temperature offset following + ! the Gade (1979) meltwater mixing line. + if (associated(fluxes%srunoff) .and. associated(fluxes%heat_content_srunoff)) then + if (fluxes%srunoff_latent_heat) then + fluxes%heat_content_srunoff(i,j) = tv%C_p*fluxes%srunoff(i,j)*T(i,1) - & + fluxes%srunoff(i,j)*fluxes%latent_heat_fusion*US%J_kg_to_Q + else + fluxes%heat_content_srunoff(i,j) = tv%C_p*fluxes%srunoff(i,j)*T(i,1) + endif + + net_heat(i) = net_heat(i) + (scale*(dt * I_Cp_Hconvert)) * fluxes%heat_content_srunoff(i,j) + + if (do_enthalpy) then + ! remove srunoff*SST here, to counteract its addition in applyBoundaryFluxesInOut + net_heat(i) = net_heat(i) - (GV%RZ_to_H * (scale * dt)) * fluxes%srunoff(i,j) * T(i,1) + + if (calculate_diags .and. associated(tv%TempxPmE)) then + tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * & + (I_Cp*fluxes%heat_content_srunoff(i,j) - fluxes%srunoff(i,j)*T(i,1)) + endif + endif + endif + ! Add explicit heat flux for calving (which is part of the ice-ocean boundary ! flux type). Calving is otherwise added with a temperature of SST. ! This should only be done when do_enthalpy is true. @@ -949,8 +992,9 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (associated(tv%TempxPmE)) then tv%TempxPmE(i,j) = (I_Cp*dt*scale) * & - ((((fluxes%heat_content_lprec(i,j) + fluxes%heat_content_fprec(i,j)) + & + (((((fluxes%heat_content_lprec(i,j) + fluxes%heat_content_fprec(i,j)) + & (fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j))) + & + fluxes%heat_content_srunoff(i,j)) + & (fluxes%heat_content_lrunoff_glc(i,j) + fluxes%heat_content_frunoff_glc(i,j))) + & (fluxes%heat_content_evap(i,j) + fluxes%heat_content_cond(i,j))) endif @@ -1363,6 +1407,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_diag)) & + call hchksum(fluxes%latent_srunoff_diag, mesg//" fluxes%latent_srunoff_diag", 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 +1449,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 +1463,9 @@ 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_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 +1569,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_diag,'latent_srunoff_diag') call locMsg(fluxes%sens,'sens') call locMsg(fluxes%evap,'evap') call locMsg(fluxes%lprec,'lprec') @@ -1532,10 +1585,12 @@ 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_lprec,'heat_content_lprec') call locMsg(fluxes%heat_content_fprec,'heat_content_fprec') call locMsg(fluxes%heat_content_vprec,'heat_content_vprec') @@ -1709,6 +1764,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', & @@ -1803,6 +1862,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 +1936,10 @@ 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, 'Effective heat content (relative to 0C) of submarine melt into ocean', & + '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', & @@ -1992,6 +2059,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 +2102,11 @@ 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) + 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 +2243,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', & @@ -2432,6 +2513,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%vprec(i,j) = wt1*fluxes%vprec(i,j) + wt2*flux_tmp%vprec(i,j) fluxes%lrunoff(i,j) = wt1*fluxes%lrunoff(i,j) + wt2*flux_tmp%lrunoff(i,j) fluxes%frunoff(i,j) = wt1*fluxes%frunoff(i,j) + wt2*flux_tmp%frunoff(i,j) + fluxes%srunoff(i,j) = wt1*fluxes%srunoff(i,j) + wt2*flux_tmp%srunoff(i,j) fluxes%lrunoff_glc(i,j) = wt1*fluxes%lrunoff_glc(i,j) + wt2*flux_tmp%lrunoff_glc(i,j) fluxes%frunoff_glc(i,j) = wt1*fluxes%frunoff_glc(i,j) + wt2*flux_tmp%frunoff_glc(i,j) fluxes%seaice_melt(i,j) = wt1*fluxes%seaice_melt(i,j) + wt2*flux_tmp%seaice_melt(i,j) @@ -2504,6 +2586,11 @@ 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_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 +2774,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 +2944,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 +2993,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 +3115,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 +3161,13 @@ 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_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 @@ -3145,6 +3252,8 @@ 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) + 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)) & @@ -3183,6 +3292,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 +3318,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 @@ -3283,16 +3394,22 @@ 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 ((handles%id_lat > 0) .or. (handles%id_total_lat > 0) .or. (handles%id_lat_ga > 0)) then + do j=js,je ; do i=is,ie + res(i,j) = 0.0 + if (associated(fluxes%latent)) res(i,j) = res(i,j) + fluxes%latent(i,j) + ! Include latent heat from srunoff since it is handled separately from fluxes%latent + if (associated(fluxes%latent_srunoff_diag)) res(i,j) = res(i,j) + fluxes%latent_srunoff_diag(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 if ((handles%id_lat_evap > 0) .and. associated(fluxes%latent_evap_diag)) then @@ -3327,6 +3444,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_diag)) then + call post_data(handles%id_lat_srunoff, fluxes%latent_srunoff_diag, diag) + endif + if (handles%id_total_lat_srunoff > 0 .and. associated(fluxes%latent_srunoff_diag)) then + total_heat_flux = global_area_integral(fluxes%latent_srunoff_diag, 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 @@ -3512,6 +3637,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, water) 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) @@ -3526,6 +3652,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%latent_fprec_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_frunoff_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_frunoff_glc_diag,isd,ied,jsd,jed, heat) + call myAlloc(fluxes%latent_srunoff_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%salt_flux,isd,ied,jsd,jed, salt) @@ -3537,6 +3664,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%heat_content_vprec,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_lrunoff,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff,isd,ied,jsd,jed, .true.) + call myAlloc(fluxes%heat_content_srunoff,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_lrunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_massout,isd,ied,jsd,jed, enthalpy_mom) @@ -3842,10 +3970,12 @@ 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_diag)) deallocate(fluxes%latent_srunoff_diag) 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_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 +3990,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) @@ -3945,6 +4076,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%vprec, turns, fluxes%vprec) call rotate_array(fluxes_in%lrunoff, turns, fluxes%lrunoff) call rotate_array(fluxes_in%frunoff, turns, fluxes%frunoff) + call rotate_array(fluxes_in%srunoff, turns, fluxes%srunoff) call rotate_array(fluxes_in%lrunoff_glc, turns, fluxes%lrunoff_glc) call rotate_array(fluxes_in%frunoff_glc, turns, fluxes%frunoff_glc) call rotate_array(fluxes_in%seaice_melt, turns, fluxes%seaice_melt) @@ -3962,6 +4094,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%latent_fprec_diag, turns, fluxes%latent_fprec_diag) call rotate_array(fluxes_in%latent_frunoff_diag, turns, fluxes%latent_frunoff_diag) call rotate_array(fluxes_in%latent_frunoff_glc_diag, turns, fluxes%latent_frunoff_glc_diag) + call rotate_array(fluxes_in%latent_srunoff_diag, turns, fluxes%latent_srunoff_diag) endif if (do_salt) then @@ -3977,6 +4110,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%heat_content_lrunoff_glc, turns, fluxes%heat_content_lrunoff_glc) call rotate_array(fluxes_in%heat_content_frunoff, turns, fluxes%heat_content_frunoff) call rotate_array(fluxes_in%heat_content_frunoff_glc, turns, fluxes%heat_content_frunoff_glc) + call rotate_array(fluxes_in%heat_content_srunoff, turns, fluxes%heat_content_srunoff) if (associated (fluxes_in%heat_content_evap)) then call rotate_array(fluxes_in%heat_content_evap, turns, fluxes%heat_content_evap) else @@ -4052,6 +4186,8 @@ 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 ! NOTE: gustless_accum_bug is set during allocation fluxes%num_msg = fluxes_in%num_msg @@ -4224,6 +4360,7 @@ 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) + 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) @@ -4245,6 +4382,7 @@ subroutine homogenize_forcing(fluxes, G, GV, US) call homogenize_field_t(fluxes%latent_fprec_diag, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%latent_frunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%latent_frunoff_glc_diag, G, tmp_scale=US%QRZ_T_to_W_m2) + call homogenize_field_t(fluxes%latent_srunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) endif if (do_salt) call homogenize_field_t(fluxes%salt_flux, G, tmp_scale=US%RZ_T_to_kg_m2s) @@ -4256,6 +4394,7 @@ 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) + call homogenize_field_t(fluxes%heat_content_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..9f103ad681 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -1058,8 +1058,8 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then do j=js,je ; do i=is,ie FW_in(i,j) = dt*G%areaT(i,j)*(fluxes%evap(i,j) + & - (((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))))) + ((((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)))) + fluxes%srunoff(i,j))) enddo ; enddo else call MOM_error(WARNING, & @@ -1100,7 +1100,8 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & (fluxes%heat_content_evap(i,j) + fluxes%heat_content_lprec(i,j) + & fluxes%heat_content_cond(i,j) + fluxes%heat_content_fprec(i,j) + & - fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j)) + fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j) + & + fluxes%heat_content_srunoff(i,j)) enddo ; enddo elseif (associated(tv%TempxPmE)) then do j=js,je ; do i=is,ie From 16f468bf3ae043380c927ee66901aa6a71497840 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Wed, 19 Aug 2026 08:47:23 +1000 Subject: [PATCH 02/11] First attempt at submarine melt depth distribution --- src/core/MOM_forcing_type.F90 | 35 +---- .../vertical/MOM_diabatic_aux.F90 | 144 +++++++++++++++++- .../vertical/MOM_diabatic_driver.F90 | 8 +- 3 files changed, 153 insertions(+), 34 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 74a17d82e7..c836799a5f 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -651,11 +651,10 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & ! net volume/mass of liquid and solid passing through surface boundary fluxes netMassInOut(i) = dt * (scale * & - ((((((((( fluxes%lprec(i,j) & + (((((((( fluxes%lprec(i,j) & + fluxes%fprec(i,j) ) & + fluxes%evap(i,j) ) & + fluxes%lrunoff(i,j) ) & - + fluxes%srunoff(i,j) ) & + fluxes%lrunoff_glc(i,j)) & + fluxes%vprec(i,j) ) & + fluxes%seaice_melt(i,j)) & @@ -664,11 +663,10 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (do_NMIOr) then ! Repeat the above code without multiplying by a timestep for legacy reasons netMassInOut_rate(i) = (scale * & - ((((((((( fluxes%lprec(i,j) & + (((((((( fluxes%lprec(i,j) & + fluxes%fprec(i,j) ) & + fluxes%evap(i,j) ) & + fluxes%lrunoff(i,j) ) & - + fluxes%srunoff(i,j) ) & + fluxes%lrunoff_glc(i,j)) & + fluxes%vprec(i,j) ) & + fluxes%seaice_melt(i,j)) & @@ -764,32 +762,6 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & endif endif - ! Add explicit heat flux for submarine melt (srunoff). Submarine melt is added at the ambient - ! temperature which is not available at depth in the coupler, so there is no option for a - ! coupler-supplied heat content. The heat content computed here optionally includes the latent - ! heat to melt the ice, accounted for by an effective temperature offset following - ! the Gade (1979) meltwater mixing line. - if (associated(fluxes%srunoff) .and. associated(fluxes%heat_content_srunoff)) then - if (fluxes%srunoff_latent_heat) then - fluxes%heat_content_srunoff(i,j) = tv%C_p*fluxes%srunoff(i,j)*T(i,1) - & - fluxes%srunoff(i,j)*fluxes%latent_heat_fusion*US%J_kg_to_Q - else - fluxes%heat_content_srunoff(i,j) = tv%C_p*fluxes%srunoff(i,j)*T(i,1) - endif - - net_heat(i) = net_heat(i) + (scale*(dt * I_Cp_Hconvert)) * fluxes%heat_content_srunoff(i,j) - - if (do_enthalpy) then - ! remove srunoff*SST here, to counteract its addition in applyBoundaryFluxesInOut - net_heat(i) = net_heat(i) - (GV%RZ_to_H * (scale * dt)) * fluxes%srunoff(i,j) * T(i,1) - - if (calculate_diags .and. associated(tv%TempxPmE)) then - tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * & - (I_Cp*fluxes%heat_content_srunoff(i,j) - fluxes%srunoff(i,j)*T(i,1)) - endif - endif - endif - ! Add explicit heat flux for calving (which is part of the ice-ocean boundary ! flux type). Calving is otherwise added with a temperature of SST. ! This should only be done when do_enthalpy is true. @@ -992,9 +964,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (associated(tv%TempxPmE)) then tv%TempxPmE(i,j) = (I_Cp*dt*scale) * & - (((((fluxes%heat_content_lprec(i,j) + fluxes%heat_content_fprec(i,j)) + & + ((((fluxes%heat_content_lprec(i,j) + fluxes%heat_content_fprec(i,j)) + & (fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j))) + & - fluxes%heat_content_srunoff(i,j)) + & (fluxes%heat_content_lrunoff_glc(i,j) + fluxes%heat_content_frunoff_glc(i,j))) + & (fluxes%heat_content_evap(i,j) + fluxes%heat_content_cond(i,j))) endif diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 25c5fdef14..c2b4038d81 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -46,6 +46,10 @@ module MOM_diabatic_aux logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff at the !! river mouths to a depth of "rivermix_depth" real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m]. + real :: srunoff_depth_distribute = 0.0 !< The depth over which submarine melt (srunoff) mass and + !! heat are distributed, spread uniformly by thickness, + !! instead of being applied entirely within the top layer + !! [H ~> m or kg m-2]. real :: dSalt_frac_max !< An upper limit on the fraction of the salt in a layer that can be !! lost to the net surface salt fluxes within a timestep [nondim] logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to @@ -1212,6 +1216,10 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t enddo ! i + ! Distribute submarine melt (srunoff) mass and heat over a range of depths, if requested. + call distribute_srunoff(CS, G, GV, US, dt, fluxes, j, h2d, T2d, tv, & + calculate_energetics, g_Hconv2, cTKE, dSV_dT, dSV_dS) + ! Step C/ in the application of fluxes ! Heat by the convergence of penetrating SW. ! SW penetrative heating uses the updated thickness from above. @@ -1240,7 +1248,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 @@ -1353,6 +1360,135 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t end subroutine applyBoundaryFluxesInOut +!> Distribute the mass and heat content of submarine melt (srunoff) uniformly by thickness +!! over the upper CS%srunoff_depth_distribute of a column. If CS%srunoff_depth_distribute +!! 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. Called once +!! per j-row from applyBoundaryFluxesInOut, alongside absorbRemainingSW, and operating on +!! the same 2-D h2d/T2d slices. +!! +!! Because srunoff is applied here rather than through netMassInOut/netMassOut above, its +!! mass is not included in the fractionOfForcing/evap_CFL_limit calculation that Loop B +!! (above) applies to the other surface mass fluxes when AGGREGATE_FW_FORCING is true. That +!! calculation is not depth-aware regardless, so including srunoff in it would not make it +!! correctly account for srunoff's eventual depth placement -- it would just add srunoff at +!! the surface for the purpose of that calculation before removing and redistributing it +!! here, which does not improve on the current behavior. This can very slightly alter how +!! much evaporation-driven forcing gets CFL-limited or spread downward at srunoff-active +!! points (i.e. near ice shelves) relative to a version of the model without this coupling, +!! and in the extreme could very slightly change the likelihood of a column grounding out +!! (netMassIn+netMassOut left over after the k-loop above) at such points. distribute_srunoff +!! itself cannot ground out or need CFL-limiting, since it only ever adds mass to a column. +subroutine distribute_srunoff(CS, G, GV, US, dt, fluxes, j, h2d, T2d, tv, & + calculate_energetics, g_Hconv2, cTKE, dSV_dT, dSV_dS) + type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: dt !< Time step [T ~> s] + type(forcing), intent(inout) :: 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) :: T2d !< Layer temperatures [C ~> degC] + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variables + logical, intent(in) :: calculate_energetics !< If true, calculate the energy + !! required to mix the srunoff added to a layer + !! throughout that layer + real, 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] + + ! 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 :: 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] + integer :: i, k, is, ie, nz + + if (.not.associated(fluxes%srunoff)) return + + is = G%isc ; ie = G%iec ; nz = GV%ke + + do i=is,ie ; if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then + + dM_tot = GV%RZ_to_H * dt * fluxes%srunoff(i,j) + + ! When CS%srunoff_depth_distribute 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, + ! as it was before this depth-distribution capability existed. Floor by GV%H_subroundoff + ! to avoid dividing by zero if the top layer has vanished; the loop below then naturally + ! rolls the deposit into the first layer down with non-negligible thickness instead. + depth = CS%srunoff_depth_distribute + 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, spread by thickness over the + ! whole column, rather than requiring a separate step to redistribute a shortfall. + total_h = 0.0 + do k=1,nz ; total_h = total_h + h2d(i,k) ; enddo + eff_depth = min(depth, total_h) + I_Habs = 1.0 / eff_depth + + h_above = 0.0 + dHeat = 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, mixing in each layer at that layer's own + ! local temperature (optionally offset by the Gade (1979) latent heat adjustment), the + ! same way Loop A above does for surface-only mass sources, just repeated over depth + ! instead of only at the top layer. + 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) + T_eff = T2d(i,k) + if (fluxes%srunoff_latent_heat) & + T_eff = T_eff - fluxes%latent_heat_fusion*US%J_kg_to_Q/tv%C_p + dHeat = dHeat + tv%C_p*dM_k*T_eff + + hOld = h2d(i,k) + if (calculate_energetics) then + ! Calculate the energy required to mix the srunoff added to this layer throughout + ! the layer, mirroring the analogous calculation in Loop A above for other incoming + ! mass sources. Srunoff is fresh water, so its salinity plays the role of Salin_in = 0 + ! there. + 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) + tv%S(i,j,k) * dSV_dS(i,j,k)) + endif + h2d(i,k) = hOld + dM_k + T2d(i,k) = (hOld*T2d(i,k) + dM_k*T_eff) / h2d(i,k) + tv%S(i,j,k) = (hOld*tv%S(i,j,k)) / h2d(i,k) ! srunoff is fresh water + endif + enddo + + ! Report the total heat added to the column, converted back to a flux, purely for + ! diagnostic purposes -- this heat has already been applied above, spread over depth, + ! not just at the surface. + if (associated(fluxes%heat_content_srunoff)) & + fluxes%heat_content_srunoff(i,j) = dHeat * GV%H_to_RZ / dt + + endif ; enddo + +end subroutine distribute_srunoff + !> This subroutine initializes the parameters and control structure of the diabatic_aux module. subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL) type(time_type), target, intent(in) :: Time !< The current model time. @@ -1436,6 +1572,12 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori CS%do_rivermix = .false. ; CS%rivermix_depth = 0.0 ; CS%ignore_fluxes_over_land = .false. endif + call get_param(param_file, mdl, "SRUNOFF_DEPTH_DISTRIBUTE", CS%srunoff_depth_distribute, & + "The depth over which submarine melt (srunoff) mass and heat are spread, "//& + "distributed uniformly by thickness, instead of being applied entirely "//& + "within the top layer. If 0, srunoff is applied entirely within the top layer.", & + units="m", default=0.0, scale=GV%m_to_H, do_not_log=.not.use_temperature) + if (GV%nkml == 0) then call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", CS%use_river_heat_content, & "If true, use the fluxes%runoff_Hflx field to set the "//& diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 5389de3e1d..3a36f2673e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1999,6 +1999,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 +3338,9 @@ 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 always distributed by distribute_srunoff (MOM_diabatic_aux.F90), "//& + "independent of this setting.", default=.true.) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & From 1bfad762089572e3298b41eb5e3aaaebb932ed5c Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Thu, 20 Aug 2026 13:58:25 +1000 Subject: [PATCH 03/11] Include dilution of other tracers from srunoff distribution --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 16 +- src/core/MOM_forcing_type.F90 | 215 ++++++++++++++++-- src/diagnostics/MOM_sum_output.F90 | 21 +- .../vertical/MOM_diabatic_aux.F90 | 190 ++++------------ .../vertical/MOM_diabatic_driver.F90 | 15 +- src/tracer/MOM_tracer_diabatic.F90 | 9 +- 6 files changed, 291 insertions(+), 175 deletions(-) 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 803c499761..2f58f81d0e 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -91,6 +91,9 @@ module MOM_surface_forcing_nuopc !! meltwater mixing line, so that melting the ice consumes latent !! heat. If false, srunoff enters the ocean at the ambient !! temperature. + real :: srunoff_depth !< The 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 [m]. 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. @@ -335,9 +338,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 = CS%srunoff_depth 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) @@ -574,10 +578,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. Instead, the latent heat - ! extraction from the ocean due to submarine melt is accounted for as a Gade (1979)-line effective - ! temperature offset via fluxes%heat_content_srunoff in MOM_forcing_type.F90, so that the same - ! code path can be extended to include depth-distribution of submarine melt in the future. + ! 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 @@ -1242,6 +1244,10 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "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", CS%srunoff_depth, & + "The 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) 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 "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index c836799a5f..e20811b892 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 @@ -235,6 +235,10 @@ module MOM_forcing_type !! (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 = 0.0 !< The 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 in surface_forcing_CS [m]. ! arrays needed in the some tracer modules, e.g., MOM_CFC_cap real, pointer, dimension(:,:) :: & @@ -1045,6 +1049,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 of a column. If fluxes%srunoff_depth 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.) .and. (fluxes%srunoff(i,j) /= 0.0)) then + + dM_tot = GV%RZ_to_H * dt * fluxes%srunoff(i,j) + + ! When fluxes%srunoff_depth 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 = GV%m_to_H * fluxes%srunoff_depth + 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 + + 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 @@ -1679,7 +1824,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') @@ -1760,11 +1905,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) !========================================================================= @@ -1965,7 +2112,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', & @@ -2484,7 +2632,6 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%vprec(i,j) = wt1*fluxes%vprec(i,j) + wt2*flux_tmp%vprec(i,j) fluxes%lrunoff(i,j) = wt1*fluxes%lrunoff(i,j) + wt2*flux_tmp%lrunoff(i,j) fluxes%frunoff(i,j) = wt1*fluxes%frunoff(i,j) + wt2*flux_tmp%frunoff(i,j) - fluxes%srunoff(i,j) = wt1*fluxes%srunoff(i,j) + wt2*flux_tmp%srunoff(i,j) fluxes%lrunoff_glc(i,j) = wt1*fluxes%lrunoff_glc(i,j) + wt2*flux_tmp%lrunoff_glc(i,j) fluxes%frunoff_glc(i,j) = wt1*fluxes%frunoff_glc(i,j) + wt2*flux_tmp%frunoff_glc(i,j) fluxes%seaice_melt(i,j) = wt1*fluxes%seaice_melt(i,j) + wt2*flux_tmp%seaice_melt(i,j) @@ -2499,6 +2646,13 @@ 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 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) .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%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) @@ -3551,7 +3705,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 @@ -3579,6 +3733,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 @@ -3608,7 +3764,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, 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) @@ -3623,7 +3779,10 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%latent_fprec_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_frunoff_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_frunoff_glc_diag,isd,ied,jsd,jed, heat) - call myAlloc(fluxes%latent_srunoff_diag,isd,ied,jsd,jed, heat) + + if (present(heat) .and. present(srunoff)) then ; if (heat .and. srunoff) then + call myAlloc(fluxes%latent_srunoff_diag,isd,ied,jsd,jed, .true.) + endif ; endif call myAlloc(fluxes%salt_flux,isd,ied,jsd,jed, salt) @@ -3635,7 +3794,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%heat_content_vprec,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_lrunoff,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff,isd,ied,jsd,jed, .true.) - call myAlloc(fluxes%heat_content_srunoff,isd,ied,jsd,jed, .true.) + call myAlloc(fluxes%heat_content_srunoff,isd,ied,jsd,jed, srunoff) call myAlloc(fluxes%heat_content_lrunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_massout,isd,ied,jsd,jed, enthalpy_mom) @@ -3728,6 +3887,13 @@ 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%latent_srunoff_diag, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%latent_srunoff_diag)) + 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, & @@ -3941,7 +4107,7 @@ 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_diag)) deallocate(fluxes%latent_srunoff_diag) + if (associated(fluxes%latent_srunoff_diag)) deallocate(fluxes%latent_srunoff_diag) 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) @@ -4047,7 +4213,6 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%vprec, turns, fluxes%vprec) call rotate_array(fluxes_in%lrunoff, turns, fluxes%lrunoff) call rotate_array(fluxes_in%frunoff, turns, fluxes%frunoff) - call rotate_array(fluxes_in%srunoff, turns, fluxes%srunoff) call rotate_array(fluxes_in%lrunoff_glc, turns, fluxes%lrunoff_glc) call rotate_array(fluxes_in%frunoff_glc, turns, fluxes%frunoff_glc) call rotate_array(fluxes_in%seaice_melt, turns, fluxes%seaice_melt) @@ -4065,7 +4230,6 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%latent_fprec_diag, turns, fluxes%latent_fprec_diag) call rotate_array(fluxes_in%latent_frunoff_diag, turns, fluxes%latent_frunoff_diag) call rotate_array(fluxes_in%latent_frunoff_glc_diag, turns, fluxes%latent_frunoff_glc_diag) - call rotate_array(fluxes_in%latent_srunoff_diag, turns, fluxes%latent_srunoff_diag) endif if (do_salt) then @@ -4081,7 +4245,6 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%heat_content_lrunoff_glc, turns, fluxes%heat_content_lrunoff_glc) call rotate_array(fluxes_in%heat_content_frunoff, turns, fluxes%heat_content_frunoff) call rotate_array(fluxes_in%heat_content_frunoff_glc, turns, fluxes%heat_content_frunoff_glc) - call rotate_array(fluxes_in%heat_content_srunoff, turns, fluxes%heat_content_srunoff) if (associated (fluxes_in%heat_content_evap)) then call rotate_array(fluxes_in%heat_content_evap, turns, fluxes%heat_content_evap) else @@ -4113,6 +4276,13 @@ 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%latent_srunoff_diag)) & + call rotate_array(fluxes_in%latent_srunoff_diag, turns, fluxes%latent_srunoff_diag) + 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)) & @@ -4159,6 +4329,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) 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 = fluxes_in%srunoff_depth ! NOTE: gustless_accum_bug is set during allocation fluxes%num_msg = fluxes_in%num_msg @@ -4331,7 +4502,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) - call homogenize_field_t(fluxes%srunoff, 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) @@ -4353,9 +4527,13 @@ subroutine homogenize_forcing(fluxes, G, GV, US) call homogenize_field_t(fluxes%latent_fprec_diag, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%latent_frunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%latent_frunoff_glc_diag, G, tmp_scale=US%QRZ_T_to_W_m2) - call homogenize_field_t(fluxes%latent_srunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) endif + ! Unlike the other latent_*_diag fields above, fluxes%latent_srunoff_diag is only associated + ! when submarine-melt coupling is actually configured, so it needs its own guard here. + if (associated(fluxes%latent_srunoff_diag)) & + call homogenize_field_t(fluxes%latent_srunoff_diag, 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 @@ -4365,7 +4543,10 @@ 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) - call homogenize_field_t(fluxes%heat_content_srunoff, G, tmp_scale=US%QRZ_T_to_W_m2) + ! Unlike the other heat_content fields above, fluxes%heat_content_srunoff is only associated + ! when submarine-melt coupling is actually configured, so it needs its own guard here. + if (associated(fluxes%heat_content_srunoff)) & + call homogenize_field_t(fluxes%heat_content_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 9f103ad681..c08644fa54 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -1058,9 +1058,16 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then do j=js,je ; do i=is,ie FW_in(i,j) = dt*G%areaT(i,j)*(fluxes%evap(i,j) + & - ((((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)))) + fluxes%srunoff(i,j))) + (((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.") @@ -1100,9 +1107,15 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & (fluxes%heat_content_evap(i,j) + fluxes%heat_content_lprec(i,j) + & fluxes%heat_content_cond(i,j) + fluxes%heat_content_fprec(i,j) + & - fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j) + & - fluxes%heat_content_srunoff(i,j)) + fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j)) enddo ; enddo + ! fluxes%heat_content_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. + if (associated(fluxes%heat_content_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_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) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index c2b4038d81..c598e9653f 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 @@ -46,10 +46,6 @@ module MOM_diabatic_aux logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff at the !! river mouths to a depth of "rivermix_depth" real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m]. - real :: srunoff_depth_distribute = 0.0 !< The depth over which submarine melt (srunoff) mass and - !! heat are distributed, spread uniformly by thickness, - !! instead of being applied entirely within the top layer - !! [H ~> m or kg m-2]. real :: dSalt_frac_max !< An upper limit on the fraction of the salt in a layer that can be !! lost to the net surface salt fluxes within a timestep [nondim] logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to @@ -730,6 +726,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings) real :: H_limit_fluxes ! Surface fluxes are scaled down fluxes when the total depth of the ocean ! drops below this value [H ~> m or kg m-2] + real :: srunoff_Teff ! The effective temperature carried by srunoff's contribution to + ! net_heat_rate [C ~> degC] real :: IforcingDepthScale ! The inverse of the layer thickness below which mass losses are ! shifted to the next deeper layer [H ~> m or kg m-2] real :: Idt ! The inverse of the timestep [T-1 ~> s-1] @@ -766,7 +764,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_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 diagnostic + ! [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] @@ -873,7 +876,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,srunoff_Teff) & !$OMP firstprivate(SurfPressure,plume_flux) do j=js,je ! Work in vertical slices for efficiency @@ -1216,9 +1220,40 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t enddo ! i - ! Distribute submarine melt (srunoff) mass and heat over a range of depths, if requested. - call distribute_srunoff(CS, G, GV, US, dt, fluxes, j, h2d, T2d, tv, & - calculate_energetics, g_Hconv2, cTKE, dSV_dT, dSV_dS) + 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 penetrating SW. + if (calculate_buoyancy) then + do i=is,ie + if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then + netmassinout_rate(i) = netmassinout_rate(i) + GV%RZ_to_H * fluxes%srunoff(i,j) + srunoff_Teff = T2d(i,1) + if (fluxes%srunoff_latent_heat) & + srunoff_Teff = srunoff_Teff - fluxes%latent_heat_fusion*US%J_kg_to_Q/tv%C_p + netheat_rate(i) = netheat_rate(i) + (GV%RZ_to_H * fluxes%srunoff(i,j)) * srunoff_Teff + endif + enddo ! i + endif + + ! Distribute submarine melt (srunoff) mass and heat over a range of depths, if requested, + ! updating temperature and salinity together. + ! Because srunoff is applied after loop B above, it's possible that groundings could be + ! reported when in fact they would not have occured had srunoff been accounted for. It + ! might be better to do this before loop A and B? + 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 contributions to relevant diagnostics. + do i=is,ie + if (associated(fluxes%heat_content_srunoff)) & + fluxes%heat_content_srunoff(i,j) = dHeat_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 ! Step C/ in the application of fluxes ! Heat by the convergence of penetrating SW. @@ -1360,135 +1395,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t end subroutine applyBoundaryFluxesInOut -!> Distribute the mass and heat content of submarine melt (srunoff) uniformly by thickness -!! over the upper CS%srunoff_depth_distribute of a column. If CS%srunoff_depth_distribute -!! 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. Called once -!! per j-row from applyBoundaryFluxesInOut, alongside absorbRemainingSW, and operating on -!! the same 2-D h2d/T2d slices. -!! -!! Because srunoff is applied here rather than through netMassInOut/netMassOut above, its -!! mass is not included in the fractionOfForcing/evap_CFL_limit calculation that Loop B -!! (above) applies to the other surface mass fluxes when AGGREGATE_FW_FORCING is true. That -!! calculation is not depth-aware regardless, so including srunoff in it would not make it -!! correctly account for srunoff's eventual depth placement -- it would just add srunoff at -!! the surface for the purpose of that calculation before removing and redistributing it -!! here, which does not improve on the current behavior. This can very slightly alter how -!! much evaporation-driven forcing gets CFL-limited or spread downward at srunoff-active -!! points (i.e. near ice shelves) relative to a version of the model without this coupling, -!! and in the extreme could very slightly change the likelihood of a column grounding out -!! (netMassIn+netMassOut left over after the k-loop above) at such points. distribute_srunoff -!! itself cannot ground out or need CFL-limiting, since it only ever adds mass to a column. -subroutine distribute_srunoff(CS, G, GV, US, dt, fluxes, j, h2d, T2d, tv, & - calculate_energetics, g_Hconv2, cTKE, dSV_dT, dSV_dS) - type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, intent(in) :: dt !< Time step [T ~> s] - type(forcing), intent(inout) :: 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) :: T2d !< Layer temperatures [C ~> degC] - type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variables - logical, intent(in) :: calculate_energetics !< If true, calculate the energy - !! required to mix the srunoff added to a layer - !! throughout that layer - real, 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] - - ! 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 :: 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] - integer :: i, k, is, ie, nz - - if (.not.associated(fluxes%srunoff)) return - - is = G%isc ; ie = G%iec ; nz = GV%ke - - do i=is,ie ; if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then - - dM_tot = GV%RZ_to_H * dt * fluxes%srunoff(i,j) - - ! When CS%srunoff_depth_distribute 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, - ! as it was before this depth-distribution capability existed. Floor by GV%H_subroundoff - ! to avoid dividing by zero if the top layer has vanished; the loop below then naturally - ! rolls the deposit into the first layer down with non-negligible thickness instead. - depth = CS%srunoff_depth_distribute - 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, spread by thickness over the - ! whole column, rather than requiring a separate step to redistribute a shortfall. - total_h = 0.0 - do k=1,nz ; total_h = total_h + h2d(i,k) ; enddo - eff_depth = min(depth, total_h) - I_Habs = 1.0 / eff_depth - - h_above = 0.0 - dHeat = 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, mixing in each layer at that layer's own - ! local temperature (optionally offset by the Gade (1979) latent heat adjustment), the - ! same way Loop A above does for surface-only mass sources, just repeated over depth - ! instead of only at the top layer. - 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) - T_eff = T2d(i,k) - if (fluxes%srunoff_latent_heat) & - T_eff = T_eff - fluxes%latent_heat_fusion*US%J_kg_to_Q/tv%C_p - dHeat = dHeat + tv%C_p*dM_k*T_eff - - hOld = h2d(i,k) - if (calculate_energetics) then - ! Calculate the energy required to mix the srunoff added to this layer throughout - ! the layer, mirroring the analogous calculation in Loop A above for other incoming - ! mass sources. Srunoff is fresh water, so its salinity plays the role of Salin_in = 0 - ! there. - 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) + tv%S(i,j,k) * dSV_dS(i,j,k)) - endif - h2d(i,k) = hOld + dM_k - T2d(i,k) = (hOld*T2d(i,k) + dM_k*T_eff) / h2d(i,k) - tv%S(i,j,k) = (hOld*tv%S(i,j,k)) / h2d(i,k) ! srunoff is fresh water - endif - enddo - - ! Report the total heat added to the column, converted back to a flux, purely for - ! diagnostic purposes -- this heat has already been applied above, spread over depth, - ! not just at the surface. - if (associated(fluxes%heat_content_srunoff)) & - fluxes%heat_content_srunoff(i,j) = dHeat * GV%H_to_RZ / dt - - endif ; enddo - -end subroutine distribute_srunoff - !> This subroutine initializes the parameters and control structure of the diabatic_aux module. subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL) type(time_type), target, intent(in) :: Time !< The current model time. @@ -1572,12 +1478,6 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori CS%do_rivermix = .false. ; CS%rivermix_depth = 0.0 ; CS%ignore_fluxes_over_land = .false. endif - call get_param(param_file, mdl, "SRUNOFF_DEPTH_DISTRIBUTE", CS%srunoff_depth_distribute, & - "The depth over which submarine melt (srunoff) mass and heat are spread, "//& - "distributed uniformly by thickness, instead of being applied entirely "//& - "within the top layer. If 0, srunoff is applied entirely within the top layer.", & - units="m", default=0.0, scale=GV%m_to_H, do_not_log=.not.use_temperature) - if (GV%nkml == 0) then call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", CS%use_river_heat_content, & "If true, use the fluxes%runoff_Hflx field to set the "//& diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 3a36f2673e..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 @@ -3338,9 +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. This does not apply to submarine melt (srunoff), "//& - "which is always distributed by distribute_srunoff (MOM_diabatic_aux.F90), "//& - "independent of this setting.", 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..6ac9459b10 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 @@ -613,6 +613,13 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, evap_CFL_lim enddo ! i + ! Distribute submarine melt (srunoff) mass over a range of depths, if requested, diluting + ! this tracer according. + ! Because srunoff is applied after loop B above, it's possible that groundings could be + ! reported when in fact they would not have occured had srunoff been accounted for. It + ! might be better to do this before loop A and B? + if (associated(fluxes%srunoff)) call distribute_srunoff(G, GV, dt, fluxes, j, h2d, Tr2d) + ! Step C/ copy updated tracer concentration from the 2d slice now back into model state. do k=1,nz ; do i=is,ie Tr(i,j,k) = Tr2d(i,k) From dbf257e7a100bc77918b639b1aa84a95d9559512 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Thu, 20 Aug 2026 21:57:55 +1000 Subject: [PATCH 04/11] Move srunoff distribution to before application of surface boundary fluxes --- .../vertical/MOM_diabatic_aux.F90 | 68 +++++++++---------- src/tracer/MOM_tracer_diabatic.F90 | 12 ++-- 2 files changed, 38 insertions(+), 42 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index c598e9653f..02c6d9c747 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -998,6 +998,39 @@ 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. + if (calculate_buoyancy) then + do i=is,ie + if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then + netmassinout_rate(i) = netmassinout_rate(i) + GV%RZ_to_H * fluxes%srunoff(i,j) + srunoff_Teff = T2d(i,1) + if (fluxes%srunoff_latent_heat) & + srunoff_Teff = srunoff_Teff - fluxes%latent_heat_fusion*US%J_kg_to_Q/tv%C_p + netheat_rate(i) = netheat_rate(i) + (GV%RZ_to_H * fluxes%srunoff(i,j)) * srunoff_Teff + 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. + do i=is,ie + if (associated(fluxes%heat_content_srunoff)) & + fluxes%heat_content_srunoff(i,j) = dHeat_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. @@ -1220,41 +1253,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t enddo ! i - 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 penetrating SW. - if (calculate_buoyancy) then - do i=is,ie - if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then - netmassinout_rate(i) = netmassinout_rate(i) + GV%RZ_to_H * fluxes%srunoff(i,j) - srunoff_Teff = T2d(i,1) - if (fluxes%srunoff_latent_heat) & - srunoff_Teff = srunoff_Teff - fluxes%latent_heat_fusion*US%J_kg_to_Q/tv%C_p - netheat_rate(i) = netheat_rate(i) + (GV%RZ_to_H * fluxes%srunoff(i,j)) * srunoff_Teff - endif - enddo ! i - endif - - ! Distribute submarine melt (srunoff) mass and heat over a range of depths, if requested, - ! updating temperature and salinity together. - ! Because srunoff is applied after loop B above, it's possible that groundings could be - ! reported when in fact they would not have occured had srunoff been accounted for. It - ! might be better to do this before loop A and B? - 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 contributions to relevant diagnostics. - do i=is,ie - if (associated(fluxes%heat_content_srunoff)) & - fluxes%heat_content_srunoff(i,j) = dHeat_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 - ! Step C/ in the application of fluxes ! Heat by the convergence of penetrating SW. ! SW penetrative heating uses the updated thickness from above. diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index 6ac9459b10..32b30784d4 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -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. @@ -613,13 +618,6 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, evap_CFL_lim enddo ! i - ! Distribute submarine melt (srunoff) mass over a range of depths, if requested, diluting - ! this tracer according. - ! Because srunoff is applied after loop B above, it's possible that groundings could be - ! reported when in fact they would not have occured had srunoff been accounted for. It - ! might be better to do this before loop A and B? - if (associated(fluxes%srunoff)) call distribute_srunoff(G, GV, dt, fluxes, j, h2d, Tr2d) - ! Step C/ copy updated tracer concentration from the 2d slice now back into model state. do k=1,nz ; do i=is,ie Tr(i,j,k) = Tr2d(i,k) From f6487e09e113ca15da84d2a0098bf17fb3ee68e0 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Tue, 29 Sep 2026 12:58:44 +1000 Subject: [PATCH 05/11] Don't post 'latent' diagnostics when neither fluxes%latent nor latent_srunoff_diag are associated --- src/core/MOM_forcing_type.F90 | 33 ++++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index e20811b892..0638c7f7eb 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1132,7 +1132,7 @@ subroutine distribute_srunoff(G, GV, dt, fluxes, j, h2d, C2d, T2d, C_p, US, & if (present(dHeat_total)) dHeat_total(i) = 0.0 if (present(dTempxPmE_total)) dTempxPmE_total(i) = 0.0 - if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then + 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) @@ -3519,20 +3519,31 @@ 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) .or. (handles%id_total_lat > 0) .or. (handles%id_lat_ga > 0)) then - do j=js,je ; do i=is,ie - res(i,j) = 0.0 - if (associated(fluxes%latent)) res(i,j) = res(i,j) + fluxes%latent(i,j) - ! Include latent heat from srunoff since it is handled separately from fluxes%latent - if (associated(fluxes%latent_srunoff_diag)) res(i,j) = res(i,j) + fluxes%latent_srunoff_diag(i,j) - enddo ; enddo - if (handles%id_lat > 0) call post_data(handles%id_lat, res, diag) + if (associated(fluxes%latent_srunoff_diag)) 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_diag + ! 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_diag(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(res, G, tmp_scale=US%QRZ_T_to_W_m2) + 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(res, G, tmp_scale=US%QRZ_T_to_W_m2) + 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 From af7f2aa40ed6d07ae7a7987526ea3a1c73782a53 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Tue, 29 Sep 2026 16:41:49 +1000 Subject: [PATCH 06/11] Handle scaling of SRUNOFF_DEPTH correctly --- .../drivers/nuopc_cap/mom_ocean_model_nuopc.F90 | 2 +- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 8 +++++--- src/core/MOM_forcing_type.F90 | 16 ++++++++-------- .../vertical/MOM_diabatic_aux.F90 | 2 +- 4 files changed, 15 insertions(+), 13 deletions(-) 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 2f58f81d0e..1a116fd9b8 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 @@ -93,7 +94,7 @@ module MOM_surface_forcing_nuopc !! temperature. real :: srunoff_depth !< The 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 [m]. + !! 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. @@ -1172,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 @@ -1247,7 +1249,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call get_param(param_file, mdl, "SRUNOFF_DEPTH", CS%srunoff_depth, & "The 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) + "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 "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 0638c7f7eb..0e4a175c7f 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -238,7 +238,7 @@ module MOM_forcing_type real :: srunoff_depth = 0.0 !< The 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 in surface_forcing_CS [m]. + !! CS%srunoff_depth 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(:,:) :: & @@ -1138,7 +1138,7 @@ subroutine distribute_srunoff(G, GV, dt, fluxes, j, h2d, C2d, T2d, C_p, US, & ! When fluxes%srunoff_depth 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 = GV%m_to_H * fluxes%srunoff_depth + depth = fluxes%srunoff_depth 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 @@ -1181,7 +1181,7 @@ subroutine distribute_srunoff(G, GV, dt, fluxes, j, h2d, C2d, T2d, C_p, US, & 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 + enddo ! k-loop for layer distribution if (present(dHeat_total)) dHeat_total(i) = dHeat if (present(dTempxPmE_total)) dTempxPmE_total(i) = dTempxPmE @@ -3791,10 +3791,6 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%latent_frunoff_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_frunoff_glc_diag,isd,ied,jsd,jed, heat) - if (present(heat) .and. present(srunoff)) then ; if (heat .and. srunoff) then - call myAlloc(fluxes%latent_srunoff_diag,isd,ied,jsd,jed, .true.) - endif ; endif - call myAlloc(fluxes%salt_flux,isd,ied,jsd,jed, salt) if (present(heat) .and. present(water)) then ; if (heat .and. water) then @@ -3805,13 +3801,17 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%heat_content_vprec,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_lrunoff,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff,isd,ied,jsd,jed, .true.) - call myAlloc(fluxes%heat_content_srunoff,isd,ied,jsd,jed, srunoff) call myAlloc(fluxes%heat_content_lrunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_massout,isd,ied,jsd,jed, enthalpy_mom) 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_diag,isd,ied,jsd,jed, .true.) + call myAlloc(fluxes%heat_content_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. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 02c6d9c747..e155813d21 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1005,7 +1005,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! analogous to how penetrating SW is handled. if (calculate_buoyancy) then do i=is,ie - if ((G%mask2dT(i,j) > 0.) .and. (fluxes%srunoff(i,j) /= 0.0)) then + 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) srunoff_Teff = T2d(i,1) if (fluxes%srunoff_latent_heat) & From 28517d3975a1ab6164f499f40be5e3698f867aee Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Tue, 27 Oct 2026 20:33:50 +1100 Subject: [PATCH 07/11] Rename SRUNOFF_DEPTH to SRUNOFF_DEPTH_MAX --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 12 ++++++------ src/core/MOM_forcing_type.F90 | 18 +++++++++--------- 2 files changed, 15 insertions(+), 15 deletions(-) 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 1a116fd9b8..cb660ea38c 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -92,9 +92,9 @@ module MOM_surface_forcing_nuopc !! meltwater mixing line, so that melting the ice consumes latent !! heat. If false, srunoff enters the ocean at the ambient !! temperature. - real :: srunoff_depth !< The 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 :: 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. @@ -342,7 +342,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, 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 = CS%srunoff_depth + 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) @@ -1246,8 +1246,8 @@ subroutine surface_forcing_init(Time, G, GV, US, param_file, diag, CS, restore_s "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", CS%srunoff_depth, & - "The depth over which the mass and heat from coupled submarine melt are spread, "//& + 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, & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 0e4a175c7f..b012a94289 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -235,10 +235,10 @@ module MOM_forcing_type !! (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 = 0.0 !< The 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 in surface_forcing_CS [H ~> m or kg m-2]. + 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(:,:) :: & @@ -1050,8 +1050,8 @@ 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 of a column. If fluxes%srunoff_depth is 0, the -!! distribution depth for a column instead falls back to that column's own top layer thickness, +!! 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. @@ -1136,9 +1136,9 @@ subroutine distribute_srunoff(G, GV, dt, fluxes, j, h2d, C2d, T2d, C_p, US, & dM_tot = GV%RZ_to_H * dt * fluxes%srunoff(i,j) - ! When fluxes%srunoff_depth is 0 (the default), fall back to this column's own top layer + ! 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 + 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 @@ -4340,7 +4340,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) 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 = fluxes_in%srunoff_depth + fluxes%srunoff_depth_max = fluxes_in%srunoff_depth_max ! NOTE: gustless_accum_bug is set during allocation fluxes%num_msg = fluxes_in%num_msg From 1ad4075cf4563b01ad9e20ba637f0077a794c766 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Wed, 28 Oct 2026 10:28:22 +1100 Subject: [PATCH 08/11] Fix double counting of srunoff latent heat diags --- src/core/MOM_forcing_type.F90 | 52 ++++++++++++++++--- src/diagnostics/MOM_sum_output.F90 | 16 ++++-- .../vertical/MOM_diabatic_aux.F90 | 14 +++-- 3 files changed, 67 insertions(+), 15 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index b012a94289..8f440c8a0a 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -154,6 +154,8 @@ module MOM_forcing_type 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] @@ -380,7 +382,7 @@ module MOM_forcing_type 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 + 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 @@ -399,7 +401,7 @@ module MOM_forcing_type 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 + 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 @@ -1582,6 +1584,9 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) 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) @@ -1707,6 +1712,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) 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') @@ -2055,9 +2061,13 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 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, 'Effective heat content (relative to 0C) of submarine melt into ocean', & + 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', & @@ -2226,6 +2236,11 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 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, & @@ -2716,6 +2731,12 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) 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) + & @@ -3293,6 +3314,13 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h 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 @@ -3377,8 +3405,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) - if (associated(fluxes%heat_content_srunoff)) & - res(i,j) = res(i,j) + fluxes%heat_content_srunoff(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)) & @@ -3810,6 +3839,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & if (present(heat) .and. present(srunoff)) then ; if (heat .and. srunoff) then call myAlloc(fluxes%latent_srunoff_diag,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) @@ -3902,6 +3932,8 @@ subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes, turns) 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_diag, G%isd, G%ied, G%jsd, G%jed, & associated(fluxes_ref%latent_srunoff_diag)) @@ -4124,6 +4156,7 @@ subroutine deallocate_forcing_type(fluxes) 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) @@ -4291,6 +4324,8 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) 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_diag)) & call rotate_array(fluxes_in%latent_srunoff_diag, turns, fluxes%latent_srunoff_diag) @@ -4554,10 +4589,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 is only associated - ! when submarine-melt coupling is actually configured, so it needs its own guard here. + ! 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 c08644fa54..c241aff51c 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -1109,11 +1109,12 @@ 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_srunoff is only associated when submarine-melt coupling is actually + ! 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. - if (associated(fluxes%heat_content_srunoff)) then + ! 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_srunoff(i,j) + 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 @@ -1126,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 e155813d21..0547803a0d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -766,10 +766,10 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t 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] dHeat_srunoff, & ! The heat added to each column by srunoff, from distribute_srunoff, for - ! the heat_content_srunoff diagnostic [Q H ~> J m-2] + ! 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 diagnostic - ! [C H ~> degC m or degC kg m-2] + ! 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] @@ -1022,10 +1022,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t 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. + ! 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) = dHeat_srunoff(i) * GV%H_to_RZ / dt + 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 From 257ed35ef7634fcf5d52c1687ea46faa10ca640f Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Wed, 28 Oct 2026 11:15:46 +1100 Subject: [PATCH 09/11] Fix srunoff contribution to net_heat_rate --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 2 +- src/core/MOM_forcing_type.F90 | 57 +++++++++++-------- .../vertical/MOM_diabatic_aux.F90 | 13 ++--- 3 files changed, 38 insertions(+), 34 deletions(-) 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 cb660ea38c..395d7510a0 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -521,7 +521,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, 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) * & + 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 diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 8f440c8a0a..23052421d8 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -114,15 +114,16 @@ 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(:,:) :: & latent_evap_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from evaporating liquid water (typically < 0) latent_fprec_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting fprec (typically < 0) latent_frunoff_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting frunoff (calving) (typically < 0) - latent_frunoff_glc_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting glacier frunoff (typically < 0) - latent_srunoff_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from submarine melt (srunoff) (typically < 0) + latent_frunoff_glc_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from melting glacier frunoff (typically < 0) ! water mass fluxes into the ocean [R Z T-1 ~> kg m-2 s-1]; these fluxes impact the ocean mass real, pointer, dimension(:,:) :: & @@ -1525,8 +1526,8 @@ 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_diag)) & - call hchksum(fluxes%latent_srunoff_diag, mesg//" fluxes%latent_srunoff_diag", G%HI, & + 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, & @@ -1690,7 +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_diag,'latent_srunoff_diag') + call locMsg(fluxes%latent_srunoff,'latent_srunoff') call locMsg(fluxes%sens,'sens') call locMsg(fluxes%evap,'evap') call locMsg(fluxes%lprec,'lprec') @@ -2661,13 +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 fields above, fluxes%srunoff is only associated when - ! submarine-melt coupling is actually configured, so it needs its own guard here. + ! 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) @@ -3548,12 +3555,12 @@ 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 (associated(fluxes%latent_srunoff_diag)) then + 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_diag + ! 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_diag(i,j) + 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 @@ -3609,11 +3616,11 @@ 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_diag)) then - call post_data(handles%id_lat_srunoff, fluxes%latent_srunoff_diag, diag) + 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_diag)) then - total_heat_flux = global_area_integral(fluxes%latent_srunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) + 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 @@ -3837,7 +3844,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & endif ; endif if (present(heat) .and. present(srunoff)) then ; if (heat .and. srunoff) then - call myAlloc(fluxes%latent_srunoff_diag,isd,ied,jsd,jed, .true.) + 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 @@ -3934,8 +3941,8 @@ subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes, turns) 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_diag, G%isd, G%ied, G%jsd, G%jed, & - associated(fluxes_ref%latent_srunoff_diag)) + 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)) @@ -4150,7 +4157,7 @@ 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_diag)) deallocate(fluxes%latent_srunoff_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) @@ -4326,8 +4333,8 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) 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_diag)) & - call rotate_array(fluxes_in%latent_srunoff_diag, turns, fluxes%latent_srunoff_diag) + 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) @@ -4575,10 +4582,10 @@ 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 latent_*_diag fields above, fluxes%latent_srunoff_diag is only associated - ! when submarine-melt coupling is actually configured, so it needs its own guard here. - if (associated(fluxes%latent_srunoff_diag)) & - call homogenize_field_t(fluxes%latent_srunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) + ! 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) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 0547803a0d..1f1eaa664c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -726,8 +726,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings) real :: H_limit_fluxes ! Surface fluxes are scaled down fluxes when the total depth of the ocean ! drops below this value [H ~> m or kg m-2] - real :: srunoff_Teff ! The effective temperature carried by srunoff's contribution to - ! net_heat_rate [C ~> degC] real :: IforcingDepthScale ! The inverse of the layer thickness below which mass losses are ! shifted to the next deeper layer [H ~> m or kg m-2] real :: Idt ! The inverse of the timestep [T-1 ~> s-1] @@ -877,7 +875,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, & !$OMP mixing_depth,A_brine,fraction_left_brine, & !$OMP plume_fraction,dK,total_h,dHeat_srunoff, & - !$OMP dTempxPmE_srunoff,srunoff_Teff) & + !$OMP dTempxPmE_srunoff) & !$OMP firstprivate(SurfPressure,plume_flux) do j=js,je ! Work in vertical slices for efficiency @@ -1002,15 +1000,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! 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. + ! 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) - srunoff_Teff = T2d(i,1) - if (fluxes%srunoff_latent_heat) & - srunoff_Teff = srunoff_Teff - fluxes%latent_heat_fusion*US%J_kg_to_Q/tv%C_p - netheat_rate(i) = netheat_rate(i) + (GV%RZ_to_H * fluxes%srunoff(i,j)) * srunoff_Teff + 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 From be7a5caa903ea06cdc6e503b7f2f518d4e37c4b1 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Wed, 28 Oct 2026 11:22:08 +1100 Subject: [PATCH 10/11] Include srunoff in net fresh water adjustment --- .../drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) 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 395d7510a0..3e62abfc44 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -694,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) From ccd68529edb2326885a21dfcde9ee1e9c0e92ea7 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Wed, 28 Oct 2026 11:27:30 +1100 Subject: [PATCH 11/11] Include srunoff latent heat in net_heat_coupler and LwLatSens --- src/core/MOM_forcing_type.F90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 23052421d8..9827f3dc05 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -3383,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) @@ -3502,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 @@ -3510,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 @@ -3519,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