From 1c8dc7e8f73aba6ba5109524fbd455c44b1ea2e7 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Wed, 14 Aug 2024 16:57:35 -0700 Subject: [PATCH] Attempt at separating sea ice melt/formation from liquid precipitation For unjustified legacy reasons, these have been bundled together in the past. They are not only bundled together in SIS2 variables, but also in the ice-ocean fluxes passed to the coupler and hence to MOM6. This commit attempts to address this Issue: https://github.com/NOAA-GFDL/SIS2/issues/213 However, I currently am getting a SEGFAULT error stating that `iobt%seaice_melt` has not been allocated memory when I try to do a checksum on it. --- src/SIS_slow_thermo.F90 | 9 ++------- src/SIS_sum_output.F90 | 1 + src/SIS_types.F90 | 5 ++++- src/combined_ice_ocean_driver.F90 | 10 ++++++---- src/ice_boundary_types.F90 | 6 +++--- src/ice_model.F90 | 4 +++- src/ice_type.F90 | 12 +++++++++--- 7 files changed, 28 insertions(+), 19 deletions(-) diff --git a/src/SIS_slow_thermo.F90 b/src/SIS_slow_thermo.F90 index d4bc6f0f..7f63788b 100644 --- a/src/SIS_slow_thermo.F90 +++ b/src/SIS_slow_thermo.F90 @@ -1322,7 +1322,8 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) ! With transmuted ice, the ice is non-conservatively changed to match the ocean properties. IOF%flux_salt(i,j) = IOF%flux_salt(i,j) + salt_to_ocn * (0.001*Idt_slow) - net_melt(i,j) = net_melt(i,j) + water_to_ocn * Idt_slow ! This goes to IOF%lprec_ocn_top. + net_melt(i,j) = net_melt(i,j) + water_to_ocn * Idt_slow + IOF%seaice_melt(i,j) = net_melt(i,j) IOF%Enth_mass_out_ocn(i,j) = IOF%Enth_mass_out_ocn(i,j) + heat_to_ocn ! With transmuted ice, the imbalances are stored to close the heat and salt budgets. @@ -1369,12 +1370,6 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) call disable_SIS_averaging(CS%diag) - ! Combine the liquid precipitation with the net melt of ice and snow for - ! passing to the ocean. These may later be kept separate. - do j=jsc,jec ; do i=isc,iec - IOF%lprec_ocn_top(i,j) = IOF%lprec_ocn_top(i,j) + net_melt(i,j) - enddo ; enddo - ! Make sure TrLay is no longer allocated if(allocated(TrLay)) deallocate(TrLay) end subroutine SIS2_thermodynamics diff --git a/src/SIS_sum_output.F90 b/src/SIS_sum_output.F90 index ecfdcfc3..575db44b 100644 --- a/src/SIS_sum_output.F90 +++ b/src/SIS_sum_output.F90 @@ -759,6 +759,7 @@ subroutine accumulate_bottom_input(IST, OSS, FIA, IOF, dt, G, US, IG, CS) do j=jsc,jec ; do i=isc,iec CS%water_in_col(i,j) = CS%water_in_col(i,j) - dt * & ( ((FIA%runoff(i,j) + FIA%calving(i,j)) + & + IOF%seaice_melt(i,j) + & (IOF%lprec_ocn_top(i,j) + IOF%fprec_ocn_top(i,j))) - IOF%evap_ocn_top(i,j) ) Flux_SW = 0.0 diff --git a/src/SIS_types.F90 b/src/SIS_types.F90 index 18d5809f..f941b6c9 100644 --- a/src/SIS_types.F90 +++ b/src/SIS_types.F90 @@ -358,6 +358,7 @@ module SIS_types flux_lh_ocn_top, & !< The upward flux of latent heat at the ocean surface [Q R Z T-1 ~> W m-2]. lprec_ocn_top, & !< The downward flux of liquid precipitation at the ocean surface [R Z T-1 ~> kg m-2 s-1]. fprec_ocn_top, & !< The downward flux of frozen precipitation at the ocean surface [R Z T-1 ~> kg m-2 s-1]. + seaice_melt, & !< The downward freshwater flux into the ocean due to sea ice melt [R Z T-1 ~> kg m-2 s-1]. flux_u_ocn, & !< The flux of x-momentum into the ocean at locations given by !! flux_uv_stagger [R Z L T-2 ~> Pa]. !! Note that regardless of the staggering, flux_u_ocn is allocated as though on an A-grid. @@ -925,6 +926,7 @@ subroutine alloc_ice_ocean_flux(IOF, HI, do_stress_mag, do_iceberg_fields, do_tr allocate(IOF%flux_sw_ocn(SZI_(HI), SZJ_(HI), NBANDS), source=0.0) allocate(IOF%lprec_ocn_top(SZI_(HI), SZJ_(HI)), source=0.0) allocate(IOF%fprec_ocn_top(SZI_(HI), SZJ_(HI)), source=0.0) + allocate(IOF%seaice_melt(SZI_(HI), SZJ_(HI)), source=0.0) allocate(IOF%flux_u_ocn(SZI_(HI), SZJ_(HI)), source=0.0) allocate(IOF%flux_v_ocn(SZI_(HI), SZJ_(HI)), source=0.0) if (alloc_stress_mag) then @@ -2171,7 +2173,7 @@ subroutine dealloc_ice_ocean_flux(IOF) deallocate(IOF%flux_sh_ocn_top, IOF%evap_ocn_top) deallocate(IOF%flux_lw_ocn_top, IOF%flux_lh_ocn_top) deallocate(IOF%flux_sw_ocn) - deallocate(IOF%lprec_ocn_top, IOF%fprec_ocn_top, IOF%flux_salt) + deallocate(IOF%lprec_ocn_top, IOF%fprec_ocn_top, IOF%seaice_melt, IOF%flux_salt) deallocate(IOF%flux_u_ocn, IOF%flux_v_ocn, IOF%pres_ocn_top, IOF%mass_ice_sn_p) if (allocated(IOF%stress_mag)) deallocate(IOF%stress_mag) if (allocated(IOF%transmutation_salt_flux)) deallocate(IOF%transmutation_salt_flux) @@ -2218,6 +2220,7 @@ subroutine IOF_chksum(mesg, IOF, G, US, mech_fluxes, thermo_fluxes) call hchksum(IOF%evap_ocn_top, trim(mesg)//" IOF%evap_ocn_top", G%HI, scale=US%RZ_T_to_kg_m2s) call hchksum(IOF%lprec_ocn_top, trim(mesg)//" IOF%lprec_ocn_top", G%HI, scale=US%RZ_T_to_kg_m2s) call hchksum(IOF%fprec_ocn_top, trim(mesg)//" IOF%fprec_ocn_top", G%HI, scale=US%RZ_T_to_kg_m2s) + call hchksum(IOF%seaice_melt, trim(mesg)//" IOF%seaice_melt", G%HI, scale=US%RZ_T_to_kg_m2s) endif if (do_mech) then call hchksum(IOF%flux_u_ocn, trim(mesg)//" IOF%flux_u_ocn", G%HI, scale=US%RZ_T_to_kg_m2s*US%L_T_to_m_s) diff --git a/src/combined_ice_ocean_driver.F90 b/src/combined_ice_ocean_driver.F90 index 3c159c84..410979b7 100644 --- a/src/combined_ice_ocean_driver.F90 +++ b/src/combined_ice_ocean_driver.F90 @@ -304,6 +304,7 @@ subroutine direct_flux_ice_to_IOB(Time, Ice, IOB, do_thermo) if (ASSOCIATED(IOB%lw_flux)) IOB%lw_flux(:,:) = Ice%flux_lw(:,:) if (ASSOCIATED(IOB%lprec)) IOB%lprec(:,:) = Ice%lprec(:,:) if (ASSOCIATED(IOB%fprec)) IOB%fprec(:,:) = Ice%fprec(:,:) + if (ASSOCIATED(IOB%seaice_melt)) IOB%seaice_melt(:,:) = Ice%seaice_melt(:,:) if (ASSOCIATED(IOB%runoff)) IOB%runoff(:,:) = Ice%runoff(:,:) if (ASSOCIATED(IOB%calving)) IOB%calving(:,:) = Ice%calving if (ASSOCIATED(IOB%runoff_hflx)) IOB%runoff_hflx(:,:) = Ice%runoff_hflx(:,:) @@ -326,10 +327,11 @@ subroutine direct_flux_ice_to_IOB(Time, Ice, IOB, do_thermo) call data_override('OCN', 'sw_flux_nir_dif', IOB%sw_flux_nir_dif, Time) call data_override('OCN', 'sw_flux_vis_dir', IOB%sw_flux_vis_dir, Time) call data_override('OCN', 'sw_flux_vis_dif', IOB%sw_flux_vis_dif, Time) - call data_override('OCN', 'lprec', IOB%lprec , Time) - call data_override('OCN', 'fprec', IOB%fprec , Time) - call data_override('OCN', 'runoff', IOB%runoff , Time) - call data_override('OCN', 'calving', IOB%calving , Time) + call data_override('OCN', 'lprec', IOB%lprec , Time) + call data_override('OCN', 'fprec', IOB%fprec , Time) + call data_override('OCN', 'seaice_melt', IOB%seaice_melt , Time) + call data_override('OCN', 'runoff', IOB%runoff , Time) + call data_override('OCN', 'calving', IOB%calving , Time) call data_override('OCN', 'runoff_hflx', IOB%runoff_hflx , Time) call data_override('OCN', 'calving_hflx', IOB%calving_hflx , Time) endif diff --git a/src/ice_boundary_types.F90 b/src/ice_boundary_types.F90 index b9ea3a82..ff8ef034 100644 --- a/src/ice_boundary_types.F90 +++ b/src/ice_boundary_types.F90 @@ -37,9 +37,9 @@ module ice_boundary_types frazil => NULL(), & !< The frazil heat rejected by the ocean [J m-2]. sea_level => NULL(), & !< The sea level after adjustment for any surface !! pressure that the ocean allows to be expressed [m]. - calving => NULL(), & !< The mass flux per unit area of the ice shelf to convert to - !!bergs [RZ_T ~> kg m-2 s-1]. - calving_hflx => NULL() !< Calving heat flux [Q R Z T-1 ~> W m-2]. + calving => NULL(), & !< The mass flux per unit area of the ice shelf to convert to + !!bergs [RZ_T ~> kg m-2 s-1]. + calving_hflx => NULL() !< Calving heat flux [Q R Z T-1 ~> W m-2]. real, dimension(:,:,:), pointer :: data =>NULL() !< S collective field for "named" fields above integer :: stagger = BGRID_NE !< A flag indicating how the velocities are staggered. integer :: xtype !< A flag indicating the exchange type, which may be set to diff --git a/src/ice_model.F90 b/src/ice_model.F90 index 1042270d..6eee20bf 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -623,6 +623,7 @@ subroutine set_ocean_top_fluxes(Ice, IST, IOF, FIA, OSS, G, US, IG, sCS) Ice%flux_sw_vis_dir(:,:) = 0.0 ; Ice%flux_sw_vis_dif(:,:) = 0.0 Ice%flux_lw(:,:) = 0.0 ; Ice%flux_lh(:,:) = 0.0 Ice%fprec(:,:) = 0.0 ; Ice%lprec(:,:) = 0.0 + Ice%seaice_melt(:,:) = 0.0 call coupler_type_rescale_data(Ice%ocean_fluxes, 0.0) i_off = LBOUND(Ice%flux_t,1) - G%isc ; j_off = LBOUND(Ice%flux_t,2) - G%jsc @@ -640,6 +641,7 @@ subroutine set_ocean_top_fluxes(Ice, IST, IOF, FIA, OSS, G, US, IG, sCS) Ice%flux_lh(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_lh_ocn_top(i,j) Ice%fprec(i2,j2) = US%RZ_T_to_kg_m2s*IOF%fprec_ocn_top(i,j) Ice%lprec(i2,j2) = US%RZ_T_to_kg_m2s*IOF%lprec_ocn_top(i,j) + Ice%seaice_melt(i2,j2) = US%RZ_T_to_kg_m2s*IOF%seaice_melt(i,j) Ice%runoff(i2,j2) = US%RZ_T_to_kg_m2s*FIA%runoff(i,j) Ice%calving(i2,j2) = US%RZ_T_to_kg_m2s*FIA%calving(i,j) Ice%runoff_hflx(i2,j2) = US%QRZ_T_to_W_m2*FIA%runoff_hflx(i,j) @@ -659,7 +661,7 @@ subroutine set_ocean_top_fluxes(Ice, IST, IOF, FIA, OSS, G, US, IG, sCS) if (allocated(IOF%melt_nudge)) then do j=jsc,jec ; do i=isc,iec i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences. - Ice%lprec(i2,j2) = Ice%lprec(i2,j2) + US%RZ_T_to_kg_m2s*IOF%melt_nudge(i,j) + Ice%seaice_melt(i2,j2) = Ice%seaice_melt(i2,j2) + US%RZ_T_to_kg_m2s*IOF%melt_nudge(i,j) enddo ; enddo endif diff --git a/src/ice_type.F90 b/src/ice_type.F90 index 1687de67..950c239e 100644 --- a/src/ice_type.F90 +++ b/src/ice_type.F90 @@ -97,9 +97,10 @@ module ice_type_mod flux_sw_vis_dif => NULL(), & !< The diffuse visible shortwave heat flux into the ocean [W m-2]. flux_sw_nir_dir => NULL(), & !< The direct near-infrared heat flux into the ocean [W m-2]. flux_sw_nir_dif => NULL(), & !< The diffuse near-infrared heat flux into the ocean [W m-2]. - flux_lh => NULL(), & !< The latent heat flux out of the ocean [W m-2]. - lprec => NULL(), & !< The liquid precipitation flux into the ocean [kg m-2]. - fprec => NULL(), & !< The frozen precipitation flux into the ocean [kg m-2]. + flux_lh => NULL(), & !< The latent heat flux out of the ocean [W m-2]. + lprec => NULL(), & !< The liquid precipitation flux into the ocean [kg m-2]. + fprec => NULL(), & !< The frozen precipitation flux into the ocean [kg m-2]. + seaice_melt => NULL(), & !< The sea ice meltwater flux into the ocean [kg m-2]. p_surf => NULL(), & !< The pressure at the ocean surface [Pa]. This may !! or may not include atmospheric pressure. runoff => NULL(), & !< Liquid runoff into the ocean [kg m-2]. @@ -198,6 +199,7 @@ subroutine ice_type_slow_reg_restarts(domain, CatIce, param_file, Ice, & call safe_alloc_ptr(Ice%flux_lh, isc, iec, jsc, jec) !NI call safe_alloc_ptr(Ice%lprec, isc, iec, jsc, jec) call safe_alloc_ptr(Ice%fprec, isc, iec, jsc, jec) + call safe_alloc_ptr(Ice%seaice_melt, isc, iec, jsc, jec) call safe_alloc_ptr(Ice%p_surf, isc, iec, jsc, jec) call safe_alloc_ptr(Ice%runoff, isc, iec, jsc, jec) call safe_alloc_ptr(Ice%calving, isc, iec, jsc, jec) @@ -233,6 +235,7 @@ subroutine ice_type_slow_reg_restarts(domain, CatIce, param_file, Ice, & call register_restart_field(Ice_restart, 'flux_lw', Ice%flux_lw) call register_restart_field(Ice_restart, 'lprec', Ice%lprec) call register_restart_field(Ice_restart, 'fprec', Ice%fprec) + call register_restart_field(Ice_restart, 'seaice_melt', Ice%seaice_melt) call register_restart_field(Ice_restart, 'runoff', Ice%runoff) call register_restart_field(Ice_restart, 'calving', Ice%calving) call register_restart_field(Ice_restart, 'runoff_hflx', Ice%runoff_hflx, mandatory=.false.) @@ -338,6 +341,7 @@ subroutine dealloc_Ice_arrays(Ice) if (associated(Ice%flux_lh)) deallocate(Ice%flux_lh) if (associated(Ice%lprec)) deallocate(Ice%lprec) if (associated(Ice%fprec)) deallocate(Ice%fprec) + if (associated(Ice%seaice_melt)) deallocate(Ice%seaice_melt) if (associated(Ice%p_surf)) deallocate(Ice%p_surf) if (associated(Ice%runoff)) deallocate(Ice%runoff) if (associated(Ice%calving)) deallocate(Ice%calving) @@ -419,6 +423,7 @@ subroutine Ice_public_type_chksum(mesg, Ice, check_fast, check_slow, check_rough call chksum(Ice%flux_lh, trim(mesg)//" Ice%flux_lh") call chksum(Ice%lprec, trim(mesg)//" Ice%lprec") call chksum(Ice%fprec, trim(mesg)//" Ice%fprec") + call chksum(Ice%seaice_melt, trim(mesg)//" Ice%seaice_melt") call chksum(Ice%p_surf, trim(mesg)//" Ice%p_surf") call chksum(Ice%calving, trim(mesg)//" Ice%calving") call chksum(Ice%runoff, trim(mesg)//" Ice%runoff") @@ -680,6 +685,7 @@ subroutine ice_data_type_chksum(mesg, timestep, Ice, init_call) chks = SIS_chksum(Ice%flux_lh ) ; if (root) write(outunit,100) 'ice_data_type%flux_lh ', chks chks = SIS_chksum(Ice%lprec ) ; if (root) write(outunit,100) 'ice_data_type%lprec ', chks chks = SIS_chksum(Ice%fprec ) ; if (root) write(outunit,100) 'ice_data_type%fprec ', chks + chks = SIS_chksum(Ice%seaice_melt ) ; if (root) write(outunit,100) 'ice_data_type%seaice_melt ', chks chks = SIS_chksum(Ice%p_surf ) ; if (root) write(outunit,100) 'ice_data_type%p_surf ', chks chks = SIS_chksum(Ice%runoff ) ; if (root) write(outunit,100) 'ice_data_type%runoff ', chks chks = SIS_chksum(Ice%calving ) ; if (root) write(outunit,100) 'ice_data_type%calving ', chks