Merge Unit testing work from RESPEC

.gitignore

@@ -34,14 +34,11 @@ __pycache__/*
 htmlcov
 .tox
 .tox/*
-.coverage
-.coverage/*
 
-# Translations
-*.mo
-
-
-# Stashes, etc.
-.hg
-*.tmp*
-tests/_LargeFileStash
+stats.dat
+.ropeproject
+.idea
+*.ech
+*.log
+*.hbnhead
+*.units.dbf

HSP2/ADCALC.py

@@ -4,6 +4,7 @@
 '''
 
 from numpy import zeros
+from numba import njit
 from HSP2.utilities import make_numba_dict
 
 # The clean way to get calculated data from adcalc() into advert() is to use a closure. 
@@ -127,7 +128,7 @@ def adcalc_(simlen, delts, nexits, crrat, ks, vol, ADFG, O, VOL, SROVOL, EROVOL,
 	return
 
 
-#@jit(nopython=True)
+#@njit(cache=True)
 def advect(imat, conc, nexits, vols, vol, srovol, erovol, sovol, eovol):
 	''' Simulate advection of constituent totally entrained in water.
 	Originally designed to be called as: advect(loop, imat, conc, omat, *ui['adcalcData'])
@@ -140,9 +141,7 @@ def advect(imat, conc, nexits, vols, vol, srovol, erovol, sovol, eovol):
 	# sovol  = SOVOL[loop,:]
 	# eovol  = EOVOL[loop,:]
 
-	omat = 0.0
-	if nexits > 1:
-		omat = zeros((nexits))
+	omat = zeros((nexits))
 	if vol > 0.0:    # reach/res contains water
 		concs = conc
 		conc = (imat + concs * (vols - srovol)) / (vol + erovol)  # material entering during interval, weighted volume of outflow based on conditions at start of ivl (srovol), and weighted volume of outflow based on conditions at end of ivl (erovol)

HSP2/GQUAL.py

@@ -1245,7 +1245,7 @@ def advqal(isqal,rsed,bsed,depscr,rosed,osed,nexits,rsqals,rbqals,errors):
 		if bsed <= 0.0:     # no bed sediments at end of interval
 			bqal = -1.0e30
 			if abs(dsqal) > 0.0 or abs(rbqals > 0.0):
-				errors[4] += 1  # ERRMSG4: error-under these conditions these values should be zero
+				errors[5] += 1  # ERRMSG4: error-under these conditions these values should be zero
 		else:     # there is bed sediment at the end of the interval
 			rbqal= dsqal + rbqals
 			bqal = rbqal / bsed

HSP2/HTRCH.py

@@ -34,12 +34,12 @@
   0.0026, 0.0026, 0.0027, 0.0027, 0.0028, 0.0028, 0.0027, 0.0026, 0.0024, 0.0023, 0.0022,
   0.0021, 0.0021, 0.0020]
 
-ERRMSG = []
+ERRMSGS =('HTRCH: Water temperature is above 66 C (150 F) -- In most cases, this indicates an instability in advection','')     #ERRMSG0
 
 def htrch(store, siminfo, uci, ts):
 	'''Simulate heat exchange and water temperature'''
 
-	errorsV = zeros(len(ERRMSG), dtype=int)
+	errorsV = zeros(len(ERRMSGS), dtype=int)
 
 	advectData = uci['advectData']
 	(nexits, vol, VOL, SROVOL, EROVOL, SOVOL, EOVOL) = advectData
@@ -172,7 +172,7 @@ def htrch(store, siminfo, uci, ts):
 
 		if tw > 66.0:
 			if adfg < 2:
-				# errormsg:  'advect:tw problem ',dtw, airtmp
+				errorsV[0] += 1  # Water temperature is above 66 C (150 F).  In most cases, this indicates an instability in advection.
 				rtw = tw
 			else:
 				tw = tws
@@ -369,7 +369,7 @@ def htrch(store, siminfo, uci, ts):
 		for i in range(nexits):
 			ts['OHEAT' + str(i+1)] = OHEAT[:, i]
 
-	return errorsV, ERRMSG
+	return errorsV, ERRMSGS
 
 def vapor(tmp):
 		'''	# define vapor function based on temperature (deg c); vapor pressure is expressed in millibars'''

HSP2/IQUAL.py

@@ -19,7 +19,7 @@
 UNDEFINED: sliqsp
 '''
 
-ERRMSG = []
+ERRMSGS =('IQUAL: A constituent must be associated with overland flow in order to receive atmospheric deposition inputs','')     #ERRMSG0
 
 def iqual(store, siminfo, uci, ts):
 	''' Simulate washoff of quality constituents (other than solids, Heat, dox, and co2)
@@ -35,7 +35,7 @@ def iqual(store, siminfo, uci, ts):
 		flags = uci['IQUAL' + iqual + '_FLAGS']
 		constituents.append(flags['QUALID'])
 
-	errorsV = zeros(len(ERRMSG), dtype=int)
+	errorsV = zeros(len(ERRMSGS), dtype=int)
 	delt60 = siminfo['delt'] / 60     # delt60 - simulation time interval in hours
 	simlen = siminfo['steps']
 	tindex = siminfo['tindex']
@@ -118,6 +118,9 @@ def iqual(store, siminfo, uci, ts):
 			elif iqadfgc == -1:
 				ts['IQADCN'] = ts['IQADCN' + str(index) + ' 1']
 
+			if QSOFG == 0 and (iqadfgf != 0 or iqadfgc != 0):
+				errorsV[0] += 1  # error - non-qualof cannot have atmospheric deposition
+
 		if 'IQADFX' not in ts:
 			ts['IQADFX'] = zeros(simlen)
 		if 'IQADCN' not in ts:
@@ -229,5 +232,5 @@ def iqual(store, siminfo, uci, ts):
 			IQADDR[loop] = adfxfx
 			IQADEP[loop] = adtot
 
-	return errorsV, ERRMSG
+	return errorsV, ERRMSGS
 

HSP2/IWTGAS.py

@@ -125,17 +125,23 @@ def iwtgas(store, siminfo, uci, ts):
 				if slico2 >= 0.0:     # there is co2 conc of surface lateral inflow
 					soco2 = slico2 * slifac + soco2 * (1.0 - slifac)
 
+			SOHT[loop]   = sotmp * suro * EFACTA  # compute outflow of heat energy in water - units are deg. c-in./ivl
+			SODOXM[loop] = sodox * suro * MFACTA  # calculate outflow mass of dox - units are mg-in./l-ivl
+			SOCO2M[loop] = soco2 * suro * MFACTA  # calculate outflow mass of co2 - units are mg-in./l-ivl
+			SOTMP[loop]  = (sotmp * 9.0 / 5.0) + 32.0
+			SODOX[loop]  = sodox
+			SOCO2[loop]  = soco2
+
 		else:
 			sotmp = -1.0e30
 			sodox = -1.0e30
 			soco2 = -1.0e30
 
-		SOHT[loop]   = sotmp * suro * EFACTA  # compute outflow of heat energy in water - units are deg. c-in./ivl
-		SODOXM[loop] = sodox * suro * MFACTA  # calculate outflow mass of dox - units are mg-in./l-ivl
-		SOCO2M[loop] = soco2 * suro * MFACTA  # calculate outflow mass of co2 - units are mg-in./l-ivl
-
-		SOTMP[loop]  = (sotmp * 9.0 / 5.0) + 32.0
-		SODOX[loop]  = sodox
-		SOCO2[loop]  = soco2
+			SOHT[loop]   = 0.0
+			SODOXM[loop] = 0.0
+			SOCO2M[loop] = 0.0
+			SOTMP[loop]  = sotmp
+			SODOX[loop]  = sodox
+			SOCO2[loop]  = soco2
 
 	return errorsV, ERRMSG

HSP2/PQUAL.py

@@ -16,7 +16,7 @@
 NEED to check all units conversions
 '''
 
-ERRMSG = []
+ERRMSGS =('PQUAL: A constituent must be associated with overland flow in order to receive atmospheric deposition inputs','')     #ERRMSG0
 
 # english system
 FACTA  = 1.0
@@ -38,7 +38,7 @@ def pqual(store, siminfo, uci, ts):
 		flags = uci['PQUAL' + pqual + '_FLAGS']
 		constituents.append(flags['QUALID'])
 
-	errorsV = zeros(len(ERRMSG), dtype=int)
+	errorsV = zeros(len(ERRMSGS), dtype=int)
 	delt60 = siminfo['delt'] / 60  # delt60 - simulation time interval in hours
 	simlen = siminfo['steps']
 	tindex = siminfo['tindex']
@@ -138,6 +138,9 @@ def pqual(store, siminfo, uci, ts):
 			elif pqadfgc == -1:
 				ts['PQADCN'] = ts['PQADCN' + str(index) + ' 1']
 
+			if QSOFG == 0 and (pqadfgf != 0 or pqadfgc != 0):
+				errorsV[0] += 1  # error - non-qualof cannot have atmospheric deposition
+
 		if 'PQADFX' not in ts:
 			ts['PQADFX'] = zeros(simlen)
 		if 'PQADCN' not in ts:
@@ -329,5 +332,5 @@ def pqual(store, siminfo, uci, ts):
 			PQADDR[loop] = adfxfx
 			PQADEP[loop] = adtot
 
-	return errorsV, ERRMSG
+	return errorsV, ERRMSGS
 

HSP2/PSTEMP.py

@@ -11,9 +11,9 @@
 from HSP2.utilities  import hoursval, initm, make_numba_dict
 
 
-ERRMSG = ['SLTMP temperature less than 100C',   # MSG0
-		'ULTMP temperature less than 100C',     # MSG1
-		'LGTMP temperature less than 100C', 	# MSG2
+ERRMSG = ['SLTMP temperature less than -100C',   # MSG0
+		'ULTMP temperature less than -100C',     # MSG1
+		'LGTMP temperature less than -100C', 	# MSG2
 		'SLTMP temperature greater than 100C', 	# MSG3
 		'ULTMP temperature greater than 100C',  # MSG4
 		'LGTMP temperature greater than 100C']  # MSG5

HSP2/PWATER.py

@@ -21,7 +21,7 @@
           'PWATER: Proute runoff did not converge',       #ERRMSG6
           'PWATER: UZI highly negative',                  #ERRMSG7
           'PWATER: Reset AGWS to zero',                   #ERRMSG8
-          'PWATER: High Water Table code not implimented', #ERRMSG9
+          'PWATER: High Water Table code not implemented', #ERRMSG9
           )
 
 def pwater(store, siminfo, uci, ts):

HSP2/SEDTRN.py

@@ -8,12 +8,18 @@
 from HSP2.ADCALC import advect
 from HSP2.utilities  import make_numba_dict
 
-ERRMSG = []
+ERRMSGS =('SEDTRN: Warning -- bed storage of sediment size fraction sand is empty',                                   #ERRMSG0
+          'SEDTRN: Warning -- bed storage of sediment size fraction silt is empty',                                   #ERRMSG1
+          'SEDTRN: Warning -- bed storage of sediment size fraction clay is empty',                                   #ERRMSG2
+          'SEDTRN: Warning -- bed depth appears excessive',                                                           #ERRMSG3
+          'SEDTRN: Fatal error ocurred in colby method- variable outside valid range- switching to toffaleti method', #ERRMSG4
+          'SEDTRN: Simulation of sediment requires all 3 "auxiliary flags" (AUX1FG, etc) in section HYDR must be turned on', #ERRMSG5
+          'SEDTRN: When specifying the initial composition of the bed, the fraction of sand, silt, and clay must sum to a value close to 1.0.')  #ERRMSG6
 
 def sedtrn(store, siminfo, uci, ts):
 	''' Simulate behavior of inorganic sediment'''
 
-	errorsV = zeros(len(ERRMSG), dtype=int)
+	errorsV = zeros(len(ERRMSGS), dtype=int)
 
 	simlen = siminfo['steps']
 	delt   = siminfo['delt']
@@ -31,6 +37,9 @@ def sedtrn(store, siminfo, uci, ts):
 	# table SANDFG
 	sandfg = ui['SANDFG']   # 1: Toffaleti method, 2:Colby method, 3:old HSPF power function
 
+	if ui['AUX3FG'] == 0:
+		errorsV[5] += 1  # error - sediment transport requires aux3fg to be on
+
 	# table SED-GENPARM
 	bedwid = ui['BEDWID']
 	bedwrn = ui['BEDWRN']
@@ -99,7 +108,7 @@ def sedtrn(store, siminfo, uci, ts):
 	clay_bedfr  = ui['CLAYFR']	
 	total_bedfr = sand_bedfr + silt_bedfr + clay_bedfr
 	if abs(total_bedfr - 1.0) > 0.01:
-		pass # error message: sum of bed sediment fractions is not close enough to 1.0
+		errorsV[6] += 1  # error message: sum of bed sediment fractions is not close enough to 1.0
 
 	# suspended sediment concentrations; table ssed-init
 	sand_ssed1  = ui['SSED1']
@@ -390,17 +399,20 @@ def sedtrn(store, siminfo, uci, ts):
 		osed4 += osed1
 		sand_rssed1 = sand_t_rsed7 = sand_rsed1 + sand_wt_rsed4  # total storage in mg.vol/l
 		if sand_wt_rsed4 == 0.0:       # warn that bed is empty
-			pass # errmsg
+			# errmsg
+			errorsV[0] += 1  # The bed storage of sediment size fraction sand is empty.
 
 		osed4 += osed2
 		silt_rssed2 = silt_t_rsed8 = silt_rsed2 + silt_wt_rsed5    # total storage in mg.vol/l
 		if silt_wt_rsed5 == 0.0:         # warn that bed is empty
-			pass # errmsg
+			# errmsg
+			errorsV[1] += 1  # The bed storage of sediment size fraction silt is empty.
 
 		osed4 += osed3
 		clay_rssed3 = clay_t_rsed9 = clay_rsed3 + clay_wt_rsed6  	# total storage in mg.vol/l
 		if clay_wt_rsed6 == 0.0:   # warn that bed is empty
-			pass # errmsg
+			# errmsg
+			errorsV[2] += 1  # The bed storage of sediment size fraction clay is empty.
 
 		# find the volume occupied by each fraction of bed sediment- ft3 or m3
 		volsed = (sand_wt_rsed4 / (sand_rho * 1.0e06)
@@ -418,7 +430,8 @@ def sedtrn(store, siminfo, uci, ts):
 		volsed = volsed / (1.0 - por)        #  allow for porosit
 		beddep = volsed / (len_ * bedwid)     # calculate thickness of bed- ft or m
 		if beddep > bedwrn:
-			pass  # Errormsg:  warn that bed depth appears excessive
+			# Errormsg:  warn that bed depth appears excessive
+			errorsV[3] += 1
 
 		svol = vol  # svol is volume at start of time step, update for next time thru
 
@@ -487,7 +500,7 @@ def sedtrn(store, siminfo, uci, ts):
 			ts['OSED3' + str(i + 1)] = OSED3[:, i]
 			ts['OSED4' + str(i + 1)] = OSED4[:, i]
 
-	return errorsV, ERRMSG
+	return errorsV, ERRMSGS
 
 
 def bdexch (avdepm, w, tau, taucd, taucs, m, vol, frcsed, susp, bed):
@@ -644,8 +657,8 @@ def colby(v, db50, fhrad, fsl, tempr):
 		if vx > v:
 			break
 	iv2 = iv1 + 1
-	yy1 = log10(VG[iv1])
-	yy2 = log10(VG[iv2])
+	yy1 = log10(VG[iv1-1])
+	yy2 = log10(VG[iv2-1])
 	yyratio = (log10(v) - yy1) / (yy2 - yy1)		
 
 	tmpr = min(100.0, max(32.0, tempr * 1.8 + 32.0))

HSP2/SNOW.py

@@ -102,7 +102,7 @@ def _snow_(ui, ts):
         icefg  = int(ui['ICEFG'])
 
     melev  = ui['MELEV']
-    packf  = ui['PACKF']     # inital df.PKSNOW += df.PKICE fixed in uciReader
+    packf  = ui['PACKF'] + ui['PACKI']
     packi  = ui['PACKI']
     packw  = ui['PACKW']
     paktmp = ui['PAKTMP']
@@ -535,7 +535,7 @@ def _snow_(ui, ts):
         PACKF[step]  = packf
         PACKI[step]  = packi
         PACKW[step]  = packw
-        PACK[step]   = packf + packi + packw
+        PACK[step]   = packf + packw # + packi
         PAKTMP[step] = paktmp
         PDEPTH[step] = pdepth
         PRAIN[step]  = prain
@@ -565,6 +565,7 @@ def vapor(SVP, temp):
 def hour6flag(siminfo, dofirst=False):
     '''timeseries with 1 at 6am and earlier each day, and zero otherwise'''
     hours24 = zeros(24)
+    hours24[0] = 1.0
     hours24[1] = 1.0
     hours24[2] = 1.0
     hours24[3] = 1.0

HSP2/main.py

@@ -88,10 +88,12 @@ def main(hdfname, saveall=False, jupyterlab=True):
                         ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData']
                         # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL']
                         ui['PARAMETERS']['HTFG'] = flags['HTRCH']
+                        ui['PARAMETERS']['AUX3FG'] = 0
                         if flags['HYDR']:
                             ui['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN']
                             ui['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH']
                             ui['PARAMETERS']['DB50'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DB50']
+                            ui['PARAMETERS']['AUX3FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX3FG']
                     if activity == 'GQUAL':
                         ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData']
                         ui['PARAMETERS']['HTFG'] = flags['HTRCH']