Merge OXRX, NUTRX, and PLANK class approach into LimnoTech develop

HSP2/GQUAL.py

@@ -643,26 +643,31 @@ def _gqual_(ui, ts):
 	if tempfg == 1:
 		if 'TW' not in ts:
 			errors[6] += 1  # ERRMSG6: timeseries not available
+			ts['TW_GQ'] = zeros(simlen)
 		else:
 			ts['TW_GQ'] = ts['TW']
 	if phflag == 1:
 		if 'PHVAL' not in ts:
 			errors[7] += 1  # ERRMSG7: timeseries not available
+			ts['PHVAL_GQ'] = zeros(simlen)
 		else:
 			ts['PHVAL_GQ'] = ts['PHVAL']
 	if roxfg == 1:
 		if 'ROC' not in ts:
 			errors[8] += 1  # ERRMSG8: timeseries not available
+			ts['ROC_GQ'] = zeros(simlen)		
 		else:
 			ts['ROC_GQ'] = ts['ROC']
 	if cldfg == 1:
 		if 'CLOUD' not in ts:
 			errors[9] += 1  # ERRMSG9: timeseries not available
+			ts['CLOUD_GQ'] = zeros(simlen)
 		else:
 			ts['CLOUD_GQ'] = ts['CLOUD']
 	if sdfg == 1:
 		if 'SSED4' not in ts:
 			errors[10] += 1  # ERRMSG10: timeseries not available
+			ts['SDCNC_GQ'] = zeros(simlen)
 		else:
 			ts['SDCNC_GQ'] = ts['SSED4']
 	if phytfg == 1:

HSP2/HTRCH.py

@@ -167,6 +167,8 @@ def _htrch_(ui, ts):
 		tgrnd = ui['TGRND']
 		kmud = ui['KMUD'] * delt60  # convert rate coefficients from kcal/m2/C/hr to kcal/m2/C/ivl
 		kgrnd = ui['KGRND'] * delt60
+	else:
+		muddep = 0.0
 
 	if uunits == 1:			# input units are deg.F, but need deg.C
 		muddep = muddep / 3.281

HSP2/NUTRX_Class.py

@@ -201,7 +201,7 @@ def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts, OXRX):
 		self.NUADFG = zeros(7, dtype=np.int32)
 
 		for j in range(1, 7):
-			self.NUADFG[j] = int(ui['NUADFG' + str(j)])
+			self.NUADFG[j] = int(ui['NUADFG(' + str(j) + ')'])
 
 		# error handling:
 		if self.TAMFG == 0 and (self.AMVFG == 1 or self.ADNHFG == 1):

HSP2/PHCARB_Class.py

@@ -0,0 +1,144 @@
+import numpy as np
+from numpy import zeros, array
+import numba as nb
+from numba.typed import Dict
+from numba.experimental import jitclass
+from math import exp, log
+
+from HSP2.ADCALC import advect
+from HSP2.RQUTIL import benth
+
+spec = [
+	('AFACT', nb.float64),
+	('alkcon', nb.int32),
+	('anaer', nb.int32),
+	('BENRFG', nb.int32),
+	('brco2', nb.float64[:]),
+	('cfcinv', nb.float64),
+	('co2', nb.float64),
+	('conv', nb.float64),
+	('delt60', nb.float64),
+	('delth', nb.float64),
+	('delts', nb.float64),
+	('errors', nb.int64[:]),
+	('ico2', nb.float64),
+	('itic', nb.float64),
+	('ncons', nb.int32),
+	('nexits', nb.int32),
+	('oco2', nb.float64[:]),
+	('otic', nb.float64[:]),
+	('ph', nb.float64),
+	('phcnt', nb.int32),
+	('roco2', nb.float64),
+	('rotic', nb.float64),
+	('uunits', nb.int32),
+	('svol', nb.float64),
+	('tic', nb.float64),
+	('vol', nb.float64),
+]
+
+@jitclass(spec)
+class PHCARB_Class:
+
+	#-------------------------------------------------------------------
+	# class initialization:
+	#-------------------------------------------------------------------
+	def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts):
+
+		''' Initialize variables for pH, carbon dioxide, total inorganic carbon,  and alkalinity '''
+
+		self.errors = zeros(int(ui['errlen']), dtype=np.int64)
+
+		delt60 = siminfo['delt'] / 60.0  # delt60 - simulation time interval in hours
+		self.delt60 = delt60
+		self.simlen = int(siminfo['steps'])
+		self.delts  = siminfo['delt'] * 60
+		self.uunits = int(siminfo['units'])
+
+		self.nexits = int(nexits)
+
+		self.vol = vol
+		self.svol = self.vol
+
+		# inflow/outflow conversion factor:
+		if self.uunits == 2:		# SI conversion: (g/m3)*(m3/ivld) --> [kg/ivld]
+			self.conv = 1.0e-3
+		else:						# Eng. conversion: (g/m3)*(ft3/ivld) --> [lb/ivld]
+			self.conv = 6.2428e-5
+
+		# required values from other modules:
+		self.BENRFG = int(ui_rq['BENRFG'])	# via table-type benth-flag
+		self.anaer = int(ui_rq['ANAER'])
+
+		# flags - table-type ph-parm1
+		self.phcnt  = int(ui['PHCNT'])    # is the maximum number of iterations to pH solution.
+		self.alkcon = int(ui['ALKCON'])   # ALKCON  is the number of the conservative substance which is alkalinity.
+
+		ncons = int(ui_rq['NCONS'])
+		if self.alkcon > ncons:
+			self.errors[0] += 1
+			# ERRMSG: Invalid CONS index specified for ALKCON (i.e., ALKCON > NCONS).
+
+		# flags - table-type ph-parm2
+		self.cfcinv = ui['CFCINV']
+
+		self.brco2 = zeros(2)
+		for i in range(2):
+			self.brco2[i] = ui['BRCO2' + str(i+1)] * self.delt60
+
+		# table-type ph-init
+		self.tic = ui['TIC']
+		self.co2 = ui['CO2']
+		self.ph   = ui['PH']
+
+		# initialize outflows:
+		self.roco2 = 0.0
+		self.rotic = 0.0
+
+		self.oco2 = zeros(nexits)
+		self.otic = zeros(nexits)
+
+		return
+
+	#-------------------------------------------------------------------
+	# simulation (single timestep):
+	#-------------------------------------------------------------------
+
+	def simulate(self, phif1, phif2, avdepe, tw, dox, advectData):
+
+		''' simulate ph, carbon dioxide, total inorganic carbon, and alkalinity'''
+
+		# hydraulics:
+		(nexits, vols, vol, srovol, erovol, sovol, eovol) = advectData
+
+		self.vol = vol
+
+		# inflows: convert from [mass/ivld] to [conc.*vol/ivld]
+		self.itic = phif1 / self.conv
+		self.ico2 = phif2 / self.conv
+
+		# advect TIC:
+		(self.tic, self.rotic, self.otic) = \
+			advect(self.itic, self.bod, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol)
+
+		# advect CO2:
+		(self.co2, self.roco2, self.oco2) = \
+			advect(self.ico2, self.co2, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol)
+
+		if self.CONSFG == 0:
+			pass
+		else:
+			# con(alkcon) is available from section cons
+			pass
+
+
+		if self.vol > 0.0:
+			pass
+
+		self.svol = self.vol  # svol is volume at start of time step, update for next time thru
+
+		return
+
+	#-------------------------------------------------------------------
+	# utility functions:
+	#-------------------------------------------------------------------

HSP2/PLANK_Class.py

@@ -321,7 +321,9 @@ def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts, OXRX, NUTRX):
 			self.numbal = self.BALFG          # single species or none
 
 		self.cfsaex = 1.0
-		if self.HTFG == 0 or 'CFSAEX' in ui:     # fraction of surface exposed - table-type surf-exposed
+		if self.HTFG > 0 and 'CFSAEX' in ui_rq:		# via heat-parm input table
+			self.cfsaex = ui_rq['CFSAEX']		
+		elif 'CFSAEX' in ui:     				# fraction of surface exposed - table-type surf-exposed
 			self.cfsaex = ui['CFSAEX']
 
 		# table-type plnk-parm1
@@ -348,10 +350,16 @@ def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts, OXRX, NUTRX):
 		self.cmmn   = ui['CMMN']
 		self.cmmnp  = ui['CMMNP']
 		self.cmmp   = ui['CMMP']
+
 		self.talgrh = ui['TALGRH']
 		self.talgrl = ui['TALGRL']
 		self.talgrm = ui['TALGRM']
 
+		if self.uunits == 1:
+			self.talgrh = (self.talgrh - 32.0) * 0.555		
+			self.talgrl = (self.talgrl - 32.0) * 0.555		
+			self.talgrm = (self.talgrm - 32.0) * 0.555		
+
 		# table-type plnk-parm3
 		self.alr20 = ui['ALR20'] * delt60   	# convert rates from 1/hr to 1/ivl
 		self.aldh  = ui['ALDH']  * delt60 
@@ -652,6 +660,9 @@ def simulate(self, tw, phval, co2, tss, OXRX, NUTRX, iphyto, izoo, iorn, iorp, i
 				else:                      # use full depth and velocity
 					self.balvel = avvele
 					self.baldep = avdepe
+			else:
+				self.balvel = 0.0
+				self.baldep = 0.0
 
 			# calculate solar radiation absorbed; solrad is the solar radiation at gage,
 			# corrected for location of reach; 0.97 accounts for surface reflection
@@ -690,7 +701,7 @@ def simulate(self, tw, phval, co2, tss, OXRX, NUTRX, iphyto, izoo, iorn, iorp, i
 			#	sestonic algae growth & respiration
 			#-----------------------------------------------------------
 			if self.PHYFG == 1:   # simulate phytoplankton, phyrx only called here
-				(po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.phyto,self.limphy,self.co2,self.phycla,
+				(po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.phyto,self.limphy,self.pyco2,self.phycla,
 				dophy,bodphy,tamphy,no3phy,po4phy,phdth,phgro,ornphy,orpphy,orcphy) \
 					= self.phyrx(self.phylit,tw,self.talgrl,self.talgrh,self.talgrm,self.malgr,self.cmmp, \
 							self.cmmnp,NUTRX.TAMFG,self.AMRFG,self.NSFG,self.cmmn,self.cmmlt,self.delt60, \
@@ -716,7 +727,7 @@ def simulate(self, tw, phval, co2, tss, OXRX, NUTRX, iphyto, izoo, iorn, iorp, i
 				#	zooplankton growth & death:
 				#-----------------------------------------------------------
 				if self.ZOOFG == 1:    # simulate zooplankton, zorx only called here
-					(dox,bod,self.zoo,self.orn,self.orp,self.orc,tam,no3,po4,zeat,zco2,dozoo,bodzoo,nitzoo,po4zoo,zgro,zdth,zorn,zorp,zorc) \
+					(dox,bod,self.zoo,self.orn,self.orp,self.orc,tam,no3,po4,zeat,self.zoco2,dozoo,bodzoo,nitzoo,po4zoo,zgro,zdth,zorn,zorp,zorc) \
 						= self.zorx(self.zfil20,self.tczfil,tw,self.phyto,self.mzoeat,self.zexdel,self.cvpb, \
 							self.zres20,self.tczres,NUTRX.anaer,self.zomass,NUTRX.TAMFG,self.refr, \
 							self.ZFOOD,self.zd,self.oxzd,self.cvbn,self.cvbp,self.cvbc,self.cvnrbo,self.cvbo, \
@@ -1253,6 +1264,11 @@ def litrch(inlit, extb, extcla, extsed, avdepe, baldep, PHYFG, BALFG):
 		''' calculate light correction factor to algal growth (cflit); 
 		determine amount of light available to phytoplankton and benthic algae'''
 		ln01 = 4.60517   # minus natural log 0.01
+
+		cflit  = 0.0
+		phylit = 0.0
+		ballit = 0.0
+
 		if inlit > 0.0:
 			# calculate extinction of light based on the base extinction
 			# coefficient of the water incremented by self-shading effects
@@ -1290,10 +1306,10 @@ def litrch(inlit, extb, extcla, extsed, avdepe, baldep, PHYFG, BALFG):
 				ballit = inlit * exp(-extco * baldep)
 				if ballit < 0.0001:
 					ballit = 0.0
+
 		else:   # there is no incident solar radiation; algal growth cannot occur
-			cflit  = 0.0
-			phylit = 0.0
-			ballit = 0.0
+			pass
+
 		return phylit, ballit, cflit
 
 	@staticmethod
@@ -1629,7 +1645,7 @@ def pksums(self, NUTRX, phyto, zoo, orn, orp, orc, no3, tam, no2, snh4, po4, spo
 		tp     = -1.0e30
 
 		if (self.vol <= 0):
-			return
+			return torn, torp, torc, potbod, tn, tp
 
 		# Calculate sums:
 		tval = 0.0