From 509e2c33e95e3a2370dc406fb2fe4d06192420a6 Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Thu, 23 Nov 2023 13:09:04 -0500 Subject: [PATCH] ice_history: refactor CMIP history variables (#906) * ice_flux: zero-initialize divu and shear in init_history_dyn 'divu' and 'shear' are accessed in 'accum_hist' when writing the initial condition before they are initialized at the start of {eap, evp, implicit_solver}. This leads to runtime error when compiling with NaN initialization. Zero-initialize 'divu' and 'shear' in init_history_dyn, where the related variable 'strength' is already zero-initialized. * ice_history_shared: disallow 'x' in history frequency variables f_*' In the current code, nothing prevents users from leaving 'x' along with active frequencies in the individual namelist history frequency variables, for example: f_aice = 'xmd' This configuration does not work correctly, however. The corresponding history fields are correctly defined in ice_history_shared::define_hist_field, but since the calls to ice_history_shared::accum_hist_field in ice_history::accum_hist are only done after checking that the first element of each frequency variable is not 'x', the corresponding variables in the history files are all zero. Prevent that behaviour by actually disallowing 'x' in history frequency variables if any other frequencies are active. To implement that, add a check in the loop in define_hist_field, which loops through vhistfreq, (corresponding to f_aice, etc. in ice_history). Since this subroutine initializes 'id(:)' to zero and then writes a (non-zero) index in 'id' for any active frequency, it suffices to check that all previous indices are non-zero. * ice_history: remove uneeded conditions around CMIP history variables In ice_history::accum_hist, after the calls to accum_hist, we loop on the different output streams, and on the history variables in the avail_hist_fields array, to mask out land points and convert units for each output variable. Since 3c99e106 (Update CICE with CMIP changes. (#191), 2018-09-27), we also use this loop to do a special treatment for some CMIP variables (namely, averaging them only for time steps where ice is present, and masking points where ice is absent). This adjustment is done if the corresponding output frequency variable (f_sithick, etc.) does not have 'x' as its first element, and if the corresponding index in avail_hist_field for that variable/frequency (n_sithick(ns)) is not zero. Both conditions are in fact uneeded since they are always true. The first condition is always true because if the variable is found in the avail_hist_field array, which is ensured by the condition on line 3645, then necessarily its corresponding namelist output frequency won't have 'x' as its first character (since this is enforced in ice_history_shared::define_hist_field). The second condition is always true because if the variable is found in the avail_hist_field array, then necessarily its index in that array, n_(ns), is non-zero (see ice_history_shared::define_hist_field). Remove these uneeded conditions. This commit is best viewed with git show --color-moved --color-moved-ws=allow-indentation-change * ice_history: use loop index directly for CMIP variables In ice_history::accum_hist, there is a special treatment for some CMIP variables where they are averaged only for time steps where ice is present, and points where there is no ice are masked. This is done on the loop on output streams (with loop index n). This special averaging is done by accessing a2D and a3Dc using the variable n_(ns), which corresponds to the index in the avail_hist_field array where this history variable/frequency is defined. By construction, this index correponds to the loop index 'n', for both the 2D and the 3D loops. Simplify the code by using 'n' directly. * ice_history_shared: add two logical components to ice_hist_field At the end of ice_history::accum_hist, we do a special processing for some CMIP variables: we average them only for time steps where ice is present, and also mask ice-free points. The code to do that is repeated for each variable to which it applies. In order to reduce code duplication, let's introduce two new logical components to our 'ice_hist_field' type, defaulting them to .false., and make them optional arguments in ice_history_shared::define_hist_field. This allows us to avoid defining them for each output variable. We'll set them for CMIP variables in a following commit. * ice_history: set avg_ice_present, mask_ice_free_points for relevant CMIP variables In the previous commit, we added two components to type ice_hist_field (avg_ice_present and mask_ice_free_points), relating to some special treatment for CMIP variables (whether to average only for time steps where the ice is present and to mask ice-free points). Set these to .true. in the call to 'define_hist_field' for the relevant 2D variables [1], and set only 'avg_ice_present' to .true. for the 3D variables siitdthick and siitdsnthick, corresponding to the code under the "Mask out land points and convert units" loop in ice_history::accum_hist. [1] sithick siage sisnthick sitemptop sitempsnic sitempbot siu siv sidmasstranx sistrxdtop sistrydtop sistrxubot sistryubot sicompstren sispeed sidir sialb sihc siflswdtop siflswutop siflswdbot sifllwdtop sifllwutop siflsenstop siflsensupbot sifllatstop siflcondtop siflcondbot sipr sifb siflsaltbot siflfwbot siflfwdrain sidragtop sirdgthick siforcetiltx siforcetilty siforcecoriolx siforcecorioly siforceintstrx siforceintstry * ice_history: use avg_ice_present, mask_ice_free_points to reduce duplication Some CMIP variables are processed differently in ice_history::accum_hist: they are averaged only for time steps when ice is present, and points where ice is absent are masked. This processing is repeated for each of these variables in the 2D and 3Dc loops. To reduce code duplication, use the new components avg_ice_present and mask_ice_free_points of ice_hist_field to perform this processing only for variables that were defined accordingly. The relevant variables already have those components defined as of the previous commit. Note that we still need a separate loop for the variable 'sialb' (sea ice albedo) to mask points below the horizon. --- cicecore/cicedyn/analysis/ice_history.F90 | 632 ++---------------- .../cicedyn/analysis/ice_history_shared.F90 | 24 +- cicecore/cicedyn/general/ice_flux.F90 | 4 +- 3 files changed, 98 insertions(+), 562 deletions(-) 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