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: #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.
NOAA-GFDL · Aug 14, 2024 · 1c8dc7e · 1c8dc7e
1 parent 2c49005
commit 1c8dc7e
Show file tree

Hide file tree

Showing 7 changed files with 28 additions and 19 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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