+Added IST%snow_to_ocn and IST%enth_snow_to_ocn

  Added two new arrays to the ice state type to facilitate the conservative
shedding of snow to the ocean due to mechanical processes, like ridging.  These
fluxes can not be directly sent to the ocean during a mechanical step, so they
are retained until the next slow thermodynamic step.  These arrays are only
allocated and used if DO_RIDGING is true.  Lines prematurely adding these fluxes
to the FIA%fprec_top have been eliminated from SIS_dynamics_trans. In addition,
there are is a new optional argument to alloc_IST_arrays and a new required
argument to ice_ridging, while an optional argument to finish_ice_transport has
been eliminated.  Because ridging is not yet properly implemented and therefore
is commented out, all answers in the MOM6-examples test cases are bitwise
identical.

src/SIS_dyn_trans.F90

@@ -308,9 +308,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
   real :: dt_adv       ! The advective timestep [s].
   real :: Idt_slow     ! The inverse of dt_slow [s-1].
   real, dimension(SZI_(G),SZJ_(G)) :: &
-   ! rdg_s2o, &  ! snow mass [kg m-2] dumped into ocean during ridging
-    rdg_rate, & ! A ridging rate [s-1], this will be calculated from the strain rates in the dynamics.
-    snow2ocn    ! Snow dumped into ocean during ridging [kg m-2]
+    rdg_rate  ! A ridging rate [s-1], this will be calculated from the strain rates in the dynamics.
   integer :: i, j, k, l, m, n, isc, iec, jsc, jec, ncat
   integer :: isd, ied, jsd, jed
   integer :: ndyn_steps, nds ! The number of dynamic steps.
@@ -587,7 +585,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
 
       if (CS%nts==0) &
         call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp, &
-                                  snow2ocn=snow2ocn, rdg_rate=rdg_rate)
+                                  rdg_rate=rdg_rate)
 
       call mpp_clock_end(iceClock8)
     endif
@@ -599,12 +597,6 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
   enddo ! nds=1,ndyn_steps
   call finish_ocean_top_stresses(IOF, G)
 
-  ! Add snow mass dumped into ocean to flux of frozen precipitation:
-  !### WARNING - rdg_s2o is never calculated!!!
-!  if (CS%do_ridging) then ; do k=1,ncat ; do j=jsc,jec ; do i=isc,iec
-!    FIA%fprec_top(i,j,k) = FIA%fprec_top(i,j,k) + rdg_s2o(i,j) * Idt_slow
-!  enddo ; enddo ; enddo ; endif
-
   ! Set appropriate surface quantities in categories with no ice.
   if (allocated(IST%t_surf)) then
     !$OMP parallel do default(shared)

src/SIS_slow_thermo.F90

@@ -806,9 +806,7 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG)
   bsnk(:,:) = 0.0
   salt_change(:,:) = 0.0
   h2o_change(:,:) = 0.0
-!$OMP parallel default(none) shared(isc,iec,jsc,jec,ncat,nb,G,IST,salt_change, &
-!$OMP                               kg_H_Nk,h2o_change,NkIce,IG,CS,IOF,FIA) &
-!$OMP                        private(part_ocn)
+  !$OMP parallel default(shared) private(part_ocn)
   if (CS%ice_rel_salin <= 0.0) then
     !$OMP do
     do j=jsc,jec ; do m=1,NkIce ; do k=1,ncat ; do i=isc,iec
@@ -836,12 +834,25 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG)
     IOF%lprec_ocn_top(i,j) = part_ocn * FIA%lprec_top(i,j,0)
     IOF%fprec_ocn_top(i,j) = part_ocn * FIA%fprec_top(i,j,0)
   enddo ; enddo
+
 ! mw/new precip will eventually be intercepted by pond eliminating need for next 3 lines
   !$OMP do
   do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
     IOF%lprec_ocn_top(i,j) = IOF%lprec_ocn_top(i,j) + &
                              IST%part_size(i,j,k) * FIA%lprec_top(i,j,k)
   enddo ; enddo ; enddo
+
+  ! Add fluxes of snow and other properties to the ocean due to recent ridging or drifting events.
+  if (allocated(IST%snow_to_ocn)) then
+    !$OMP do
+    do j=jsc,jec ; do i=isc,iec ; if (IST%snow_to_ocn(i,j) > 0.0) then
+      IOF%fprec_ocn_top(i,j) = IOF%fprec_ocn_top(i,j) + IST%snow_to_ocn(i,j) * Idt_slow
+      IOF%Enth_Mass_out_ocn(i,j) = IOF%Enth_Mass_out_ocn(i,j) - &
+              IST%snow_to_ocn(i,j) * IST%enth_snow_to_ocn(i,j)
+      ! h2o_change(i,j) = h2o_change(i,j) - IST%snow_to_ocn(i,j)
+      IST%snow_to_ocn(i,j) = 0.0 ;  IST%enth_snow_to_ocn(i,j) = 0.0
+    endif ; enddo ; enddo
+  endif
 !$OMP end parallel
 
   ! Set up temporary tracer array

src/SIS_sum_output.F90

@@ -445,6 +445,11 @@ subroutine write_ice_statistics(IST, day, n, G, IG, CS, message, check_column, t
                   ((0.001*IST%mH_ice(i,j,k)*kg_H_nlay) * IST%sal_ice(i,j,k,L))
       enddo
     endif ; enddo
+    if (allocated(IST%snow_to_ocn)) then ; if (IST%snow_to_ocn(i,j) > 0.0) then
+      area_pt = G%areaT(i,j) * G%mask2dT(i,j)
+      col_mass(i,j,hem) = col_mass(i,j,hem) + area_pt * IST%snow_to_ocn(i,j)
+      col_heat(i,j,hem) = col_heat(i,j,hem) + area_pt * (IST%snow_to_ocn(i,j) * IST%enth_snow_to_ocn(i,j))
+    endif ; endif
     if (ice_area(i,j,hem) > 0.1*G%AreaT(i,j)) ice_extent(i,j,hem) = G%AreaT(i,j)
 
   enddo ; enddo

src/SIS_transport.F90

@@ -211,15 +211,13 @@ end subroutine ice_cat_transport
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 !> finish_ice_transport completes the ice transport and thickness class redistribution
-subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, snow2ocn, rdg_rate)
+subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, rdg_rate)
   type(cell_average_state_type),     intent(inout) :: CAS !< A structure with ocean-cell averaged masses.
   type(ice_state_type),              intent(inout) :: IST !< A type describing the state of the sea ice
   type(SIS_hor_grid_type),           intent(inout) :: G   !< The horizontal grid type
   type(ice_grid_type),               intent(inout) :: IG  !< The sea-ice specific grid type
   type(SIS_tracer_registry_type),    pointer       :: TrReg !< The registry of SIS ice and snow tracers.
   type(SIS_transport_CS),            pointer       :: CS  !< A pointer to the control structure for this module
-  real, dimension(SZI_(G),SZJ_(G)), &
-                           optional, intent(inout) :: snow2ocn !< Snow dumped into ocean during ridging [kg m-2]
   real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: rdg_rate !< The ice ridging rate [s-1].
 
   ! Local variables
@@ -231,6 +229,8 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, snow2ocn, rdg_rate)
     mca0_ice, &  ! The initial mass of ice per unit ocean area in a cell [H ~> kg m-2].
     mca0_snow    ! The initial mass of snow per unit ocean area in a cell [H ~> kg m-2].
 !### These will be needed when the ice ridging is properly implemented.
+!  real :: snow2ocn !< Snow dumped into ocean during ridging [kg m-2]
+!  real :: enth_snow2ocn !< Mass-averaged enthalpy of the now dumped into ocean during ridging [J kg-1]
 !  real, dimension(SZI_(G),SZJ_(G)) :: &
 !    rdg_open, & ! formation rate of open water due to ridging
 !    rdg_vosh    ! rate of ice mass shifted from level to ridged ice
@@ -278,14 +278,21 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, snow2ocn, rdg_rate)
 !   !  ### heat_snow AND OTHER TRACERS ARE OMITTED.
 !   if (CS%do_ridging) then
 !     do j=jsc,jec ; do i=isc,iec
-!       snow2ocn(i,j) = 0.0 !TOM> initializing snow2ocean
 !       if (sum(IST%mH_ice(i,j,:)) > 1.e-10*CS%Rho_ice .and. &
-!           sum(IST%part_size(i,j,1:nCat)) > 0.01) &
+!           sum(IST%part_size(i,j,1:nCat)) > 0.01) then
 !         call ice_ridging(nCat, IST%part_size(i,j,:), IST%mH_ice(i,j,:), &
 !             IST%mH_snow(i,j,:), &
 !             heat_ice(i,j,:,1), heat_ice(i,j,:,2), & !Niki: Is this correct? Bob: No, 2-layers hard-coded.
-!             age_ice(i,j,:), snow2ocn(i,j), rdg_rate(i,j), IST%rgd_mice(i,j,:), &
+!             age_ice(i,j,:), snow2ocn, enth_snow2ocn, rdg_rate(i,j), IST%rgd_mice(i,j,:), &
 !             CAS%dt_sum, IG%mH_cat_bound, rdg_open(i,j), rdg_vosh(i,j))
+!         ! Store the snow mass (and related properties?) that will be passed to the ocean at the
+!         ! next opportunity.
+!         if (snow2ocn > 0.0) then
+!           IST%enth_snow_to_ocn(i,j) = (IST%enth_snow_to_ocn(i,j) * IST%snow_to_ocn(i,j) + enth_snow2ocn * snow2ocn) / &
+!                                       (IST%snow_to_ocn(i,j) + snow2ocn)
+!           IST%snow_to_ocn(i,j) = IST%snow_to_ocn(i,j) + snow2ocn
+!         endif
+!       endif
 !     enddo ; enddo
 !   endif   ! do_ridging
 

src/SIS_types.F90

@@ -77,6 +77,12 @@ module SIS_types
     mH_ice, &   !< The mass per unit area of the ice in each category [H ~> kg m-2].
     t_surf      !< The surface temperature [Kelvin].
 
+  real, allocatable, dimension(:,:) :: &
+    snow_to_ocn, & !< The mass per unit ocean area of snow that will be dumped into the
+                   !! ocean due to recent mechanical activities like ridging or drifting [kg m-2].
+    enth_snow_to_ocn !< The average enthalpy of the snow that will be dumped into the
+                   !! ocean due to recent mechanical activities like ridging or drifting [Enth ~> J kg-1].
+
   real, allocatable, dimension(:,:,:,:) :: sal_ice  !< The salinity of the sea ice
                 !! in each category and fractional thickness layer [gSalt kg-1].
   real, allocatable, dimension(:,:,:,:) :: enth_ice !< The enthalpy of the sea ice
@@ -403,12 +409,13 @@ module SIS_types
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 !> alloc_IST_arrays allocates the arrays in an ice_state_type.
-subroutine alloc_IST_arrays(HI, IG, IST, omit_velocities, omit_Tsurf)
+subroutine alloc_IST_arrays(HI, IG, IST, omit_velocities, omit_Tsurf, do_ridging)
   type(hor_index_type), intent(in)    :: HI  !< The horizontal index type describing the domain
   type(ice_grid_type),  intent(in)    :: IG  !< The sea-ice specific grid type
   type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice
   logical,    optional, intent(in)    :: omit_velocities !< If true, do not allocate velocity arrays
   logical,    optional, intent(in)    :: omit_Tsurf !< If true, do not allocate the surface temperature array
+  logical,    optional, intent(in)    :: do_ridging !< If true, allocate arrays related to ridging
 
   integer :: isd, ied, jsd, jed, CatIce, NkIce, idr
   logical :: do_vel, do_Tsurf
@@ -427,6 +434,11 @@ subroutine alloc_IST_arrays(HI, IG, IST, omit_velocities, omit_Tsurf)
   allocate(IST%enth_ice( isd:ied, jsd:jed, CatIce, NkIce)) ; IST%enth_ice(:,:,:,:) = 0.0
   allocate(IST%sal_ice(  isd:ied, jsd:jed, CatIce, NkIce)) ; IST%sal_ice(:,:,:,:) = 0.0
 
+  if (present(do_ridging)) then ; if (do_ridging) then
+    allocate(IST%snow_to_ocn(isd:ied, jsd:jed)) ; IST%snow_to_ocn(:,:) = 0.0
+    allocate(IST%enth_snow_to_ocn(isd:ied, jsd:jed)) ; IST%enth_snow_to_ocn(:,:) = 0.0
+  endif ; endif
+
   if (do_vel) then
     ! These velocities are only required for the slow ice processes, and hence
     ! can use the memory macros.
@@ -487,6 +499,13 @@ subroutine ice_state_register_restarts(IST, G, IG, Ice_restart, restart_file)
     idr = register_restart_field(Ice_restart, restart_file, 'sal_ice', IST%sal_ice, &
                                  domain=mpp_domain, mandatory=.false., units="kg/kg")
 
+    if (allocated(IST%snow_to_ocn)) then
+      idr = register_restart_field(Ice_restart, restart_file, 'snow_to_ocn', IST%snow_to_ocn, &
+                                   domain=mpp_domain, mandatory=.false., units="kg m-2")
+      idr = register_restart_field(Ice_restart, restart_file, 'enth_snow_to_ocn', IST%enth_snow_to_ocn, &
+                                   domain=mpp_domain, mandatory=.false., units="J kg-1")
+    endif
+
     if (IST%Cgrid_dyn) then
       if (G%symmetric) then
         idr = register_restart_field(Ice_restart, restart_file, 'sym_u_ice_C', IST%u_ice_C, &
@@ -1979,6 +1998,8 @@ subroutine dealloc_IST_arrays(IST)
   deallocate(IST%part_size, IST%mH_snow, IST%mH_ice)
   deallocate(IST%mH_pond) ! mw/new
   deallocate(IST%enth_snow, IST%enth_ice, IST%sal_ice)
+  if (allocated(IST%snow_to_ocn)) deallocate(IST%snow_to_ocn)
+  if (allocated(IST%enth_snow_to_ocn)) deallocate(IST%enth_snow_to_ocn)
   if (allocated(IST%t_surf)) deallocate(IST%t_surf)
 
   if (allocated(IST%u_ice_C)) deallocate(IST%u_ice_C)

src/ice_model.F90

@@ -2129,7 +2129,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow,
     call ice_type_slow_reg_restarts(sGD%mpp_domain, CatIce, &
                       param_file, Ice, Ice%Ice_restart, restart_file)
 
-    call alloc_IST_arrays(sHI, sIG, sIST, omit_tsurf=Eulerian_tsurf)
+    call alloc_IST_arrays(sHI, sIG, sIST, omit_tsurf=Eulerian_tsurf, do_ridging=do_ridging)
     call ice_state_register_restarts(sIST, sG, sIG, Ice%Ice_restart, restart_file)
 
     call alloc_ocean_sfc_state(Ice%sCS%OSS, sHI, sIST%Cgrid_dyn, gas_fields_ocn)
@@ -2486,7 +2486,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow,
       ! However, in some model this causes answers to change after a restart
       ! because these state variables are changed only after a restart and
       ! not at each timestep. Provide a switch to turn this option off.
-      if(do_mask_restart) then
+      if (do_mask_restart) then
         do k=1,CatIce
           sIST%mH_snow(:,:,k) = sIST%mH_snow(:,:,k) * H_rescale_snow * sG%mask2dT(:,:)
           sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:)

src/ice_ridge.F90

@@ -196,7 +196,7 @@ end function ridge_rate
 !TOM>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 !> ice_ridging parameterizes mechanical redistribution of thin (undeformed) ice
 !! into thicker (deformed/ridged) ice categories
-subroutine ice_ridging(km, cn, hi, hs, t1, t2, age, snow_to_ocn, rdg_rate, hi_rdg, &
+subroutine ice_ridging(km, cn, hi, hs, t1, t2, age, snow_to_ocn, enth_snow_to_ocn, rdg_rate, hi_rdg, &
                        dt, hlim_in, rdg_open, vlev)
   !  Subroutine written by T. Martin, 2008
   integer,             intent(in)    :: km  !< The number of ice thickness categories
@@ -208,7 +208,10 @@ subroutine ice_ridging(km, cn, hi, hs, t1, t2, age, snow_to_ocn, rdg_rate, hi_rd
   real, dimension(1:), intent(inout) :: t1  !< Volume integrated upper layer temperature [degC m3]?
   real, dimension(1:), intent(inout) :: t2  !< Volume integrated upper layer temperature [degC m3]?
   real, dimension(1:), intent(inout) :: age !< Volume integrated ice age [m3 years]?
-  real,                intent(out)   :: snow_to_ocn !< total snow volume dumped into ocean during ridging
+  real,                intent(out)   :: enth_snow_to_ocn !< average of enthalpy of the snow dumped into
+                                            !! ocean due to this ridging event [J kg-1]
+  real,                intent(out)   :: snow_to_ocn !< total snow mass dumped into ocean due to this
+                                            !! ridging event [kg m-2]
   real,                intent(in)    :: rdg_rate    !< Ridging rate from subroutine ridge_rate
   real, dimension(1:), intent(inout) :: hi_rdg      !< A diagnostic of the ridged ice volume in each category.
   real,                intent(in)    :: dt          !< time step dt has units seconds
@@ -232,6 +235,7 @@ subroutine ice_ridging(km, cn, hi, hs, t1, t2, age, snow_to_ocn, rdg_rate, hi_rd
   real                  :: ardg, vrdg  ! area and volume of newly formed rdiged (vlev=vrdg!!!)
   real, dimension(1:km) :: hmin, hmax, efold, rdg_ratio, hlim
   real                  :: hl, hr
+  real                  :: snow_dump, enth_dump
   real                  :: cn_tot, part_undef_sum
   real                  :: div_adv, Rnet, Rdiv, Rtot, rdg_area, rdgtmp, hlimtmp
   real                  :: area_frac
@@ -244,7 +248,7 @@ subroutine ice_ridging(km, cn, hi, hs, t1, t2, age, snow_to_ocn, rdg_rate, hi_rd
   hlimtmp = hlim_in(km)
   hlim(km) = hlim_unlim   ! ensures all ridged ice is smaller than thickest ice allowed
   frac_hs_rdg = 1.0-s2o_frac
-  !snow_to_ocn = 0.0 -> done in subroutine transport
+  snow_to_ocn = 0.0 ; enth_snow_to_ocn = 0.0
   alev=0.0; ardg=0.0; vlev=0.0; vrdg=0.0
   !
   call ice_ridging_init(km, cn, hi, part_undef, part_undef_sum, &
@@ -346,7 +350,12 @@ subroutine ice_ridging(km, cn, hi, hs, t1, t2, age, snow_to_ocn, rdg_rate, hi_rd
         hi_rdg(kd) = max(hi_rdg(kd),0.0)      ! ensure hi_rdg >= 0
 
         ! dump part of the snow in ocean (here, sum volume, transformed to flux in update_ice_model_slow)
-        snow_to_ocn = snow_to_ocn + frac_hs(kd)*(1.0-frac_hs_rdg)
+        snow_dump = frac_hs(kd)*(1.0-frac_hs_rdg)
+        if (snow_to_ocn > 0.0) then
+          enth_dump = t1(kd)  !### THIS IS WRONG, BUT IS A PLACEHOLDER FOR NOW.
+          enth_snow_to_ocn = (enth_snow_to_ocn*snow_to_ocn + enth_dump*snow_dump) / (snow_to_ocn + snow_dump)
+          snow_to_ocn = snow_to_ocn + snow_dump
+        endif
 
       enddo