Sea ice mass budget diagnostics for CMIP7

src/SIS2_ice_thm.F90

@@ -1008,8 +1008,9 @@ end subroutine ice_check
 subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
                            snow, rain, evap, tmlt, bmlt, NkIce, npassive, TrLay, &
                            heat_to_ocn, h2o_ice_to_ocn, h2o_ocn_to_ice, evap_from_ocn, &
-                           snow_to_ice, salt_to_ice, ITV, US, CS, ablation, &
-                           enthalpy_evap, enthalpy_melt, enthalpy_freeze)
+                           snow_to_ice, s2i_i, s2i_s, salt_to_ice, ITV, US, CS, ablation, &
+                           ablation_i, ablation_s, enthalpy_evap, enthalpy_melt, enthalpy_freeze, &
+                           evap_i, evap_s)
   ! mw/new - melt pond - added first two arguments & rain
   real, intent(in   ) :: a_ice       !< area of ice (1-open_water_frac) for pond retention [nondim]
   real, intent(inout) :: m_pond      !< melt pond mass [R Z ~> kg m-2]
@@ -1023,7 +1024,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
         intent(inout) :: Salin       !< Conserved ice bulk salinity by layer [S ~> gSalt kg-1]
   real, intent(in   ) :: snow        !< new snow [R Z ~> kg m-2]
   real, intent(in   ) :: rain        !< rain for pond source [R Z ~> kg m-2] - not yet active
-  real, intent(in   ) :: evap        !< ice evaporation/sublimation [R Z ~> kg m-2]
+  real, intent(in   ) :: evap        !< ice and snow evaporation/sublimation [R Z ~> kg m-2]
   real, intent(in   ) :: tmlt        !< top melting energy [Q R Z ~> J m-2]
   real, intent(in   ) :: bmlt        !< bottom melting energy [Q R Z ~> J m-2]
   integer, intent(in) :: NkIce       !< The number of ice layers.
@@ -1035,18 +1036,24 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
   real, intent(  out) :: h2o_ocn_to_ice !< liquid water flux from ocean [R Z ~> kg m-2]
   real, intent(  out) :: evap_from_ocn  !< evaporation flux from ocean [R Z ~> kg m-2]
   real, intent(  out) :: snow_to_ice !< snow below waterline becomes ice [R Z ~> kg m-2]
+  real, intent(  out) :: s2i_i       !< snow below waterline becomes ice (ice flux) [R Z ~> kg m-2]
+  real, intent(  out) :: s2i_s       !< snow below waterline becomes ice (snow flux) [R Z ~> kg m-2]
   real, intent(  out) :: salt_to_ice !< Net flux of salt to the ice [R Z S ~> gSalt m-2].
   type(ice_thermo_type), intent(in) :: ITV !< The ice thermodynamic parameter structure.
   type(unit_scale_type), intent(in) :: US  !< A structure with unit conversion factors
   type(SIS2_ice_thm_CS), intent(in) :: CS  !< The SIS2_ice_thm control structure.
 
-  real, intent(  out) :: ablation      !< The mass loss from bottom melt [R Z ~> kg m-2].
+  real, intent(  out) :: ablation      !< The mass loss from bottom melt of ice and snow [R Z ~> kg m-2].
+  real, intent(  out) :: ablation_i    !< The mass loss from bottom melt of ice [R Z ~> kg m-2].
+  real, intent(  out) :: ablation_s    !< The mass loss from bottom melt of snow [R Z ~> kg m-2].
   real, intent(  out) :: enthalpy_evap !< The enthalpy loss due to the mass loss
                                        !! by evaporation / sublimation. [Q R Z ~> J m-2]
   real, intent(  out) :: enthalpy_melt !< The enthalpy loss due to the mass loss
                                        !! by melting [Q R Z ~> J m-2].
   real, intent(  out) :: enthalpy_freeze !< The enthalpy gain due to the mass gain
                                        !! by freezing [Q R Z ~> J m-2].
+  real, intent(  out) :: evap_i        !< ice evaporation/sublimation [R Z ~> kg m-2]
+  real, intent(  out) :: evap_s        !< snow evaporation/sublimation [R Z ~> kg m-2]
 
   real :: top_melt, bot_melt, melt_left ! Heating amounts, all in [Q R Z ~> J m-2]
   real :: mtot_ice    ! The summed ice mass [R Z ~> kg m-2].
@@ -1097,9 +1104,11 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
 
   evap_from_ocn = 0.0 ! for excess evap-melt
   h2o_ocn_to_ice = 0.0 ; h2o_ice_to_ocn = 0.0 ; snow_to_ice  = 0.0
+  s2i_i = 0.0 ; s2i_s = 0.0
   h2o_to_pond = 0.0
   h2o_from_pond = 0.0
   salt_to_ice = 0.0
+  evap_s = 0.0 ; evap_i = 0.0
   enthM_freezing = 0.0 ; enthM_melt = 0.0 ; enthM_evap = 0.0 ; enthM_snowfall = 0.0
 
   ! raining on cold ice led to unphysical temperature oscillations
@@ -1117,6 +1126,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
   if (evap < 0.0) then
     m_lay(0) = m_lay(0) - evap ! Treat frost formation like snow.
     enthM_snowfall = enthM_snowfall - evap*enthalpy(0)
+    evap_s = evap_s + evap
   endif
 
   if (top_melt < 0.0 .and. CS%do_pond) then ! mw/new: add fresh/0C ice to top layer
@@ -1185,6 +1195,11 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
       evap_here = min(evap_left, m_lay(k))
       evap_left = evap_left - evap_here
       m_lay(k) = m_lay(k) - evap_here
+      if (k == 0) then
+        evap_s = evap_s + evap_here
+      else
+        evap_i = evap_i + evap_here
+      endif
       ! Assume that evaporation does not make ice salty?
       if (k>0) salt_to_ice = salt_to_ice - Salin(k) * evap_here
       enthM_evap = enthM_evap + evap_here * enthalpy(k)
@@ -1229,6 +1244,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
   ! apply bottom melt heat flux
 
   melt_left = bot_melt ; ablation = 0.0
+  ablation_i = 0.0 ; ablation_s = 0.0
   if (melt_left > 0.0) then ! melt ice and snow from below
     do k=NkIce,0,-1
       if (melt_left < m_lay(k) * (enth_fr(k) - Enthalpy(k))) then
@@ -1243,7 +1259,11 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
       h2o_ice_to_ocn = h2o_ice_to_ocn + M_melt
       enthM_melt = enthM_melt + M_melt*enth_fr(k)
       ablation = ablation + M_melt
-
+      if (k > 0) then
+        ablation_i = ablation_i + M_melt
+      else
+        ablation_s = ablation_s + M_melt
+      endif
       if (melt_left <= 0.0) exit ! All melt energy has been used.
     enddo
 
@@ -1287,6 +1307,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
     snow_to_ice = min(m_submerged - mtot_ice, m_lay(0)) ! need ice from snow
 
     m_lay(0) = m_lay(0) - snow_to_ice
+    s2i_s = sn2i_s - snow_to_ice
 
     ! Add ice to the topmost layer and dilute its salinity.
     Enthalpy(1) = (m_lay(1)*Enthalpy(1) + snow_to_ice*Enthalpy(0)) / &
@@ -1296,6 +1317,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, &
       TrLay(1,tr) = TrLay(1,tr) * m_lay(1) / (m_lay(1) + snow_to_ice)
     enddo
     m_lay(1) = m_lay(1) + snow_to_ice
+    s2i_i = s2i_i + snow_to_ice
   else
     snow_to_ice = 0.0
   endif

src/SIS_ctrl_types.F90

@@ -238,6 +238,17 @@ subroutine ice_diagnostics_init(IOF, OSS, FIA, G, US, IG, diag, Time, Cgrid)
                'frozen runoff sensible heat flux', 'W/m^2', conversion=US%QRZ_T_to_W_m2, missing_value=missing)
   FIA%id_evap     = register_SIS_diag_field('ice_model', 'EVAP',diag%axesT1, Time, &
                'evaporation', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, missing_value=missing)
+
+  !CMOR evaporation diagnostics
+  FIA%id_evap_i   = register_SIS_diag_field('ice_model', 'sidmassevapsubl',diag%axesT1, Time, &
+               'Sea-Ice Mass Change Through Evaporation and Sublimation', 'kg/(m^2*s)', &
+               conversion=US%RZ_T_to_kg_m2s, missing_value=missing, &
+               cmor_standard_name='water_evapotranspiration_flux')
+  FIA%id_evap_s   = register_SIS_diag_field('ice_model', 'sisndmasssubl',diag%axesT1, Time, &
+               'Snow Mass Rate of Change Through Evaporation or Sublimation', 'kg/(m^2*s)', &
+               conversion=US%RZ_T_to_kg_m2s, missing_value=missing, &
+               cmor_standard_name='tendency_of_atmosphere_mass_content_of_water_vapor_due_to_sublimation_of_surface_snow_and_ice')
+
   IOF%id_saltf    = register_SIS_diag_field('ice_model', 'SALTF', diag%axesT1, Time, &
                'ice to ocean salt flux', 'kg/(m^2*s)', conversion=US%S_to_ppt*US%RZ_T_to_kg_m2s, missing_value=missing)
   FIA%id_tmelt    = register_SIS_diag_field('ice_model', 'TMELT', diag%axesT1, Time, &
@@ -337,6 +348,12 @@ subroutine ice_diagnostics_init(IOF, OSS, FIA, G, US, IG, diag, Time, Cgrid)
 
   OSS%id_frazil   = register_SIS_diag_field('ice_model', 'FRAZIL', diag%axesT1, Time, &
                'energy flux of frazil formation', 'W/m^2', conversion=US%QRZ_T_to_W_m2, missing_value=missing)
+
+  !CMOR frazil diagnostic
+  OSS%id_frazilmass   = register_SIS_diag_field('ice_model', 'sidmassgrowthwat', diag%axesT1, Time, &
+               'Sea-Ice Mass Change Through Growth in Supercooled Open Water (Frazil)', 'kg m-2 s-1', & 
+               conversion=US%RZ_to_kg_m2, missing_value=missing, & 
+               cmor_standard_name='tendency_of_sea_ice_amount_due_to_frazil_ice_accumulation_in_leads')
 
   if (coupler_type_initialized(OSS%tr_fields)) &
     call coupler_type_set_diags(OSS%tr_fields, 'ice_model', diag%axesT1%handles, Time)

src/SIS_ice_diags.F90

@@ -222,7 +222,7 @@ subroutine post_ice_state_diagnostics(IDs, IST, OSS, IOF, dt_slow, Time, G, US,
 
   ! Write out diagnostics of the ocean surface state, as seen by the slow sea ice.
   ! These fields do not change over the course of the sea-ice time stepping.
-  call post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag)
+  call post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag, IST)
 
   if (IDs%id_e2m>0) then
     tmp2d(:,:) = 0.0
@@ -264,15 +264,18 @@ end subroutine post_ice_state_diagnostics
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 !> Offer diagnostics of the ocean surface field, as seen by the sea ice.
-subroutine post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag)
+subroutine post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag, IST)
   type(ocean_sfc_state_type), intent(in)    :: OSS  !< A structure containing the arrays that describe
                                                     !! the ocean's surface state for the ice model.
   real,                       intent(in)    :: dt_slow  !< The time interval of these diagnostics [T ~> s]
   type(time_type),            intent(in)    :: Time     !< The ending time of these diagnostics
   type(SIS_hor_grid_type),    intent(inout) :: G    !< The horizontal grid type
   type(SIS_diag_ctrl),        pointer       :: diag !< A structure that is used to regulate diagnostic output
+  type(ice_state_type),       optional,intent(in) :: IST !< A type describing the state of the sea ice
 
-  real :: Idt_slow ! The inverse of the thermodynamic step [T-1 ~> s-1].
+  real :: Idt_slow  ! The inverse of the thermodynamic step [T-1 ~> s-1].
+  real :: LatHtFus  ! The latent heat of fusion of ice [Q ~> J kg-1].
+  real :: ILatHtFus ! The inverse of latent heat of fusion of ice [Q ~> kg J-1].
   Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow
 
   ! Write out diagnostics of the ocean surface state, as seen by the slow sea ice.
@@ -289,6 +292,11 @@ subroutine post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag)
   endif
   if (OSS%id_frazil>0) &
     call post_data(OSS%id_frazil, OSS%frazil*Idt_slow, diag)
+  if (OSS%id_frazilmass>0) then
+    call get_SIS2_thermo_coefs(IST%ITV, Latent_fusion=LatHtFus)
+    ILatHtFus = 1.0 / LatHtFus
+    call post_data(OSS%id_frazilmass, ILatHtFus*OSS%frazil*Idt_slow, diag)
+  endif
 
   if (coupler_type_initialized(OSS%tr_fields)) &
     call coupler_type_send_data(OSS%tr_fields, Time)