diff --git a/cicecore/cicedyn/analysis/ice_history.F90 b/cicecore/cicedyn/analysis/ice_history.F90
index 6c440cc86..34e5a9131 100644
--- a/cicecore/cicedyn/analysis/ice_history.F90
+++ b/cicecore/cicedyn/analysis/ice_history.F90
@@ -1496,42 +1496,42 @@ subroutine init_hist (dt)
          call define_hist_field(n_sithick,"sithick","m",tstr2D, tcstr, &
              "sea ice thickness",                             &
              "volume divided by area", c1, c0, &
-             ns1, f_sithick)
+             ns1, f_sithick, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siage,"siage","s",tstr2D, tcstr,    &
              "sea ice age",                             &
              "none", c1, c0,                                      &
-             ns1, f_siage)
+             ns1, f_siage, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sisnthick,"sisnthick","m",tstr2D, tcstr,    &
              "sea ice snow thickness",                             &
              "snow volume divided by area", c1, c0, &
-             ns1, f_sisnthick)
+             ns1, f_sisnthick, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sitemptop,"sitemptop","K",tstr2D, tcstr,    &
              "sea ice surface temperature", &
              "none", c1, c0,           &
-             ns1, f_sitemptop)
+             ns1, f_sitemptop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sitempsnic,"sitempsnic","K",tstr2D, tcstr,    &
              "snow ice interface temperature", &
              "surface temperature when no snow present", c1, c0, &
-             ns1, f_sitempsnic)
+             ns1, f_sitempsnic, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sitempbot,"sitempbot","K",tstr2D, tcstr,    &
              "sea ice bottom temperature",                             &
              "none", c1, c0,           &
-             ns1, f_sitempbot)
+             ns1, f_sitempbot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siu,"siu","m/s",ustr2D, ucstr,  &
              "ice x velocity component", &
              "none", c1, c0,         &
-             ns1, f_siu)
+             ns1, f_siu, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siv,"siv","m/s",ustr2D, ucstr,  &
              "ice y velocity component", &
              "none", c1, c0,         &
-             ns1, f_siv)
+             ns1, f_siv, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sidmasstranx,"sidmasstranx","kg/s",ustr2D, ucstr,  &
              "x component of snow and sea ice mass transport", &
@@ -1546,32 +1546,32 @@ subroutine init_hist (dt)
          call define_hist_field(n_sistrxdtop,"sistrxdtop","N m-2",ustr2D, ucstr,  &
              "x component of atmospheric stress on sea ice", &
              "none", c1, c0,         &
-             ns1, f_sistrxdtop)
+             ns1, f_sistrxdtop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sistrydtop,"sistrydtop","N m-2",ustr2D, ucstr,  &
              "y component of atmospheric stress on sea ice", &
              "none", c1, c0,         &
-             ns1, f_sistrydtop)
+             ns1, f_sistrydtop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sistrxubot,"sistrxubot","N m-2",ustr2D, ucstr,  &
              "x component of ocean stress on sea ice", &
              "none", c1, c0,         &
-             ns1, f_sistrxubot)
+             ns1, f_sistrxubot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sistryubot,"sistryubot","N m-2",ustr2D, ucstr,  &
              "y component of ocean stress on sea ice", &
              "none", c1, c0,         &
-             ns1, f_sistryubot)
+             ns1, f_sistryubot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sicompstren,"sicompstren","N m-1",tstr2D, tcstr,  &
              "compressive sea ice strength",                      &
              "none", c1, c0,         &
-             ns1, f_sicompstren)
+             ns1, f_sicompstren, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sispeed,"sispeed","m/s",ustr2D, ucstr, &
              "ice speed",                                  &
              "none", c1, c0,         &
-             ns1, f_sispeed)
+             ns1, f_sispeed, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sidir,"sidir","deg",ustr2D, ucstr,  &
              "ice direction",                                         &
@@ -1581,7 +1581,7 @@ subroutine init_hist (dt)
          call define_hist_field(n_sialb,"sialb","1",tstr2D, tcstr,  &
              "sea ice albedo",                                  &
              "none", c1, c0,         &
-             ns1, f_sialb)
+             ns1, f_sialb, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sihc,"sihc","J m-2",tstr2D, tcstr,  &
              "sea ice heat content",                                  &
@@ -1666,117 +1666,117 @@ subroutine init_hist (dt)
          call define_hist_field(n_siflswdtop,"siflswdtop","W/m2",tstr2D, tcstr, &
              "down shortwave flux over sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflswdtop)
+             ns1, f_siflswdtop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflswutop,"siflswutop","W/m2",tstr2D, tcstr, &
              "upward shortwave flux over sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflswutop)
+             ns1, f_siflswutop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflswdbot,"siflswdbot","W/m2",tstr2D, tcstr, &
              "down shortwave flux at bottom of ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflswdbot)
+             ns1, f_siflswdbot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sifllwdtop,"sifllwdtop","W/m2",tstr2D, tcstr, &
              "down longwave flux over sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_sifllwdtop)
+             ns1, f_sifllwdtop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sifllwutop,"sifllwutop","W/m2",tstr2D, tcstr, &
              "upward longwave flux over sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_sifllwutop)
+             ns1, f_sifllwutop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflsenstop,"siflsenstop","W/m2",tstr2D, tcstr, &
              "sensible heat flux over sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflsenstop)
+             ns1, f_siflsenstop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflsensupbot,"siflsensupbot","W/m2",tstr2D, tcstr, &
              "sensible heat flux at bottom of sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflsensupbot)
+             ns1, f_siflsensupbot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sifllatstop,"sifllatstop","W/m2",tstr2D, tcstr, &
              "latent heat flux over sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_sifllatstop)
+             ns1, f_sifllatstop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflcondtop,"siflcondtop","W/m2",tstr2D, tcstr, &
              "conductive heat flux at top of sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflcondtop)
+             ns1, f_siflcondtop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflcondbot,"siflcondbot","W/m2",tstr2D, tcstr, &
              "conductive heat flux at bottom of sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflcondbot)
+             ns1, f_siflcondbot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sipr,"sipr","kg m-2 s-1",tstr2D, tcstr, &
              "rainfall over sea ice", &
              "none", c1, c0,                            &
-             ns1, f_sipr)
+             ns1, f_sipr, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sifb,"sifb","m",tstr2D, tcstr, &
              "sea ice freeboard above sea level",                &
              "none", c1, c0,                            &
-             ns1, f_sifb)
+             ns1, f_sifb, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflsaltbot,"siflsaltbot","kg m-2 s-1",tstr2D, tcstr, &
              "salt flux from sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflsaltbot)
+             ns1, f_siflsaltbot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflfwbot,"siflfwbot","kg m-2 s-1",tstr2D, tcstr, &
              "fresh water flux from sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflfwbot)
+             ns1, f_siflfwbot, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siflfwdrain,"siflfwdrain","kg m-2 s-1",tstr2D, tcstr, &
              "fresh water drainage through sea ice", &
              "positive downward", c1, c0,                            &
-             ns1, f_siflfwdrain)
+             ns1, f_siflfwdrain, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sidragtop,"sidragtop","1",tstr2D, tcstr, &
              "atmospheric drag over sea ice", &
              "none", c1, c0,                            &
-             ns1, f_sidragtop)
+             ns1, f_sidragtop, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sirdgthick,"sirdgthick","m",tstr2D, tcstr, &
              "sea ice ridge thickness", &
              "vrdg divided by ardg", c1, c0, &
-             ns1, f_sirdgthick)
+             ns1, f_sirdgthick, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siforcetiltx,"siforcetiltx","N m-2",tstr2D, tcstr, &
              "sea surface tilt term", &
              "none", c1, c0,                            &
-             ns1, f_siforcetiltx)
+             ns1, f_siforcetiltx, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siforcetilty,"siforcetilty","N m-2",tstr2D, tcstr, &
              "sea surface tile term", &
              "none", c1, c0,                            &
-             ns1, f_siforcetilty)
+             ns1, f_siforcetilty, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siforcecoriolx,"siforcecoriolx","N m-2",tstr2D, tcstr, &
              "coriolis term",                                      &
              "none", c1, c0,                            &
-             ns1, f_siforcecoriolx)
+             ns1, f_siforcecoriolx, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siforcecorioly,"siforcecorioly","N m-2",tstr2D, tcstr, &
              "coriolis term",                                      &
              "none", c1, c0,                            &
-             ns1, f_siforcecorioly)
+             ns1, f_siforcecorioly, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siforceintstrx,"siforceintstrx","N m-2",tstr2D, tcstr, &
              "internal stress term", &
              "none", c1, c0,                            &
-             ns1, f_siforceintstrx)
+             ns1, f_siforceintstrx, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_siforceintstry,"siforceintstry","N m-2",tstr2D, tcstr, &
              "internal stress term", &
              "none", c1, c0,                            &
-             ns1, f_siforceintstry)
+             ns1, f_siforceintstry, avg_ice_present=.true., mask_ice_free_points=.true.)
 
          call define_hist_field(n_sistreave,"sistreave","N m-1",ustr2D, ucstr, &
              "average normal stress",                       &
@@ -1866,11 +1866,11 @@ subroutine init_hist (dt)
 
            call define_hist_field(n_siitdthick,"siitdthick","m",tstr3Dc, tcstr, &
               "ice thickness, categories","none", c1, c0, &
-              ns1, f_siitdthick)
+              ns1, f_siitdthick, avg_ice_present=.true.)
 
            call define_hist_field(n_siitdsnthick,"siitdsnthick","m",tstr3Dc, tcstr, &
               "snow thickness, categories","none", c1, c0, &
-              ns1, f_siitdsnthick)
+              ns1, f_siitdsnthick, avg_ice_present=.true.)
 
       endif ! if (histfreq(ns1) /= 'x') then
       enddo ! ns1
@@ -3654,501 +3654,29 @@ subroutine accum_hist (dt)
               enddo             ! j
 
               ! Only average for timesteps when ice present
-              if (index(avail_hist_fields(n)%vname,'sithick') /= 0) then
-                 if (f_sithick(1:1) /= 'x' .and. n_sithick(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sithick(ns),iblk) = &
-                             a2D(i,j,n_sithick(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sithick(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siage') /= 0) then
-                 if (f_siage(1:1) /= 'x' .and. n_siage(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siage(ns),iblk) = &
-                             a2D(i,j,n_siage(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siage(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sisnthick') /= 0) then
-                 if (f_sisnthick(1:1) /= 'x' .and. n_sisnthick(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sisnthick(ns),iblk) = &
-                             a2D(i,j,n_sisnthick(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sisnthick(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sitemptop') /= 0) then
-                 if (f_sitemptop(1:1) /= 'x' .and. n_sitemptop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sitemptop(ns),iblk) = &
-                             a2D(i,j,n_sitemptop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sitemptop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sitempsnic') /= 0) then
-                 if (f_sitempsnic(1:1) /= 'x' .and. n_sitempsnic(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sitempsnic(ns),iblk) = &
-                             a2D(i,j,n_sitempsnic(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sitempsnic(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sitempbot') /= 0) then
-                 if (f_sitempbot(1:1) /= 'x' .and. n_sitempbot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sitempbot(ns),iblk) = &
-                             a2D(i,j,n_sitempbot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sitempbot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siu') /= 0) then
-                 if (f_siu(1:1) /= 'x' .and. n_siu(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siu(ns),iblk) = &
-                             a2D(i,j,n_siu(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siu(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siv') /= 0) then
-                 if (f_siv(1:1) /= 'x' .and. n_siv(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siv(ns),iblk) = &
-                             a2D(i,j,n_siv(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siv(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sistrxdtop') /= 0) then
-                 if (f_sistrxdtop(1:1) /= 'x' .and. n_sistrxdtop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sistrxdtop(ns),iblk) = &
-                             a2D(i,j,n_sistrxdtop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sistrxdtop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sistrydtop') /= 0) then
-                 if (f_sistrydtop(1:1) /= 'x' .and. n_sistrydtop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sistrydtop(ns),iblk) = &
-                             a2D(i,j,n_sistrydtop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sistrydtop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sistrxubot') /= 0) then
-                 if (f_sistrxubot(1:1) /= 'x' .and. n_sistrxubot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sistrxubot(ns),iblk) = &
-                             a2D(i,j,n_sistrxubot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sistrxubot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sistryubot') /= 0) then
-                 if (f_sistryubot(1:1) /= 'x' .and. n_sistryubot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sistryubot(ns),iblk) = &
-                             a2D(i,j,n_sistryubot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sistryubot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sicompstren') /= 0) then
-                 if (f_sicompstren(1:1) /= 'x' .and. n_sicompstren(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sicompstren(ns),iblk) = &
-                             a2D(i,j,n_sicompstren(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sicompstren(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sispeed') /= 0) then
-                 if (f_sispeed(1:1) /= 'x' .and. n_sispeed(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sispeed(ns),iblk) = &
-                             a2D(i,j,n_sispeed(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sispeed(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
+              if (avail_hist_fields(n)%avg_ice_present) then
+                 do j = jlo, jhi
+                 do i = ilo, ihi
+                    if (tmask(i,j,iblk)) then
+                          a2D(i,j,n,iblk) = &
+                          a2D(i,j,n,iblk)*avgct(ns)*ravgip(i,j)
+                    endif
+                    ! Mask ice-free points
+                    if (avail_hist_fields(n)%mask_ice_free_points) then
+                       if (ravgip(i,j) == c0) a2D(i,j,n,iblk) = spval_dbl
+                    endif
+                 enddo             ! i
+                 enddo             ! j
               endif
+
+              ! CMIP albedo: also mask points below horizon
               if (index(avail_hist_fields(n)%vname,'sialb') /= 0) then
-                 if (f_sialb(1:1) /= 'x' .and. n_sialb(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sialb(ns),iblk) = &
-                             a2D(i,j,n_sialb(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sialb(ns),iblk) = spval_dbl
-                       if (albcnt(i,j,iblk,ns) <= puny) a2D(i,j,n_sialb(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflswdtop') /= 0) then
-                 if (f_siflswdtop(1:1) /= 'x' .and. n_siflswdtop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflswdtop(ns),iblk) = &
-                             a2D(i,j,n_siflswdtop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflswdtop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflswutop') /= 0) then
-                 if (f_siflswutop(1:1) /= 'x' .and. n_siflswutop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflswutop(ns),iblk) = &
-                             a2D(i,j,n_siflswutop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflswutop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflswdbot') /= 0) then
-                 if (f_siflswdbot(1:1) /= 'x' .and. n_siflswdbot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflswdbot(ns),iblk) = &
-                             a2D(i,j,n_siflswdbot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflswdbot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sifllwdtop') /= 0) then
-                 if (f_sifllwdtop(1:1) /= 'x' .and. n_sifllwdtop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sifllwdtop(ns),iblk) = &
-                             a2D(i,j,n_sifllwdtop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sifllwdtop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sifllwutop') /= 0) then
-                 if (f_sifllwutop(1:1) /= 'x' .and. n_sifllwutop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sifllwutop(ns),iblk) = &
-                             a2D(i,j,n_sifllwutop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sifllwutop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflsenstop') /= 0) then
-                 if (f_siflsenstop(1:1) /= 'x' .and. n_siflsenstop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflsenstop(ns),iblk) = &
-                             a2D(i,j,n_siflsenstop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflsenstop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflsensupbot') /= 0) then
-                 if (f_siflsensupbot(1:1) /= 'x' .and.  n_siflsensupbot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflsensupbot(ns),iblk) = &
-                             a2D(i,j,n_siflsensupbot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflsensupbot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sifllatstop') /= 0) then
-                 if (f_sifllatstop(1:1) /= 'x' .and. n_sifllatstop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sifllatstop(ns),iblk) = &
-                             a2D(i,j,n_sifllatstop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sifllatstop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sipr') /= 0) then
-                 if (f_sipr(1:1) /= 'x' .and. n_sipr(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sipr(ns),iblk) = &
-                             a2D(i,j,n_sipr(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sipr(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sifb') /= 0) then
-                 if (f_sifb(1:1) /= 'x' .and. n_sifb(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sifb(ns),iblk) = &
-                             a2D(i,j,n_sifb(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sifb(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflcondtop') /= 0) then
-                 if (f_siflcondtop(1:1) /= 'x' .and. n_siflcondtop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflcondtop(ns),iblk) = &
-                             a2D(i,j,n_siflcondtop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflcondtop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflcondbot') /= 0) then
-                 if (f_siflcondbot(1:1) /= 'x' .and. n_siflcondbot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflcondbot(ns),iblk) = &
-                             a2D(i,j,n_siflcondbot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflcondbot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflsaltbot') /= 0) then
-                 if (f_siflsaltbot(1:1) /= 'x' .and. n_siflsaltbot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflsaltbot(ns),iblk) = &
-                             a2D(i,j,n_siflsaltbot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflsaltbot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflfwbot') /= 0) then
-                 if (f_siflfwbot(1:1) /= 'x' .and. n_siflfwbot(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflfwbot(ns),iblk) = &
-                             a2D(i,j,n_siflfwbot(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflfwbot(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siflfwdrain') /= 0) then
-                 if (f_siflfwdrain(1:1) /= 'x' .and. n_siflfwdrain(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siflfwdrain(ns),iblk) = &
-                             a2D(i,j,n_siflfwdrain(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siflfwdrain(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-             if (index(avail_hist_fields(n)%vname,'sidragtop') /= 0) then
-                 if (f_sidragtop(1:1) /= 'x' .and. n_sidragtop(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sidragtop(ns),iblk) = &
-                             a2D(i,j,n_sidragtop(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sidragtop(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'sirdgthick') /= 0) then
-                 if (f_sirdgthick(1:1) /= 'x' .and. n_sirdgthick(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_sirdgthick(ns),iblk) = &
-                             a2D(i,j,n_sirdgthick(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_sirdgthick(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siforcetiltx') /= 0) then
-                 if (f_siforcetiltx(1:1) /= 'x' .and. n_siforcetiltx(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siforcetiltx(ns),iblk) = &
-                             a2D(i,j,n_siforcetiltx(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siforcetiltx(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siforcetilty') /= 0) then
-                 if (f_siforcetilty(1:1) /= 'x' .and. n_siforcetilty(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siforcetilty(ns),iblk) = &
-                             a2D(i,j,n_siforcetilty(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siforcetilty(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siforcecoriolx') /= 0) then
-                 if (f_siforcecoriolx(1:1) /= 'x' .and.  n_siforcecoriolx(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siforcecoriolx(ns),iblk) = &
-                             a2D(i,j,n_siforcecoriolx(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siforcecoriolx(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siforcecorioly') /= 0) then
-                 if (f_siforcecorioly(1:1) /= 'x' .and.  n_siforcecorioly(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siforcecorioly(ns),iblk) = &
-                             a2D(i,j,n_siforcecorioly(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siforcecorioly(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siforceintstrx') /= 0) then
-                 if (f_siforceintstrx(1:1) /= 'x' .and.  n_siforceintstrx(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siforceintstrx(ns),iblk) = &
-                             a2D(i,j,n_siforceintstrx(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siforceintstrx(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
-              if (index(avail_hist_fields(n)%vname,'siforceintstry') /= 0) then
-                 if (f_siforceintstry(1:1) /= 'x' .and.  n_siforceintstry(ns) /= 0) then
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a2D(i,j,n_siforceintstry(ns),iblk) = &
-                             a2D(i,j,n_siforceintstry(ns),iblk)*avgct(ns)*ravgip(i,j)
-                       endif
-                       if (ravgip(i,j) == c0) a2D(i,j,n_siforceintstry(ns),iblk) = spval_dbl
-                    enddo             ! i
-                    enddo             ! j
-                 endif
-              endif
+               do j = jlo, jhi
+               do i = ilo, ihi
+                  if (albcnt(i,j,iblk,ns) <= puny) a2D(i,j,n,iblk) = spval_dbl
+               enddo             ! i
+               enddo             ! j
+            endif
 
               ! back out albedo/zenith angle dependence
               if (avail_hist_fields(n)%vname(1:6) == 'albice') then
@@ -4259,33 +3787,17 @@ subroutine accum_hist (dt)
               enddo             ! i
               enddo             ! j
               enddo             ! k
-              if (index(avail_hist_fields(nn)%vname,'siitdthick') /= 0) then
-                 if (f_siitdthick(1:1) /= 'x' .and. n_siitdthick(ns)-n2D /= 0) then
-                    do k = 1, ncat_hist
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk) = &
-                             a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk)*avgct(ns)*ravgipn(i,j,k)
-                       endif
-                    enddo             ! i
-                    enddo             ! j
-                    enddo             ! k
-                 endif
-              endif
-              if (index(avail_hist_fields(nn)%vname,'siitdsnthick') /= 0) then
-                 if (f_siitdsnthick(1:1) /= 'x' .and.  n_siitdsnthick(ns)-n2D /= 0) then
-                    do k = 1, ncat_hist
-                    do j = jlo, jhi
-                    do i = ilo, ihi
-                       if (tmask(i,j,iblk)) then
-                             a3Dc(i,j,k,n_siitdsnthick(ns)-n2D,iblk) = &
-                             a3Dc(i,j,k,n_siitdsnthick(ns)-n2D,iblk)*avgct(ns)*ravgipn(i,j,k)
-                       endif
-                    enddo             ! i
-                    enddo             ! j
-                    enddo             ! k
-                 endif
+              if (avail_hist_fields(nn)%avg_ice_present) then
+                 do k = 1, ncat_hist
+                 do j = jlo, jhi
+                 do i = ilo, ihi
+                    if (tmask(i,j,iblk)) then
+                          a3Dc(i,j,k,n,iblk) = &
+                          a3Dc(i,j,k,n,iblk)*avgct(ns)*ravgipn(i,j,k)
+                    endif
+                 enddo             ! i
+                 enddo             ! j
+                 enddo             ! k
               endif
 
               endif
diff --git a/cicecore/cicedyn/analysis/ice_history_shared.F90 b/cicecore/cicedyn/analysis/ice_history_shared.F90
index 3c31f23ca..6d4850119 100644
--- a/cicecore/cicedyn/analysis/ice_history_shared.F90
+++ b/cicecore/cicedyn/analysis/ice_history_shared.F90
@@ -81,6 +81,8 @@ module ice_history_shared
           real (kind=dbl_kind) :: conb    ! additive conversion factor
           character (len=1) :: vhistfreq  ! frequency of history output
           integer (kind=int_kind) :: vhistfreq_n ! number of vhistfreq intervals
+          logical (kind=log_kind) :: avg_ice_present ! only average where ice is present
+          logical (kind=log_kind) :: mask_ice_free_points ! mask ice-free points
       end type
 
       integer (kind=int_kind), parameter, public :: &
@@ -811,7 +813,7 @@ end subroutine construct_filename
 
       subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
                                    vdesc, vcomment, cona, conb, &
-                                   ns, vhistfreq)
+                                   ns, vhistfreq, avg_ice_present, mask_ice_free_points)
 
       use ice_calendar, only: histfreq, histfreq_n
 
@@ -837,14 +839,28 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
       integer (kind=int_kind), intent(in) :: &
          ns             ! history file stream index
 
+      logical (kind=log_kind), optional, intent(in) :: &
+         avg_ice_present       , & ! compute average only when ice is present
+         mask_ice_free_points      ! mask ice-free points
+
       integer (kind=int_kind) :: &
          ns1        , & ! variable stream loop index
          lenf           ! length of namelist string
 
       character (len=40) :: stmp
 
+      logical (kind=log_kind) :: &
+         l_avg_ice_present       , & ! compute average only when ice is present
+         l_mask_ice_free_points      ! mask ice-free points
+
       character(len=*), parameter :: subname = '(define_hist_field)'
 
+      l_avg_ice_present = .false.
+      l_mask_ice_free_points = .false.
+
+      if(present(avg_ice_present)) l_avg_ice_present = avg_ice_present
+      if(present(mask_ice_free_points)) l_mask_ice_free_points = mask_ice_free_points
+
       if (histfreq(ns) == 'x') then
          call abort_ice(subname//'ERROR: define_hist_fields has histfreq x')
       endif
@@ -855,6 +871,10 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
       do ns1 = 1, lenf
          if (vhistfreq(ns1:ns1) == histfreq(ns)) then
 
+            if (ns1 > 1 .and. index(vhistfreq(1:ns1-1),'x') /= 0) then
+               call abort_ice(subname//'ERROR: history frequency variable f_' // vname // ' can''t contain ''x'' along with active frequencies')
+            endif
+
             num_avail_hist_fields_tot = num_avail_hist_fields_tot + 1
 
             if (vcoord(11:14) == 'time') then
@@ -917,6 +937,8 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
             avail_hist_fields(id(ns))%conb = conb
             avail_hist_fields(id(ns))%vhistfreq = vhistfreq(ns1:ns1)
             avail_hist_fields(id(ns))%vhistfreq_n = histfreq_n(ns)
+            avail_hist_fields(id(ns))%avg_ice_present = l_avg_ice_present
+            avail_hist_fields(id(ns))%mask_ice_free_points = l_mask_ice_free_points
 
          endif
       enddo
diff --git a/cicecore/cicedyn/general/ice_flux.F90 b/cicecore/cicedyn/general/ice_flux.F90
index 0fffa06b3..4c37a0696 100644
--- a/cicecore/cicedyn/general/ice_flux.F90
+++ b/cicecore/cicedyn/general/ice_flux.F90
@@ -1022,7 +1022,7 @@ end subroutine init_history_therm
 
       subroutine init_history_dyn
 
-      use ice_state, only: aice, vice, trcr, strength
+      use ice_state, only: aice, vice, trcr, strength, divu, shear
       use ice_grid,  only: grid_ice
 
       logical (kind=log_kind) :: &
@@ -1041,6 +1041,8 @@ subroutine init_history_dyn
 
       sig1    (:,:,:) = c0
       sig2    (:,:,:) = c0
+      divu    (:,:,:) = c0
+      shear   (:,:,:) = c0
       taubxU  (:,:,:) = c0
       taubyU  (:,:,:) = c0
       strength (:,:,:) = c0