SEDMNT.py and SOLIDS.py -- refinements and testing

HSP2/SEDMNT.py

@@ -4,9 +4,9 @@
 
 Conversion of HSPF HPERSED.FOR module into Python''' 
 
-from numpy import zeros, where
-from numba import jit
-from HSP2  import initm
+from numpy import zeros, where, int64
+from numba import njit
+from HSP2.utilities  import initm, make_numba_dict, hourflag
 
 ERRMSG = []
 
@@ -18,39 +18,42 @@
 PDETS= DETS*MFACTA # convert dimensional variables to external units
 '''
 
-ERRMSG = []
-
-def sedmnt(store, general, ui, ts):
+def sedmnt(store, siminfo, uci, ts):
 	''' Produce and remove sediment from the land surface'''
 
 	errorsV = zeros(len(ERRMSG), dtype=int)
-	delt   = general['sim_delt'] 
-	simlen = general['sim_len']
-	tindex = general['tindex']
-
+
+	delt   = siminfo['delt']
+	simlen = siminfo['steps']
+	tindex = siminfo['tindex']
+
+	ui = make_numba_dict(uci)  # Note: all values converted to float automatically
 	VSIVFG = ui['VSIVFG']
 	SDOPFG = ui['SDOPFG']
+	CSNOFG = int(ui['CSNOFG'])
 	smpf   = ui['SMPF']
 	krer   = ui['KRER']
 	jrer   = ui['JRER']
 	affix  = ui['AFFIX']
-	nvsi   = ui['NVSI'] * DELT / 1440.
+	nvsi   = ui['NVSI'] * delt / 1440.
 	kser   = ui['KSER']
 	jser   = ui['JSER']
 	kger   = ui['KGER']
 	jger   = ui['JGER']
-	nvsim  = ui['NVSIM']
-
-	initm(general, ui, ts, 'CRVFG', 'COVERM', 'COVER')
-	COVER = ts['COVER']
-	cover = COVER[0]   # for numba
 
-	if RAINF in ts:
+	u = uci['PARAMETERS']
+	ts['COVERI'] = initm(siminfo, uci, u['CRVFG'], 'MONTHLY_COVER', u['COVER'])
+	COVERI = ts['COVERI']
+	cover = COVERI[0]   # for numba
+	ts['NVSI'] = initm(siminfo, uci, u['VSIVFG'], 'MONTHLY_NVSI', u['NVSI'])
+	NVSI = ts['NVSI'] * delt / 1440.
+
+	if 'RAINF' in ts:
 		RAIN = ts['RAINF']
 	else:
 		RAIN = ts['PREC']
 
-	for name in ['PREC', 'SLSED', 'SNOCOV', 'SURO', 'SURS',]:
+	for name in ['PREC', 'SLSED', 'SNOCOV', 'SURO', 'SURS']:
 		if name not in ts:
 			ts[name] = zeros(simlen)
 	PREC   = ts['PREC']
@@ -60,48 +63,54 @@ def sedmnt(store, general, ui, ts):
 	SURS   = ts['SURS']
 
 	# preallocate output arrays
-	DETS  = ts['DETS']  = zeros(sim_len)
-	WSSD  = ts['WSSD']  = zeros(sim_len)
-	SCRSD = ts['SCRSD'] = zeros(sim_len)
-	SOSED = ts['SOSED'] = zeros(sim_len)
+	DETS  = ts['DETS']  = zeros(simlen)
+	WSSD  = ts['WSSD']  = zeros(simlen)
+	SCRSD = ts['SCRSD'] = zeros(simlen)
+	SOSED = ts['SOSED'] = zeros(simlen)
+	COVER = ts['COVER'] = zeros(simlen)
+
+        # HSPF 12.5 has only one sediment block
+	dets = ui['DETS']
 
 	# BLOCK SPECIFIC VALUES
-	nblks = ui['NBLKS']
+	# nblks = ui['NBLKS']
 
-	''' the initial storages are repeated for each block'''
-	sursb = ui['SURSB']
-	uzsb  = ui['UZSB']
-	ifwsb = ui['IFWSB']
-	detsb  = ui['DETSB']
+	# ''' the initial storages are repeated for each block'''
+	# sursb = ui['SURSB']
+	# uzsb  = ui['UZSB']
+	# ifwsb = ui['IFWSB']
+	# detsb  = ui['DETSB']
 
-	if nblks > 1:
-		SUROB  = ts['SIROB']  # if NBLKS > 1
-		SURB   = ts['SURSB']  # if NBLKS > 1
-
-		DETSB  = ts['DETSB']  = zeros(NBLKS, sim_len)
-		WSSDB  = ts['WSSDB']  = zeros(NBLKS, sim_len)
-		SCRSDB = ts['SCRSDB'] = zeros(NBLKS, sim_len)
-		SOSDB  = ts['SOSDB']  = zeros(NBLKS, sim_len)
-
+	# if nblks > 1:
+	#	SUROB  = ts['SIROB']  # if NBLKS > 1
+	#	SURB   = ts['SURSB']  # if NBLKS > 1
+
+	#	DETSB  = ts['DETSB']  = zeros(nblks, simlen)
+	#	WSSDB  = ts['WSSDB']  = zeros(nblks, simlen)
+	#	SCRSDB = ts['SCRSDB'] = zeros(nblks, simlen)
+	#	SOSDB  = ts['SOSDB']  = zeros(nblks, simlen)
 
 	# initialize sediment transport capacity
+	surs = SURS[0]
+	delt60 = delt / 60.0  # simulation interval in hours
 	stcap = delt60 * kser * (surs / delt60)**jser if SDOPFG else 0.0
 
-	DAYFG = where(tindex.hour==1, True, False)   # ??? need to check if minute == 0	
-	DAYFG[0] = True
+	# DAYFG = where(tindex.hour==1, True, False)   # ??? need to check if minute == 0
+	DAYFG = hourflag(siminfo, 0, dofirst=True).astype(bool)
 
 	DRYDFG = 1
 
-	for loop in range(sim_len):
+	for loop in range(simlen):
 		dayfg  = DAYFG[loop]
 		rain   = RAIN[loop]
 		prec   = PREC[loop]
 		suro   = SURO[loop]
 		surs   = SURS[loop]
 		snocov = SNOCOV[loop]
+		slsed  = SLSED[loop]
 
 		if dayfg:
-			cover = COVER[loop]
+			cover = COVERI[loop]
 
 		# estimate the quantity of sediment particles detached from the soil surface by rainfall and augment the detached sediment storage
 		# detach()
@@ -119,8 +128,8 @@ def sedmnt(store, general, ui, ts):
 			det = 0.0
 		# end detach()
 
-		if DAYFG:     # it is the first interval of the day
-			nvsi = ???  # update nvsi
+		if dayfg:     # it is the first interval of the day
+			nvsi = NVSI[loop]
 			if VSIVFG == 2 and DRYDFG:  # last day was dry, add a whole days load in first interval detailed output will show load added over the whole day.;
 				dummy = nvsi * (1440. / delt)
 				dets  = dets + dummy
@@ -200,5 +209,6 @@ def sedmnt(store, general, ui, ts):
 		WSSD[loop]  = wssd
 		SCRSD[loop] = scrsd
 		SOSED[loop] = sosed
+		COVER[loop] = cover
 
 	return errorsV, ERRMSG

HSP2/SOLIDS.py

@@ -5,23 +5,24 @@
 Conversion of HSPF HIMPSLD.FOR module into Python''' 
 
 from numpy import zeros, where
-from numba import jit
-from HSP2  import initm
+from numba import njit
+from HSP2.utilities  import initm, make_numba_dict, hourflag
 
 
 MFACTA = 1.0  # english units
 
 ERRMSG = []
 
 
-def solids(store, general, ui, ts):
+def solids(store, siminfo, uci, ts):
 	'''Accumulate and remove solids from the impervious land segment'''
 
 	errorsV = zeros(len(ERRMSG), dtype=int)
-	delt60 = general['sim_delt'] / 60     # delt60 - simulation time interval in hours
-	simlen = general['sim_len']
-	tindex = general['tindex']
+	delt60 = siminfo['delt'] / 60     # delt60 - simulation time interval in hours
+	simlen = siminfo['steps']
+	tindex = siminfo['tindex']
 
+	ui = make_numba_dict(uci)  # Note: all values converted to float automatically
 	SDOPFG = ui['SDOPFG']
 	keim   = ui['KEIM']
 	jeim   = ui['JEIM']
@@ -40,15 +41,17 @@ def solids(store, general, ui, ts):
 	SLDS  = ts['SLDS']  = zeros(simlen)
 
 	drydfg = 1  # assume day is dry
-	DAYFG = where(tindex.hour==1, True, False)   # ??? need to check if minute == 0
-
+	DAYFG = hourflag(siminfo, 0, dofirst=True).astype(bool)
+
+	u = uci['PARAMETERS']
 	# process optional monthly arrays to return interpolated data or constant array
-	initm(general, ui, ts, 'VASDFG', 'ACCSDM', 'ACCSDP')
-	initm(general, ui, ts, 'VRSDFG', 'REMSDM', 'REMSDP')
+	ts['ACCSDP'] = initm(siminfo, uci, u['VASDFG'], 'ACCSDM', u['ACCSDP'])
+	ts['REMSDP'] = initm(siminfo, uci, u['VRSDFG'], 'REMSDM', u['REMSDP'])
 	ACCSDP = ts['ACCSDP']
 	REMSDP = ts['REMSDP']
-
-	for loop in range(sim_length):
+
+	slds = ui['SLDS']
+	for loop in range(simlen):
 		suro   = SURO[loop]
 		surs   = SURS[loop]
 		prec   = PREC[loop]
@@ -91,9 +94,12 @@ def solids(store, general, ui, ts):
 			if drydfg == 1:       # precipitation did not occur during the previous day
 				# update storage due to accumulation and removal which  occurs independent of runoff - units are lbs/acre
 				slds = accsdp + slds * (1.0 - remsdp)
-
+
+			if prec > 0.0:
 				# there is precipitation on the first interval of the new day
-				drydfg = 0 if prec > 0.0 else drydfg = 1
+				drydfg = 0
+			else:
+				drydfg = 1
 
 		SOSLD[loop] = sosld  # * MFACTA
 		SLDS[loop]  = slds   # * MFACTA

HSP2/configuration.py

@@ -10,7 +10,9 @@
 from HSP2.ATEMP  import atemp
 from HSP2.SNOW   import snow
 from HSP2.PWATER import pwater
+from HSP2.SEDMNT import sedmnt
 from HSP2.IWATER import iwater
+from HSP2.SOLIDS import solids
 from HSP2.HYDR   import hydr
 
 def noop (store, siminfo, ui, ts):
@@ -20,10 +22,10 @@ def noop (store, siminfo, ui, ts):
 
 # Note: This is the ONLY place in HSP2 that defines activity execution order
 activities = {
-  'PERLND': {'ATEMP':atemp, 'SNOW':snow, 'PWATER':pwater, 'SEDMNT':noop,
+  'PERLND': {'ATEMP':atemp, 'SNOW':snow, 'PWATER':pwater, 'SEDMNT':sedmnt,
      'PSTEMP':noop, 'PWTGAS':noop, 'PQUAL':noop, 'MSTLAY':noop, 'PEST':noop,
      'NITR':noop, 'PHOS':noop, 'TRACER':noop},
-  'IMPLND': {'ATEMP':atemp, 'SNOW':snow, 'IWATER':iwater, 'SOLIDS':noop,
+  'IMPLND': {'ATEMP':atemp, 'SNOW':snow, 'IWATER':iwater, 'SOLIDS':solids,
      'IWTGAS':noop, 'IQUAL':noop},
   'RCHRES': {'HYDR':hydr, 'ADCALC':noop, 'CONS':noop, 'HTRCH':noop,
      'SEDTRN':noop, 'GQUAL':noop, 'OXRX':noop, 'NUTRX':noop, 'PLANK':noop,

HSP2/main.py

@@ -58,6 +58,9 @@ def main(hdfname, saveall=False):
                 if operation == 'RCHRES':
                     get_flows(store,ts,activity,segment,ddlinks,ddmasslinks,siminfo['steps'], msg)
                 ui = uci[(operation, activity, segment)]   # ui is a dictionary
+                if operation == 'PERLND' and activity == 'SEDMNT':
+                    # special exception here to make CSNOFG available
+                    ui['PARAMETERS']['CSNOFG'] = uci[(operation, 'PWATER', segment)]['PARAMETERS']['CSNOFG']
 
                 ############ calls activity function like snow() ##############
                 errors, errmessages = function(store, siminfo, ui, ts)