Merge PR MPAS-Dev#282 'xylar/ocean/time_varying_land_ice_pressure' in…

…to ocean/develop

This merge adds time-varying land-ice forcing via the following fields:

landIceDraft
landIcePressure
landIceFraction

A test case isomip_plus/2km/time_varying_Ocean0 is added as a
(relatively trivial) test of this forcing. In this test case,
landIcePressure and landIceDraft are initially 10% of their "correct"
values and increase linearly over 1 year to their full values, then
remain constant for an additional year. However, the test case crashes
sometime during the second year because layers become too thin and/or
tilted, or velocities become too large. This is an intentional feature
of the test case, showing why further work is needed to handle coupled
ice dynamics.

Modifications are made to the ISOMIP+ viz (MPAS-Dev#275) to correctly handle
repeated running of the viz while the simulation is in progress. A few
other tweaks (e.g. the format of dates in the figure titles) are also
included.

Modifications were made to how sea-level is maintained in ISOMIP+ test
cases despite freshwater fluxes and changes in land-ice pressure. A
fictitious evaporation was already used at the northern boundary, but
this has been modified to adjust to sea-level itself rather than melt
fluxes to accommodate changing land-ice pressure.

Finally, a division-by-zero error in timeSeriesStatsGlobal was fixed. If
melt fluxes are turned on but there are no locations where the melt mask
nonzero (which happened in my testing), global stats was trying to
compute mean fluxes and other stats with zero area in the denominator.
This has been fixed.

Registry.xml

@@ -758,6 +758,30 @@
 			description="Number of days to ramp up time varying forcing"
 			possible_values="Any positive real number"
 		/>
+		<nml_option name="config_use_time_varying_land_ice_forcing" type="logical" default_value=".false." units="unitless"
+			description="If true calculate input forcing fields."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_time_varying_land_ice_forcing_start_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+			description="Forcing time to use at the simulation start time"
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_land_ice_forcing_reference_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+			description="Reference time for the forcing"
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_land_ice_forcing_cycle_start" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+			description="Start time for the forcing cycle."
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_land_ice_forcing_cycle_duration" type="character" default_value="2-00-00_00:00:00" units="unitless"
+			description="Duration of the forcing cycle."
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_land_ice_forcing_interval" type="character" default_value="01:00:00" units="unitless"
+			description="Time between forcing inputs"
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
 	</nml_record>
 	<nml_record name="coupling" mode="init;forward">
 		<nml_option name="config_ssh_grad_relax_timescale" type="real" default_value="0.0" units="seconds"
@@ -1306,6 +1330,7 @@
 	</nml_record>
 	<packages>
 		<package name="timeVaryingAtmosphericForcingPKG" description="This package includes variables required for time varying atmospheric forcing"/>
+		<package name="timeVaryingLandIceForcingPKG" description="This package includes variavles required for time varying land-ice forcing"/>
 		<package name="variableShortwave" description="This package includes variables required to compute spatially variable shortwave extinction coefficients"/>
 
 		<package name="splitTimeIntegrator" description="This package includes variables required for either the split or unsplit explicit time integrators."/>
@@ -1713,6 +1738,17 @@
 			<var name="windSpeedU"/>
 			<var name="windSpeedV"/>
 		</stream>
+		<stream name="land_ice_forcing"
+				type="input"
+				filename_template="land_ice_forcing.nc"
+				filename_interval="none"
+				input_interval="none"
+				packages="timeVaryingLandIceForcingPKG"
+				immutable="true">
+			<var name="landIcePressureForcing"/>
+			<var name="landIceFractionForcing"/>
+			<var name="landIceDraftForcing"/>
+		</stream>
 
 		<!-- Diagnostics streams for forward mode -->
 		<stream name="additional_output"
@@ -3125,6 +3161,18 @@
 			description="Pressure at the sea surface due to the atmosphere."
 			packages="timeVaryingAtmosphericForcingPKG"
 		/>
+		<var name="landIceFractionForcing" type="real" dimensions="nCells Time" units="unitless"
+			description="The fraction of each cell covered by land ice"
+			packages="timeVaryingLandIceForcingPKG"
+		/>
+		<var name="landIcePressureForcing" type="real" dimensions="nCells Time" units="Pa"
+			description="Pressure defined at the sea surface due to land ice."
+			packages="timeVaryingLandIceForcingPKG"
+		/>
+		<var name="landIceDraftForcing" type="real" dimensions="nCells Time" units="m"
+			description="The elevation of the interface between land ice and the ocean."
+			packages="timeVaryingLandIceForcingPKG"
+		/>
 	</var_struct>
 	<var_struct name="scratch" time_levs="1">
 		<var name="normalThicknessFlux" persistence="scratch" type="real" dimensions="nVertLevels nEdges Time"

analysis_members/mpas_ocn_global_stats.F

@@ -1315,7 +1315,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
 
       ! landIceFreshwaterFlux
       variableIndex = variableIndex + 1
-      if (associated(landIceFreshwaterFlux)) then
+      if (associated(landIceFreshwaterFlux) .and. landIceFloatingAreaSum > 0.0_RKIND) then
          averages(variableIndex) = sums(variableIndex)/(landIceFloatingAreaSum)
          rms(variableIndex) = sqrt(sumSquares(variableIndex)/(landIceFloatingAreaSum))
       else
@@ -1328,7 +1328,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
 
       ! accumulatedLandIceMass
       variableIndex = variableIndex + 1
-      if (associated(accumulatedLandIceMass)) then
+      if (associated(accumulatedLandIceMass) .and. landIceFloatingAreaSum > 0.0_RKIND) then
          averages(variableIndex) = sums(variableIndex)/(landIceFloatingAreaSum)
          rms(variableIndex) = sqrt(sumSquares(variableIndex)/(landIceFloatingAreaSum))
       else
@@ -1338,7 +1338,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
 
       ! accumulatedLandIceHeat
       variableIndex = variableIndex + 1
-      if (associated(accumulatedLandIceHeat)) then
+      if (associated(accumulatedLandIceHeat) .and. landIceFloatingAreaSum > 0.0_RKIND) then
          averages(variableIndex) = sums(variableIndex)/(landIceFloatingAreaSum)
          rms(variableIndex) = sqrt(sumSquares(variableIndex)/(landIceFloatingAreaSum))
       else
@@ -1348,7 +1348,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
 
       ! accumulatedLandIceFrazilMass
       variableIndex = variableIndex + 1
-      if (landIceFluxesPkgActive .and. frazilIcePkgActive) then
+      if (landIceFluxesPkgActive .and. frazilIcePkgActive .and. landIceFloatingAreaSum > 0.0_RKIND) then
          averages(variableIndex) = sums(variableIndex)/(landIceFloatingAreaSum)
          rms(variableIndex) = sqrt(sumSquares(variableIndex)/(landIceFloatingAreaSum))
       else

driver/mpas_ocn_core_interface.F

@@ -124,6 +124,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: variableShortwaveActive
       logical, pointer :: gmActive
       logical, pointer :: timeVaryingAtmosphericForcingPKGActive
+      logical, pointer :: timeVaryingLandIceForcingPKGActive
 
       type (mpas_pool_iterator_type) :: pkgItr
       logical, pointer :: packageActive
@@ -149,6 +150,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: config_use_tidal_forcing
       logical, pointer :: config_use_standardGM
       logical, pointer :: config_use_time_varying_atmospheric_forcing
+      logical, pointer :: config_use_time_varying_land_ice_forcing
 
       character (len=StrKIND), pointer :: config_time_integrator
       character (len=StrKIND), pointer :: config_ocean_run_mode
@@ -312,6 +314,12 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
          timeVaryingAtmosphericForcingPKGActive = .true.
       endif
 
+      call mpas_pool_get_package(packagePool, 'timeVaryingLandIceForcingPKGActive', timeVaryingLandIceForcingPKGActive)
+      call mpas_pool_get_config(configPool, 'config_use_time_varying_land_ice_forcing', config_use_time_varying_land_ice_forcing)
+      if (config_use_time_varying_land_ice_forcing) then
+         timeVaryingLandIceForcingPKGActive = .true.
+      endif
+
       !
       ! call into analysis member driver to set analysis member packages
       !

mode_forward/mpas_ocn_forward_mode.F

@@ -489,6 +489,7 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
       logical, pointer :: config_use_activeTracers_surface_restoring
       logical, pointer :: config_use_surface_salinity_monthly_restoring
       logical, pointer :: config_use_time_varying_atmospheric_forcing
+      logical, pointer :: config_use_time_varying_land_ice_forcing
 
       ierr = 0
 
@@ -502,6 +503,8 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
                                                    config_use_surface_salinity_monthly_restoring)
       call MPAS_pool_get_config(domain % configs, 'config_use_time_varying_atmospheric_forcing', &
                                                    config_use_time_varying_atmospheric_forcing)
+      call MPAS_pool_get_config(domain % configs, 'config_use_time_varying_land_ice_forcing', &
+                                                   config_use_time_varying_land_ice_forcing)
 
       ! Eventually, dt should be domain specific
       timeStep = mpas_get_clock_timestep(domain % clock, ierr=ierr)
@@ -644,7 +647,8 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
             if (trim(config_sw_absorption_type)=='ohlmann00' .or. &
                 config_use_ecosysTracers .or. &
                 (config_use_activeTracers_surface_restoring .and. config_use_surface_salinity_monthly_restoring) .or. &
-                config_use_time_varying_atmospheric_forcing) then
+                config_use_time_varying_atmospheric_forcing .or. &
+                config_use_time_varying_land_ice_forcing) then
 
                 call ocn_time_varying_forcing_write_restart_times(domain)
             end if