merge in latest dev/gfdl updates (#36)

* TC4 integration into test suite

This patch renames the tc4 test to activate it in the test suite.  It
also modifies the Makefile to build the input field test scripts.

It also modifies the Python build scripts to be PEP8-conformant.

We temporarily disable tc4 in the restart tests, since they currently
fail.  This needs to be addressed before we can merge this into the main
branch.

The patch does not enable the necessary Python modules for running on
Travis, that will also be addressed later.

* Travis python support; tc4 Makefile

The custom TC4 Makefile has been added (oops), and the presumed Python
Ubuntu packages have been added for Travis.

* Verify ENABLE_THERMODYNAMICS is True before posting C_p diagnostic

* Make tc4 faster

* remove trailing whitespace

* add unit scaling

* fix restart fail for tc4 and some cleanup

* remove trailiny ws

* Enable tc4.restart test

* +Pass timeesteps to tracer diagnostics in [T]

  Pass timeesteps to the tracer diagnistics routines post_tracer_diagnostics and
postALE_tracer_diagnostics and to adiabatic in units of [T}.  All answers are
bitwise identical.

* +Rescaled tracer advective flux diagnostics

  Rescaled the internal units of the tracer advective flux diagnostics to units
of [conc H L2 T-1] for code simplicity and dimensional consistency testing.
Also corrected the units of some tracer fluxes as documented in comments and
commented out unused elements of the tracer_type.  All answers are bitwise
identical.

* +Pass timesteps to ALE_main in [T]

  Pass the timesteps to ALE_main, ALE_main_offline, and ALE_main_accelerated in
units of [T] for code simplicity and dimensional consistency testing.  This also
includes the rescaling of remapping-driven tracer tendencies.  All answers and
diagnostics are bitwise identical.

* +Pass timesteps to tracer column_physics in [T]

  Pass timesteps to the various tracer column_physics routines in [T] for
dimensional consistency testing.  Also added a new unit_scale_type argument to
these routines.  All answers are bitwise identical, but there are minor
interface changes to 13 subroutines.

* +Pass timesteps to applyTracerBoundaryFluxesInOut in [T]

  Pass timesteps to applyTracerBoundaryFluxesInOut in [T], and use units of
[T-1] for internal source and decay rates for the oil tracer and in fluxes of
CFCs.  Also modified extract_offline_main to return timesteps as real values
with units of [T].  Also there is a new unit_scale_type argument to
register_oil_tracer.   All answers in the MOM6_examples test cases and
regression tests are bitwise identical.

* Simplified expressions in MOM_PointAccel

  Simplified expressions inside of MOM_PointAccel, taking into account that all
velocities use the same units of [L T-1].  All answers are bitwise identical.

* Corrected dimensional epsilons in downscaling

  Added distinct negligible volumes, face areas, horizonal areas and lengths
with proper dimensional rescaling in the downsample field routines.  With these
changes, downscaled diagnostics should now pass the dimensional rescaling tests,
whereas previously there would have been a problem when the numbers used to
represent lengths are smaller than about 1e-8 times their MKS values.  All
answers are bitwise identical without dimensional rescaling.

* Simplified expressions in MOM_offline_aux

  Simplified expressions in distribute_residual_uh_barotropic.  All answers are
bitwise identical.

* Revised wave_speed to return speed in [L T-1]

  Revised wave_speed to return the internal wave speed in units of [L T-1] and
to use mono_N2_depth in units of [Z] for code simplification and expanded
dimensional consistency testing.  Also revised the internal units of some
related diagnostics in calculate_diagnostic_fields.  All answers and diagnostics
are bitwise identical.

* Rescaled internal variables in wave_speed

  Rescale internal calculations in wave_speed and wave_speeds for greater
robustness and dimensional consistency testing.  All answers are bitwise
identical and pass dimensional scaling tests.

* +Changed the units of minimum_forcing_depth to [H]

  Changed the units of minimum_forcing_depth passed to applyBoundaryFluxesInOut
and applyTracerBoundaryFluxesInOut to [H].  All answers are bitwise identical.

* Correction of documented units in comments

  Corrected some units in comments and eliminated some unused variables.
All answers are bitwise identical.

* Adiabatic clock ID bugfix

This patch fixes an initialization bug of the diabatic timer, which was
being used to measure adiabatic time but was never initialized if the
experiment was configured as adiabatic.

We fix this by introducing a separate timer for the adiabatic solver.
Although we could have reused the diabatic timer, the addition of a new
variable should not add any overhead on modern compilers.

* Corrected an OMP declaration

  Added a variable to an OMP declaration.  All answers are bitwise identical,
and a recently added compile-time error with openMP was fixed.

* Update MOM.F90

Fixed Alistair's embarrassing error.

* Dimensional rescaling in MOM_open_boundary.F90

  Added rescaling for dimensional consistency testing in MOM_open_boundary.F90,
including splitting variables with different units that had previously shared
the same variable and adding more extensive documentation of variables.  Also
changed the dimensions of the timesteps passed to radiation_open_bdry_conds and
update_segment_tracer_reservoirs to [T] and added vertical_grid_type and
unit_scale_type arguments to open_boundary_init and open_boundary_test_extern_h.
All answers are bitwise identical, although some probably bugs have been noted
in comments and there are new or altered arguments to several routines.

* (*)Fixed invariance bugs in MOM_open_boundary.F90

  Corrected dimensional consistency bugs in update_segment_tracer_reservoirs and
horizontal indexing and related bugs in gradient_at_q_points with oblique_grad
OBCs.  These will both change answers in test cases that use some open boundary
condition options, but not in any of the MOM6-examples test cases.

.testing/Makefile

@@ -251,6 +251,7 @@ results/%/ocean.stats.$(1): ../build/$(2)/MOM6
 	if [ $(3) ]; then find ../build/$(2) -name *.gcda -exec rm -f '{}' \; ; fi
 	mkdir -p work/$$*/$(1)
 	cp -rL $$*/* work/$$*/$(1)
+	cd work/$$*/$(1) && if [ -f Makefile ]; then make; fi
 	mkdir -p work/$$*/$(1)/RESTART
 	echo $(4) > work/$$*/$(1)/MOM_override
 	cd work/$$*/$(1) && $$(call MPIRUN_CMD,$(5)) -n $(6) ../../../$$< 2> debug.out > std.out \
@@ -285,6 +286,7 @@ results/%/ocean.stats.restart: ../build/symmetric/MOM6
 	rm -rf work/$*/restart
 	mkdir -p work/$*/restart
 	cp -rL $*/* work/$*/restart
+	cd work/$*/restart && if [ -f Makefile ]; then make; fi
 	mkdir -p work/$*/restart/RESTART
 	# Generate the half-period input namelist
 	# TODO: Assumes runtime set by DAYMAX, will fail if set by input.nml
@@ -294,20 +296,19 @@ results/%/ocean.stats.restart: ../build/symmetric/MOM6
 		&& if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \
 		&& printf -v timeunit_int "%.f" "$${timeunit}" \
 		&& halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \
-		&& printf "\n&ocean_solo_nml\n    seconds = $${halfperiod}\n/\n" >> input.nml \
-		&& echo $${daymax} $${timeunit}
+		&& printf "\n&ocean_solo_nml\n    seconds = $${halfperiod}\n/\n" >> input.nml
 	# Run the first half-period
-	cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug.out > std.out \
-		|| ! sed 's/^/$*.restart1: /' std.out debug.out \
-		&& sed 's/^/$*.restart1: /' std.out
+	cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug1.out > std1.out \
+		|| ! sed 's/^/$*.restart1: /' std1.out debug1.out \
+		&& sed 's/^/$*.restart1: /' std1.out
 	# Setup the next inputs
 	cd work/$*/restart && rm -rf INPUT && mv RESTART INPUT
 	mkdir work/$*/restart/RESTART
 	cd work/$*/restart && sed -i -e "s/input_filename *= *'n'/input_filename = 'r'/g" input.nml
 	# Run the second half-period
-	cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug.out > std.out \
-		|| ! sed 's/^/$*.restart2: /' std.out debug.out \
-		&& sed 's/^/$*.restart2: /' std.out
+	cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug2.out > std2.out \
+		|| ! sed 's/^/$*.restart2: /' std2.out debug2.out \
+		&& sed 's/^/$*.restart2: /' std2.out
 	# Archive the results and cleanup
 	mkdir -p $(@D)
 	cp work/$*/restart/ocean.stats $@

.testing/_tc4/MOM_input → .testing/tc4/MOM_input

@@ -7,10 +7,15 @@
 USE_REGRIDDING = True           !   [Boolean] default = False
                                 ! If True, use the ALE algorithm (regridding/remapping). If False, use the
                                 ! layered isopycnal algorithm.
-DT = 300.0                      !   [s]
+DT = 1200.0                      !   [s]
                                 ! The (baroclinic) dynamics time step.  The time-step that is actually used will
                                 ! be an integer fraction of the forcing time-step (DT_FORCING in ocean-only mode
                                 ! or the coupling timestep in coupled mode.)
+DT_THERM = 3600.0                !   [s] default = 300.0
+                                ! The thermodynamic and tracer advection time step. Ideally DT_THERM should be
+                                ! an integer multiple of DT and less than the forcing or coupling time-step,
+                                ! unless THERMO_SPANS_COUPLING is true, in which case DT_THERM can be an integer
+                                ! multiple of the coupling timestep.  By default DT_THERM is set to DT.				
 C_P = 3925.0                    !   [J kg-1 K-1] default = 3991.86795711963
                                 ! The heat capacity of sea water, approximated as a constant. This is only used
                                 ! if ENABLE_THERMODYNAMICS is true. The default value is from the TEOS-10
@@ -377,10 +382,15 @@ WIND_CONFIG = "zero"            !
 ! === module MOM_restart ===
 
 ! === module MOM_main (MOM_driver) ===
-DAYMAX = 1.0                    !   [days]
+DAYMAX = 0.25                    !   [days]
                                 ! The final time of the whole simulation, in units of TIMEUNIT seconds.  This
                                 ! also sets the potential end time of the present run segment if the end time is
                                 ! not set via ocean_solo_nml in input.nml.
+
+ENERGYSAVEDAYS = 0.125          !   [days] default = 1.44E+04
+                                ! The interval in units of TIMEUNIT between saves of the
+                                ! energies of the run and other globally summed diagnostics.
+
 RESTART_CONTROL = 3             ! default = 1
                                 ! An integer whose bits encode which restart files are written. Add 2 (bit 1)
                                 ! for a time-stamped file, and odd (bit 0) for a non-time-stamped file. A
@@ -397,3 +407,6 @@ MAXCPU = 2.88E+04               !   [wall-clock seconds] default = -1.0
                                 ! processors used.
 
 ! === module MOM_file_parser ===
+
+DIAG_AS_CHKSUM = True
+DEBUG = True

.testing/tc4/Makefile

@@ -0,0 +1,3 @@
+all:
+	python build_grid.py
+	python build_data.py

.testing/tc4/build_data.py

@@ -0,0 +1,80 @@
+import netCDF4 as nc
+import numpy as np
+
+x = nc.Dataset('ocean_hgrid.nc').variables['x'][1::2, 1::2]
+y = nc.Dataset('ocean_hgrid.nc').variables['y'][1::2, 1::2]
+zbot = nc.Dataset('topog.nc').variables['depth'][:]
+zbot0 = zbot.max()
+
+
+def t_fc(x, y, z, radius=5.0, tmag=1.0):
+    """a radially symmetric anomaly in the center of the domain.
+    units are meters and degC.
+    """
+    ny, nx = x.shape
+    nz = z.shape[0]
+
+    x0 = x[int(ny/2), int(nx/2)]
+    y0 = y[int(ny/2), int(nx/2)]
+
+    tl = np.zeros((nz, ny, nx))
+    zb = z[-1]
+    if len(z) > 1:
+        zd = z / zb
+    else:
+        zd = [0.]
+    for k in np.arange(len(zd)):
+        r = np.sqrt((x - x0)**2 + (y - y0)**2)
+        tl[k, :] += (1.0 - np.minimum(r / radius, 1.0)) * tmag * (1.0 - zd[k])
+    return tl
+
+
+ny, nx = x.shape
+nz = 3
+z = (np.arange(nz) * zbot0) / nz
+
+temp = t_fc(x, y, z)
+salt = np.zeros(temp.shape)+35.0
+fl = nc.Dataset('temp_salt_ic.nc', 'w', format='NETCDF3_CLASSIC')
+fl.createDimension('lon', nx)
+fl.createDimension('lat', ny)
+fl.createDimension('depth', nz)
+fl.createDimension('Time', None)
+zv = fl.createVariable('depth', 'f8', ('depth'))
+lonv = fl.createVariable('lon', 'f8', ('lon'))
+latv = fl.createVariable('lat', 'f8', ('lat'))
+timev = fl.createVariable('Time', 'f8', ('Time'))
+timev.calendar = 'noleap'
+timev.units = 'days since 0001-01-01 00:00:00.0'
+timev.modulo = ' '
+tv = fl.createVariable('ptemp', 'f8', ('Time', 'depth', 'lat', 'lon'),
+                       fill_value=-1.e20)
+sv = fl.createVariable('salt', 'f8', ('Time', 'depth', 'lat', 'lon'),
+                       fill_value=-1.e20)
+tv[:] = temp[np.newaxis, :]
+sv[:] = salt[np.newaxis, :]
+zv[:] = z
+lonv[:] = x[0, :]
+latv[:] = y[:, 0]
+timev[0] = 0.
+fl.sync()
+fl.close()
+
+
+# Make Sponge forcing file
+dampTime = 20.0     # days
+secDays = 8.64e4
+fl = nc.Dataset('sponge.nc', 'w', format='NETCDF3_CLASSIC')
+fl.createDimension('lon', nx)
+fl.createDimension('lat', ny)
+lonv = fl.createVariable('lon', 'f8', ('lon'))
+latv = fl.createVariable('lat', 'f8', ('lat'))
+spv = fl.createVariable('Idamp', 'f8', ('lat', 'lon'), fill_value=-1.e20)
+Idamp = np.zeros((ny, nx))
+if dampTime > 0.:
+    Idamp = 0.0 + 1.0 / (dampTime * secDays)
+spv[:] = Idamp
+lonv[:] = x[0, :]
+latv[:] = y[:, 0]
+fl.sync()
+fl.close()