diff --git a/HSP2/ADCALC.py b/HSP2/ADCALC.py index 96bdc7a0..75c7c376 100644 --- a/HSP2/ADCALC.py +++ b/HSP2/ADCALC.py @@ -3,7 +3,7 @@ License: LGPL2 ''' -from numpy import zeros +from numpy import zeros, array from numba import njit from HSP2.utilities import make_numba_dict @@ -29,6 +29,7 @@ def adcalc(store, siminfo, uci, ts): nexits = int(ui['NEXITS']) # table type GEN-INFO ui['simlen'] = siminfo['steps'] ui['delts'] = siminfo['delt'] * 60.0 # delts is the simulation interval in seconds + ui['uunits'] = siminfo['units'] # calculated timeseries for advect() if 'SROVOL' not in ts: @@ -67,13 +68,23 @@ def _adcalc_(ui, ts): simlen = int(ui['simlen']) nexits = int(ui['NEXITS']) delts = ui['delts'] + uunits = ui['uunits'] + + # units conversion constants, 1 ACRE is 43560 sq ft. assumes input in acre-ft + VFACT = 43560.0 + AFACT = 43560.0 + if uunits == 2: + # si units conversion constants, 1 hectare is 10000 sq m, assumes area input in hectares, vol in Mm3 + VFACT = 1.0e6 + AFACT = 10000.0 + # table ADCALC-DATA if 'CRRAT' in ui: crrat = ui['CRRAT'] else: crrat = 1.5 if 'VOL' in ui: - vol = ui['VOL'] + vol = ui['VOL'] * VFACT else: vol = 0.0 @@ -98,6 +109,8 @@ def _adcalc_(ui, ts): for index in range(nexits): ts['EOVOL' + str(index + 1)] = EOVOL[:, index] + ROS = ui['ROS'] + # external time series O = zeros((simlen, nexits)) if nexits > 1: @@ -107,10 +120,10 @@ def _adcalc_(ui, ts): O[:, 0] = ts['RO'] for loop in range(simlen): - vols = VOL[loop-1] * 43560 if loop > 0 else vol + vols = VOL[loop-1] * VFACT if loop > 0 else vol o = O[loop] - os = O[loop-1] if loop > 0 else O[loop] + os = O[loop-1] if loop > 0 else array([ROS]) ro = 0.0 ros= 0.0 for index in range(nexits): diff --git a/HSP2/COPY.py b/HSP2/COPY.py new file mode 100644 index 00000000..58be5cae --- /dev/null +++ b/HSP2/COPY.py @@ -0,0 +1,48 @@ +from HSP2.utilities import get_timeseries +import pandas as pd +from typing import List, Dict + +class Copy(): + """ + Partial implementation of the COPY module. + In HSPF, COPY module supports ability able to 'copy' in this case output timeseries to + locations specified in NETWORK block or to EXTERNAL SOURCES. + This functionality is not currently implemented, presently only loading from EXT SOURCES + """ + + _ts = {} + _ts['MEAN'] = {} + _ts['POINT'] = {} + + def __init__(self, store: pd.HDFStore, sim_info: Dict, ext_sources: List) -> None: + ts = get_timeseries(store, ext_sources, sim_info) + for source in ext_sources: + themn = source.TMEMN + themsb = source.TMEMSB + self.set_ts(ts[f'{themn}{themsb}'], themn, themsb) + + def set_ts(self, ts: pd.Series, themn: str, themsb: str) -> None: + """Set the provided timeseries to ts dictionary + + ts: pd.Series + pandas Series class instance corresponding to a timeseries + tmemn: str, {'MEAN', 'POINT'} + Target member name, specifies if target timeseries is in mean-valued + or point-valued dictionaries + tmemsb: str, + Target member name subscripts, acts as key for mean-valued and point-valued dictionaries + Original HSPF restricts this to 0-20 but no limit enforced in HSP2 + """ + self._ts[themn][themsb] = ts + + def get_ts(self, tmemn: str, tmemsb: str) -> pd.Series: + """Gets the specified timeseries from copy class instance based + + tmemn: str, {'MEAN', 'POINT'} + Target member name, specifies if target timeseries is in mean-valued + or point-valued dictionaries + tmemsb: str, + Target member name subscripts, acts as key for mean-valued and point-valued dictionaries + Original HSPF restricts this to 0-20 but no limit enforced in HSP2 + """ + return self._ts[tmemn][tmemsb] diff --git a/HSP2/GENER.py b/HSP2/GENER.py new file mode 100644 index 00000000..12a10bdb --- /dev/null +++ b/HSP2/GENER.py @@ -0,0 +1,147 @@ +from numba import njit +import numpy as np +import pandas as pd +from typing import Dict, List + +from HSP2.utilities import get_timeseries + +class Gener(): + """ + Partial implementation of the GENER module. + Currently supports OPCODES 1-7, 9-23, & 25-26 + """ + + ts_input_1 = pd.Series() # type: pd.Series + ts_input_2 = pd.Series() # type: pd.Series + ts_output = pd.Series() # type: pd.Series + opcode = 0 # type: int + k = -1.0E30 # type: float + + def __init__(self, segment: str, copies: Dict, geners: Dict, ddlinks: Dict, ddgener: Dict) -> None: + self.opcode = ddgener['OPCODE'][segment] + if self.opcode in [9,10,11,24,25,26]: + self.k = ddgener['PARM'][segment] + + for link in ddlinks[segment]: + ts = pd.Series() + if link.SVOL == 'COPY': + copy = copies[link.SVOLNO] + ts = copy.get_ts(link.SMEMN,link.SMEMSB1) + elif link.SVOL == 'GENER': + gener = geners[link.SVOLNO] + ts = gener.get_ts() + else: + raise NotImplementedError(f"Invalid SVOL. GENER module does not currently support reading TimeSeries for '{link.SVOL}'") + + if link.TGRPN == 'INPUT' and link.TMEMN == 'ONE': + self.ts_input_1 = ts + elif link.TGRPN == 'INPUT' and link.TMEMN == 'TWO': + self.ts_input_2 = ts + else: + raise AttributeError(f"No attribute {link.TGRPN}{link.THEMN} to assign TimeSeries. Should be either 'INPUTONE' or 'INPUTWO'") + + self._execute_gener() + + def get_ts(self) -> pd.Series: + """ + Returns the result TimeSeries generated from executing the operation specified by the OPCODE. + """ + return self.ts_output + + def _execute_gener(self) -> None: + gener_op = getattr(self, f'_opcode{self.opcode}') + ts_result = gener_op() + #May need additional logic here to set default of 1.0E30 to be consistent with FORTRAN code + self.ts_output = ts_result + + def _opcode1(self) -> pd.Series: + return np.abs(self.ts_input_1) + + def _opcode2(self) -> pd.Series: + return np.sqrt(self.ts_input_1) + + def _opcode3(self) -> pd.Series: + return np.trunc(self.ts_input_1) + + def _opcode4(self) -> pd.Series: + return np.ceil(self.ts_input_1) + + def _opcode5(self) -> pd.Series: + return np.floor(self.ts_input_1) + + def _opcode6(self) -> pd.Series: + return np.log(self.ts_input_1) + + def _opcode7(self) -> pd.Series: + return np.log10(self.ts_input_1) + + def _opcode8(self) -> pd.Series: + #Not presently implemented, read UCI would need to modify to + #process NTERMS and COEFFS sub blocks of GENER block + raise NotImplementedError("GENER OPCODE 8 is not currently supported") + + def _opcode9(self) -> pd.Series: + return np.power(self.k , self.ts_input_1) + + def _opcode10(self) -> pd.Series: + return np.power(self.ts_input_1, self.k) + + def _opcode11(self) -> pd.Series: + return np.add(self.ts_input_1, self.k) + + def _opcode12(self) -> pd.Series: + return np.sin(self.ts_input_1) + + def _opcode13(self) -> pd.Series: + return np.cos(self.ts_input_1) + + def _opcode14(self) -> pd.Series: + return np.tan(self.ts_input_1) + + def _opcode15(self) -> pd.Series: + return np.cumsum(self.ts_input_1) + + def _opcode16(self) -> pd.Series: + return np.add(self.ts_input_1, self.ts_input_2) + + def _opcode17(self) -> pd.Series: + return np.subtract(self.ts_input_1, self.ts_input_2) + + def _opcode18(self) -> pd.Series: + return np.multiply(self.ts_input_1, self.ts_input_2) + + def _opcode19(self) -> pd.Series: + return np.divide(self.ts_input_1, self.ts_input_2) + + def _opcode20(self) -> pd.Series: + return np.maximum(self.ts_input_1, self.ts_input_2) + + def _opcode21(self) -> pd.Series: + return np.minimum(self.ts_input_1, self.ts_input_2) + + def _opcode22(self) -> pd.Series: + return np.power(self.ts_input_1, self.ts_input_2) + + def _opcode23(self) -> pd.Series: + ts_out = pd.Series(index=self.ts_input_1.index) + s = 0 + for idx in self.ts_input_1.index: + result = self.ts_input_1[idx] - self.ts_input_2[idx] - s + if result < 0: + s = s - self.ts_input_1[idx] + self.ts_input_2[idx] + result = 0 + else: + s = 0 + ts_out[idx] = result + return ts_out + + def _opcode24(self) -> pd.Series: + #skip for now + #would need to figure out timeseries length component + raise NotImplementedError("GENER OPCODE 24 is not currently supported") + + def _opcode25(self) -> pd.Series: + return np.maximum(self.ts_input_1, self.k) + + def _opcode26(self) -> pd.Series: + return np.minimum(self.ts_input_1, self.k) diff --git a/HSP2/GQUAL.py b/HSP2/GQUAL.py index 48c0bf85..67632fd1 100644 --- a/HSP2/GQUAL.py +++ b/HSP2/GQUAL.py @@ -239,7 +239,7 @@ def gqual(store, siminfo, uci, ts): return errors, ERRMSGS -@njit(cache=True) +#@njit(cache=True) def _gqual_(ui, ts): ''' Simulate the behavior of a generalized quality constituent''' errors = zeros(int(ui['errlen'])).astype(int64) @@ -668,6 +668,7 @@ def _gqual_(ui, ts): if phytfg == 1: if 'PHYTO' not in ts: errors[11] += 1 # ERRMSG11: timeseries not available + ts['PHYTO_GQ'] = zeros(simlen) else: ts['PHYTO_GQ'] = ts['PHYTO'] @@ -713,35 +714,28 @@ def _gqual_(ui, ts): ts['GQUAL' + str(index) + '_ISQAL3'] = zeros(simlen) ISQAL3 = ts['GQUAL' + str(index) + '_ISQAL3'] - DEPSCR1 = ts['DEPSCR1'] - DEPSCR2 = ts['DEPSCR2'] - DEPSCR3 = ts['DEPSCR3'] - ROSED1 = ts['ROSED1'] - ROSED2 = ts['ROSED2'] - ROSED3 = ts['ROSED3'] - - # OSED1 = zeros((simlen, nexits)) - # OSED2 = zeros((simlen, nexits)) - # OSED3 = zeros((simlen, nexits)) - OSED1 = [] - OSED2 = [] - OSED3 = [] - for timeindex in range(simlen): - tarray1 = [] - tarray2 = [] - tarray3 = [] - if nexits > 1: - for xindex in range(nexits): - tarray1.append(ts['OSED1' + str(xindex + 1)][timeindex]) - tarray2.append(ts['OSED2' + str(xindex + 1)][timeindex]) - tarray3.append(ts['OSED3' + str(xindex + 1)][timeindex]) - else: - tarray1.append(ts['ROSED1'][timeindex]) - tarray2.append(ts['ROSED2'][timeindex]) - tarray3.append(ts['ROSED3'][timeindex]) - OSED1.append(tarray1) - OSED2.append(tarray2) - OSED3.append(tarray3) + if qalfg[7] == 1: # this constituent is associated with sediment + DEPSCR1 = ts['DEPSCR1'] + DEPSCR2 = ts['DEPSCR2'] + DEPSCR3 = ts['DEPSCR3'] + ROSED1 = ts['ROSED1'] + ROSED2 = ts['ROSED2'] + ROSED3 = ts['ROSED3'] + + OSED1 = zeros((simlen, nexits)) + OSED2 = zeros((simlen, nexits)) + OSED3 = zeros((simlen, nexits)) + + for timeindex in range(simlen): + if nexits > 1: + for xindex in range(nexits): + OSED1[timeindex, xindex] = ts['OSED1' + str(xindex + 1)][timeindex] + OSED2[timeindex, xindex] = ts['OSED2' + str(xindex + 1)][timeindex] + OSED3[timeindex, xindex] = ts['OSED3' + str(xindex + 1)][timeindex] + else: + OSED1[timeindex, 0] = ts['ROSED1'][timeindex] + OSED2[timeindex, 0] = ts['ROSED2'][timeindex] + OSED3[timeindex, 0] = ts['ROSED3'][timeindex] # this number is used to adjust reaction rates for temperature # TW20 = TW - 20.0 @@ -838,47 +832,40 @@ def _gqual_(ui, ts): vol = VOL[loop] * AFACT toqal = TOQAL[loop] tosqal = TOSQAL[loop] - if uunits == 1: - depscr1 = DEPSCR1[loop] / 3.121E-08 - depscr2 = DEPSCR2[loop] / 3.121E-08 - depscr3 = DEPSCR3[loop] / 3.121E-08 - rosed1 = ROSED1[loop] / 3.121E-08 - rosed2 = ROSED2[loop] / 3.121E-08 - rosed3 = ROSED3[loop] / 3.121E-08 - osed1 = array(OSED1[loop]) - osed2 = array(OSED2[loop]) - osed3 = array(OSED3[loop]) - osed1 = osed1 / 3.121E-08 - osed2 = osed2 / 3.121E-08 - osed3 = osed3 / 3.121E-08 - # osed1 = [x / 3.121E-08 for x in osed1] - # osed2 = [x / 3.121E-08 for x in osed2] - # osed3 = [x / 3.121E-08 for x in osed3] - rsed[1] = RSED1[loop] / 3.121E-08 - rsed[2] = RSED2[loop] / 3.121E-08 - rsed[3] = RSED3[loop] / 3.121E-08 - rsed[4] = RSED4[loop] / 3.121E-08 - rsed[5] = RSED5[loop] / 3.121E-08 - rsed[6] = RSED6[loop] / 3.121E-08 - else: - depscr1 = DEPSCR1[loop] / 1E-06 - depscr2 = DEPSCR2[loop] / 1E-06 - depscr3 = DEPSCR3[loop] / 1E-06 # 2.83E-08 - rosed1 = ROSED1[loop] / 1E-06 - rosed2 = ROSED2[loop] / 1E-06 - rosed3 = ROSED3[loop] / 1E-06 - osed1 = array(OSED1[loop]) - osed2 = array(OSED2[loop]) - osed3 = array(OSED3[loop]) - osed1 = osed1 / 1E-06 - osed2 = osed2 / 1E-06 - osed3 = osed3 / 1E-06 - rsed[1] = RSED1[loop] / 1E-06 - rsed[2] = RSED2[loop] / 1E-06 - rsed[3] = RSED3[loop] / 1E-06 - rsed[4] = RSED4[loop] / 1E-06 - rsed[5] = RSED5[loop] / 1E-06 - rsed[6] = RSED6[loop] / 1E-06 + + # initialize sediment-related variables: + if qalfg[7] == 1: # constituent is sediment-associated + osed1 = zeros(nexits) + osed2 = zeros(nexits) + osed3 = zeros(nexits) + + # define conversion factor: + cf = 1.0 + if uunits == 1: + cf = 3.121e-8 + else: + cf = 1.0e-6 + + depscr1 = DEPSCR1[loop] / cf + depscr2 = DEPSCR2[loop] / cf + depscr3 = DEPSCR3[loop] / cf + rosed1 = ROSED1[loop] / cf + rosed2 = ROSED2[loop] / cf + rosed3 = ROSED3[loop] / cf + + for i in range(nexits): + osed1[i] = OSED1[loop, i] / cf + osed2[i] = OSED2[loop, i] / cf + osed3[i] = OSED3[loop, i] / cf + + rsed = zeros(7) + rsed[1] = RSED1[loop] / cf + rsed[2] = RSED2[loop] / cf + rsed[3] = RSED3[loop] / cf + rsed[4] = RSED4[loop] / cf + rsed[5] = RSED5[loop] / cf + rsed[6] = RSED6[loop] / cf + isqal1 = ISQAL1[loop] isqal2 = ISQAL2[loop] isqal3 = ISQAL3[loop] diff --git a/HSP2/HTRCH.py b/HSP2/HTRCH.py index fb49c86c..eced0105 100644 --- a/HSP2/HTRCH.py +++ b/HSP2/HTRCH.py @@ -75,12 +75,12 @@ def htrch(store, siminfo, uci, ts): if 'DELH' in ui: delh = ui['DELH'] else: - delh = zeros(tstop) + delh = zeros(int(tstop)) ts['DELH'] = delh if 'DELTT' in ui: deltt = ui['DELTT'] else: - deltt = zeros(tstop) + deltt = zeros(int(tstop)) ts['DELTT'] = deltt u = uci['PARAMETERS'] @@ -168,6 +168,11 @@ def _htrch_(ui, ts): kmud = ui['KMUD'] * delt60 # convert rate coefficients from kcal/m2/C/hr to kcal/m2/C/ivl kgrnd = ui['KGRND'] * delt60 + if uunits == 1: # input units are deg.F, but need deg.C + muddep = muddep / 3.281 + tgrnd = (tgrnd - 32.0) * 0.555 + tstop = (tstop - 32.0) * 0.555 + # if bedflg == 3: delh = ts['DELH'] deltt = ts['DELTT'] @@ -349,6 +354,9 @@ def _htrch_(ui, ts): # print('bedflg') if bedflg == 1: # one-layer bed conduction model tgrnd = TGRND[loop] + if uunits == 1: + tgrnd = (tgrnd - 32.0) * 0.555 + qbed = kmud * (tgrnd - tw) elif bedflg == 2: # two-layer bed conduction model # Following is subrouting #$BEDHT2 diff --git a/HSP2/HYDR.py b/HSP2/HYDR.py index 2ae78e38..54190543 100644 --- a/HSP2/HYDR.py +++ b/HSP2/HYDR.py @@ -124,6 +124,9 @@ def hydr(store, siminfo, uci, ts): if 'O' in ts: del ts['O'] if 'OVOL' in ts: del ts['OVOL'] + # save initial outflow(s) from reach: + uci['PARAMETERS']['ROS'] = ui['ROS'] + return errors, ERRMSGS @@ -249,6 +252,9 @@ def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels): #rirsht = 0.0 irrdem = 0.0 + # store initial outflow from reach: + ui['ROS'] = ro + # HYDR (except where noted) for step in range(steps): convf = CONVF[step] diff --git a/HSP2/NUTRX.py b/HSP2/NUTRX.py index 19017225..6056da8e 100644 --- a/HSP2/NUTRX.py +++ b/HSP2/NUTRX.py @@ -55,22 +55,20 @@ def nutrx(store, siminfo, uci, ts): # ERRMSG: error - sediment associated nh4 and/or po4 is being simulated,but sediment is not being simulated in section sedtrn uafxm = zeros((13,4)) + sf = 1.0 + + if uunits == 1: # convert from lb/ac.day to mg.ft3/l.ft2.ivl + sf = 0.3677 * delt60 / 24.0 + else: # convert from kg/ha.day to mg.m3/l.m2.ivl + sf = 0.1 * delt60 / 24.0 + + # convert units to internal: if NUADFG[1] > 0: - uafxm[:,1] = ui['NUAFXM1'] + uafxm[:,1] = ui['NUAFXM1'] * sf if NUADFG[2] > 0: - uafxm[:,2] = ui['NUAFXM2'] + uafxm[:,2] = ui['NUAFXM2'] * sf if NUADFG[3] > 0: - uafxm[:,3] = ui['NUAFXM3'] - - # convert units to internal - if uunits == 1: # convert from lb/ac.day to mg.ft3/l.ft2.ivl - uafxm[:,1] *= 0.3677 * delt60 / 24.0 - uafxm[:,2] *= 0.3677 * delt60 / 24.0 - uafxm[:,3] *= 0.3677 * delt60 / 24.0 - else: # convert from kg/ha.day to mg.m3/l.m2.ivl - uafxm[:,1] *= 0.1 * delt60 / 24.0 - uafxm[:,2] *= 0.1 * delt60 / 24.0 - uafxm[:,3] *= 0.1 * delt60 / 24.0 + uafxm[:,3] = ui['NUAFXM3'] * sf # conversion factors - table-type conv-val1 cvbo = ui['CVBO'] diff --git a/HSP2/NUTRX_Class.py b/HSP2/NUTRX_Class.py new file mode 100644 index 00000000..293b3927 --- /dev/null +++ b/HSP2/NUTRX_Class.py @@ -0,0 +1,1232 @@ +import numpy as np +from numpy import zeros, array +from math import log +import numba as nb +from numba.experimental import jitclass + +from HSP2.ADCALC import advect +from HSP2.RQUTIL import sink, decbal, benth +from HSP2.OXRX_Class import OXRX_Class +from HSP2.utilities import make_numba_dict, initm + +spec = [ + ('adnh4', nb.float64[:]), + ('ADNHFG', nb.int32), + ('adnhpm', nb.float64[:]), + ('adpo4', nb.float64[:]), + ('ADPOFG', nb.int32), + ('adpopm', nb.float64[:]), + ('AFACT', nb.float64), + ('AMVFG', nb.int32), + ('anaer', nb.float64), + ('benpo4', nb.float64), + ('BENRFG', nb.int32), + ('bentam', nb.float64), + ('bnh4', nb.float64[:]), + ('bnrpo4', nb.float64), + ('bnrpo4', nb.float64), + ('bnrtam', nb.float64), + ('bnrtam', nb.float64), + ('bodno3', nb.float64), + ('bodpo4', nb.float64), + ('bodtam', nb.float64), + ('bpcntc', nb.float64), + ('bpo4', nb.float64[:]), + ('brpo4', nb.float64[:]), + ('brtam', nb.float64[:]), + ('conv', nb.float64), + ('cvbn', nb.float64), + ('cvbo', nb.float64), + ('cvbp', nb.float64), + ('cvbpc', nb.float64), + ('cvbpn', nb.float64), + ('cvoc', nb.float64), + ('cvon', nb.float64), + ('cvop', nb.float64), + ('decco2', nb.float64), + ('decnit', nb.float64), + ('decpo4', nb.float64), + ('delt60', nb.float64), + ('delts', nb.float64), + ('denbod', nb.float64), + ('DENFG', nb.int32), + ('denno3', nb.float64), + ('denoxt', nb.float64), + ('dnust', nb.float64[:]), + ('dnust2', nb.float64[:]), + ('dsnh4', nb.float64[:]), + ('dspo4', nb.float64[:]), + ('errors', nb.int64[:]), + ('expnvg', nb.float64), + ('expnvl', nb.float64), + ('ino2', nb.float64), + ('ino3', nb.float64), + ('ipo4', nb.float64), + ('isnh4', nb.float64[:]), + ('ispo4', nb.float64[:]), + ('itam', nb.float64), + ('kno220', nb.float64), + ('kno320', nb.float64), + ('ktam20', nb.float64), + ('nexits', nb.int32), + ('nh3', nb.float64), + ('nh3vlt', nb.float64), + ('nh4', nb.float64), + ('nitdox', nb.float64), + ('nitno2', nb.float64), + ('nitno3', nb.float64), + ('nittam', nb.float64), + ('no2', nb.float64), + ('NO2FG', nb.int32), + ('no3', nb.float64), + ('nuaddr', nb.float64[:]), + ('nuadep', nb.float64[:]), + ('NUADCN', nb.float64[:,:]), + ('NUADFG', nb.int32[:]), + ('NUADFX', nb.float64[:,:]), + ('nuadpm', nb.float64[:]), + ('nuadwt', nb.float64[:]), + ('nucf1', nb.float64[:]), + ('nucf2', nb.float64[:,:]), + ('nucf3', nb.float64[:,:]), + ('nucf4', nb.float64[:]), + ('nucf5', nb.float64[:]), + ('nucf6', nb.float64[:]), + ('nucf7', nb.float64[:]), + ('nucf8', nb.float64[:,:]), + ('nuecnt', nb.float64[:]), + ('nust', nb.float64[:,:]), + ('ono2', nb.float64[:]), + ('ono3', nb.float64[:]), + ('opo4', nb.float64[:]), + ('osnh4', nb.float64[:,:]), + ('ospo4', nb.float64[:,:]), + ('otam', nb.float64[:]), + ('PHFLAG', nb.int32), + ('phval', nb.float64), + ('phvalm', nb.float64), + ('PLKFG', nb.int32), + ('po4', nb.float64), + ('PO4FG', nb.int32), + ('rnh3', nb.float64), + ('rnh4', nb.float64), + ('rno2', nb.float64), + ('rno3', nb.float64), + ('rono2', nb.float64), + ('rono3', nb.float64), + ('ropo4', nb.float64), + ('rosnh4', nb.float64[:]), + ('rospo4', nb.float64[:]), + ('rotam', nb.float64), + ('rpo4', nb.float64), + ('rrno2', nb.float64), + ('rrno3', nb.float64), + ('rrpo4', nb.float64), + ('rrtam', nb.float64), + ('rsed', nb.float64[:]), + ('RSED1', nb.float64[:]), + ('RSED2', nb.float64[:]), + ('RSED3', nb.float64[:]), + ('RSED4', nb.float64[:]), + ('RSED5', nb.float64[:]), + ('RSED6', nb.float64[:]), + ('RSED7', nb.float64[:]), + ('rsnh4', nb.float64[:]), + ('rspo4', nb.float64[:]), + ('rtam', nb.float64), + ('SEDFG', nb.int32), + ('simlen', nb.int32), + ('snh4', nb.float64[:]), + ('spo4', nb.float64[:]), + ('svol', nb.float64), + ('tam', nb.float64), + ('TAMFG', nb.int32), + ('tcden', nb.float64), + ('tcnit', nb.float64), + ('tnucf1', nb.float64[:]), + ('tnucf2', nb.float64[:,:]), + ('tnuif', nb.float64[:]), + ('totno3', nb.float64), + ('totpo4', nb.float64), + ('tottam', nb.float64), + ('uunits', nb.int32), + ('vol', nb.float64), + ('volnh3', nb.float64) +] + +@jitclass(spec) +class NUTRX_Class: + + #------------------------------------------------------------------- + # class initialization: + #------------------------------------------------------------------- + def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts, OXRX): + + ''' Initialize instance variables for nutrient simulation ''' + + 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 + + # table-type nut-flags + self.TAMFG = int(ui['NH3FG']) + self.NO2FG = int(ui['NO2FG']) + self.PO4FG = int(ui['PO4FG']) + self.AMVFG = int(ui['AMVFG']) + self.DENFG = int(ui['DENFG']) + self.ADNHFG = int(ui['ADNHFG']) + self.ADPOFG = int(ui['ADPOFG']) + self.PHFLAG = int(ui['PHFLAG']) + + self.PLKFG = int(ui_rq['PLKFG']) + self.SEDFG = int(ui_rq['SEDFG']) + self.BENRFG = int(ui_rq['BENRFG']) + + # table-type nut-ad-flags + self.NUADFG = zeros(7, dtype=np.int32) + + for j in range(1, 7): + self.NUADFG[j] = int(ui['NUADFG' + str(j)]) + + # error handling: + if self.TAMFG == 0 and (self.AMVFG == 1 or self.ADNHFG == 1): + self.errors[0] += 1 + # ERRMSG: tam is not being simulated and nh3 volat. or + # nh4 adsorption is being simulated + + if (self.PO4FG == 0 and self.ADPOFG == 1): + self.errors[1] += 1 + # ERRMSG: po4 is not being simulated, and + # po4 adsorption is being simulated + + if (self.ADNHFG == 1 or self.ADPOFG == 1) and self.SEDFG == 0: + self.errors[2] += 1 + # ERRMSG: sediment associated nh4 and/or po4 is being simulated,but sediment is not being simulated in section sedtrn + + # atmospheric deposition - initialize time series: + self.NUADFX = zeros((self.simlen,4)) + self.NUADCN = zeros((self.simlen,4)) + + for j in range(1, 4): + if 'NUADFX' + str(j) in ts: self.NUADFX[:,j] = ts['NUADFX' + str(j)] + if 'NUADCN' + str(j) in ts: self.NUADCN[:,j] = ts['NUADCN' + str(j)] + + ''' + NUADFX = zeros(self.simlen) + NUADCN = zeros(self.simlen) + + for j in range (1, 4): + n = (2 * j) - 1 + nuadfg1 = self.NUADFG[n] + nuadfg2 = self.NUADFG[n+1] + + # dry deposition: + if nuadfg1 > 0: + NUADFX = initm(siminfo, ui, nuadfg1, 'NUTRX_MONTHLY/NUADFX', 0.0) + elif nuadfg1 == -1: + NUADFX = ts['NUADFX' + str(j)] + + # wet deposition: + if nuadfg2 > 0: + NUADCN = initm(siminfo, ui, nuadfg2, 'NUTRX_MONTHLY/NUADCN', 0.0) + elif nuadfg2 == -1: + NUADCN = ts['NUADCN' + str(j)] + ''' + + ''' + sf = 1.0 + + if self.uunits == 1: # convert from lb/ac.day to mg.ft3/l.ft2.ivl + sf = 0.3677 * delt60 / 24.0 + else: # convert from kg/ha.day to mg.m3/l.m2.ivl + sf = 0.1 * delt60 / 24.0 + ''' + + ''' + # convert units to internal: + if self.NUADFG[1] > 0: + self.uafxm[:,1] = ui['NUADFXM1'] * sf + if self.NUADFG[2] > 0: + self.uafxm[:,2] = ui['NUADFXM2'] * sf + if self.NUADFG[3] > 0: + self.uafxm[:,3] = ui['NUADFXM3'] * sf + ''' + + # conversion factors - table-type conv-val1 + self.cvbo = ui['CVBO'] + self.cvbpc = ui['CVBPC'] + self.cvbpn = ui['CVBPN'] + self.bpcntc = ui['BPCNTC'] + + # calculate derived values + self.cvbp = (31.0 * self.bpcntc) / (1200.0 * self.cvbpc) + self.cvbn = 14.0 * self.cvbpn * self.cvbp / 31.0 + + self.cvoc = self.bpcntc / (100.0 * self.cvbo) + self.cvon = self.cvbn / self.cvbo + self.cvop = self.cvbp / self.cvbo + + # benthic release parameters - table-type nut-benparm + self.anaer = ui['ANAER'] + self.brtam = zeros(2) + self.brpo4 = zeros(2) + + if self.BENRFG == 1 or self.PLKFG == 1: # benthal release parms - table-type nut-benparm + self.brtam[0] = ui['BRNIT1'] * delt60 # convert units from 1/hr to 1/ivl + self.brtam[1] = ui['BRNIT2'] * delt60 # convert units from 1/hr to 1/ivl + self.brpo4[0] = ui['BRPO41'] * delt60 # convert units from 1/hr to 1/ivl + self.brpo4[1] = ui['BRPO42'] * delt60 # convert units from 1/hr to 1/ivl + + self.bnrtam = 0.0 + self.bnrpo4 = 0.0 + + # nitrification parameters - table-type nut-nitdenit + self.ktam20 = ui['KTAM20'] * delt60 # convert units from 1/hr to 1/ivl + self.kno220 = ui['KNO220'] * delt60 # convert units from 1/hr to 1/ivl + self.tcnit = ui['TCNIT'] + self.kno320 = ui['KNO320'] * delt60 # convert units from 1/hr to 1/ivl + self.tcden = ui['TCDEN'] + self.denoxt = ui['DENOXT'] + + if self.TAMFG == 1 and self.AMVFG == 1: # ammonia volatilization parameters table nut-nh3volat + self.expnvg = ui['EXPNVG'] + self.expnvl = ui['EXPNVL'] + + if self.TAMFG == 1 and self.PHFLAG == 3: # monthly ph values table mon-phval, not in RCHRES.SEQ + self.phvalm = ui['PHVALM'] + + #self.nupm3 = zeros(7) + self.nuadpm = zeros(7) + self.rsnh4 = zeros(13) + self.rspo4 = zeros(13) + + # sediment mass storages: + self.RSED1 = zeros(self.simlen); self.RSED2 = zeros(self.simlen) + self.RSED3 = zeros(self.simlen); self.RSED4 = zeros(self.simlen) + self.RSED5 = zeros(self.simlen); self.RSED6 = zeros(self.simlen) + + if 'RSED1' in ts: self.RSED1 = ts['RSED1'] + if 'RSED2' in ts: self.RSED2 = ts['RSED2'] + if 'RSED3' in ts: self.RSED3 = ts['RSED3'] + if 'RSED4' in ts: self.RSED4 = ts['RSED4'] + if 'RSED5' in ts: self.RSED5 = ts['RSED5'] + if 'RSED6' in ts: self.RSED6 = ts['RSED6'] + if 'RSED7' in ts: self.RSED7 = ts['RSED7'] + + self.rsed = zeros(8) + cf = 3.121e-8 if self.uunits == 1 else 1.00e-6 + + if self.SEDFG: + self.rsed[1] = ui_rq['SSED1'] * self.vol + self.rsed[2] = ui_rq['SSED2'] * self.vol + self.rsed[3] = ui_rq['SSED3'] * self.vol + self.rsed[4] = self.rsed[1] + self.rsed[2] + self.rsed[3] + + self.rsed[5] = self.RSED5[0] / cf + self.rsed[6] = self.RSED6[0] / cf + self.rsed[7] = self.RSED7[0] / cf + + # bed sediment concentrations of nh4 and po4 - table nut-bedconc, not in RCHRES.SEQ + # initialize constant bed concentrations (NH4, PO4) + # (convert concentrations from mg/kg to internal units of mg/mg) + self.bnh4 = zeros(4); self.bpo4 = zeros(4) + + for i in range(1, 4): + key = 'BNH4' + str(i) + if key in ui: self.bnh4[i] = ui[key] / 1.0e6 + + key = 'BPO4' + str(i) + if key in ui: self.bpo4[i] = ui[key] / 1.0e6 + + self.adnhpm = zeros(4); self.adpopm = zeros(4) + + if (self.TAMFG == 1 and self.ADNHFG == 1) or (self.PO4FG == 1 and self.ADPOFG == 1): + #self.nupm3[:] = ui['NUPM3'] / 1.0E6 # convert concentrations from mg/kg to internal units of mg/mg + + # initialize adsorbed nutrient mass storages in bed + self.rsnh4[8] = 0.0 + self.rspo4[8] = 0.0 + + for i in range(5, 8): + self.rsnh4[i] = self.bnh4[i-4] * self.rsed[i] + self.rspo4[i] = self.bpo4[i-4] * self.rsed[i] + self.rsnh4[8] += self.rsnh4[i] + self.rspo4[8] += self.rspo4[i] + + + # adsorption parameters - table-type nut-adsparm + for i in range(1, 4): + self.adnhpm[i] = ui['ADNHPM(' + str(i) + ')'] + self.adpopm[i] = ui['ADPOPM(' + str(i) + ')'] + + # initial conditions - table-type nut-dinit + self.dnust = zeros(7); self.dnust2 = zeros(7) + self.dnust[1] = ui['NO3']; self.dnust2[1] = self.dnust[1] * self.vol + self.dnust[2] = ui['TAM']; self.dnust2[2] = self.dnust[2] * self.vol + self.dnust[3] = ui['NO2']; self.dnust2[3] = self.dnust[3] * self.vol + self.dnust[4] = ui['PO4']; self.dnust2[4] = self.dnust[4] * self.vol + + self.phval = 0.0 + + if self.TAMFG == 1: # do the tam-associated initial values (nh4 nh3 phval) + self.phval = ui['PHVAL'] + + # assume nh4 and nh3 are 0.99 x tam and 0.01 x tam respectively + self.dnust[5] = 0.99 * self.dnust[2] + self.dnust2[5] = self.dnust[5] * self.vol + self.dnust[6] = 0.01 * self.dnust[2] + self.dnust2[6] = self.dnust[6] * self.vol + + self.snh4 = zeros(4) + self.spo4 = zeros(4) + + if (self.TAMFG == 1 and self.ADNHFG == 1) or (self.PO4FG == 1 and self.ADPOFG == 1): + # suspended sediment concentrations of nh4 and po4 - table nut-adsinit + # (input concentrations are mg/kg - these are converted to mg/mg for + # internal computations) + for i in range(1, 4): + self.snh4[i] = ui['SNH4' + str(i)] / 1.0e6 # suspended nh4 (sand, silt, clay) + self.spo4[i] = ui['SPO4' + str(i)] / 1.0e6 # suspended po4 (sand, silt, clay) + + # initialize adsorbed nutrient mass storages in suspension + self.rsnh4[4] = 0.0 + self.rspo4[4] = 0.0 + + for i in range(1, 4): + self.rsnh4[i] = self.snh4[i] * self.rsed[i] + self.rspo4[i] = self.spo4[i] * self.rsed[i] + self.rsnh4[4] += self.rsnh4[i] + self.rspo4[4] += self.rspo4[i] + + # initialize totals on sand, silt, clay, and grand total + self.rsnh4[9] = self.rsnh4[1] + self.rsnh4[5] + self.rsnh4[10] = self.rsnh4[2] + self.rsnh4[6] + self.rsnh4[11] = self.rsnh4[3] + self.rsnh4[7] + self.rsnh4[12] = self.rsnh4[4] + self.rsnh4[8] + self.rspo4[9] = self.rspo4[1] + self.rspo4[5] + self.rspo4[10] = self.rspo4[2] + self.rspo4[6] + self.rspo4[11] = self.rspo4[3] + self.rspo4[7] + self.rspo4[12] = self.rspo4[4] + self.rspo4[8] + + # initialize total storages of nutrients in reach + self.nust = zeros((5,2)) + + self.nust[1,1] = self.dnust2[1] + self.nust[2,1] = self.dnust2[2] + if self.ADNHFG == 1: + self.nust[2,1] += self.rsnh4[4] + + self.nust[3,1] = self.dnust2[3] + self.nust[4,1] = self.dnust2[4] + if self.ADPOFG == 1: + self.nust[4,1] += self.rspo4[4] + + # initialize nutrient flux if nutrient is not simulated + self.otam = zeros(nexits); self.ono2 = zeros(nexits); self.opo4 = zeros(nexits) + self.rosnh4 = zeros(5); self.rospo4 = zeros(5) + self.dspo4 = zeros(5); self.dsnh4 = zeros(5) + self.adpo4 = zeros(5); self.adnh4 = zeros(5) + self.ospo4 = zeros((nexits, 5)); self.osnh4 = zeros((nexits, 5)) + self.nucf1 = zeros(5) + self.nucf4 = zeros(8) + self.nucf5 = zeros(9) + self.nucf6 = zeros(2) + self.nucf7 = zeros(7) + + self.nucf2 = zeros((5,3)) + self.nucf3 = zeros((5,3)) + self.nucf8 = zeros((5,3)) + + self.tnucf1 = zeros(5) + self.tnucf2 = zeros((nexits,5)) + self.tnuif = zeros(5) + + # initialize outflow variables: + self.rono3 = 0.0; self.ono3 = zeros(nexits) + self.rotam = 0.0; self.otam = zeros(nexits) + self.rono2 = 0.0; self.ono2 = zeros(nexits) + self.ropo4 = 0.0; self.opo4 = zeros(nexits) + + # initialize process variables: + self.decnit = 0.0; self.decpo4 = 0.0; self.decco2 = 0.0 + self.nitdox = 0.0; self.denbod = 0.0; self.nittam = 0.0 + self.bnrtam = 0.0; self.volnh3 = 0.0; self.bodtam = 0.0 + self.nitno2 = 0.0; self.nitno3 = 0.0; self.denno3 = 0.0 + self.bodno3 = 0.0 + self.bnrpo4 = self.bodpo4 = 0.0 + self.nh3vlt = 0.0 + + self.nuecnt = zeros(4) + + # initialize nutrient states: + self.no3 = ui['NO3'] + self.tam = ui['TAM'] + self.no2 = ui['NO2'] + self.po4 = ui['PO4'] + + self.nh3 = 0.01 * self.tam + self.nh4 = 0.99 * self.tam + + # initialize total nutrient masses: + self.update_mass() + + return + + def simulate(self, loop, tw, wind, phval, OXRX, ino3, itam, ino2, ipo4, isnh4, ispo4, + prec, sarea, scrfac, avdepe, depcor, sed_depscr, sed_rosed, sed_osed, advectData): + ''' Determine primary inorganic nitrogen and phosphorus balances''' + + # hydraulics: + (nexits, vols, vol, srovol, erovol, sovol, eovol) = advectData + + self.vol = vol + + # inflows: convert from [mass/ivld] to [conc.*vol/ivld] + self.ino3 = ino3 / self.conv + self.itam = itam / self.conv + self.ino2 = ino2 / self.conv + self.ipo4 = ipo4 / self.conv + + self.isnh4 = isnh4 / self.conv + self.ispo4 = ispo4 / self.conv + + #compute atmospheric deposition influx (! - to be implemented) + self.nuadep = zeros(7) + self.nuaddr = zeros(7) + self.nuadwt = zeros(7) + + ''' LTI - temporarily removed + for i in range(1,4): + n = 2*(i-1)+ 1 + # dry deposition + if self.NUADFG[n] <= -1: + nuadfx = ts['nuadfx'] + self.nuaddr[i] = sarea*nuadfx + elif self.NUADFG[n] >= 1: + nuaddr[i] = sarea*dayval(nuafxm[mon,i],nuafxm[nxtmon,i],day, ndays) + else: + self.nuaddr[i] = 0.0 + # wet deposition + if self.NUADFG[n+1] <= -1: + nuadcn = ts['nuadcn'] + nuadwt[n]= prec*sarea*nuadcn + elif self.NUADFG[n+1] >= 1: + nuadwt[i] = prec*sarea*dayval(nuacnm[mon,i],nuacnm[nxtmon,i],day,ndays) + else: + nuadwt[i] = 0.0 + + nuadep[i] = nuaddr[i] + nuadwt[i] + ''' + + # advect nitrate + self.tnuif[1] = self.ino3 + inno3 = self.ino3 + self.nuadep[1] + + self.no3, self.rono3, self.ono3 = \ + advect(inno3, self.no3, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + + self.tnucf1[1] = self.rono3 + if self.nexits > 1: + self.tnucf2[:,1] = self.ono3[:] # nexits + + # advect total ammonia: + if self.TAMFG == 1: + intam = self.itam + self.nuadep[2] + + self.tam, self.rotam, self.otam = \ + advect(intam, self.tam, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + + self.tnucf1[2] = self.rotam + if self.nexits > 1: + self.tnucf2[:,2] = self.otam[:] + + # advect nitrite: + if self.NO2FG == 1: + self.tnuif[3] = self.ino2 + + self.no2, self.rono2, self.ono2 = \ + advect(self.ino2, self.no2, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + + self.tnucf1[3] = self.rono2 + if self.nexits > 1: + self.tnucf2[:,3] = self.ono2[:] # nexits + + # advect dissolved PO4: + if self.PO4FG == 1: + inpo4 = self.ipo4 + self.nuadep[3] + self.po4, self.ropo4, self.opo4 = advect(inpo4,self.po4, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + + # sediment variables (require unit conversion): + self.rsed = zeros(8) + rosed = zeros(4) + osed = zeros((self.nexits,4)) + depscr = zeros(4) + + cf = 3.121e-8 if self.uunits == 1 else 1.00e-6 + + if self.SEDFG: + self.rsed[1] = self.RSED1[loop] / cf + self.rsed[2] = self.RSED2[loop] / cf + self.rsed[3] = self.RSED3[loop] / cf + self.rsed[4] = self.RSED4[loop] / cf + self.rsed[5] = self.RSED5[loop] / cf + self.rsed[6] = self.RSED6[loop] / cf + self.rsed[7] = self.RSED7[loop] / cf + + for j in range(1, 4): + rosed[j] = sed_rosed[j] / cf + depscr[j] = sed_depscr[j] / cf + + for i in range(nexits): + osed[i,j] = sed_osed[i,j] / cf + + # advect adsorbed PO4: + if self.PO4FG == 1 and self.ADPOFG == 1: + # zero the accumulators + self.ispo4[4] = 0.0 + self.dspo4[4] = 0.0 + self.rospo4[4] = 0.0 + if self.nexits > 1: + self.ospo4[:,4] = 0.0 # nexits + + # repeat for each sediment fraction (LTI) + for j in range(1, 4): # get data on sediment-associated phosphate + osed_ = osed[:,j] # all exits for sed class "j" + ospo4_ = zeros(nexits) + + (self.spo4[j], self.dspo4[j], self.rospo4[j], ospo4_) \ + = self.advnut(self.ispo4[j],self.rsed[j],self.rsed[j+4],depscr[j],rosed[j],osed_,self.nexits, + self.rspo4[j],self.rspo4[j + 4],self.bpo4[j]) + + self.ospo4[:,j] = ospo4_ + self.ispo4[4] += self.ispo4[j] + self.dspo4[4] += self.dspo4[j] + self.rospo4[4] += self.rospo4[j] + + if self.nexits > 1: + self.ospo4[:,4] += self.ospo4[:,j] # nexits + + self.tnuif[4] = self.ipo4 + self.ispo4[4] + self.tnucf1[4] = self.ropo4 + self.rospo4[4] + if self.nexits > 1: + self.tnucf2[:,4] = self.opo4[:]+ self.ospo4[:,4] # nexits + else: # no adsorbed fraction + self.tnuif[4] = self.ipo4 + self.tnucf1[4] = self.ropo4 + if self.nexits > 1: + self.tnucf2[:,4] = self.opo4[:] + + # advect adsorbed ammonium + if self.TAMFG == 1 and self.ADNHFG == 1: # advect adsorbed ammonium + + # zero the accumulators + self.isnh4[4] = 0.0; + self.dsnh4[4] = 0.0 + self.rosnh4[4] = 0.0 + if self.nexits > 1: + self.osnh4[:,4] = 0.0 # nexits + + # repeat for each sediment fraction + for j in range(1, 4): + osed_ = osed[:,j] # all exits for sed class "j" + osnh4_ = zeros(nexits) + + (self.snh4[j],self.dsnh4[j],self.rosnh4[j],osnh4_) \ + = self.advnut(self.isnh4[j],self.rsed[j],self.rsed[j + 3],depscr[j],rosed[j],osed_,self.nexits, + self.rsnh4[j],self.rsnh4[j + 4],self.bnh4[j]) + + self.osnh4[:,j] = osnh4_ + self.isnh4[4] += self.isnh4[j] + self.dsnh4[4] += self.dsnh4[j] + self.rosnh4[4] += self.rosnh4[j] + + if self.nexits > 1: + self.osnh4[:,4] += self.osnh4[:,j] # nexits + + self.tnuif[2] = self.itam + self.isnh4[4] + self.tnucf1[2] = self.rotam + self.rosnh4[4] + if self.nexits > 1: + self.tnucf2[:,2] = self.otam[:] + self.osnh4[:,4] # nexits + + else: # no adsorbed fraction + self.tnuif[2] = self.itam + self.tnucf1[2] = self.rotam + if self.nexits > 1: + self.tnucf2[:,2] = self.otam[:] # nexits + + + # calculate ammonia ionization in water column + if self.TAMFG == 1: + # get ph values + # assign last computed value from RQUAL (i.e., via time series, monthly inputs, or constant): + if (phval >= 0.0): + self.phval = phval + + # compute ammonia ionization + (self.nh3, self.nh4) = self.ammion(tw, self.phval, self.tam) + + if avdepe > 0.17: + if self.BENRFG == 1: + # simulate benthal release of inorganic nitrogen and + # ortho-phosphorus; and compute associated fluxes + if self.TAMFG == 1: + (self.tam, self.bentam) = benth(OXRX.dox,self.anaer,self.brtam,scrfac,depcor,self.tam) + self.bnrtam = self.bentam * self.vol + + if self.PO4FG == 1: + self.po4, self.benpo4 = benth(OXRX.dox,self.anaer,self.brpo4,scrfac,depcor, self.po4) + self.bnrpo4 = self.benpo4 * self.vol + + if self.TAMFG == 1: + if self.AMVFG == 1: # compute ammonia volatilization + twkelv = tw + 273.16 # convert water temperature to degrees kelvin + avdepm = avdepe * 0.3048 # convert depth to meters + (self.tam, self.nh3vlt) = self.nh3vol(self.expnvg,self.expnvl,OXRX.korea,wind,self.delt60,self.delts,avdepm,twkelv,tw,self.phval,self.tam) + self.volnh3 = -self.nh3vlt * self.vol + else: + self.volnh3 = 0.0 + + # calculate amount of nitrification; nitrification does not + # take place if the do concentration is less than 2.0 mg/l + (self.tam,self.no2,self.no3,OXRX.dox,dodemd,tamnit,no2ntc,no3nit) = \ + self.nitrif(self.ktam20,self.tcnit,tw,self.NO2FG,self.kno220,self.tam,self.no2,self.no3,OXRX.dox) + + # compute nitrification fluxes + self.nitdox = -dodemd * self.vol + self.nittam = -tamnit * self.vol + self.nitno2 = no2ntc * self.vol + self.nitno3 = no3nit * self.vol + + if self.DENFG == 1: # consider denitrification processes, and compute associated fluxes + no3de = 0.0 + (self.no3, no3de) = self.denit(self.kno320, self.tcden, tw, OXRX.dox, self.denoxt, self.no3) + self.denno3 = -no3de * self.vol + + # calculate amount of inorganic constituents released by bod decay in reach water + self.decnit = OXRX.bodox * self.cvon + self.decpo4 = OXRX.bodox * self.cvop + self.decco2 = OXRX.bodox * self.cvoc + + # update state variables of inorganic constituents which + # are end products of bod decay; and compute associated fluxes + (self.tam, self.no3, self.po4) = decbal(self.TAMFG, self.PO4FG, self.decnit, self.decpo4, + self.tam, self.no3, self.po4) + if self.TAMFG == 1: + self.bodtam = self.decnit * self.vol + else: + self.bodno3 = self.decnit * self.vol + + if self.PO4FG == 1: + self.bodpo4 = self.decpo4 * self.vol + + if self.PO4FG == 1 and self.SEDFG == 1 and self.ADPOFG == 1: # compute adsorption/desorption of phosphate + dumxxx = 0.0 + + (self.po4, self.spo4, dumxxx, self.adpo4) \ + = self.addsnu(self.vol, self.rsed, self.adpopm, self.po4, self.spo4, dumxxx, self.adpo4) + + if self.TAMFG == 1 and self.SEDFG == 1 and self.ADNHFG == 1: # compute adsorption/desorption of ammonium + # first compute ammonia ionization + (self.nh3, self.nh4) = self.ammion(tw, self.phval, self.tam) + (self.nh4, self.snh4, self.tam, self.adnh4) \ + = self.addsnu(self.vol, self.rsed, self.adnhpm, self.nh4, self.snh4, self.tam, self.adnh4) + # then re-compute ammonia ionization + (self.nh3, self.nh4) = self.ammion(tw, self.phval, self.tam) + else: + # too little water is in reach to warrant simulation of quality processes + self.decnit = 0.0; self.decpo4 = 0.0; self.decco2 = 0.0 + self.nitdox = 0.0; self.denbod = 0.0; self.nittam = 0.0 + self.bnrtam = 0.0; self.volnh3 = 0.0; self.bodtam = 0.0 + self.nitno2 = 0.0; self.nitno3 = 0.0; self.denno3 = 0.0 + self.bodno3 = 0.0 + self.bnrpo4 = self.bodpo4 = 0.0 + + self.adnh4[1:5] = 0.0 + self.adpo4[1:5] = 0.0 + + #self.totdox = self.readox + self.boddox + self.bendox + self.nitdox + #self.totbod = self.decbod + self.bnrbod + self.snkbod + self.denbod + self.totno3 = self.nitno3 + self.denno3 + self.bodno3 + self.tottam = self.nittam + self.volnh3 + self.bnrtam + self.bodtam + self.totpo4 = self.bnrpo4 + self.bodpo4 + + if self.PO4FG == 1 and self.SEDFG == 1 and self.ADPOFG == 1: # find total quantity of phosphate on various forms of sediment + totpm1 = 0.0; totpm2 = 0.0; totpm3 = 0.0 + + for j in range(1, 4): + self.rspo4[j] = self.spo4[j] * self.rsed[j] # compute mass of phosphate adsorbed to each suspended fraction + self.rspo4[j + 4] = self.bpo4[j] * self.rsed[j + 3] # compute mass of phosphate adsorbed to each bed fraction + self.rspo4[j + 8] = self.rspo4[j] + self.rspo4[j + 4] # compute total mass of phosphate on each sediment fraction + + totpm1 += self.rspo4[j] + totpm2 += self.rspo4[j + 4] + totpm3 += self.rspo4[j + 8] + + self.rspo4[4] = totpm1 # compute total suspended phosphate + self.rspo4[8] = totpm2 # compute total bed phosphate + self.rspo4[12] = totpm3 # compute total sediment-associated phosphate + + # calculate total amount of ammonium on various forms of sediment + if self.TAMFG == 1 and self.SEDFG == 1 and self.ADNHFG == 1: + totnm1 = 0.0; totnm2 = 0.0; totnm3 = 0.0 + + for j in range(1, 4): + self.rsnh4[j] = self.snh4[j] * self.rsed[j] # compute mass of ammonium adsorbed to each suspended fraction + self.rsnh4[j + 4] = self.bnh4[j] * self.rsed[j + 3] # compute mass of ammonium adsorbed to each bed fraction + self.rsnh4[j + 8] = self.rsnh4[j] + self.rsnh4[j + 4] # compute total mass of ammonium on each sediment fraction + + totnm1 += self.rsnh4[j] + totnm2 += self.rsnh4[j + 4] + totnm3 += self.rsnh4[j + 8] + + self.rsnh4[4] = totnm1 # compute total suspended ammonium + self.rsnh4[8] = totnm2 # compute total bed ammonium + self.rsnh4[12] = totnm3 # compute total sediment-associated ammonium + + + self.svol = self.vol # svol is volume at start of time step, update for next time thru + + # update total nutrient masses (handle after PLANK): + #self.update_mass() + + + return OXRX + + def update_mass(self): + # calculate total resident mass of nutrient species + + self.rno3 = self.no3 * self.vol + self.rtam = self.tam * self.vol + self.rno2 = self.no2 * self.vol + self.rpo4 = self.po4 * self.vol + self.rnh4 = self.nh4 * self.vol + self.rnh3 = self.nh3 * self.vol + + self.rrno3 = self.no3 * self.vol + self.rrtam = self.tam * self.vol + + if self.ADNHFG == 1: + self.rrtam += self.rsnh4[4] # add adsorbed suspended nh4 to dissolved + + self.rrno2 = self.no2 * self.vol + self.rrpo4 = self.po4 * self.vol + + if self.ADPOFG == 1: + self.rrpo4 += self.rspo4[4] # add adsorbed suspended po4 to dissolved + + return + + #-------------------------------------------------------------- + # static methods + #-------------------------------------------------------------- + @staticmethod + def addsnu(vol, rsed, adpm, dnut, snut, dnutxx, adnut): + ''' simulate exchange of nutrient (phosphate or ammonium) between the + dissolved state and adsorption on suspended sediment- 3 adsorption + sites are considered: 1- suspended sand 2- susp. silt + 3- susp. clay + assumes instantaneous linear equilibrium''' + + if vol > 0.0: # adsorption/desorption can take place + # establish nutrient equilibrium between reach water and suspended sediment; first find the new dissolved nutrient conc. in reach water + dnutin = dnut + num = vol * dnut + denom = vol + + for j in range(1, 4): + if rsed[j] > 0.0: # accumulate terms for numerator and denominator in dnut equation + num += snut[j] * rsed[j] + denom += adpm[j] * rsed[j] + + dnut = num / denom # calculate new dissolved concentration-units are mg/l + dnutxx= dnutxx - (dnutin - dnut) # also calculate new tam conc if doing nh4 adsorption + + # calculate new conc on each sed class and the corresponding adsorption/desorption flux + adnut[4] = 0.0 + + for j in range(1, 4): + if rsed[j] > 0.0: # this sediment class is present-calculate data pertaining to it + temp = dnut * adpm[j] # new concentration + + # quantity of material transferred + # adnut[j]= (temp - snut[j])*rsed[j] + snut[j] = temp + + # accumulate total adsorption/desorption flux above bed + adnut[4] += adnut[j] + + else: # this sediment class is absent + adnut[j] = 0.0 + # snut[j] is unchanged-"undefined" + + else: # no water, no adsorption/desorption + adnut[1:5] = 0.0 + # snut(1 thru 3) and dnut should already have been set to undefined values + + return dnut, snut, dnutxx, adnut + + + def advnut(self,isnut,rsed,bsed,depscr,rosed,osed,nexits,rsnuts,rbnuts,bnut): + + ''' simulate the advective processes, including deposition and scour for the + inorganic nutrient adsorbed to one sediment size fraction''' + + if depscr < 0.0: # there was sediment scour during the interval + # compute flux of nutrient mass into water column with scoured sediment fraction + dsnut = bnut * depscr + + # calculate concentration in suspension-under these conditions, denominator should never be zero + snut = (isnut + rsnuts - dsnut) / (rsed + rosed) + rosnut = rosed * snut + else: # there was deposition or no scour/deposition during the interval + denom = rsed + depscr + rosed + if denom == 0.0: # there was no sediment in suspension during the interval + snut = -1.0e30 + rosnut = 0.0 + dsnut = 0.0 + + # fix sed-nut problem caused by very small sediment loads that are stored in + # wdm file as zero (due to wdm attribute tolr > 0.0) when adsorbed nut load + # is not zero; changed comparison from 0.0 to 1.0e-3; this should not cause + # any mass balance errors since the condition is not likely to exist over a + # long period and will be insignificant compared to + # the total mass over a printout period; note that 1.0e-3 mg*ft3/l is 0.028 mg + # (a very, very small mass) + if abs(isnut) > 1.0e-3 or abs(rsnuts) > 1.0e-3: + self.errors[3] += 1 + # errmsg: error-under these conditions these values should be zero + else: # there was some suspended sediment during the interval + # calculate conc on suspended sed + snut = (isnut + rsnuts) / denom + rosnut= rosed * snut + dsnut = depscr * snut + + if rsed == 0.0: + # rchres ended up without any suspended sediment-revise + # value for snut, but values obtained for rosnut, and dsnut are still ok + snut = -1.0e30 + + # calculate conditions on the bed + if bsed == 0.0: + # no bed sediments at end of interval + if abs(dsnut) > 0.0 or abs(rbnuts) > 0.0: + self.errors[4] += 1 + # errmsg: error-under this condition these values should be zero + + osnut = zeros(nexits) + if nexits > 1: + # compute outflow through each individual exit + if rosed == 0.0: # all zero + osnut[:] = 0.0 + else: + osnut[:] = rosnut * osed[:] / rosed + + return snut, dsnut, rosnut, osnut + + @staticmethod + def ammion(tw, ph, tam): + ''' simulate ionization of ammonia to ammonium using empirical relationships developed by loehr, 1973''' + + if tam >= 0.0: # tam is defined, compute fractions + # adjust very low or high values of water temperature to fit limits of dat used to develop empirical relationship + if tw < 5.0: twx = 5.0 + elif tw > 35.0: twx = 35.0 + else: twx = tw + + if ph < 4.0: phx = 4.0 + elif ph > 10.0: phx = 10.0 + else: phx = ph + + # compute ratio of ionization constant values for aqueous ammonia and water at current water temperatue + ratio = (-3.39753 * log(0.02409 * twx)) * 1.0e9 + + # compute fraction of total ammonia that is un-ionized + frac = 10.0**(phx) / (10.0**phx + ratio) + + # update nh3 and nh4 state variables to account for ionization + nh3 = frac * tam + nh4 = tam - nh3 + else: # tam conc undefined + nh3 = -1.0e30 + nh4 = -1.0e30 + return nh3, nh4 + + @staticmethod + def denit(kno320, tcden, tw, dox, denoxt, no3): + ''' calculate amount of denitrification; denitrification does not take place + if the do concentration is above user-specified threshold do value (denoxt)''' + denno3 = 0.0 + + if dox <= denoxt: # calculate amount of no3 denitirified to nitrogen gas + denno3 = 0.0 + if no3 > 0.001: + denno3 = kno320 * (tcden**(tw - 20.0)) * no3 + no3 = no3 - denno3 + if no3 < 0.001: # adjust amount of no3 denitrified so that no3 state variable is not a negative number; set no3 to a value of .001 mg/l + denno3 = denno3 + no3 - 0.001 + no3 = 0.001 + else: + pass # denitrification does not occur + + return no3, denno3 + + + def hcintp (self, phval, tw): + ''' calculate henry's constant for ammonia based on ph and water temperature''' + + xtw = array([4.44, 15.56, 26.67, 37.78]) + xhplus = array([1.0, 10.0, 100.0, 1000.0, 10000.0]) + yhenc = array([[0.000266, 0.000754, 0.00198, 0.00486], + [0.00266, 0.00753, 0.0197, 0.0480], + [0.0263, 0.0734, 0.186, 0.428], + [0.238, 0.586, 1.20, 2.05], + [1.2, 1.94, 2.65, 3.31]]) # dimensions: fortran 4,5 + yhenc = np.transpose(yhenc) + + # adjust very low or very high values of water temperature to fit limits of henry's contant data range + if tw < 4.44: # use low temperature range values for henry's constant (4.4 degrees c or 40 degrees f) + twx = 4.44 + elif tw > 37.78: # use high temperature range values for henry's constant (37.78 degrees c or 100 degrees f) + twx = 37.78 + else: # use unmodified water temperature value in interpolation + twx = tw + + # convert ph value to a modified version of hydrogen ion concentration + # because our interpolation routine cant seem to work with small numbers + hplus = 10.0**(phval) * 1.0e-6 + + # adjust very low or very high values of hydrogen ion concentration to fit limits of henry's constant data range + if hplus > 10000.0: # use low hydrogen ion concentration range values for henry's constant + hplus = 10000.0 + elif hplus < 1.0: # use high hydrogen ion concentration range values for henry's constant + hplus = 1.0 + + # perform two-dimensional interpolation of henry's constant values to estimate henry's + # constant for water temperature and ph conditions in water column (based on p. 97 of numerical recipes) + i4 = 4 + i5 = 5 + yhtmp = zeros(5) + ytwtmp = zeros(4) + + for i in range(4): # do 10 i= 1, 4 + for j in range(5): # do 20 j= 1, 5 + yhtmp[j] = yhenc[i,j] # copy row into temporary storage + # 20 continue + # perform linear interpolation within row of values + ytwtmp[i] = self.intrp1(xhplus, yhtmp, i5, hplus) + # 10 continue + + # do final interpolation in remaining dimension + hcmf = self.intrp1(xtw, ytwtmp, i4, twx) + + # convert henry's constant from molar fraction form to units of atm.m3/mole: assume + # 1) dilute air and water solutions + # 2) ideal gas law + # 3) stp i.e., 1 atm total pressure + # 4) 1 gram water = 1 cm3 + + # xa(air) 1 + # --------- * ----------------------------------------- + # xa(water) (1.e+6 m3/g water)/(18.01 g/mole water) + + hcnh3 = hcmf * (18.01 * 1.e-6) + + return hcnh3 + + + @staticmethod + def intrp1(xarr0, yarr0, len_, xval): + ''' perform one-dimensional interpolation of henry's constant values for ammonia (based on p. 82 of numerical recipes)''' + + c = zeros(11); d = zeros(11) + + # modify array indexing (1-based): + cnt = len(xarr0) + 1 + xarr = zeros(cnt) + yarr = zeros(cnt) + + xarr[1:cnt] = xarr0[0:cnt-1] + yarr[1:cnt] = yarr0[0:cnt-1] + + # interpolate: + ns = 1 + dif = abs(xval-xarr[1]) + # find the index ns of the closest array entry + for i in range(1, len_+1): # do 10 i= 1, len + dift = abs(xval - xarr[i]) + if dift < dif: + ns = i + dif = dift + + # initialize correction array values + c[i] = yarr[i] + d[i] = yarr[i] + + # select intial approximation of yval + yval = yarr[ns] + ns = ns - 1 + # loop over the current values in correction value arrays (c & d) to update them + for j in range(1, len_): # do 30 j = 1, len -1 + + for i in range (1, len_ - j + 1): # do 20 i = 1, len - j + ho = xarr[i] - xval + hp = xarr[i + j] - xval + w = c[i + 1] - d[i] + den = ho - hp + den = w / den + + # update correction array values + d[i] = hp * den + c[i] = ho * den + + # select correction to yval + if 2 * ns < len_-j: + dyval = c[ns + 1] + else: + dyval= d[ns] + ns = ns - 1 + + # compute yval + yval = yval + dyval + + return yval + + def nh3vol(self, expnvg, expnvl, korea, wind, delt60, delts, avdepm, twkelv, tw, phval, tam): + ''' calculate ammonia volatilization using two-layer theory''' + + if tam > 0.0: + # convert reaeration coefficient into units needed for computatuion + # of bulk liquid film gas transfer coefficient (cm/hr) based on + # average depth of water + dokl = korea * (avdepm * 100.0) / delt60 + + # compute bulk liquid film gas transfer coefficient for ammonia using + # equation 183 of mccutcheon; 1.8789 equals the ratio of oxygen + # molecule molecular weight to ammonia molecular weight + nh3kl = dokl * 1.8789**(expnvl / 2.0) + + # convert wind speed from meters/ivl (wind) to meters/sec (windsp) + windsp = wind / delts + + # compute bulk gas film gas transfer coefficient (cm/hr) for ammonia + # using equation 184 of mccutcheon; the product of the expression + # (700.*windsp) is expressed in cm/hr; 1.0578 equals the ratio of water + # molecule molecular weight to ammonia molecular weight + if windsp <= 0.0: + windsp = 0.001 + nh3kg = 700.0 * windsp * 1.0578**(expnvg / 2.0) + + # compute henry's constant for ammonia as a function of temperature + # hcinp() called only here + hcnh3 = self.hcintp(phval, tw) + + # avoid divide by zero errors + chk = nh3kl * hcnh3 + if chk > 0.0: + # compute overall mass transfer coefficient for ammonia (kr) in cm/hr + # using equation 177 of mccutcheon; first calculate the inverse of kr + # (krinv); 8.21e-05 equals ideal gas constant value expressed as + # atm/degrees k mole + krinv = (1.0 / nh3kl) + ((8.21e-05) * twkelv) / (hcnh3 * nh3kg) + kr = (1.0 / krinv) + + # compute reach-specific gas transfer coefficient (units are /interval) + knvol = (kr / (avdepm * 100.0)) * delt60 + else: # korea or hcnh3 was zero (or less) + knvol = 0.0 + + # compute ammonia flux out of reach due to volatilization; assumes that + # equilibrium concentration of ammonia is sufficiently small to be considered zero + nh3vlt = knvol * tam + if nh3vlt >= tam: + nh3vlt = 0.99 * tam + tam = 0.01 * tam + else: + tam = tam - nh3vlt + else: # no ammonia present; hence, no volatilization occurs + nh3vlt = 0.0 + return tam, nh3vlt + + @staticmethod + def nitrif(ktam20, tcnit, tw, no2fg, kno220, tam, no2, no3, dox): + ''' calculate amount of nitrification; nitrification does not take place if the do concentration is less than 2.0 mg/l''' + + if dox >= 2.0: + # calculate amount of tam oxidized to no2; tamnit is expressed as mg tam-n/l + tamnit = 0.0 + if tam > 0.001: + tamnit = ktam20 * (tcnit**(tw - 20.0)) * tam + tam = tam - tamnit + if tam < 0.001: # adjust amount of tam oxidized so that tam state variable is not a negative number; set tam to a value of .001 mg/l + tamnit = tamnit + tam - .001 + tam = .001 + if no2fg == 1: # calculate amount of no2 oxidized to no3; no2nit is expressed as mg no2-n/l + no2nit = 0.0 + if no2 > 0.001: + no2nit = kno220 * (tcnit**(tw - 20.0)) * no2 + + # update no2 state variable to account for nitrification + if no2nit > 0.0: + if no2 + tamnit - no2nit <= 0.0: + no2nit = 0.9 * (no2 + tamnit) + no2 = 0.1 * (no2 + tamnit) + else: + no2 = no2 + tamnit - no2nit + else: + no2 = no2 + tamnit + no2ntc = tamnit - no2nit + else: # no2 is not simulated; tam oxidized is fully oxidized to no3 + no2nit = tamnit + no2ntc = 0.0 + + # update no3 state variable to account for nitrification and compute concentration flux of no3 + no3 = no3 + no2nit + no3nit = no2nit + + # find oxygen demand due to nitrification + dodemd = 3.22 * tamnit + 1.11 * no2nit + + if dox < dodemd: + # adjust nitrification demands on oxygen so that dox will not be zero; + # routine proportionally reduces tam oxidation to no2 and no2 oxidation to no3 + rho = dox / dodemd + if rho < 0.001: + rho = 0.0 + rhoc3 = (1.0 - rho) * tamnit + rhoc2 = (1.0 - rho) * no2nit + tam = tam + rhoc3 + if no2fg == 1: + no2 = no2 - rhoc3 + rhoc2 + no3 = no3 - rhoc2 + dodemd = dox + dox = 0.0 + tamnit = tamnit - rhoc3 + no2nit = no2nit - rhoc2 + no3nit = no3nit - rhoc2 + if no2fg == 1: + no2ntc = no2ntc - rhoc3 + rhoc2 + else: # projected do value is acceptable + dox = dox - dodemd + else: # nitrification does not occur + tamnit = 0.0 + no2nit = 0.0 + dodemd = 0.0 + no2ntc = 0.0 + no3nit = 0.0 + + return tam, no2, no3, dox, dodemd, tamnit, no2ntc, no3nit \ No newline at end of file diff --git a/HSP2/OXRX_Class.py b/HSP2/OXRX_Class.py new file mode 100644 index 00000000..47c4bc14 --- /dev/null +++ b/HSP2/OXRX_Class.py @@ -0,0 +1,346 @@ +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 + +from HSP2.ADCALC import advect, oxrea +from HSP2.RQUTIL import sink +from HSP2.utilities import make_numba_dict + +spec = [ + ('AFACT', nb.float64), + ('benod', nb.float64), + ('bendox', nb.float64), + ('benox', nb.float64), + ('BENRFG', nb.int32), + ('bnrbod', nb.float64), + ('bod', nb.float64), + ('bodbnr', nb.float64), + ('boddox', nb.float64), + ('bodox', nb.float64), + ('BRBOD', nb.float64[:]), + ('cforea', nb.float64), + ('cfpres', nb.float64), + ('conv', nb.float64), + ('decbod', nb.float64), + ('delt60', nb.float64), + ('delth', nb.float64), + ('delts', nb.float64), + ('doben', nb.float64), + ('dorea', nb.float64), + ('dox', nb.float64), + ('errors', nb.int64[:]), + ('expod', nb.float64), + ('expred', nb.float64), + ('exprel', nb.float64), + ('exprev', nb.float64), + ('GQFG', nb.int32), + ('GQALFG4', nb.int32), + ('idox', nb.float64), + ('ibod', nb.float64), + ('kbod20', nb.float64), + ('kodset', nb.float64), + ('korea', nb.float64), + ('len_', nb.float64), + ('LKFG', nb.int32), + ('nexits', nb.int32), + ('obod', nb.float64[:]), + ('odox', nb.float64[:]), + ('rbod', nb.float64), + ('rdox', nb.float64), + ('readox', nb.float64), + ('reak', nb.float64), + ('reakt', nb.float64), + ('REAMFG', nb.int32), + ('rdox', nb.float64), + ('rbod', nb.float64), + ('relbod', nb.float64), + ('robod', nb.float64), + ('rodox', nb.float64), + ('satdo', nb.float64), + ('simlen', nb.int32), + ('snkbod', nb.float64), + ('supsat', nb.float64), + ('svol', nb.float64), + ('tcben', nb.float64), + ('tcbod', nb.float64), + ('tcginv', nb.float64), + ('totdox', nb.float64), + ('totbod', nb.float64), + ('uunits', nb.int32), + ('vol', nb.float64) +] + +@jitclass(spec) +class OXRX_Class: + + #------------------------------------------------------------------- + # class initialization: + #------------------------------------------------------------------- + def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts): + + ''' Initialize variables for primary DO, BOD balances ''' + + 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 + + # gqual flags + self.GQFG = int(ui_rq['GQFG']) + self.GQALFG4 = int(ui_rq['GQALFG4']) + + # table-type ox-genparm + self.kbod20 = ui['KBOD20'] * delt60 # convert units from 1/hr to 1/ivl + self.tcbod = ui['TCBOD'] + self.kodset = ui['KODSET'] * delt60 # convert units from 1/hr to 1/ivl + self.supsat = ui['SUPSAT'] + + # table-type ox-init + self.dox = ui['DOX'] + self.bod = ui['BOD'] + self.satdo = ui['SATDO'] + + # other required values + self.BENRFG = int(ui_rq['BENRFG']) # via table-type benth-flag + self.REAMFG = int(ui['REAMFG']) # table-type ox-flags + elev = ui['ELEV'] # table-type elev + + self.cfpres = ((288.0 - 0.001981 * elev) / 288.0)**5.256 # pressure correction factor - + ui['CFPRES'] = self.cfpres + + self.LKFG = int(ui_rq['LKFG']) + + self.cforea = 0.0 + self.delth = 0.0 + self.reak = 0.0; self.reakt = 1.0 + self.expred = 0.0; self.exprev = 0.0 + self.expod = 0.0; self.exprel = 0.0 + + if self.LKFG == 1: + self.cforea = ui['CFOREA'] # reaeration parameter from table-type ox-cforea + + elif self.REAMFG == 1: # tsivoglou method; table-type ox-tsivoglou + self.reakt = ui['REAKT'] + self.tcginv = ui['TCGINV'] + + self.len_ = ui['LEN'] * 5280.0 # mi to feet + self.delth = ui['DELTH'] + if self.uunits == 2: + self.len_ = ui['LEN'] * 1000.0 # length of reach, in meters + self.delth = ui['DELTH'] * 1000.0 # convert to meters + + elif self.REAMFG == 2: # owen/churchill/o'connor-dobbins; table-type ox-tcginv + self.tcginv = ui['TCGINV'] + self.reak = ui['REAK'] + self.expred = ui['EXPRED'] + + elif self.REAMFG == 3: # user formula - table-type ox-reaparm + self.tcginv = ui['TCGINV'] + self.reak = ui['REAK'] + self.expred = ui['EXPRED'] + + if self.BENRFG == 1: # benthic release parms - table-type ox-benparm + self.benod = ui['BENOD'] * self.delt60 # convert units from 1/hr to 1/ivl + self.tcben = ui['TCBEN'] + self.expod = ui['EXPOD'] + self.exprel = ui['EXPREL'] + + self.BRBOD = zeros(2) + self.BRBOD[0] = ui['BRBOD1'] * self.delt60 # convert units from 1/hr to 1/ivl + self.BRBOD[1] = ui['BRBOD2'] * self.delt60 # convert units from 1/hr to 1/ivl + + #self.BRBOD = array([ui['BRBOD1'] , ui['BRBOD2']]) * self.delt60 # convert units from 1/hr to 1/ivl + + self.snkbod = 0.0 + + self.rdox = self.dox * self.vol + self.rbod = self.bod * self.vol + + self.odox = zeros(nexits) + self.obod = zeros(nexits) + + self.korea = 0.0 + + return + + #------------------------------------------------------------------- + # simulation (single timestep): + #------------------------------------------------------------------- + + def simulate(self, oxif1, oxif2, wind, scrfac, avdepe, avvele, depcor, tw, advectData): + + # hydraulics: + (nexits, vols, vol, srovol, erovol, sovol, eovol) = advectData + + self.vol = vol + + # inflows: convert from [mass/ivld] to [conc.*vol/ivld] + self.idox = oxif1 / self.conv + self.ibod = oxif2 / self.conv + + # advect dissolved oxygen + (self.dox, self.rodox, self.odox) = \ + advect(self.idox, self.dox, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + + # advect bod + (self.bod, self.robod, self.obod) = \ + advect(self.ibod, self.bod, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + + # initialize variables: + self.bodox = 0.0 + self.readox = 0.0 + self.boddox = 0.0 + self.bendox = 0.0 + self.decbod = 0.0 + + if avdepe > 0.17: # benthal influences are considered + # sink bod + self.bod, self.snkbod = sink(self.vol, avdepe, self.kodset, self.bod) + self.snkbod = -self.snkbod + + if self.BENRFG == 1: + #$OXBEN # simulate benthal oxygen demand and benthal release of bod, and compute associated fluxes + # calculate amount of dissolved oxygen required to satisfy benthal oygen demand (mg/m2.ivl) + self.benox = self.benod * (self.tcben**(tw -20.0)) * (1.0 -exp(-self.expod * self.dox)) + + # adjust dissolved oxygen state variable to acount for oxygen lost to benthos, and compute concentration flux + self.doben = self.dox + self.dox = self.dox - (self.benox * depcor) + self.doben = self.benox * depcor + if self.dox >= 0.001: + self.doben = self.benox * depcor + else: + self.dox = 0.0 + + # calculate benthal release of bod; release is a function of dissolved oxygen + # (dox) and a step function of stream velocity; brbod(1) is the aerobic benthal + # release rate; brbod(2) is the base increment to benthal release under + # decreasing do concentration; relbod is expressed as mg bod/m2.ivl + self.relbod = (self.BRBOD[0] + self.BRBOD[1] * exp(-self.exprel * self.dox)) * scrfac + + # add release to bod state variable and compute concentration flux + self.bod = self.bod + self.relbod * depcor + self.bodbnr = self.relbod * depcor + + # end #$OXBEN + self.bendox = -self.doben * self.vol + self.bnrbod = self.bodbnr * self.vol + + if self.LKFG == 1: + wind = 0.0 + + # calculate oxygen reaeration + if not (self.GQFG == 1 and self.GQALFG4 == 1): + self.korea = oxrea( + self.LKFG,wind,self.cforea,avvele,avdepe,self.tcginv, + self.REAMFG,self.reak,self.reakt,self.expred,self.exprev,self.len_, + self.delth,tw,self.delts,self.delt60,self.uunits) + + # calculate oxygen saturation level for current water + # temperature; satdo is expressed as mg oxygen per liter + self.satdo = 14.652 + tw * (-0.41022 + tw * (0.007991 - 0.7777e-4 * tw)) + + # adjust satdo to conform to prevalent atmospheric pressure + # conditions; cfpres is the ratio of site pressure to sea level pressure + self.satdo = self.cfpres * self.satdo + + if self.satdo < 0.0: + self.errors[0] += 1 + # warning - this occurs only when water temperature is very high - above + # about 66 c. usually means error in input gatmp (or tw if htrch is not being simulated). + self.satdo = 0.0 # reset saturation level + + # compute dissolved oxygen value after reaeration,and the reaeration flux + dorea = self.korea * (self.satdo - self.dox) + self.dox = self.dox + dorea + self.readox = dorea * self.vol + + #$BODDEC + '''calculate concentration of oxygen required to satisfy computed bod decay''' + self.bodox = (self.kbod20 * (self.tcbod**(tw -20.0))) * self.bod # bodox is expressed as mg oxygen/liter.ivl + if self.bodox > self.bod: + self.bodox = self.bod + + # adjust dissolved oxygen state variable to acount for oxygen lost to bod decay, and compute concentration flux + if self.bodox >= self.dox: + self.bodox = self.dox + self.dox = 0.0 + else: + self.dox = self.dox - self.bodox + + # adjust bod state variable to account for bod decayed + self.bod -= self.bodox + if self.bod < 0.0001: + self.bod = 0.0 + # end #$BODDEC + + self.boddox = -self.bodox * self.vol + self.decbod = -self.bodox * self.vol + + else: # there is too little water to warrant simulation of quality processes + self.bodox = 0.0 + + self.readox = 0.0 + self.boddox = 0.0 + self.bendox = 0.0 + self.decbod = 0.0 + self.bnrbod = 0.0 + self.snkbod = 0.0 + + self.totdox = self.readox + self.boddox + self.bendox + self.totbod = self.decbod + self.bnrbod + self.snkbod + + self.rdox = self.dox * self.vol + self.rbod = self.bod * self.vol + + self.svol = self.vol # svol is volume at start of time step, update for next time thru + + return + + #------------------------------------------------------------------- + # utility functions: + #------------------------------------------------------------------- + def adjust_dox(self, nitdox, denbod, phydox, zoodox, baldox): + # if dox exceeds user specified level of supersaturation, then release excess do to the atmosphere + + doxs = self.dox + + if self.dox > self.supsat * self.satdo: + self.dox = self.supsat * self.satdo + + self.readox = self.readox + (self.dox - doxs) * self.vol + self.totdox = self.readox + self.boddox + self.bendox \ + + nitdox + phydox + zoodox + baldox + + + # update dissolved totals and totals of nutrients + self.rdox = self.dox * self.vol + self.rbod = self.bod * self.vol + + return + + ''' + def update_totals(self, nitdox, denbod): + + self.totdox = self.readox + self.boddox + self.bendox + nitdox + self.totbod = self.decbod + self.bnrbod + self.snkbod + denbod + + return + ''' \ No newline at end of file diff --git a/HSP2/PLANK.py b/HSP2/PLANK.py index 332f409f..83de46d0 100644 --- a/HSP2/PLANK.py +++ b/HSP2/PLANK.py @@ -25,23 +25,28 @@ def plank(store, siminfo, uci, ts): ui = make_numba_dict(uci) # flags - table-type PLNK-FLAGS - PHYFG = ui['PHYFG'] - ZOOFG = ui['ZOOFG'] - BALFG = ui['BALFG'] - SDLTFG = ui['SDLTFG'] - AMRFG = ui['AMRFG'] - DECFG = ui['DECFG'] - NSFG = ui['NSFG'] - ZFOOD = ui['ZFOOD'] - BNPFG = ui['BNPFG'] - - HTFG = ui['HTFG'] - TAMFG = ui['NH3FG'] - PO4FG = ui['PO4FG'] + PHYFG = int(ui['PHYFG']) + ZOOFG = int(ui['ZOOFG']) + BALFG = int(ui['BALFG']) + SDLTFG = int(ui['SDLTFG']) + AMRFG = int(ui['AMRFG']) + DECFG = int(ui['DECFG']) + NSFG = int(ui['NSFG']) + ZFOOD = int(ui['ZFOOD']) + BNPFG = int(ui['BNPFG']) + + HTFG = int(ui['HTFG']) + TAMFG = int(ui['NH3FG']) + NO2FG = int(ui['NO2FG']) + PO4FG = int(ui['PO4FG']) + + ADNHFG = int(ui['ADNHFG']) + ADPOFG = int(ui['ADPOFG']) bpcntc = ui['BPCNTC'] cvbo = ui['CVBO'] cvbpc = ui['CVBPC'] + cvbpn = ui['CVBPN'] if ZOOFG == 1 and PHYFG == 0: pass # ERRMSG: error - zooplankton cannot be simulated without phytoplankton @@ -52,12 +57,12 @@ def plank(store, siminfo, uci, ts): if BALFG == 2: # user has selected multiple species with more complex kinetics # additional benthic algae flags - table-type BENAL-FLAG - numbal = ui['NUMBAL'] - BINVFG = ui['BINVFG'] - BFIXFG1 = ui['BFIXFG1'] - BFIXFG2 = ui['BFIXFG2'] - BFIXFG3 = ui['BFIXFG3'] - BFIXFG4 = ui['BFIXFG4'] + numbal = int(ui['NUMBAL']) + BINVFG = int(ui['BINVFG']) + BFIXFG1 = int(ui['BFIXFG1']) + BFIXFG2 = int(ui['BFIXFG2']) + BFIXFG3 = int(ui['BFIXFG3']) + BFIXFG4 = int(ui['BFIXFG4']) else: numbal = BALFG # single species or none @@ -79,6 +84,7 @@ def plank(store, siminfo, uci, ts): cvnrbo = nonref * cvbo cvbp = (31.0 * bpcntc) / (1200.0 * cvbpc) + cvbn = 14.0 * cvbpn * cvbp / 31.0 cvpb = 31.0 / (1000.0 * cvbp) cvbcl = 31.0 * ratclp / cvpb @@ -210,48 +216,21 @@ def plank(store, siminfo, uci, ts): orc = ui['ORC'] # atmospheric deposition flags - PLADFG = array(ui['PLADFG']) # six PLADFG1 to PLADFG6; table-type PLNK-AD-FLAGS - PLAFXM = zeros(4) - - for j in range(1, 4): - n = 2*(j - 1) + 1 - if PLADFG[n] > 0: # monthly flux must be read - PLAFXM[j] = array(ui[('PLAFXM',j)]) - if uunits == 1: # convert from lb/ac.day to mg.ft3/l.ft2.ivl - PLAFXM[j] *= 0.3677 * DELT60 / 24.0 - elif uunits == 2: # convert from kg/ha.day to mg.m3/l.m2.ivl - PLAFXM[j] *= 0.1 * DELT60 / 24.0 - # Note: get PLAFX array from monthly (above), constant, or time series - # PLAFX is dimension (simlen, 3) - # same with PLADCN - - # compute atmospheric deposition influx; [N, P, C] - pladdr = zeros(4) - pladwt = zeros(4) - pladep = zeros(4) - flxbal = zeros((4,numbal)) - - for i in range(1,4): - n = 2 * (i - 1) + 1 - pladdr[i] = SAREA * PLADFX[i] # dry deposition - pladwt[i] = PREC * SAREA * PLADCN[i] # wet deposition - pladep[i] = pladdr[i] + pladwt[i] - + PLADFG = zeros(7) + for j in range(1,7): + PLADFG[j] = ui['PLADFG(' + str(j) + ')'] + + # variable initialization: if PHYFG == 0: # initialize fluxes of inactive constituent rophyt = 0.0 ophyt[:] = 0.0 #nexits - phydox = 0.0 - phybod = 0.0 - phytam = 0.0 - phyno3 = 0.0 - phypo4 = 0.0 - phyorn = 0.0 - phyorp = 0.0 - phyorc = 0.0 + phydox = phybod = 0.0 + phytam = phyno3 = phypo4 = 0.0 + phyorn = phyorp = phyorc = 0.0 pyco2 = 0.0 - dthphy = 0.0 - grophy = 0.0 - totphy = 0.0 + dthphy = grophy = totphy = 0.0 + + ozoo = zeros(nexits) if ZOOFG == 1: # convert zoo to mg/l zoo *= zomass @@ -259,21 +238,16 @@ def plank(store, siminfo, uci, ts): # initialize fluxes of inactive constituent rozoo = 0.0 ozoo[:] = 0.0 #nexits - zoodox = 0.0 - zoobod = 0.0 - zootam = 0.0 - zoono3 = 0.0 - zoopo4 = 0.0 - zooorn = 0.0 - zooorp = 0.0 - zooorc = 0.0 + zoodox = zoobod = 0.0 + zootam = zoono3 = zoopo4 = 0.0 + zooorn = zooorp = zooorc = 0.0 zoophy = 0.0 zoco2 = 0.0 - grozoo = 0.0 - dthzoo = 0.0 - totzoo = 0.0 + grozoo = dthzoo = totzoo = 0.0 benal = zeros(numbal) + flxbal = zeros((4,5)) + if numbal == 1: # single species benal[1] = ui['BENAL'] # points to table-type plnk-init above for rvals elif numbal >= 2: # multiple species - table-type benal-init @@ -281,32 +255,30 @@ def plank(store, siminfo, uci, ts): str_ba = 'BENAL' + str(n+1) benal[n] = ui[str_ba] # how to get multiples??? else: # no benthic algae simulated - baldox = 0.0 - balbod = 0.0 - baltam = 0.0 - balno3 = 0.0 - balpo4 = 0.0 - balorn = 0.0 - balorp = 0.0 - balorc = 0.0 + baldox = balbod = 0.0 + baltam = balno3 = balpo4 = 0.0 + balorn = balorp = balorc = 0.0 baco2 = 0.0 flxbal[1:4,1:5] = 0.0 # compute derived quantities phycla = phyto * cvbcl - balcla = zeros(numbal) + balcla = zeros(4) for i in range(numbal): balcla[i] = benal[i] * cvbcl + lsnh4 = lspo4 = zeros(4) + if vol > 0.0: # compute initial summary concentrations for i in range(1, 4): - lsnh4[i] = rsnh4[i] / vol - lspo4[i] = rspo4[i] / vol - - torn,torp,torc,potbod,tn,tp = \ - pksums (phyfg,zoofg,tamfg,no2fg,po4fg,adnhfg,adpofg, + #lsnh4[i] = rsnh4[i] / vol #LTI-need to address this NUTRX var. + #lspo4[i] = rspo4[i] / vol #LTI-need to address this NUTRX var. + pass + + (torn,torp,torc,potbod,tn,tp) = \ + pksums (PHYFG,ZOOFG,TAMFG,NO2FG,PO4FG,ADNHFG,ADPOFG, cvbn,cvbp,cvbc,cvbo,cvnrbo,phyto,zoo,orn,orp,orc,no3,tam,no2,lsnh4[1],lsnh4[2], - lsnh4[3],po4,lspo4[1],lspo4[2],lspo4[3],bod,torn,torp,torc,potbod,tn,tp) + lsnh4[3],po4,lspo4[1],lspo4[2],lspo4[3],bod) #,torn,torp,torc,potbod,tn,tp) else: # undefined summary concentrations torn = -1.0e30 torp = -1.0e30 @@ -318,8 +290,34 @@ def plank(store, siminfo, uci, ts): return errors, ERRMSGS @jit(nopython = True) - def plank(dox, bod, iphyto, izoo, iorn, iorp, iorc, tw, wash, solrad, prec, sarea, advData): + def plank_run(dox, bod, iphyto, izoo, iorn, iorp, iorc, rsnh4, rspo4, tw, wash, solrad, prec, sarea, advData): '''Simulate behavior of plankton populations and associated reactions''' + + # update atmospheric deposition: + PLAFXM = zeros(4) + + for j in range(1, 4): + n = 2*(j - 1) + 1 + if PLADFG[n] > 0: # monthly flux must be read + PLAFXM[j] = array(ui[('PLAFXM',j)]) #LTI-needs testing + if uunits == 1: # convert from lb/ac.day to mg.ft3/l.ft2.ivl + PLAFXM[j] *= 0.3677 * DELT60 / 24.0 + elif uunits == 2: # convert from kg/ha.day to mg.m3/l.m2.ivl + PLAFXM[j] *= 0.1 * DELT60 / 24.0 + # Note: get PLAFX array from monthly (above), constant, or time series + # PLAFX is dimension (simlen, 3) + # same with PLADCN + + # compute atmospheric deposition influx; [N, P, C] + pladdr = pladwt = pladep = zeros(4) + + for i in range(1,4): + n = 2 * (i - 1) + 1 + pladdr[i] = SAREA * PLADFX[i] # dry deposition + pladwt[i] = PREC * SAREA * PLADCN[i] # wet deposition + pladep[i] = pladdr[i] + pladwt[i] + + if PHYFG == 1: # phytoplankton simulated # advecvt phytoplankton phyto, rophyt,ophyt = advplk(iphyto,vols,srovol,vol,erovol,sovol, eovol,nexits,oref,mxstay,seed,delts,phyto,rophyt,ophyt) @@ -616,15 +614,15 @@ def plank(dox, bod, iphyto, izoo, iorn, iorp, iorc, tw, wash, solrad, prec, sare adnhfg,adpofg,cvbn,cvbp,cvbc,cvbo,cvnrbo,rophyt,rozoo,roorn,roorp, roorc,rono3, rotam,rono2,rosnh4[1],rosnh4[2],rosnh4[3],ropo4,rospo4[1],rospo4[2],rospo4[3],robod, rotorn,rotorp,rotorc,dumval,rototn,rototp) - + if nexits > 1: # outflows by exit for i in range(nexits): otorn[i],otorp[i],otorc[i],dumval,ototn[i], ototp[i] \ = pksums(phyfg,zoofg, tamfg,no2fg,po4fg,adnhfg,adpofg,cvbn,cvbp,cvbc,cvbo,cvnrbo,ophyt[i],ozoo[i], oorn[i],oorp[i],oorc[i],ono3[i],otam[i],ono2[i], osnh4[i,1],osnh4[i,2], - osnh4[i,3],opo4[i], ospo4[i,1],ospo4[i,2],ospo4[i,3],obod[i],otorn[i], - otorp[i],otorc[i],dumval,ototn[i],ototp[i]) + osnh4[i,3],opo4[i], ospo4[i,1],ospo4[i,2],ospo4[i,3],obod[i]) + #otorn[i],otorp[i],otorc[i],dumval,ototn[i],ototp[i]) # cbrb added call to ammion to redistribute tam after algal influence nh3,nh4 = ammion(tw, phval, tam, nh3, nh4) diff --git a/HSP2/PLANK_Class.py b/HSP2/PLANK_Class.py new file mode 100644 index 00000000..7a0fc1aa --- /dev/null +++ b/HSP2/PLANK_Class.py @@ -0,0 +1,2040 @@ +import numpy as np +from numpy import zeros, array +from math import log, exp +import numba as nb +from numba.experimental import jitclass + +from HSP2.ADCALC import advect +from HSP2.OXRX_Class import OXRX_Class +from HSP2.NUTRX_Class import NUTRX_Class +from HSP2.RQUTIL import sink, decbal +from HSP2.utilities import make_numba_dict, initm + +spec = [ + ('OXRX', OXRX_Class.class_type.instance_type), + ('NUTRX', NUTRX_Class.class_type.instance_type), + ('aldh', nb.float64), + ('aldl', nb.float64), + ('alnpr', nb.float64), + ('alr20', nb.float64), + ('AMRFG', nb.int32), + ('baco2', nb.float64), + ('balbod', nb.float64), + ('balcla', nb.float64[:]), + ('baldep', nb.float64), + ('baldox', nb.float64), + ('BALFG', nb.int32), + ('ballit', nb.float64), + ('balno3', nb.float64), + ('balorc', nb.float64), + ('balorn', nb.float64), + ('balorp', nb.float64), + ('balpo4', nb.float64), + ('balr20', nb.float64[:]), + ('baltam', nb.float64), + ('balvel', nb.float64), + ('benal', nb.float64[:]), + ('BFIXFG', nb.int32[:]), + ('binv', nb.float64), + ('BINVFG', nb.int32), + ('binvm', nb.float64), + ('BNPFG', nb.int32), + ('bpcntc', nb.float64), + ('campr', nb.float64), + ('cbnrbo', nb.float64), + ('cfbalg', nb.float64), + ('cfbalr', nb.float64), + ('cflit', nb.float64), + ('cfsaex', nb.float64), + ('cktrb1', nb.float64), + ('cktrb2', nb.float64), + ('claldh', nb.float64), + ('cmingr', nb.float64), + ('cmmbi', nb.float64), + ('cmmd1', nb.float64[:]), + ('cmmd2', nb.float64[:]), + ('cmmlt', nb.float64), + ('cmmn', nb.float64), + ('cmmnb', nb.float64[:]), + ('cmmnp', nb.float64), + ('cmmp', nb.float64), + ('cmmpb', nb.float64[:]), + ('cmmv', nb.float64), + ('co2', nb.float64), + ('conv', nb.float64), + ('cremvl', nb.float64), + ('cslit', nb.float64[:]), + ('cslof1', nb.float64[:]), + ('cslof2', nb.float64[:]), + ('ctrbq1', nb.float64), + ('ctrbq2', nb.float64), + ('cvbc', nb.float64), + ('cvbcl', nb.float64), + ('cvbn', nb.float64), + ('cvbo', nb.float64), + ('cvbp', nb.float64), + ('cvbpc', nb.float64), + ('cvbpn', nb.float64), + ('cvnrbo', nb.float64), + ('cvpb', nb.float64), + ('DECFG', nb.int32), + ('delt', nb.float64), + ('delt60', nb.float64), + ('delts', nb.float64), + ('dthbal', nb.float64[:]), + ('dthphy', nb.float64), + ('dthtba', nb.float64), + ('dthzoo', nb.float64), + ('errors', nb.int64[:]), + ('extb', nb.float64), + ('flxbal', nb.float64[:,:]), + ('fravl', nb.float64), + ('frrif', nb.float64), + ('grobal', nb.float64[:]), + ('grophy', nb.float64), + ('grores', nb.float64[:]), + ('grotba', nb.float64), + ('grozoo', nb.float64), + ('HTFG', nb.int32), + ('iorc', nb.float64), + ('iorn', nb.float64), + ('iorp', nb.float64), + ('iphyto', nb.float64), + ('itorc', nb.float64), + ('itorn', nb.float64), + ('itorp', nb.float64), + ('itotn', nb.float64), + ('itotp', nb.float64), + ('izoo', nb.float64), + ('limbal', nb.int32[:]), + ('limphy', nb.int32), + ('litsed', nb.float64), + ('lmingr', nb.float64), + ('lsnh4', nb.float64[:]), + ('lspo4', nb.float64[:]), + ('malgr', nb.float64), + ('mbal', nb.float64), + ('mbalgr', nb.float64[:]), + ('minbal', nb.float64), + ('mxstay', nb.float64), + ('mzoeat', nb.float64), + ('naldh', nb.float64), + ('nexits', nb.int32), + ('nmaxfx', nb.float64), + ('nminc', nb.float64), + ('nmingr', nb.float64), + ('nonref', nb.float64), + ('NSFG', nb.int32), + ('numbal', nb.int32), + ('oorc', nb.float64[:]), + ('oorn', nb.float64[:]), + ('oorp', nb.float64[:]), + ('ophyt', nb.float64[:]), + ('orc', nb.float64), + ('oref', nb.float64), + ('orn', nb.float64), + ('orp', nb.float64), + ('otorc', nb.float64[:]), + ('otorn', nb.float64[:]), + ('otorp', nb.float64[:]), + ('ototn', nb.float64[:]), + ('ototp', nb.float64[:]), + ('oxald', nb.float64), + ('oxzd', nb.float64), + ('ozoo', nb.float64[:]), + ('paldh', nb.float64), + ('paradf', nb.float64), + ('PHFG', nb.int32), + ('phybod', nb.float64), + ('phycla', nb.float64), + ('phydox', nb.float64), + ('PHYFG', nb.int32), + ('phylit', nb.float64), + ('phyno3', nb.float64), + ('phyorc', nb.float64), + ('phyorn', nb.float64), + ('phyorp', nb.float64), + ('phypo4', nb.float64), + ('physet', nb.float64), + ('phytam', nb.float64), + ('phyto', nb.float64), + ('PLADCN', nb.float64[:]), + ('PLADFG', nb.int32[:]), + ('PLADFX', nb.float64[:]), + ('pmingr', nb.float64), + ('potbod', nb.float64), + ('pyco2', nb.float64), + ('ratclp', nb.float64), + ('refr', nb.float64), + ('refset', nb.float64), + ('rifcq1', nb.float64), + ('rifcq2', nb.float64), + ('rifcq3', nb.float64), + ('rifdf', nb.float64[:]), + ('rifvf', nb.float64[:]), + ('roorc', nb.float64), + ('roorn', nb.float64), + ('roorp', nb.float64), + ('rophyt', nb.float64), + ('rotorc', nb.float64), + ('rotorn', nb.float64), + ('rotorp', nb.float64), + ('rototn', nb.float64), + ('rototp', nb.float64), + ('rozoo', nb.float64), + ('SDLTFG', nb.int32), + ('seed', nb.float64), + ('simlen', nb.int32), + ('snkphy', nb.float64), + ('snkorc', nb.float64), + ('snkorn', nb.float64), + ('snkorp', nb.float64), + ('svol', nb.float64), + ('talgrh', nb.float64), + ('talgrl', nb.float64), + ('talgrm', nb.float64), + ('tbenal', nb.float64[:]), + ('tcbalg', nb.float64[:]), + ('tcbalr', nb.float64[:]), + ('tcgraz', nb.float64), + ('tczfil', nb.float64), + ('tczres', nb.float64), + ('tn', nb.float64), + ('torc', nb.float64), + ('torn', nb.float64), + ('torp', nb.float64), + ('totbod', nb.float64), + ('totdox', nb.float64), + ('totno3', nb.float64), + ('totorc', nb.float64), + ('totorn', nb.float64), + ('totorp', nb.float64), + ('totphy', nb.float64), + ('totpo4', nb.float64), + ('tottam', nb.float64), + ('tottba', nb.float64), + ('totzoo', nb.float64), + ('tp', nb.float64), + ('uunits', nb.int32), + ('vol', nb.float64), + ('zd', nb.float64), + ('zexdel', nb.float64), + ('zfil20', nb.float64), + ('ZFOOD', nb.int32), + ('zoco2', nb.float64), + ('zomass', nb.float64), + ('zoo', nb.float64), + ('zoobod', nb.float64), + ('zoodox', nb.float64), + ('ZOOFG', nb.int32), + ('zoono3', nb.float64), + ('zooorc', nb.float64), + ('zooorn', nb.float64), + ('zooorp', nb.float64), + ('zoophy', nb.float64), + ('zoopo4', nb.float64), + ('zootam', nb.float64), + ('zres20', nb.float64), +] + +@jitclass(spec) +class PLANK_Class: + + #------------------------------------------------------------------- + # class initialization: + #------------------------------------------------------------------- + def __init__(self, siminfo, nexits, vol, ui_rq, ui, ts, OXRX, NUTRX): + + ''' Initialize instance variables for lower food web simulation ''' + + self.errors = zeros(int(ui['errlen']), dtype=np.int64) + + ''' + self.limit = array(8, type=nb.char) #np.chararray(8, itemsize=4) + self.limit[1] = 'LIT' + self.limit[2] = 'NON' + self.limit[3] = 'TEM' + self.limit[4] = 'NIT' + self.limit[5] = 'PO4' + self.limit[6] = 'NONE' + self.limit[7] = 'WAT' + ''' + + self.delt = siminfo['delt'] + 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 + + # flags - table-type PLNK-FLAGS + self.PHYFG = int(ui['PHYFG']) + self.ZOOFG = int(ui['ZOOFG']) + self.BALFG = int(ui['BALFG']) + self.SDLTFG = int(ui['SDLTFG']) + self.AMRFG = int(ui['AMRFG']) + self.DECFG = int(ui['DECFG']) + self.NSFG = int(ui['NSFG']) + self.ZFOOD = int(ui['ZFOOD']) + self.BNPFG = int(ui['BNPFG']) + + self.HTFG = int(ui_rq['HTFG']) + self.PHFG = int(ui_rq['PHFG']) + + self.bpcntc = NUTRX.bpcntc + self.cvbo = NUTRX.cvbo + self.cvbpc = NUTRX.cvbpc + self.cvbpn = NUTRX.cvbpn + + if self.ZOOFG == 1 and self.PHYFG == 0: + self.errors[0] += 1 + # ERRMSG: error - zooplankton cannot be simulated without phytoplankton + if self.NSFG == 1 and NUTRX.TAMFG == 0: + self.errors[1] += 1 + # ERRMSG: error - ammonia cannot be included in n supply if it is not + if NUTRX.PO4FG == 0: + self.errors[2] += 1 + # ERRMSG: error - phosphate must be simulated if plankton are being + + self.numbal = 0 + self.BFIXFG = zeros(5, dtype=np.int32) + + if self.BALFG == 2: # user has selected multiple species with more complex kinetics + # additional benthic algae flags - table-type BENAL-FLAG + self.numbal = int(ui['NUMBAL']) + self.BINVFG = int(ui['BINVFG']) + + for i in range(1,5): + self.BFIXFG[i] = int(ui['BFIXFG' + str(i)]) + else: + 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 + self.cfsaex = ui['CFSAEX'] + + # table-type plnk-parm1 + self.ratclp = ui['RATCLP'] + self.nonref = ui['NONREF'] + self.litsed = ui['LITSED'] + self.alnpr = ui['ALNPR'] + self.extb = ui['EXTB'] + self.malgr = ui['MALGR'] * self.delt60 + self.paradf = ui['PARADF'] + self.refr = 1.0 - self.nonref # define fraction of biomass which is refractory material + + # compute derived conversion factors + self.cvbc = self.bpcntc / 100.0 + self.cvnrbo = self.nonref * self.cvbo + + self.cvbp = (31.0 * self.bpcntc) / (1200.0 * self.cvbpc) + self.cvbn = 14.0 * self.cvbpn * self.cvbp / 31.0 + self.cvpb = 31.0 / (1000.0 * self.cvbp) + self.cvbcl = 31.0 * self.ratclp / self.cvpb + + # table-type plnk-parm2 + self.cmmlt = ui['CMMLT'] + self.cmmn = ui['CMMN'] + self.cmmnp = ui['CMMNP'] + self.cmmp = ui['CMMP'] + self.talgrh = ui['TALGRH'] + self.talgrl = ui['TALGRL'] + self.talgrm = ui['TALGRM'] + + # table-type plnk-parm3 + self.alr20 = ui['ALR20'] * delt60 # convert rates from 1/hr to 1/ivl + self.aldh = ui['ALDH'] * delt60 + self.aldl = ui['ALDL'] * delt60 + self.oxald = ui['OXALD'] * delt60 + self.naldh = ui['NALDH'] * delt60 + self.paldh = ui['PALDH'] + + # table-type plnk-parm4 + self.nmingr = ui['NMINGR'] + self.pmingr = ui['PMINGR'] + self.cmingr = ui['CMINGR'] + self.lmingr = ui['LMINGR'] + self.nminc = ui['NMINC'] + + # phytoplankton-specific parms - table-type phyto-parm + # this table must always be input so that REFSET is read + self.seed = ui['SEED'] + self.mxstay = ui['MXSTAY'] + self.oref = ui['OREF'] + self.claldh = ui['CLALDH'] + self.physet = ui['PHYSET'] * delt60 # change settling rates to units of 1/ivl + self.refset = ui['REFSET'] * delt60 # change settling rates to units of 1/ivl + + if self.PHYFG == 1 and self.ZOOFG == 1: # zooplankton-specific parameters + # table-type zoo-parm1 + self.mzoeat = ui['MZOEAT'] * delt60 # convert rates from 1/hr to 1/ivl + self.zfil20 = ui['ZFIL20'] * delt60 # convert rates from 1/hr to 1/ivl + self.zres20 = ui['ZRES20'] * delt60 # convert rates from 1/hr to 1/ivl + self.zd = ui['ZD'] * delt60 # convert rates from 1/hr to 1/ivl + self.oxzd = ui['OXZD'] * delt60 # convert rates from 1/hr to 1/ivl + # table-type zoo-parm2 + self.tczfil = ui['TCZFIL'] + self.tczres = ui['TCZRES'] + self.zexdel = ui['ZEXDEL'] + self.zomass = ui['ZOMASS'] + + if self.BALFG >= 1: # benthic algae-specific parms; table-type benal-parm + self.mbal = ui['MBAL'] / self.cvpb # convert maximum benthic algae to micromoles of phosphorus + self.cfbalr = ui['CFBALR'] + self.cfbalg = ui['CFBALG'] + self.minbal = ui['MINBAL'] / self.cvpb # convert maximum benthic algae to micromoles of phosphorus + self.campr = ui['CAMPR'] + self.fravl = ui['FRAVL'] + self.nmaxfx = ui['NMAXFX'] + + self.mbalgr = zeros(self.numbal) + self.tcbalg = zeros(self.numbal) + self.cmmnb = zeros(self.numbal) + self.cmmpb = zeros(self.numbal) + self.cmmd1 = zeros(self.numbal) + self.cmmd2 = zeros(self.numbal) + self.cslit = zeros(self.numbal) + + self.balr20 = zeros(self.numbal) + self.tcbalr = zeros(self.numbal) + self.cslof1 = zeros(self.numbal) + self.cslof2 = zeros(self.numbal) + self.grores = zeros(self.numbal) + + if self.BALFG == 2: # user has selected multiple species with more complex kinetics + for i in range(self.numbal): + # species-specific growth parms - table type benal-grow + self.mbalgr[i] = ui['MBALGR'] * self.delt60 + self.tcbalg[i] = ui['TCBALG'] + self.cmmnb[i] = ui['CMMNB'] + self.cmmpb[i] = ui['CMMPB'] + self.cmmd1[i] = ui['CMMD1'] + self.cmmd2[i] = ui['CMMD2'] + self.cslit[i] = ui['CSLIT'] + # species-specific resp and scour parms - table type benal-resscr + self.balr20[i] = ui['BALR20'] * self.delt60 + self.tcbalr[i] = ui['TCBALR'] + self.cslof1[i] = ui['CSLOF1'] * self.delt60 + self.cslof2[i] = ui['CSLOF2'] + self.grores[i] = ui['GRORES'] + + # grazing and disturbance parms - table-type benal-graze + self.cremvl = ui['CREMVL'] + self.cmmbi = ui['CMMBI'] + self.binv = ui['BINV'] + self.tcgraz = ui['TCGRAZ'] + + hrpyr = 8760.0 #constant + self.cremvl = (self.cremvl / self.cvpb) / hrpyr * self.delt60 + + if self.SDLTFG == 2: # turbidity regression parms - table-type benal-light + self.ctrbq1 = ui['CTRBQ1'] + self.ctrbq2 = ui['CTRBQ2'] + self.cktrb1 = ui['CKTRB1'] + self.cktrb2 = ui['CKTRB2'] + + if self.BINVFG == 3: # monthly benthic invertebrate density - table-type mon-binv + self.binvm = ui['BINVM'] + + # table-type benal-riff1 + self.frrif = ui['FRRIF'] + self.cmmv = ui['CMMV'] + self.rifcq1 = ui['RIFCQ1'] + self.rifcq2 = ui['RIFCQ2'] + self.rifcq3 = ui['RIFCQ3'] + + # table-type benal-riff2 + self.rifvf = zeros(5); self.rifdf = zeros(5) + + for i in range(1,5): + self.rifvf[i] = ui['RIFVF' + str(i)] + self.rifdf[i] = ui['RIFDF' + str(i)] + + # table-type plnk-init + self.phyto = ui['PHYTO'] + self.zoo = ui['ZOO'] + #benal = ui['BENAL'] + self.orn = ui['ORN'] + self.orp = ui['ORP'] + self.orc = ui['ORC'] + + # atmospheric deposition flags + self.PLADFG = zeros(7, dtype=np.int32) + for j in range(1,7): + self.PLADFG[j] = int(ui['PLADFG(' + str(j) + ')']) + + # variable initialization: + if self.PHYFG == 0: # initialize fluxes of inactive constituent + self.rophyt = 0.0 + self.ophyt[:] = 0.0 #nexits + self.phydox = self.phybod = 0.0 + self.phytam = 0.0; self.phyno3 = 0.0; self.phypo4 = 0.0 + self.phyorn = 0.0; self.phyorp = 0.0; self.phyorc = 0.0 + self.pyco2 = 0.0 + self.dthphy = 0.0; self.grophy = 0.0; self.totphy = 0.0 + + self.rozoo = 0.0 + self.ozoo = zeros(nexits) + + if self.ZOOFG == 1: # convert zoo to mg/l + self.zoo *= self.zomass + + else: # zooplankton not simulated, but use default values + # initialize fluxes of inactive constituent + self.rozoo = 0.0 + self.ozoo[:] = 0.0 #nexits + self.zoodox = 0.0; self.zoobod = 0.0 + self.zootam = 0.0; self.zoono3 = 0.0; self.zoopo4 = 0.0 + self.zooorn = 0.0; self.zooorp = 0.0; self.zooorc = 0.0 + self.zoophy = 0.0 + self.zoco2 = 0.0 + self.grozoo = 0.0; self.dthzoo = 0.0; self.totzoo = 0.0 + + # benthic algae initialization: + self.benal = zeros(self.numbal) + self.flxbal = zeros((4,5)) + + if self.numbal == 1: # single species + self.benal[0] = ui['BENAL'] # points to table-type plnk-init above for rvals + elif self.numbal >= 2: # multiple species - table-type benal-init + for n in range(self.numbal): + self.benal[n] = ui['BENAL' + str(n+1)] + else: # no benthic algae simulated + self.baldox = 0.0; self.balbod = 0.0 + self.baltam = 0.0; self.balno3 = 0.0; self.balpo4 = 0.0 + self.balorn = 0.0; self.balorp = 0.0; self.balorc = 0.0 + self.baco2 = 0.0 + + # compute derived quantities + self.phycla = self.phyto * self.cvbcl + + self.balcla = zeros(self.numbal) + for i in range(self.numbal): + self.balcla[i] = self.benal[i] * self.cvbcl + + self.lsnh4 = zeros(4) + self.lspo4 = zeros(4) + + if self.vol > 0.0: # compute initial summary concentrations + for i in range(1, 4): + self.lsnh4[i] = NUTRX.rsnh4[i] / self.vol + self.lspo4[i] = NUTRX.rspo4[i] / self.vol + + # calculate summary concentrations: + (self.torn, self.torp, self.torc, self.potbod, self.tn, self.tp) \ + = self.pksums(NUTRX,self.phyto,self.zoo,self.orn,self.orp,self.orc, + NUTRX.no3,NUTRX.tam,NUTRX.no2,self.lsnh4,NUTRX.po4,self.lspo4,OXRX.bod) + + return + + def simulate(self, tw, phval, co2, tss, OXRX, NUTRX, iphyto, izoo, iorn, iorp, iorc, + wash, solrad, prec, sarea, avdepe, avvele, depcor, ro, advData): + + ''' ''' + + # hydraulics: + (nexits, vols, vol, srovol, erovol, sovol, eovol) = advData + self.vol = vol + + # initialize temp. vars for DO/BOD and NUTRX states: + po4 = NUTRX.po4 + no3 = NUTRX.no3 + no2 = NUTRX.no2 + tam = NUTRX.tam + + dox = OXRX.dox + bod = OXRX.bod + + self.co2 = co2 + + # inflows: convert from [mass/ivld] to [conc.*vol/ivld] + self.iphyto = iphyto / self.conv + self.izoo = izoo / self.conv + self.iorn = iorn / self.conv + self.iorp = iorp / self.conv + self.iorc = iorc / self.conv + + # update atmospheric deposition (TO-DO! - best approach is likely to handle in RQUAL and pass in PLAD* values for current step): + self.PLADFX = zeros(6) + self.PLADCN = zeros(6) + + ''' + for j in range(1, 4): + n = 2*(j - 1) + 1 + + if self.PLADFG[n] > 0: # monthly flux must be read + #self.PLAFXM[j] = array(ui[('PLAFXM',j)]) + + if uunits == 1: # convert from lb/ac.day to mg.ft3/l.ft2.ivl + self.PLAFXM[j] *= 0.3677 * self.delt60 / 24.0 + elif uunits == 2: # convert from kg/ha.day to mg.m3/l.m2.ivl + self.PLAFXM[j] *= 0.1 * self.delt60 / 24.0 + + # Note: get PLAFX array from monthly (above), constant, or time series + # PLAFX is dimension (simlen, 3) + # same with PLADCN + ''' + + pladdr = zeros(4); pladwt = zeros(4); pladep = zeros(4) + + for i in range(1,4): + n = 2 * (i - 1) + 1 + pladdr[i] = sarea * self.PLADFX[i] # dry deposition + pladwt[i] = prec * sarea * self.PLADCN[i] # wet deposition + pladep[i] = pladdr[i] + pladwt[i] + + + #----------------------------------------------------------- + # advection: + #----------------------------------------------------------- + if self.PHYFG == 1: + + # advect phytoplankton + (self.phyto, self.rophyt, self.ophyt) \ + = self.advplk(self.iphyto,self.svol,srovol,self.vol,erovol,sovol,eovol,nexits, + self.oref,self.mxstay,self.seed,self.delts,self.phyto) + + (self.phyto, snkphy) = sink(self.vol,avdepe,self.physet,self.phyto) + self.snkphy = -snkphy + + if self.ZOOFG == 1: # zooplankton on; advect zooplankton + (self.zoo, self.rozoo, self.ozoo) \ + = self.advplk(self.izoo,self.svol,srovol,vol,erovol,sovol,eovol,nexits, + self.oref,self.mxstay,self.seed,self.delts,self.zoo) + + # advect organic nitrogen + inorn = self.iorn + pladep[1] + (self.orn,self.roorn,self.oorn) = advect (inorn, self.orn, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + (self.orn, snkorn) = sink(self.vol, avdepe, self.refset, self.orn) + self.snkorn = -snkorn + + # advect organic phosphorus + inorp = self.iorp + pladep[2] + (self.orp, self.roorp, self.oorp) = advect(inorp, self.orp, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + (self.orp, snkorp) = sink(self.vol, avdepe, self.refset, self.orp) + self.snkorp = -snkorp + + # advect total organic carbon + inorc = self.iorc + pladep[3] + (self.orc, self.roorc, self.oorc) = advect (inorc, self.orc, nexits, self.svol, self.vol, srovol, erovol, sovol, eovol) + (self.orc, snkorc) = sink(self.vol, avdepe, self.refset, self.orc) + self.snkorc = -snkorc + + if avdepe > 0.17: # enough water to warrant computation of water quality reactions + + if self.BALFG > 0: + if self.frrif < 1.0: + # make adjustments to average water velocity and depth for the + # portion of the reach that consists of riffle areas. + if ro < self.rifcq1: # below first cutoff flow + i= 1 + elif ro < self.rifcq2: # below second cutoff flow + i= 2 + elif ro < self.rifcq3: # below third cutoff flow + i= 3 + else: # above third cutoff flow + i= 4 + + # calculate the adjusted velocity and depth for riffle sections + self.balvel = self.rifvf[i] * avvele + self.baldep = self.rifdf[i] * avdepe + else: # use full depth and velocity + self.balvel = avvele + self.baldep = avdepe + + # calculate solar radiation absorbed; solrad is the solar radiation at gage, + # corrected for location of reach; 0.97 accounts for surface reflection + # (assumed 3 per cent); cfsaex is the ratio of radiation incident to water + # surface to gage radiation values (accounts for regional differences, shading + # of water surface, etc); inlit is a measure of light intensity immediately below + # surface of reach/res and is expressed as ly/min, adjusted for fraction that is + # photosynthetically active. + + inlit = 0.97 * self.cfsaex * solrad / self.delt * self.paradf + extsed = 0.0 + + if self.SDLTFG == 1: # estimate contribution of sediment to light extinction + extsed = self.litsed * tss + elif self.SDLTFG == 2: # equations from dssamt for estimating the extinction coefficient based on discharge and turbidity + # estimate turbidity based on linear regression on flow + turb = self.ctrbq1 * ro**self.ctrbq2 + + # estimate the portion of the extinction coefficient due to + # sediment based upon a system-wide regression of total + # extinction to turbidity, and then subtracting the + # base extinction coefficient + extsed = (self.cktrb1 * turb**self.cktrb2) - self.extb + if extsed < 0.0: # no effective sediment shading + extsed = 0.0 + else: # sediment light extinction not considered + extsed = 0.0 + + # calculate contribution of phytoplankton to light extinction (self-shading) + extcla = 0.00452 * self.phyto * self.cvbcl + + # calculate light available for algal growth, litrch only called here + (self.phylit, self.ballit, self.cflit) = self.litrch (inlit,self.extb,extcla,extsed,avdepe,self.baldep,self.PHYFG,self.BALFG) + + #----------------------------------------------------------- + # 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, + 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, \ + self.cflit,self.alr20,self.cvbpn,self.PHFG,self.DECFG,self.cvbpc,self.paldh, \ + self.naldh,self.claldh,self.aldl,self.aldh,NUTRX.anaer,self.oxald,self.alnpr, \ + self.cvbo,self.refr,self.cvnrbo,self.cvpb,self.cvbcl,self.co2, \ + self.nmingr,self.pmingr,self.cmingr,self.lmingr,self.nminc, \ + po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.phyto) + + # compute associated fluxes + self.phydox = dophy * self.vol + self.phybod = bodphy * self.vol + self.phytam = tamphy * self.vol + self.phyno3 = no3phy * self.vol + self.phypo4 = po4phy * self.vol + self.dthphy = -phdth * self.vol + self.grophy = phgro * self.vol + self.phyorn = ornphy * self.vol + self.phyorp = orpphy * self.vol + self.phyorc = orcphy * self.vol + + #----------------------------------------------------------- + # 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) \ + = 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, \ + dox,bod,self.zoo,self.orn,self.orp,self.orc,tam,no3,po4) + + # compute associated fluxes + self.zoodox = -dozoo * self.vol + self.zoobod = bodzoo * self.vol + + if NUTRX.TAMFG != 0: # ammonia on, so nitrogen excretion goes to ammonia + self.zootam = nitzoo * self.vol + self.zoono3 = 0.0 + else: # ammonia off, so nitrogen excretion goes to nitrate + self.zoono3 = nitzoo * self.vol + self.zootam = 0.0 + + self.zoopo4 = po4zoo * self.vol + self.zoophy = -zeat * self.vol + self.zooorn = zorn * self.vol + self.zooorp = zorp * self.vol + self.zooorc = zorc * self.vol + self.grozoo = zgro * self.vol + self.dthzoo = -zdth * self.vol + self.totzoo = self.grozoo + self.dthzoo + + # update phytoplankton state variable to account for zooplankton predation + self.phyto = self.phyto - zeat + + # convert phytoplankton expressed as mg biomass/l to chlorophyll a expressed as ug/l + self.phycla = self.phyto * self.cvbcl + + self.totphy = self.snkphy + self.zoophy + self.dthphy + self.grophy + + #----------------------------------------------------------- + # benthic algae growth & respiration + #----------------------------------------------------------- + self.limbal = zeros(self.numbal, dtype=np.int32) + + if self.BALFG > 0: + bgro = zeros(self.numbal) + bdth = zeros(self.numbal) + + if self.BALFG == 1: # simulate benthic algae + (po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.benal[0],self.limbal[0],self.baco2,self.balcla[0], + dobalg,bodbal,tambal,no3bal,po4bal,bgro[0],bdth[0],ornbal,orpbal,orcbal) \ + = self.balrx(self.ballit,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, \ + self.cflit,self.alr20,self.cvbpn,self.PHFG,self.DECFG,self.cvbpc,self.paldh, \ + self.naldh,self.aldl,self.aldh,NUTRX.anaer,self.oxald,self.cfbalg,self.cfbalr, \ + self.alnpr,self.cvbo,self.refr,self.cvnrbo,self.cvpb,self.mbal,depcor, \ + self.cvbcl,self.co2,self.nmingr,self.pmingr,self.cmingr,self.lmingr,self.nminc, \ + po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.benal[0]) + + elif self.BALFG == 2: # simulate enhanced benthic algae equations from dssamt (!) + # then perform reactions, balrx2 only called here + (po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.benal,self.limbal,self.baco2,self.balcla, + dobalg,bodbal,tambal,no3bal,po4bal,bgro,bdth,ornbal,orpbal,orcbal) \ + = self.balrx2 (self.ballit,tw,NUTRX.TAMFG,self.NSFG,self.delt60,self.cvbpn,self.PHFG,self.DECFG, \ + self.cvbpc,self.alnpr,self.cvbo,self.refr,self.cvnrbo,self.cvpb,depcor, \ + self.cvbcl,co2,self.numbal,self.mbalgr,self.cmmpb,self.cmmnb, \ + self.balr20,self.tcbalg,self.balvel,self.cmmv,self.BFIXFG,self.cslit,self.cmmd1, \ + self.cmmd2,self.tcbalr,self.frrif,self.cremvl,self.cmmbi,self.binv,self.tcgraz, \ + self.cslof1,self.cslof2,self.minbal,self.fravl,self.BNPFG,self.campr,self.nmingr, \ + self.pmingr,self.cmingr,self.lmingr,self.nminc,self.nmaxfx,self.grores, \ + po4,no3,tam,dox,self.orn,self.orp,self.orc,bod,self.benal) + + #compute associated fluxes + self.baldox = dobalg * self.vol + self.balbod = bodbal * self.vol + self.baltam = tambal * self.vol + self.balno3 = no3bal * self.vol + self.balpo4 = po4bal * self.vol + self.balorn = ornbal * self.vol + self.balorp = orpbal * self.vol + self.balorc = orcbal * self.vol + + self.grobal = zeros(self.numbal) + self.dthbal = zeros(self.numbal) + + for i in range(self.numbal): + self.grobal[i] = bgro[i] + self.dthbal[i] = -bdth[i] + + else: # not enough water in reach/res to warrant simulation of quality processes + self.phyorn = 0.0 + self.balorn = 0.0 + self.zooorn = 0.0 + self.phyorp = 0.0 + self.balorp = 0.0 + self.zooorp = 0.0 + self.phyorc = 0.0 + self.balorc = 0.0 + self.zooorc = 0.0 + self.pyco2 = 0.0 + self.baco2 = 0.0 + self.zoco2 = 0.0 + self.phydox = 0.0 + self.zoodox = 0.0 + self.baldox = 0.0 + self.phybod = 0.0 + self.zoobod = 0.0 + self.balbod = 0.0 + self.phytam = 0.0 + self.zootam = 0.0 + self.baltam = 0.0 + self.phyno3 = 0.0 + self.zoono3 = 0.0 + self.balno3 = 0.0 + self.phypo4 = 0.0 + self.zoopo4 = 0.0 + self.balpo4 = 0.0 + + if self.PHYFG == 1: # water scarcity limits phytoplankton growth + self.limphy = 7 #'WAT' + self.phycla = self.phyto * self.cvbcl + self.grophy = 0.0 + self.dthphy = 0.0 + self.zoophy = 0.0 + self.totphy = self.snkphy + + if self.BALFG > 0: # water scarcity limits benthic algae growth + limc = 7 #'WAT' + self.limbal = zeros(self.numbal, dtype=np.int32) + self.balcla = zeros(self.numbal) + self.grobal = zeros(self.numbal) + self.dthbal = zeros(self.numbal) + + for i in range(self.numbal): + self.balcla[i] = self.benal[i] * self.cvbcl + self.limbal[i] = limc + + if self.ZOOFG == 1: # water scarcity limits zooplankton growth + self.grozoo= 0.0 + self.dthzoo= 0.0 + self.totzoo= 0.0 + + #----------------------------------------------------------- + # store final benthic sums and fluxes + #----------------------------------------------------------- + if self.BALFG > 0: + self.tbenal = zeros(3) + self.grotba = 0.0 + self.dthtba = 0.0 + + for i in range(self.numbal): + self.flxbal[1,i] = self.grobal[i] + self.flxbal[2,i] = self.dthbal[i] + self.flxbal[3,i] = self.grobal[i] + self.dthbal[i] + self.tbenal[1] += self.benal[i] + self.grotba += self.grobal[i] + self.dthtba += self.dthbal[i] + + self.tbenal[2] = self.tbenal[1] * self.cvbcl + self.tottba = self.grotba + self.dthtba + + #----------------------------------------------------------- + # compute final process fluxes for oxygen, nutrients and organics + #----------------------------------------------------------- + self.totdox = OXRX.readox + OXRX.boddox + OXRX.bendox + NUTRX.nitdox + self.phydox + self.zoodox + self.baldox + self.totbod = OXRX.decbod + OXRX.bnrbod + OXRX.snkbod + NUTRX.denbod + self.phybod + self.zoobod + self.balbod + self.totno3 = NUTRX.nitno3 + NUTRX.denno3 + NUTRX.bodno3 + self.phyno3 + self.zoono3 + self.balno3 + self.tottam = NUTRX.nittam + NUTRX.volnh3 + NUTRX.bnrtam + NUTRX.bodtam + self.phytam + self.zootam + self.baltam + self.totpo4 = NUTRX.bnrpo4 + NUTRX.bodpo4 + self.phypo4 + self.zoopo4 + self.balpo4 + + self.totorn = self.snkorn + self.phyorn + self.zooorn + self.balorn + self.totorp = self.snkorp + self.phyorp + self.zooorp + self.balorp + self.totorc = self.snkorc + self.phyorc + self.zooorc + self.balorc + + #----------------------------------------------------------- + # compute summaries of total organics, total n and p, and potbod + #----------------------------------------------------------- + + # concentrations: + if self.vol > 0.0: + for i in range(1, 4): + self.lsnh4[i] = NUTRX.rsnh4[i] / self.vol + self.lspo4[i] = NUTRX.rspo4[i] / self.vol + + (self.torn, self.torp, self.torc, self.potbod, self.tn, self.tp) \ + = self.pksums(NUTRX,self.phyto,self.zoo,self.orn,self.orp,self.orc, + no3,tam,no2,self.lsnh4,po4,self.lspo4,bod) + else: + self.torn = -1.0e30 + self.torp = -1.0e30 + self.torc = -1.0e30 + self.potbod = -1.0e30 + self.tn = -1.0e30 + self.tp = -1.0e30 + + # total inflows: + inno3 = NUTRX.ino3 + NUTRX.nuadep[1] + intam = NUTRX.itam + NUTRX.nuadep[2] + inpo4 = NUTRX.ipo4 + NUTRX.nuadep[3] + + (self.itorn, self.itorp, self.itorc, dumval, self.itotn, self.itotp) \ + = self.pksums(NUTRX,self.iphyto,self.izoo,self.iorn,self.iorp,self.iorc, + NUTRX.ino3,NUTRX.itam,NUTRX.ino2,NUTRX.isnh4,NUTRX.ipo4,NUTRX.ispo4,OXRX.ibod) + + # total outflows: + (self.rotorn, self.rotorp, self.rotorc, dumval, self.rototn, self.rototp) \ + = self.pksums(NUTRX,self.rophyt,self.rozoo,self.roorn,self.roorp,self.roorc, + NUTRX.rono3,NUTRX.rotam,NUTRX.rono2,NUTRX.rosnh4,NUTRX.ropo4,NUTRX.rospo4,OXRX.robod) + + # outflows by exit: + self.otorn = zeros(nexits); self.otorp = zeros(nexits); self.otorc = zeros(nexits) + self.ototn = zeros(nexits); self.ototp = zeros(nexits) + + if nexits > 1: + for i in range(nexits): + (self.otorn[i], self.otorp[i], self.otorc[i], dumval, self.ototn[i], self.ototp[i]) \ + = self.pksums(NUTRX,self.ophyt[i],self.ozoo[i],self.oorn[i],self.oorp[i],self.oorc[i], + NUTRX.ono3[i],NUTRX.otam[i],NUTRX.ono2[i],NUTRX.osnh4[i],NUTRX.opo4[i],NUTRX.ospo4[i],OXRX.obod[i]) + + #----------------------------------------------------------- + # update DO/BOD and nutrient states (for OXRX/NUTRX): + #----------------------------------------------------------- + OXRX.dox = dox + OXRX.bod = bod + + NUTRX.po4 = po4 + NUTRX.no3 = no3 + NUTRX.no2 = no2 + NUTRX.tam = tam + + # redistribute TAM after algal influence: + (NUTRX.nh3,NUTRX.nh4) = NUTRX.ammion(tw, phval, tam) + + self.svol = self.vol # svol is volume at start of time step, update for next time thru + + return OXRX, NUTRX + + + @staticmethod + def advplk(iplank,vols,srovol,vol,erovol,sovol,eovol,nexits,oref,mxstay,seed,delts,plank): + ''' advect plankton''' + + # calculate concentration of plankton not subject to advection during interval + oflo = (srovol + erovol) / delts + if oref > 0.0 and oflo / oref <= 100.0: + stay = (mxstay - seed) * (2.0**(-oflo / oref)) + seed + else: + stay = seed + + if plank > stay: + # convert stay to units of mass; this mass will be converted + # back to units of concentration based on the volume of the + # reach/res at the end of the interval + mstay = stay * vols + + # determine concentration of plankton subject to advection; + # this value is passed into subroutine advect + plnkad = plank - stay + + # advect plankton + (plnkad, roplk, oplk) = advect(iplank,plnkad,nexits,vols,vol,srovol,erovol,sovol,eovol) + + # determine final concentration of plankton in reach/res after advection + plank = plnkad + mstay / vol if vol > 0.0 else plnkad + else: # no plankton leaves the reach/res + roplk = 0.0 + oplk = zeros(nexits) + mstay = plank * vols + plank = (mstay + iplank) / vol if vol > 0.0 else -1.0e30 + + return plank, roplk, oplk + + @staticmethod + def algro(light,po4,no3,tw,talgrl,talgrh,talgrm,malgr, cmmp,cmmnp, + TAMFG,AMRFG,tam,NSFG,cmmn,cmmlt,alr20,cflit,delt60, + nmingr,pmingr,lmingr): + + ''' calculate unit growth and respiration rates for algae + population; both are expressed in units of per interval''' + lim = -999 + + if light > lmingr: # sufficient light to support growth + if po4 > pmingr and no3 > nmingr: # sufficient nutrients to support growth + if talgrh > tw > talgrl: # water temperature allows growth + if tw < talgrm: # calculate temperature correction fraction + tcmalg = (tw - talgrl) / (talgrm - talgrl) + else: + # no temperature correction to maximum unit growth rate + # is necessary; water temperature is in the optimum + # range for phytoplankton growth + tcmalg = 1.0 + + # perform temperature correction to maximum unit growth + # rate; units of malgrt are per interval + malgrt = malgr * tcmalg + + # calculate maximum phosphorus limited unit growth rate + grop = malgrt * po4 * no3 / ((po4 + cmmp) * (no3 + cmmnp)) + + # calculate the maximum nitrogen limited unit growth rate + if TAMFG: # consider influence of tam on growth rate + if AMRFG: # calculate tam retardation to nitrogen limited growth rate + malgn = malgrt - 0.757 * tam + 0.051 * no3 + + # check that calculated unit growth rate does not + # exceed maximum allowable growth rate + if malgn > malgrt: + malgn = malgrt + else: + # check that calculated unit growth rate is not + # less than .001 of the maximum unit growth rate; + # if it is, set the unit growth rate equal to .001 + # of the maximum unit growth rate + lolim = 0.001 * malgrt + if malgn < lolim: + malgn = lolim + else: # ammonia retardation is not considered + malgn = malgrt + + if NSFG: # include tam in nitrogen pool for calculation of nitrogen limited growth rate + mmn = no3 + tam + else: # tam is not included in nitrogen pool for calculation of nitrogen limited growth + mmn = no3 + else: # tam is not simulated + malgn = malgrt + mmn = no3 + + # calculate the maximum nitrogen limited unit growth rate + gron = (malgn * mmn) / (cmmn + mmn) + + # calculate the maximum light limited unit growth rate + grol = (malgrt * light) / (cmmlt + light) + if grop < gron and grop < grol: + gro = grop + lim = 5 #'po4' + elif gron < grol: + gro = gron + lim = 4 #'nit' + else: + gro = grol + lim = 1 #'lit' + if gro < 0.000001 * delt60: + gro = 0.0 + + if gro > 0.95 * malgrt: + # there is no limiting factor to cause less than maximum growth rate + lim = 6 #'none' + + # adjust growth rate if control volume is not entirely + # contained within the euphotic zone; e.g. if only one + # half of the control volume is in the euphotic zone, gro + # would be reduced to one half of its specified value + gro = gro * cflit + else: # water temperature does not allow algal growth + gro = 0.0 + lim = 3 #'tem' + else: # no algal growth occurs; necessary nutrients are not available + gro = 0.0 + lim = 2 #'non' + else: # no algal growth occurs; necessary light is not available + gro = 0.0 + lim = 1 #'lit' + + # calculate unit algal respiration rate; res is expressed in + # units of per interval; alr20 is the respiration rate at 20 degrees c + res = alr20 * tw / 20.0 + return lim, gro, res + + + @staticmethod + def baldth(nsfg,no3,tam,po4,paldh,naldh,aldl,aldh,mbal,dox,anaer,oxald,bal,depcor): + ''' calculate benthic algae death''' + + slof = 0.0 + # determine whether to use high or low unit death rate; all + # unit death rates are expressed in units of per interval + + # determine available inorganic nitrogen pool for test of nutrient scarcity + nit = no3 + tam if nsfg > 0 else no3 + if po4 > paldh and nit > naldh: # unit death rate is not incremented by nutrient scarcity check for benthic algae overcrowding + balmax = mbal * depcor + if bal < balmax: # unit death rate is not incremented by benthic algae overcrowding + ald = aldl + slof = 0.0 + else: + # augment unit death rate to account for benthic algae + # overcrowding; set bal value equal to maximum; calculate + # amount of benthic algae in excess of mbal; these benthic + # algae are added to aldth + ald = aldh + slof = bal - balmax + bal = balmax + else: # augment unit death rate to account for nutrient scarcity + ald = aldh + + if dox < anaer: # conditions are anaerobic, augment unit death rate + ald += oxald + + # use unit death rate to compute death rate; dthbal is expressed + # as umoles of phosphorus per liter per interval + return (ald * bal) + slof # dthbal + + + def balrx(self, ballit,tw,talgrl,talgrh,talgrm,malgr,cmmp, cmmnp,tamfg,amrfg,nsfg,cmmn,cmmlt,delt60, + cflit,alr20,cvbpn,phfg,decfg,cvbpc,paldh, naldh,aldl,aldh,anaer,oxald,cfbalg,cfbalr, + alnpr,cvbo,refr,cvnrbo,cvpb,mbal,depcor,cvbcl,baco2,nmingr,pmingr,cmingr,lmingr, + nminc, po4,no3,tam,dox,orn,orp,orc,bod,benal): + ''' simulate behavior of benthic algae in units of umoles p per + liter; these units are used internally within balrx so that + algal subroutines may be shared by phyto and balrx; externally, + the benthic algae population is expressed in terms of areal + mass, since the population is resident entirely on the + bottom surface''' + + # convert benal to units of umoles phosphorus/l (bal) for internal calculations + i0 = 0 + bal = (benal / cvpb) * depcor + + # compute unit growth and respiration rates for benthic algae; determine growth limiting factor + (limbal,gro,res) \ + = self.algro(ballit,po4,no3,tw,talgrl,talgrh,talgrm,malgr,cmmp, cmmnp,tamfg,amrfg, + tam,nsfg,cmmn,cmmlt,alr20, cflit,delt60,nmingr,pmingr,lmingr) + + # calculate net growth rate of algae; grobal is expressed as + # umoles phosphorus per liter per interval; benthic algae growth + # will be expressed in terms of volume rather than area for the + # duration of the subroutines subordinate to balrx; the output + # values for benthic algae are converted to either mg biomass per + # sq meter or mg chla per sq meter, whichever the user + # specifies; cfbalg and cfbalr are the specified ratio of benthic + # algae growth rate to phytoplankton growth rate and ratio of + # benthic algae respiration rate to phytoplankton respiration + # rate, respectively + + grobal = (gro * cfbalg - res * cfbalr) * bal + if grobal > 0.0: # adjust growth rate to account for limitations imposed by availability of required nutrients + grtotn = grobal + grobal = self.grochk (po4,no3,tam,phfg,decfg,baco2,cvbpc,cvbpn,nsfg,nmingr,pmingr,cmingr,i0,grtotn,grobal) + + # calculate benthic algae death, baldth only called here + dthbal = self.baldth(nsfg,no3,tam,po4,paldh,naldh,aldl,aldh,mbal,dox,anaer,oxald,bal,depcor) + + bal += grobal # determine the new benthic algae population + + # adjust net growth rate, if necessary, so that population does not fall below minimum level + minbal = 0.0001 * depcor + if bal < minbal: + grobal += minbal - bal + bal = minbal + bal -= dthbal + + # adjust death rate, if necessary, so that population does not drop below minimum level + if bal < minbal: + dthbal = dthbal - (minbal - bal) + bal = minbal + + # update do state variable to account for net effect of benthic algae photosynthesis and respiration + dobalg = cvpb * cvbo * grobal + # dox = dox + (cvpb*cvbo*grobal) + if dox > -dobalg: + dox = dox + dobalg + else: + dobalg = -dox + dox = 0.0 + + # calculate amount of refractory organic constituents which result from benthic algae death + balorn = refr * dthbal * cvbpn * 0.014 + balorp = refr * dthbal * 0.031 + balorc = refr * dthbal * cvbpc * 0.012 + + # calculate amount of nonrefractory organics (bod) which result from benthic algae death + bodbal = cvnrbo * cvpb * dthbal + + # perform materials balance resulting from benthic algae death + orn,orp,orc,bod = self.orgbal(balorn,balorp,balorc,bodbal,orn,orp,orc,bod) + + # perform materials balance resulting from uptake of nutrients by benthic algae + po4,tam,no3,baco2,tambal,no3bal,po4bal = self.nutrup(grobal,nsfg,cvbpn,alnpr,cvbpc,phfg,decfg,nminc,po4,tam,no3,baco2) + baco2 = -baco2 + + # convert bal back to external units; benal is expressed as + # mg biomass/m2 and balcla is expressed as ug chlorophyll a/m2 + benal = (bal * cvpb) / depcor + balgro = (grobal * cvpb) / depcor + bdth = (dthbal * cvpb) / depcor + balcla = benal * cvbcl + + return po4,no3,tam,dox,orn,orp,orc,bod,benal,limbal,baco2,balcla,dobalg,bodbal,tambal,no3bal,po4bal,balgro,bdth,balorn,balorp,balorc + + @staticmethod + def grochk (po4,no3,tam,phfg,decfg,co2,cvbpc,cvbpn,nsfg,nmingr,pmingr,cmingr,nfixfg,grtotn,grow): + ''' check whether computed growth rate demands too much of any + nutrient; adjust growth rate, if necessary, so that at least + minimum allowed level of each nutrient remains after growth''' + + # calculate growth rate which results in minimum free po4 + # remaining after growth; uplimp is expressed as umoles + # phosphorus per liter per interval + uplimp = (po4 - pmingr) * 32.29 + + # calculate growth rate which results in minimum free + # inorganic nitrogen remaining after growth; uplimn is expressed + # as umoles phosphorus per interval + if nsfg == 0: # tam is not considered as a possible nutrient + uplimn = (no3 - nmingr) * 71.43 / cvbpn + else: + uplimn = (no3 + tam - nmingr) * 71.43 / cvbpn + + uplimc = 1.0e30 + if phfg > 0 and phfg != 2 and decfg == 0: + # phcarb is on, and co2 is being considered as a possible + # limiting nutrient to algal growth + if co2 >= 0.0: + # calculate growth rate which results in minimum free + # carbon dioxide remaining after growth; uplimc is expressed + # as umoles phosphorus per liter per interval + uplimc = (co2 - cmingr) * 83.33 / cvbpc + + # calculate difference between available nutrient concentrations and + # nutrients needed for modeled growth; amount needed for growth may differ + # for n if nitrogen-fixation occurs in any of the algal types + uplimp -= grow + uplimn -= grtotn + uplimc -= grow + + # check that calculated growth does not result in less than + # minimum allowed concentrations of orthophosphate, inorganic + # nitrogen, or carbon dioxide; if it does, adjust growth + if nfixfg: # n-fixation is not occurring for this algal type + uplim = min(uplimp, uplimn, uplimc) + else: # n-fixation is occurring, nitrogen does not limit growth + uplim = min(uplimp, uplimc) + if uplim < 0.0: # reduce growth rate to limit + grow += uplim + return grow + + + @staticmethod + 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 + if inlit > 0.0: + # calculate extinction of light based on the base extinction + # coefficient of the water incremented by self-shading effects + # of phytoplankton and light extinction due to total sediment + # suspension + extco = extb + extcla + extsed + + # calculate euphotic depth; euphotic depth is the distance, + # in feet, below the surface of the water body at which one + # percent of incident light is present + eudep = ln01 / extco + if eudep < avdepe: + # calculate fraction of layer which is contained in the + # euphotic zone; this fraction, cflit, will be multiplied + # by calculated growth in algro to assure that growth only + # occurs in the euphotic zone + cflit = eudep / avdepe + if cflit < 0.0001: + cflit = 0.0 + else: + cflit = 1.0 + + if PHYFG: + # calculate amount of light available to phytoplankton; all + # phytoplankton are assumed to be at mid-depth of the reach; + # light is expressed as langleys per minute + phylit = inlit * exp(-extco * (0.5 * min(eudep, avdepe))) + if phylit < 0.0001: + phylit = 0.0 + + if BALFG: + # calculate amount of light available to benthic algae; all + # benthic algae are assumed to be at the bottom depth of the + # reach;light is expressed as langleys per minute + 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 + return phylit, ballit, cflit + + @staticmethod + def nutrup(grow, NSFG, cvbpn, alnpr, cvbpc, PHFG, DECFG, nminc, po4, tam, no3, alco2): + ''' perform materials balance for transformation from inorganic to + organic material; uptake of po4, no3, tam, and co2 are considered''' + + # calculate po4 balance subsequent to algal uptake or release; + # 0.031 is the conversion from umoles p per liter to mg of p per liter + po4 = po4 - 0.031 * grow + po4alg = -0.031 * grow + tamalg = 0.0 + if NSFG: + # calculate tam balance subsequent to algal uptake or release + # express calculated growth rate in terms of equivalent + # nitrogen; grown is expressed as umoles nitrogen per interval + grown = grow * cvbpn + if grow < 0.0: # algal respiration exceeds growth + # nitrogen released by respiration is released in the form of tam + # no uptake or release of no3 occurs + altam = grown + alno3 = 0.0 + else: + # calculate amount of n uptake which is no3 and amount which is tam + alno3 = alnpr * grown + altam = grown - alno3 + # check that computed uptake of no3 does not consume more + # than 99 percent of available free no3; if it does, satisfy + # excess demand with free tam; no3lim is expressed as umoles + # n per liter per interval + no3lim = 70.72 * no3 + if alno3 > no3lim: + altam += alno3 - no3lim + alno3 = no3lim + else: + # check that calculated uptake of tam does not consume + # more than 99 percent of available free tam; if it does, + # satisfy excess demand with free no3; tamlim is expressed + # as umoles n per liter per interval + tamlim = 70.72 * tam + if altam > tamlim: + alno3 += altam - tamlim + altam = tamlim + + # calculate net uptake or release of tam by algae; .014 is + # the conversion from umoles of n per liter per interval to + # mg n per liter per interval + tams = tam + tamalg = -0.014 * altam + tam = tam - 0.014 * altam + if tam < nminc: + tamalg = -tams + tam = 0.0 + + else: # all inorganic n is in the form of no3 + alno3 = grow * cvbpn + + # calculate no3 balance subsequent to algal uptake or release; + # eliminate insignificant values of no3 + no3s = no3 + no3alg = -0.014 * alno3 + no3 = no3 - 0.014 * alno3 + if no3 < nminc: + no3alg = -no3s + no3 = 0.0 + + if PHFG and DECFG == 0: # calculate amount of algal uptake of co2 in mg co2-c/liter + alco2 = grow * cvbpc * 0.012 + else: + alco2 = 0.0 + return po4, tam, no3, alco2, tamalg, no3alg, po4alg + + @staticmethod + def orgbal(dthorn, dthorp, dthorc, dthbod, orn, orp, orc, bod): + ''' perform materials balance for transformation from living to dead organic material''' + + # calculate dead refractory organic nitrogen balance subsequent to plankton death; + # plankton death may be either algal death, zooplankton death, or phytoplankton + # filtered by zooplankton but not assimilated + orn += dthorn + + # calculate dead refractory organic phosphorus balance subsequent to plankton death + orp += dthorp + + # calculate dead refractory organic carbon balance subsequent to plankton death + orc += dthorc + + # calculate bod balance subsequent to plankton death + bod += dthbod + + return orn, orp, orc, bod + + @staticmethod + def phydth(nsfg,no3,tam,po4,paldh,naldh,phycla,claldh,aldl,aldh,dox,anaer,oxald,stc): + ''' calculate phytoplankton death''' + + # determine whether to use high or low unit death rate; all unit + # death rates are expressed in units of per interval + # determine available inorganic nitrogen pool for test of nutrient scarcity + nit = no3 + tam if nsfg else no3 + if po4 > paldh and nit > naldh: + # unit death rate is not incremented by nutrient scarcity + # check for phytoplankton overcrowding + if phycla < claldh: # unit death rate is not incremented by phytoplankton overcrowding + ald = aldl + else: + ald = aldh # augment unit death rate to account for overcrowding + else: # augment unit death rate to account for nutrient scarcity + ald = aldh + + # augment unit death rate if conditions are anaerobic + if dox < anaer: + ald += oxald + + # use unit death rate to compute death rate; aldth is expressed + # as umoles of phosphorus per liter per interval + return ald * stc # dthphy + + + def phyrx(self,phylit,tw,talgrl,talgrh,talgrm,malgr,cmmp,cmmnp,tamfg,amrfg,nsfg,cmmn,cmmlt,delt60,cflit,alr20,cvbpn,phfg,decfg,cvbpc,paldh, + naldh,claldh,aldl,aldh,anaer,oxald,alnpr,cvbo,refr,cvnrbo,cvpb,cvbcl,co2,nmingr,pmingr,cmingr,lmingr,nminc, + po4,no3,tam,dox,orn,orp,orc,bod,phyto): + ''' simulate behavior of phytoplankton, as standing crop, in units of umoles p per liter''' + + # convert phyto to units of umoles phosphorus (stc) and ug chlorophyll a/l (phycla) for internal calculations + stc = phyto / cvpb + phycla = phyto * cvbcl + + # compute unit growth and respiration rates for phytoplankton; + # determine growth limiting factor + (limphy,gro,res) \ + = self.algro(phylit,po4,no3,tw,talgrl,talgrh,talgrm,malgr,cmmp,cmmnp,tamfg,amrfg, + tam,nsfg,cmmn,cmmlt,alr20,cflit,delt60,nmingr,pmingr,lmingr) + + # calculate net growth rate of phytoplankton; grophy is + # expressed as umol phosphorus per liter per interval + grophy = (gro - res) * stc + + if grophy > 0.0: + # adjust growth rate to account for limitations imposed by + # availability of required nutrients + grtotn = grophy + i0 = 0 + + grophy = self.grochk (po4,no3,tam,phfg,decfg,co2,cvbpc,cvbpn,nsfg,nmingr,pmingr,cmingr,i0,grtotn,grophy) + + # calculate phytoplankton death + dthphy = self.phydth(nsfg,no3,tam,po4,paldh,naldh,phycla,claldh,aldl,aldh,dox,anaer,oxald,stc) + + # determine the new phytoplankton population + stc += grophy + + # adjust net growth rate, if necessary, so population does not fall below minimum level + if stc < .0025: + grophy -= (0.0025 - stc) + stc = 0.0025 + stc -= dthphy + + # adjust death rate, if necessary, so that population does not drop below minimum level + if stc < .0025: + dthphy -= (0.0025 - stc) + stc = 0.0025 + + # update do state variable to account for net effect of phytoplankton photosynthesis and respiration + dophy = cvpb * cvbo * grophy + if dox > -dophy: + dox += dophy + else: + dophy = -dox + dox = 0.0 + + # calculate amount of refractory organic constituents which result from phytoplankton death + phyorn = refr * dthphy * cvbpn * 0.014 + phyorp = refr * dthphy * 0.031 + phyorc = refr * dthphy * cvbpc * 0.012 + + # calculate amount of nonrefractory organics (bod) which result from phytoplankton death + phybd = cvnrbo * cvpb * dthphy + bodphy = phybd + + # perform materials balance resulting from phytoplankton death + orn,orp,orc,bod = self.orgbal(phyorn,phyorp,phyorc,phybd, orn,orp,orc,bod) + + # perform materials balance resulting from uptake of nutrients by phytoplankton + (po4,tam,no3,pyco2,tamphy,no3phy,po4phy) = self.nutrup(grophy,nsfg,cvbpn,alnpr,cvbpc,phfg,decfg,nminc,po4,tam,no3,co2) + pyco2 = -pyco2 + + # convert stc to units of mg biomass/l (phyto) and ug chlorophyll a/l (phycla) + phyto = stc * cvpb + phgro = grophy * cvpb + phdth = dthphy * cvpb + phycla = phyto * cvbcl + + return po4,no3,tam,dox,orn,orp,orc,bod,phyto,limphy,pyco2,phycla,dophy,bodphy,tamphy,no3phy,po4phy,phdth,phgro,phyorn,phyorp,phyorc + + + def zorx(self,zfil20,tczfil,tw,phyto,mzoeat,zexdel,cvpb,zres20,tczres,anaer,zomass,tamfg,refr, + zfood,zd,oxzd,cvbn,cvbp,cvbc,cvnrbo,cvbo,dox,bod,zoo,orn,orp,orc,tam,no3,po4): + ''' calculate zooplankton population balance''' + + # calculate zooplankton unit grazing rate expressed as liters + # of water filtered per mg zooplankton per interval + zfil = zfil20 * (tczfil**(tw - 20.0)) + + # calculate mass of phytoplankton biomass ingested per mg + # zooplankton per interval + zoeat = zfil * phyto + + # check that calculated unit ingestion rate does not exceed + # maximum allowable unit ingestion rate (mzoeat); if it does, + # set unit ingestion rate equal to mzoeat + if zoeat >= mzoeat: + zoeat = mzoeat + # nonrefractory portion of excretion is only partially + # decomposed; zexdec is the fraction of nonrefractory + # material which is decomposed + zexdec = zexdel + else: + # calculated unit ingestion rate is acceptable + # all nonrefractory excretion is decomposed + zexdec = 1.0 + + # calculate phytoplankton consumed by zooplankton; zeat is + # expressed as mg biomass per interval + zeat = zoeat * zoo + + # check that calculated ingestion does not reduce phytoplankton + # concentration to less than .0025 umoles p; if it does, adjust + # ingestion so that .0025 umoles of phytoplankton (expressed as p) remain + if phyto - zeat < 0.0025 * cvpb: + zeat = phyto - (0.0025 * cvpb) + + # calculate zooplankton assimilation efficiency + if zfood == 1: + # calculate assimilation efficiency resulting from ingestion + # of high quality food; zeff is dimensionless + zeff = -0.06 * phyto + 1.03 + if zeff > 0.99: # set upper limit on efficiency at 99 percent + zeff = 0.99 + elif zfood == 2: + # calculate assimilation efficiency resulting from ingestion of medium quality food + zeff = -0.03 * phyto + 0.47 + if zeff < 0.20: # set lower limit on efficiency at 20 percent + zeff = 0.20 + else: + # calculate assimilation efficiency resulting from ingestion of low quality food + zeff = -0.013 * phyto + 0.17 + if zeff < 0.03: + zeff= 0.03 # # set lower limit on efficiency at 3 percent + + # calculate zooplankton growth; zogr is expressed as mg biomass per liter per interval + zogr = zeff * zeat + + # calculate total zooplankton excretion (zexmas),excretion + # decomposed to inorganic constituents (zingex), excretion + # released as dead refractory constituents (zrefex), and + # excretion released as dead nonrefractory material (znrfex) + zexmas = zeat - zogr + zrefex = refr * zexmas + zingex = zexdec * (zexmas - zrefex) + znrfex = zexmas - zrefex - zingex + + # calculate zooplankton respiration; zres is expressed as mg + # biomass per liter per interval + zres = zres20 * (tczres**(tw - 20.0)) * zoo + + # calculate zooplankton death; zdth is expressed as mg biomass per liter per interval + if dox > anaer: + # calculate death using aerobic death rate + zdth = zd * zoo + else: + # calculate death using sum of aerobic death rate and anaerobic increment + zdth = (zd + oxzd) * zoo + + # calculate zooplankton population after growth, respiration, + # and death; adjust respiration and death, if necessary, to + # assure minimum population of zooplankton + # first, account for net growth (growth - respiration) + zoo += zogr - zres + + # maintain minimum population of .03 organisms per liter; zomass + # is a user specified conversion factor from organisms/l to mg biomass/l + lolim = 0.03 * zomass + if zoo < lolim: + zres = zres + zoo - lolim + zoo = lolim + + # subtract oxygen required to satisfy zooplankton respiration from do state variable + dozoo = 1.1 * zres + dox -= dozoo + zbod = 0.0 + if dox < 0.0: # include oxygen deficit in bod value + zbod = -dox + dox = 0.0 + + # subtract computed zooplankton death from zooplankton state variable + zoo = zoo - zdth + if zoo < lolim: + zdth = zdth + zoo - lolim + zoo = lolim + + # calculate amount of inorganic constituents which are released + # by zooplankton respiration and inorganic excretion + znit = (zingex + zres) * cvbn + zpo4 = (zingex + zres) * cvbp + zco2 = (zingex + zres) * cvbc + + # update state variables for inorganic constituents to account + # for additions from zooplankton respiration and inorganic excretion + i1 = 1 + tam, no3, po4 = decbal(tamfg, i1, znit, zpo4, tam, no3, po4) + + # calculate amount of refractory organic constituents which result from zooplankton death and excretion + zorn = ((refr * zdth) + zrefex) * cvbn + zorp = ((refr * zdth) + zrefex) * cvbp + zorc = ((refr * zdth) + zrefex) * cvbc + + # calculate amount of nonrefractory organics (bod) which result from zooplankton death and excretion + zbod = zbod + (zdth * cvnrbo) + (znrfex * cvbo) + orn, orp, orc, bod = self.orgbal(zorn, zorp, zorc, zbod, orn, orp, orc, bod) + + return dox,bod,zoo,orn,orp,orc,tam,no3,po4,zeat,zco2,dozoo,zbod,znit,zpo4,zogr,zdth,zorn,zorp,zorc + + def pksums(self, NUTRX, phyto, zoo, orn, orp, orc, no3, tam, no2, snh4, po4, spo4, bod): + ''' computes summaries of: total organic n, p, c; total n, p; potbod''' + + # undefined summary concentrations: + torn = -1.0e30 + torp = -1.0e30 + torc = -1.0e30 + potbod = -1.0e30 + tn = -1.0e30 + tp = -1.0e30 + + if (self.vol <= 0): + return + + # Calculate sums: + tval = 0.0 + + if self.PHYFG == 1: + tval += phyto + if self.ZOOFG == 1: + tval += zoo + + torn = orn + self.cvbn * tval + torp = orp + self.cvbp * tval + + tval += bod / self.cvbo + torc = orc + self.cvbc * tval + potbod = bod + + # total N + tn = torn + no3 + if NUTRX.TAMFG == 1: + tn += tam + if NUTRX.NO2FG == 1: + tn += no2 + if NUTRX.ADNHFG == 1: + for i in range(1, 4): + tn += snh4[i] + + # total P + tp = torp + if NUTRX.PO4FG == 1: + tp += po4 + if NUTRX.ADPOFG == 1: + for i in range(1, 4): + tp += spo4[i] + + # potential BOD + if self.PHYFG == 1: + potbod += (self.cvnrbo * phyto) + + if self.ZOOFG == 1: + potbod += (self.cvnrbo * zoo) + + return torn, torp, torc, potbod, tn, tp + + @staticmethod + def algro2 (ballit,po4,no3,tw,mbalgr,cmmp,TAMFG,tam,NSFG,cmmn,balr20,delt60,tcbalg,balvel, + cmmv,BFIXFG,cslit,cmmd1,cmmd2,sumba,tcbalr,nmingr,pmingr,lmingr,nmaxfx,grores): + + ''' calculate unit growth and respiration rates for benthic algae using more + complex kinetics; both are expressed in units of per interval''' + + groba = 0.0 + NFIXFG = 0 + lim = -999 + + if ballit > lmingr: # sufficient light to support growth + if po4 > pmingr and no3 > nmingr: # sufficient nutrients to support growth + + # calculate temperature correction fraction + tcmbag = tcbalg**(tw - 20.0) + + # calculate velocity limitation on nutrient availability + grofv = balvel / (cmmv + balvel) + + # calculate phosphorus limited unit growth factor + grofp = po4 * grofv / (cmmp + po4 * grofv) + + # calculate the nitrogen limited unit growth factor + if TAMFG: # consider influence of tam on growth rate + if NSFG: + # include tam in nitrogen pool for calculation of + # nitrogen limited growth rate + mmn = no3 + tam + else: + # tam is not included in nitrogen pool for calculation + # of nitrogen limited growth + mmn = no3 + else: + mmn = no3 # tam is not simulated + + if BFIXFG != 1: # calculate the maximum nitrogen limited unit growth rate + NFIXFG = 0 + grofn = (mmn * grofv) / (cmmn + mmn * grofv) + else: + # n-fixing blue-green algae; determine if available nitrogen + # concentrations are high enough to suppress nitrogen fixation + if mmn >= nmaxfx: + NFIXFG = 0 + grofn = (mmn * grofv) / (cmmn + mmn * grofv) + else: + NFIXFG = 1 # available nitrogen concentrations do not suppress n-fixation + grofn = 1.0 # nitrogen does not limit growth rate + + # calculate the maximum light limited unit growth rate + grofl = (ballit / cslit) * exp(1.0 - (ballit / cslit)) + + # calculate density limitation on growth rate + grofd = (cmmd1 * sumba + cmmd2) / (sumba + cmmd2) + if grofp < grofn and grofp < grofl: # phosphorus limited + gromin = grofp + lim = 5 #'po4' + elif grofn < grofl: # nitrogen limited + gromin = grofn + lim = 4 #'nit' + else: # light limited + gromin = grofl + lim = 1 #'lit' + if gromin > 0.95: # there is no limiting factor to cause less than maximum growth rate + lim = 6 #'none' + # calculate overall growth rate in units of per interval + groba = mbalgr * tcmbag * gromin * grofd + if groba < 1.0e-06 * delt60: + groba = 0.0 + else: # no algal growth occurs; necessary nutrients are not available + groba = 0.0 + lim = 2 #'non' + else: # no algal growth occurs; necessary light is not available + groba = 0.0 + lim = 1 #'lit' + + # calculate unit algal respiration rate in units of per interval + # balr20 is the benthic algal respiration rate at 20 degrees c + resba = balr20 * tcbalr**(tw - 20.0) + grores * groba + return NFIXFG, lim, groba, resba + + @staticmethod + def balrem(crem,sumba,sumbal,cmmbi,tcgraz,tw,binv,bal,cslof1,cslof2,balvel): + ''' calculate benthic algae removal due to macroinvertebrates and scouring, based + upon equations from dssamt. This subroutine provides similar functionality for the + enhanced periphyton kinetics (balrx2) that baldth provides for the original hspf benthic algae (balrx).''' + + # calculate the total benthic algae removed per interval due to grazing and disturbance of macroinvertebrates + reminv = crem * (sumba / (sumba + cmmbi)) * tcgraz**(tw - 20.0) * binv + + # calculate portion of total algae removed represented by the present benthic algal group + remba = reminv * bal / sumbal + + # calculate scour loss rate of benthic algae in per interval + maxvel = log(1.0 / cslof1) / cslof2 + slof = cslof1 * exp(cslof2 * balvel) if balvel < maxvel else 1.0 + + # calculate amount of present benthic algae group removed through grazing, disturbance, and scouring + remba += slof * bal + return remba + + @staticmethod + def nutup2 (grow,NSFG,cvbpn,alnpr,cvbpc,PHFG,DECFG,BNPFG,campr,sumgrn,nminc,po4,tam,no3): + ''' perform materials balance for transformation from inorganic to + organic material; uptake of po4, no3, tam, and co2 are considered. + used instead of nutrup by the balrx2 subroutine; adds a + calculated nitrogen preference function and adjustments for + n-fixation to the nutrup code.''' + + nmgpum = 0.0140067 + cmgpum = 0.012011 + pmgpum = 0.0309738 + prct99 = 0.99 + + # calculate po4 balance subsequent to algal uptake or release; + # pmgpum is the conversion from umoles p per liter to mg of p per liter + po4 = po4 - pmgpum * grow + po4alg = -pmgpum * grow + tamalg = 0.0 + if NSFG: + # calculate tam balance subsequent to algal uptake or release + # express calculated growth rate in terms of equivalent + # nitrogen; grown is expressed as umoles nitrogen per interval; + # use accumulated growth that affects available n (sumgrn) + grown = sumgrn * cvbpn + if sumgrn < 0.0: # algal respiration exceeds growth + # nitrogen released by respiration is released in the form of tam; no uptake or release of no3 occurs + altam = grown + alno3 = 0.0 + else: # calculate amount of n uptake which is no3 and amount which is tam + # check nitrogen preference flag; if bnpfg.ne.1 use original + # approach of %no3, if bnpfg.eq.1 use preference function + # from dssamt/ssamiv + if BNPFG == 1: + # use dssamt ammonia preference function, equal to 1-alnpr + alnpr2 = 1.0 - (campr * tam) / (campr * tam + no3) + else: + alnpr2 = alnpr + alno3 = alnpr2 * grown + altam = grown - alno3 + # check that computed uptake of no3 does not consume more + # than 99 percent of available free no3; if it does, satisfy + # excess demand with free tam; no3lim is expressed as umoles + # n per liter per interval + no3lim= (prct99/nmgpum)*no3 + + if alno3 > no3lim: + altam = altam + alno3 - no3lim + alno3 = no3lim + else: + # check that calculated uptake of tam does not consume + # more than 99 percent of available free tam; if it does, + # satisfy excess demand with free no3; tamlim is expressed + # as umoles n per liter per interval + tamlim = (prct99 / nmgpum) * tam + if altam > tamlim: + alno3 = alno3 + altam - tamlim + altam = tamlim + # calculate net uptake or release of tam by algae; .014 is + # the conversion from umoles of n per liter per interval to + # mg n per liter per interval + tams = tam + tamalg = -nmgpum * altam + tam = tam - nmgpum * altam + if tam < nminc: + tamalg = -tams + tam = 0.0 + else: # all inorganic n is in the form of no3 + alno3 = sumgrn * cvbpn + + # calculate no3 balance subsequent to algal uptake or release; + # eliminate insignificant values of no3 + no3s = no3 + no3alg = -nmgpum * alno3 + no3 = no3 - nmgpum * alno3 + if no3 < nminc: + no3alg = -no3s + no3 = 0.0 + + # calculate amount of algal uptake of co2; alco2 is expressed as mg co2-c/liter + alco2 = grow * cvbpc * cmgpum if PHFG and DECFG == 0 else 0.0 + + return po4, tam, no3, alco2, tamalg, no3alg, po4alg + + + def balrx2(self,ballit,tw,TAMFG,NSFG,delt60,cvbpn,PHFG,DECFG,cvbpc,alnpr,cvbo,refr,cvnrbo,cvpb, + depcor,cvbcl,co2,numbal,mbalgr,cmmpb,cmmnb,balr20,tcbalg,balvel,cmmv,BFIXFG, + cslit,cmmd1,cmmd2,tcbalr,frrif,cremvl,cmmbi,binv,tcgraz,cslof1,cslof2,minbal, + fravl,BNPFG,campr,nmingr,pmingr,cmingr,lmingr,nminc,nmaxfx,grores, + po4,no3,tam,dox,orn,orp,orc,bod,benal): + + ''' simulate behavior of up to four types of benthic algae using + algorithms adapted from the dssamt model. this subroutine was + adapted from balrx''' + nmgpum = 0.0140067 + cmgpum = 0.012011 + pmgpum = 0.0309738 + + # compute total biomass of all benthic algal types + sumba = 0.0 + for i in range(numbal): + sumba = sumba + benal[i] + + # convert to umoles p/l + sumbal = (sumba / cvpb) * depcor * frrif + + # initialize variables for cumulative net growth/removal + # of benthic algae types, and for growth affecting available n conc. + grotot = 0.0 + grtotn = 0.0 + sumgro = 0.0 + sumgrn = 0.0 + sumdth = 0.0 + + bal = zeros(numbal) + grobal = zeros(numbal) + dthbal = zeros(numbal) + limbal = zeros(numbal, dtype=np.int32) + + for i in range(numbal): # do 20 i = 1, numbal + # convert benal to units of umoles phosphorus/l (bal) for + # internal calculations, and adjust for % riffle to which + # periphyton are limited + bal[i] = (benal[i] / cvpb) * depcor * frrif + + # compute unit growth and respiration rates + (NFIXFG,limbal[i],groba,resba) \ + = self.algro2 (ballit,po4,no3,tw,mbalgr[i],cmmpb[i],TAMFG,tam,NSFG,cmmnb[i],balr20[i], + delt60,tcbalg[i],balvel,cmmv,BFIXFG[i],cslit[i],cmmd1[i],cmmd2[i],sumba,tcbalr[i], + nmingr,pmingr,lmingr,nmaxfx,grores[i]) + + # calculate net growth rate of algae; grobal is expressed as + # umoles phosphorus per liter per interval; benthic algae growth + # will be expressed in terms of volume rather than area for the + # duration of the subroutines subordinate to balrx; the output + # values for benthic algae are converted to either mg biomass per + # sq meter or mg chla per sq meter, whichever the user specifies + grobal[i] = (groba - resba) * bal[i] + + # track growth that affects water column available n concentrations + if NFIXFG != 1: # algae are not fixing n, net growth affects concentrations + groban = grobal[i] + else: # algae that are fixing n affect water column n through respiration only + groban = -resba * bal[i] + + # calculate cumulative net growth of algal types simulated so far + grotot = grotot + grobal[i] + + # calculate cumulative algal growth that affects available n + grtotn = grtotn + groban + + if grobal[i] > 0.0 and grotot > 0.0: + # check that cumulative growth rate of algal types does not exceed + # limitations imposed by the availability of required nutrients; + # if so, reduce growth rate of last algal type + + grotmp = grotot # set temporary variable for comparison purposes + grotot = self.grochk (po4,no3,tam,PHFG,DECFG,co2,cvbpc,cvbpn,NSFG,nmingr,pmingr,cmingr,NFIXFG,grtotn,grotot) + if grotot < 0.0: # this should never happen + grotot = 0.0 + + # compare nutrient-checked growth to original cumulative total, + if grotot < grotmp: + # adjust growth rate of last algal type + grobal[i] = grobal[i] - (grotmp - grotot) + # track changes in growth that affect available n concentrations + if NFIXFG != 1: # n-fixation not occurring, all growth affects n + groban = grobal[i] + else: # n-fixation is occurring, proportionately adjust respiration (groban) + groban = groban * grotot / grotmp + + # calculate benthic algae removal + crem = cremvl * depcor * frrif + dthbal[i] = self.balrem (crem,sumba,sumbal,cmmbi,tcgraz,tw,binv,bal[i],cslof1[i],cslof2[i],balvel) + + # add the net growth + bal[i] += grobal[i] + balmin = minbal * depcor * frrif + + if bal[i] < balmin and grobal[i] < 0.0: # adjust net growth rate so that population does not fall below minimum level + grotm2 = grobal[i] # set temporary variable for growth + grobal[i] = grobal[i] + (balmin- bal[i]) + bal[i] = balmin + # adjust growth that affects available n concentrations + if NFIXFG != 1: # n-fixation not occurring, all growth affects n + groban = grobal[i] + else: # n-fixation is occurring, proportionately adjust respiration (groban) + groban = groban * grobal[i] / grotm2 + + # subtract death/removal + bal[i] -= dthbal[i] + if bal[i] < balmin: # adjust death rate so that population does not drop below minimum level + dthbal[i] = dthbal[i] - (balmin- bal[i]) + if dthbal[i] < 0.0: + dthbal[i] = 0.0 + bal[i] = balmin + + # calculate total net growth and removal of all benthic algae types + sumgro = sumgro + grobal[i] + sumdth = sumdth + dthbal[i] + sumgrn = sumgrn + groban + # update internal loop tracking variables for cumulative growth + # to account for grochk and minimum biomass adjustments + grotot = sumgro + grtotn = sumgrn + # 20 continue + + # update do state variable to account for net effect of benthic algae photosynthesis and respiration + dobalg = cvpb * cvbo * sumgro + if dox > -dobalg: # enough oxygen available to satisfy demand + dox += dobalg + else: # take only available oxygen + dobalg = -dox + dox = 0.0 + + # calculate amount of refractory organic constituents which result from benthic algae death + balorn = refr * sumdth * cvbpn * nmgpum + balorp = refr * sumdth * pmgpum + balorc = refr * sumdth * cvbpc * cmgpum + + # add to orc the carbon associated with nutrients immediately + # released to the available pool from removed benthic algal biomass + balorc = balorc + fravl * (1.0 - refr) * sumdth * cvbpc * cmgpum + + # calculate amount of nonrefractory organics (bod) which result from benthic algae death + bodbal = cvnrbo * (1.0 - fravl) * cvpb * sumdth + + # perform materials balance resulting from benthic algae death + (orn,orp,orc,bod) = self.orgbal (balorn,balorp,balorc,bodbal,orn,orp,orc,bod) + + # perform materials balance resulting from uptake of nutrients by benthic algae + (po4,tam,no3,baco2,tambal,no3bal,po4bal) \ + = self.nutup2 (sumgro,NSFG,cvbpn,alnpr,cvbpc,PHFG,DECFG,BNPFG, + campr,sumgrn,nminc,po4,tam,no3) + baco2 = -baco2 + + # update available nutrient pools with nutrients immediately + # released from benthic algal biomass removal processes; nutrients + # immediately cycled are calculated as a fraction (fravl) of the + # nonrefractory biomass + avlpho = fravl * (1.0 - refr) * sumdth * pmgpum + po4 = po4 + avlpho + po4bal = po4bal + avlpho + + # nitrogen is released as tam if tam is simulated otherwise as no3 + avlnit = fravl * (1.0 - refr) * sumdth * cvbpn * nmgpum + if TAMFG: + tam = tam + avlnit + tambal = tambal + avlnit + else: + no3 = no3 + avlnit + no3bal = no3bal + avlnit + + # convert biomass to units of mg biomass/m2 and ug chlorophyll a/m2 + balgro = zeros(numbal) + bdth = zeros(numbal) + balcla = zeros(numbal) + + for i in range(numbal): + benal[i] = (bal[i] * cvpb) / (depcor * frrif) + balgro[i] = (grobal[i] *cvpb) / (depcor * frrif) + bdth[i] = (dthbal[i] *cvpb) / (depcor * frrif) + balcla[i] = benal[i] * cvbcl + + return po4,no3,tam,dox,orn,orp,orc,bod,benal,limbal,baco2,balcla,dobalg,bodbal,tambal,no3bal,po4bal,balgro,bdth,balorn,balorp,balorc \ No newline at end of file diff --git a/HSP2/RQUAL.py b/HSP2/RQUAL.py index 99991c37..7fbcd3ab 100644 --- a/HSP2/RQUAL.py +++ b/HSP2/RQUAL.py @@ -3,262 +3,356 @@ License: LGPL2 ''' -from numpy import where -from HSP2.utilities import make_numba_dict -from HSP2.OXRX import oxrx -#from HSP2.NUTRX import nutrx -#from HSP2.PLANK import plank -#from HSP2.PHCARB import phcarb +import logging +import numpy as np +from numpy import where, zeros, array, float64 +from numba import types +from numba.typed import Dict -UUNITS = 1 +from HSP2.utilities import make_numba_dict, initm +from HSP2.RQUAL_Class import RQUAL_Class +ERRMSGS_oxrx = ('OXRX: Warning -- SATDO is less than zero. This usually occurs when water temperature is very high (above ~66 deg. C). This usually indicates an error in input GATMP (or TW, if HTRCH is not being simulated).',) +ERRMSGS_nutrx = ('NUTRX: Error -- Inconsistent flags for NH4; TAM is not being simulated, but NH3 volatilization and/or NH4 adsorption are being simulated.', + 'NUTRX: Error -- Inconsistent flags for PO4; PO4 is not being simulated, but PO4 adsorption is being simualted.', + 'NUTRX: Error -- Sediment-associated NH4 and/or PO4 is being simulated, but sediment is not being simulated in module SEDTRN.', + 'NUTRX: Error -- Inorganic nutrient mass stored in or leaving the reach non-zero, but is expected to be non-zero due to lack of suspended sediment mass.', + 'NUTRX: Error -- Inorganic nutrient mass in bed is expected to be zero (due to the lack of bed sediments).') +ERRMSGS_plank = ('PLANK: Error -- Zooplankton cannot be simulated without phytoplankton.', + 'PLANK: Error -- Ammonia cannot be included in the N supply if it is not being simulated.', + 'PLANK: Error -- Phosphate must be simulated if plankton are being simulated.') +ERRMSGS_phcarb = ('PHCARB: Error -- Invalid CONS index specified for ALKCON (i.e., ALKCON > NCONS).', + 'PHCARB: Error -- A satisfactory solution for pH has not been reached.') - -def rqual(store, siminfo, uci, ts): +def rqual(store, siminfo, uci, uci_oxrx, uci_nutrx, uci_plank, uci_phcarb, ts): ''' Simulate constituents involved in biochemical transformations''' + #ERRMSGS =('') + #errors = zeros(len(ERRMSGS), dtype=np.int32) + + # simulation information: + delt60 = siminfo['delt'] / 60 # delt60 - simulation time interval in hours simlen = siminfo['steps'] + delts = siminfo['delt'] * 60 + uunits = siminfo['units'] + + siminfo_ = Dict.empty(key_type=types.unicode_type, value_type=types.float64) + for key in set(siminfo.keys()): + value = siminfo[key] + + if type(value) in {int, float}: + siminfo_[key] = float(value) + # module flags: ui = make_numba_dict(uci) - BENRFG = ui['BENFGX'] # table-type benth-flag + NUTFG = int(ui['NUTFG']) + PLKFG = int(ui['PLKFG']) + PHFG = int(ui['PHFG']) - # table type ACTIVITY - NUTFG = ui['NUTFG'] - PLKFG = ui['PLKFG'] - PHFG = ui['PHFG'] - - # get external time series - PREC = ts['PREC'] - SAREA = ts['SAREA'] - AVDEP = ts['AVDEP'] - AVVEL = ts['AVVEL'] - TW = ts['TW'] - if LKFG == 1: - WIND = ts['WIND'] - - nexits, vol, VOL, SROVOL, EROVOL, SOVOL, EOVOL = ui['advectData'] - - TW = where(TW < -100.0, 20.0, TW) # fix undefined temps if present - AVDEPE = where(UUNITS == 2, AVDEP * 3.28, AVDEP) # convert to english units) in feet - AVVELE = where(UUNITS == 2, AVVEL * 3.28, AVVEL) # convert to english units) - DEPCOR = where(avdepe > 0.0, 3.28084e-3 / AVDEPE, -1.e30) # # define conversion factor from mg/m2 to mg/l - if BENRFG == 1: - scrvel = ui['SCRVEL'] # table-type scour-parms - scrmul = ui['SCRMUL'] # table-type scour-parms - SCRFAC = where(AVVELE > scrvel, scrmul, 1.0) # calculate scouring factor - - # use Closures to capture 'ui' data to minimize calling arguments. - #### OXRX #### - - IDOX = ts['IDOX'] # optional, input flow - IBOD = ts['IBOD'] # optional, input flow - - # preallocate storage for OXRX calculated results - DOX = ts['DOX'] = zeros(simlen) # concentration, state variable - BOD = ts['BOD'] = zeros(simlen) # concentration, state variable - SATDO = ts['SATDO'] = zeros(simlen) # concentration, state variable - RODOX = ts['RODOX'] = zeros(simlen) # reach outflow of DOX - ROBOD = ts['ROBOD'] = zeros(simlen) # reach outflow of BOD - ODOX = ts['ODOX'] = zeros((simlen, nexits)) # reach outflow per gate of DOX - OBOD = ts['OBOD'] = zeros((simlen, nexits)) # reach outflow per gate of BOD - - oxrx = poxrx() # returns Numba accelerated function in closure - - if NUTFG: - # get NUTRX specific time series - INO3 = ts['INO3'] # optional, input - INH3 = ts['INH3'] # optional, input - INO2 = ts['INO2'] # optional, input - IPO4 = ts['IPO4'] # optional, input - - NUAFX = setit() # NUAFXM monthly, constant or time series - NUACN = setit() # NUACNM monthly, constant or time series - - # preallocate storage for computed time series - NO3 = ts['NO3'] = zeros(simlen) # concentration, state variable - NO2 = ts['NO2'] = zeros(simlen) # concentration, state variable - NH3 = ts['NH3'] = zeros(simlen) # concentration, state variable - PO4 = ts['PO4'] = zeros(simlen) # concentration, state variable - TAM = ts['TAM'] = zeros(simlen) # concentration, state variable - RONO3 = ts['RONO3'] = zeros(simlen) # outflow - RONO2 = ts['RONO2'] = zeros(simlen) # outflow - RONH3 = ts['RONH3'] = zeros(simlen) # outflow - ROPO4 = ts['ROPO4'] = zeros(simlen) # outflow - ONO3 = ts['ONO3'] = zeros((simlen, NEXITS)) # outflow - ONO2 = ts['ONO2'] = zeros((simlen, NEXITS)) # outflow - ONH3 = ts['ONH3'] = zeros((simlen, NEXITS)) # outflow - OPO4 = ts['OPO4'] = zeros((simlen, NEXITS)) # outflow - - nutrx = pnutrx() # returns Numba accelerated function in closure - - if PLKFG: - # get PLANK specific time series - IPHYTO = ts['IPHYTO'] # optional - IZOO = ts['IZOO'] # optional - IORN = ts['IORN'] # optional - IORP = ts['IORP'] # optional - IORC = ts['IORC'] # optional - WASH = ts['WASH'] - SOLRAD = ts['SOLRAD'] - - # preallocate arrays for better performance - ORN = ts['PKST3_ORN'] = zeros(simlen) # state variable - ORP = ts['PKST3_ORP'] = zeros(simlen) # state variable - ORC = ts['PKST3_ORC'] = zeros(simlen) # state variable - TORN = ts['PKST3_TORN'] = zeros(simlen) # state variable - TORP = ts['PKST3_TORP'] = zeros(simlen) # state variable - TORC = ts['PKST3_TORC'] = zeros(simlen) # state variable - POTBOD = ts['PKST3_POTBOD'] = zeros(simlen) # state variable - - PHYTO = ts['PHYTO'] = zeros(simlen) # concentration - ZOO = ts['ZOO'] = zeros(simlen) # concentration - BENAL = ts['BENAL'] = zeros(simlen) # concentration - PHYCLA = ts['PHYCLA'] = zeros(simlen) # concentration - BALCLA = ts['BALCLA'] = zeros(simlen) # concentration - - ROPHYTO = ts['ROPHYTO'] = zeros(simlen) # total outflow - ROZOO = ts['ROZOO'] = zeros(simlen) # total outflow - ROBENAL = ts['ROBENAL'] = zeros(simlen) # total outflow - ROPHYCLA = ts['ROPHYCLA'] = zeros(simlen) # total outflow - ROBALCLA = ts['ROBALCLA'] = zeros(simlen) # total outflow - - OPHYTO = ts['OPHYTO'] = zeros((simlen, nexits)) # outflow by gate - OZOO = ts['OZOO'] = zeros((simlen, nexits)) # outflow by gate - OBENAL = ts['OBENAL'] = zeros((simlen, nexits)) # outflow by gate - OPHYCLA = ts['OPHYCLA'] = zeros((simlen, nexits)) # outflow by gate - OBALCLA = ts['OBALCLA'] = zeros((simlen, nexits)) # outflow by gate - - BINV = setit() # ts (BINVFG==1), monthly (BINVFG) - PLADFX = setit() # time series, monthly(PLAFXM) - PLADCN = setit() # time series, monthly(PLAFXM) - - plank = pplank() # returns Numba accelerated function in closure - - if PHFG: - # get PHCARB() specific external time series - ALK = ts['CON'] # ALCON only - ITIC = ts['ITIC'] # or zero if not present - ICO2 = ts['ICO2'] # or zero if not present + # create numba dictionaries (empty if not simulated): + ui_oxrx = make_numba_dict(uci_oxrx) + ui_oxrx['errlen'] = len(ERRMSGS_oxrx) - - # preallocate output arrays for speed - PH = ts['PH'] = zeros(simlen) # state variable - TIC = ts['TIC'] = zeros(simlen) # state variable - CO2 = ts['CO2'] = zeros(simlen) # state variable - ROTIC = ts['ROTIC'] = zeros(simlen) # reach total outflow - ROCO2 = ts['ROCO2'] = zeros(simlen) # reach total outflow - OTIC = ts['OTIC'] = zeros((simlen, nexits)) # outflow by exit - OCO2 = ts['OCO2'] = zeros((simlen, nexits)) # outflow by exit - TOTCO2 = ts['TOTCO2'] = zeros(simlen) # ??? computed, but not returned??? - - phcarb = pphcarb() # returns Numba accelerated function in closure - - ############## master simulation loop ##################### + ui_nutrx = Dict.empty(key_type=types.unicode_type, value_type=types.float64) + if NUTFG == 1: + ui_nutrx = make_numba_dict(uci_nutrx) + ui_nutrx['errlen'] = len(ERRMSGS_nutrx) + ui_plank = Dict.empty(key_type=types.unicode_type, value_type=types.float64) + if PLKFG == 1: + ui_plank = make_numba_dict(uci_plank) + ui_plank['errlen'] = len(ERRMSGS_plank) - for loop in range(simlen): - avdepe = AVDEPE[loop] - avvele = AVVELE[loop] - tw = TW[loop] - depcor = DEPCOR[loop] - - advData = nexits, vol, VOL[loop], SROVOL[loop], EROVOL[loop], SOVOL[loop], EOVOL[loop] - - # simulate primary do and bod balances - (dox, bod, satdo, rodo, robod, odo, - obod) = oxrx( - IDOX[loop], IBOD[loop], WIND[loop], avdepe, avvele, tw, depcor, BENRFG, advData) - - - - if NUTFG == 1: # simulate primary inorganic nitrogen and phosphorus balances - (dox, bod, NO3[loop], NO2[loop], NH3[loop], PO4[loop], TAM[loop], RONO3[loop], - RONO2[loop], RONH3[loop], ROPO4[loop], ONO3[loop], ONO2[loop], ONH3[loop], - OPO4[loop]) = nutrx( - dox, bod, tw, INO3[loop], INH3[loop], INO2[loop], IPO4[loop], NUAFX[loop], - NUACN[loop], PREC[loop], SAREA[loop], advData) - - if PLKFG == 1: # simulate plankton populations and associated reactions - (dox, bod, ORN[loop], ORP[loop], ORC[loop], TORN[loop], TORP[loop], - TORC[loop], POTBOD[loop], PHYTO[loop], ZOO[loop], BENAL[loop], - PHYCLA[loop], BALCLA[loop], ROPHYTO[loop], ROZOO[loop], ROBENAL[loop], - ROPHYCLA[loop], ROBALCLA[loop], OPHYTO[loop], OZOO[loop], - OBENAL[loop], OPHYCLA[loop], OBALCLA[loop], BINV[loop], PLADFX[loop], - PLADCN[loop]) = plank( - dox, bod, IPHYTO[loop], IZOO[loop], IORN[loop], IORP[loop], - IORC[loop], tw, WASH[loop], SOLRAD[loop], PREC[loop], SAREA[loop], advData) - - if PHFG == 1: # simulate ph and carbon species + ui_phcarb = Dict.empty(key_type=types.unicode_type, value_type=types.float64) + if PHFG == 1: + ui_phcarb = make_numba_dict(uci_phcarb) + ui_phcarb['errlen'] = len(ERRMSGS_phcarb) + + # hydraulic results: + advectData = uci['advectData'] + (nexits, vol, VOL, SROVOL, EROVOL, SOVOL, EOVOL) = advectData + + ui['nexits'] = nexits + ui['vol'] = vol + + ts['VOL'] = VOL + ts['SROVOL'] = SROVOL + ts['EROVOL'] = EROVOL + + for i in range(nexits): + ts['SOVOL' + str(i + 1)] = SOVOL[:, i] + ts['EOVOL' + str(i + 1)] = EOVOL[:, i] + + #--------------------------------------------------------------------- + # input time series processing (atm. deposition, benthic inverts, etc.) + # TO-DO! - needs testing + #--------------------------------------------------------------------- + # NUTRX atmospheric deposition - initialize time series: + if NUTFG == 1: + for j in range(1, 4): + n = (2 * j) - 1 + + # dry deposition: + nuadfg_dd = int(ui_nutrx['NUADFG' + str(j)]) + NUADFX = zeros(simlen) + + if nuadfg_dd > 0: + NUADFX = initm(siminfo, ui, nuadfg_dd, 'NUTRX_MONTHLY/NUADFX', 0.0) + elif nuadfg_dd == -1: + if 'NUADFX' + str(j) in ts: + NUADFX = ts['NUADFX' + str(j)] + else: + pass #ERRMSG? + ts['NUADFX' + str(j)] = NUADFX + + # wet deposition: + nuadfg_wd = int(ui_nutrx['NUADFG' + str(j+1)]) + NUADCN = zeros(simlen) + + if nuadfg_wd > 0: + NUADCN = initm(siminfo, ui, nuadfg_wd, 'NUTRX_MONTHLY/NUADCN', 0.0) + elif nuadfg_wd == -1: + if 'NUADCN' + str(j) in ts: + NUADCN = ts['NUADCN' + str(j)] + else: + pass #ERRMSG? + ts['NUADCN' + str(j)] = NUADCN + + if PLKFG == 1: + # PLANK atmospheric deposition - initialize time series: + for j in range(1, 4): + n = (2 * j) - 1 + + # dry deposition: + PLADFX = zeros(simlen) + pladfg_dd = int(ui_plank['PLADFG' + str(j)]) + + if pladfg_dd > 0: + PLADFX = initm(siminfo, ui, pladfg_dd, 'PLANK_MONTHLY/PLADFX', 0.0) + elif pladfg_dd == -1: + if 'PLADFX' + str(j) in ts: + PLADFX = ts['PLADFX' + str(j)] + else: + pass #ERRMSG? + ts['PLADFX' + str(j)] = PLADFX + + # wet deposition: + PLADCN = zeros(simlen) + pladfg_wd = int(ui_plank['PLADFG' + str(j+1)]) + + if pladfg_wd > 0: + PLADCN = initm(siminfo, ui, pladfg_wd, 'PLANK_MONTHLY/PLADCN', 0.0) + elif pladfg_wd == -1: + if 'PLADCN' + str(j) in ts: + PLADCN = ts['PLADCN' + str(j)] + else: + pass #ERRMSG? + ts['PLADCN' + str(j)] = PLADCN + + # PLANK - benthic invertebrates: + ts['BINV'] = zeros(simlen) # TO-DO! - needs implementation + + #--------------------------------------------------------------------- + # initialize & run integerated WQ simulation: + #--------------------------------------------------------------------- + + # initialize WQ simulation: + RQUAL = RQUAL_Class(siminfo_, ui, ui_oxrx, ui_nutrx, ui_plank, ui_phcarb, ts) + + # run WQ simulation: + RQUAL.simulate(ts) + + #--------------------------------------------------------------------- + # compile errors & return: + #--------------------------------------------------------------------- + errlen_oxr = len(RQUAL.OXRX.errors) + errlen_ntr = 0; errlen_plk = 0; errlen_phc = 0 + + if NUTFG == 1: + errlen_ntr = len(RQUAL.NUTRX.errors) + if PLKFG == 1: + errlen_plk += len(RQUAL.PLANK.errors) + if PHFG == 1: + #errlen_phc += len(RQUAL.PHCARB.errors) + pass + + errlen = errlen_oxr + errlen_ntr + errlen_plk + errlen_phc + + errors = zeros(errlen, dtype=np.int64) + ERRMSGS = () + + ierr = -1 + for i in range(errlen_oxr): + ierr += 1 + errors[ierr] = RQUAL.OXRX.errors[i] + ERRMSGS += (ERRMSGS_oxrx[i],) + + for i in range(errlen_ntr): + ierr += 1 + errors[ierr] = RQUAL.NUTRX.errors[i] + ERRMSGS += (ERRMSGS_nutrx[i],) + + for i in range(errlen_plk): + ierr += 1 + errors[ierr] = RQUAL.PLANK.errors[i] + ERRMSGS += (ERRMSGS_plank[i],) + + for i in range(errlen_phc): + ierr += 1 + errors[ierr] = RQUAL.PHCARB.errors[i] + ERRMSGS += (ERRMSGS_phcarb[i],) + + return errors, ERRMSGS + + +#------------------------------------------------------------------- +# mass links: +#------------------------------------------------------------------- + +def expand_OXRX_masslinks(flags, uci, dat, recs): + if flags['OXRX']: + for i in range(1,3): + rec = {} + rec['MFACTOR'] = dat.MFACTOR + rec['SGRPN'] = 'OXRX' + + if dat.SGRPN == "ROFLOW": + rec['SMEMN'] = 'OXCF1' + rec['SMEMSB1'] = str(i) # species index + rec['SMEMSB2'] = '' + else: + rec['SMEMN'] = 'OXCF2' + rec['SMEMSB1'] = dat.SMEMSB1 # first sub is exit number + rec['SMEMSB2'] = str(i) # species index - (dox, bod, PH[loop], TIC[loop], CO2[loop], ROTIC[loop], - ROCO2[loop], OTIC[loop], OCO2[loop], TOTCO2[loop], - ) = phcarb( - dox, bod, ALK[loop], ITIC[loop], ICO2[loop], tw, avdepe, SCRFAC[loop], advData) + rec['TMEMN'] = 'OXIF' + rec['TMEMSB1'] = str(i) # species index + rec['TMEMSB2'] = '1' + rec['SVOL'] = dat.SVOL + + recs.append(rec) + + return + +def expand_NUTRX_masslinks(flags, uci, dat, recs): + if flags['NUTRX']: + # dissolved species: + for i in range(1,5): + rec = {} + rec['MFACTOR'] = dat.MFACTOR + rec['SGRPN'] = 'NUTRX' + + if dat.SGRPN == "ROFLOW": + rec['SMEMN'] = 'NUCF1' + rec['SMEMSB1'] = str(i) # species index + rec['SMEMSB2'] = '' + else: + rec['SMEMN'] = 'NUCF9' + rec['SMEMSB2'] = dat.SMEMSB1 # exit number + rec['SMEMSB2'] = str(i) # species index + + rec['TMEMN'] = 'NUIF1' + rec['TMEMSB1'] = str(i) # species index + rec['TMEMSB2'] = '' + rec['SVOL'] = dat.SVOL + recs.append(rec) + + # particulate species (NH4, PO4): + for j in range(1,5): # sediment type - # update totals of nutrients - rno3 = no3 * vol - rtam = tam * vol - rno2 = no2 * vol - rpo4 = po4 * vol - rnh4 = nh4 * vol - rnh3 = nh3 * vol - rrno3 = no3 * vol - rrtam = tam * vol - if ADNHFG == 1: - rrtam += rsnh4(4) # add adsorbed suspended nh4 to dissolved - - rrno2 = no2 * vol - rrpo4 = po4 * vol - if ADPOFG == 1: - rrpo4 += rspo4(4) # add adsorbed suspended po4 to dissolved - - # check do level; if dox exceeds user specified level of supersaturation, then release excess do to the atmosphere - doxs = dox - if dox > supsat * satdo: - dox = supsat * satdo - readox = readox + (dox - doxs) * vol - totdox = readox + boddox + bendox + nitdox + phydox + zoodox + baldox + # adsorbed NH4: + if flags['TAMFG'] and flags['ADNHFG']: + rec = {} + rec['MFACTOR'] = dat.MFACTOR + rec['SGRPN'] = 'NUTRX' + + if dat.SGRPN == "ROFLOW": + rec['SMEMN'] = 'NUCF2' + rec['SMEMSB1'] = str(j) # sediment type + rec['SMEMSB2'] = '1' # NH4 index + else: + rec['SMEMN'] = 'OSNH4' + rec['SMEMSB2'] = dat.SMEMSB1 # exit number + rec['SMEMSB2'] = str(j) # sediment type + + rec['TMEMN'] = 'NUIF2' + rec['TMEMSB1'] = str(j) # sediment type + rec['TMEMSB2'] = '1' # NH4 index + rec['SVOL'] = dat.SVOL + recs.append(rec) + + # adsorbed PO4: + if flags['PO4FG'] and flags['ADPOFG']: + rec = {} + rec['MFACTOR'] = dat.MFACTOR + rec['SGRPN'] = 'NUTRX' + + if dat.SGRPN == "ROFLOW": + rec['SMEMN'] = 'NUCF2' + rec['SMEMSB1'] = str(j) # sediment type + rec['SMEMSB2'] = '2' # PO4 index + else: + rec['SMEMN'] = 'OSPO4' + rec['SMEMSB2'] = dat.SMEMSB1 # exit number + rec['SMEMSB2'] = str(j) # sediment type + + rec['TMEMN'] = 'NUIF2' + rec['TMEMSB1'] = str(j) # sediment type + rec['TMEMSB2'] = '2' + rec['SVOL'] = dat.SVOL + recs.append(rec) + + return + +def expand_PLANK_masslinks(flags, uci, dat, recs): + if flags['PLANK']: - # update dissolved totals and totals of nutrients - rdox = dox * vol - rbod = bod * vol + for i in range(1,6): + rec = {} + rec['MFACTOR'] = dat.MFACTOR + rec['SGRPN'] = 'PLANK' + + if dat.SGRPN == "ROFLOW": + rec['SMEMN'] = 'PKCF1' + rec['SMEMSB1'] = str(i) # species index + rec['SMEMSB2'] = '' + else: + rec['SMEMN'] = 'TPKCF2' + rec['SMEMSB2'] = dat.SMEMSB1 # exit number + rec['SMEMSB2'] = str(i) # species index + + rec['TMEMN'] = 'PKIF' + rec['TMEMSB1'] = str(i) #dat.TMEMSB1 + rec['TMEMSB2'] = '' + rec['SVOL'] = dat.SVOL + recs.append(rec) + return +def expand_PHCARB_masslinks(flags, uci, dat, recs): + + if flags['PHCARB']: + + for i in range(1,3): + rec = {} + rec['MFACTOR'] = dat.MFACTOR + rec['SGRPN'] = 'PHCARB' + + if dat.SGRPN == "ROFLOW": + rec['SMEMN'] = 'PHCF1' + rec['SMEMSB1'] = str(i) # species index + rec['SMEMSB2'] = '' + else: + rec['SMEMN'] = 'PHCF2' + rec['SMEMSB2'] = dat.SMEMSB1 # exit number + rec['SMEMSB2'] = str(i) # species index + + rec['TMEMN'] = 'PHIF' + rec['TMEMSB1'] = str(i) # species index + rec['TMEMSB2'] = '' + rec['SVOL'] = dat.SVOL + recs.append(rec) -# #@jit(nopython=True) -# def benth (dox, anaer, BRCON, scrfac, depcor, conc): -# ''' simulate benthal release of constituent''' -# # calculate benthal release of constituent; release is a step function of aerobic/anaerobic conditions, and stream velocity; -# # scrfac, the scouring factor dependent on stream velocity and depcor, the conversion factor from mg/m2 to mg/l, -# # both calculated in rqual; releas is expressed in mg/m2.ivl -# releas = BRCON[0] * scrfac * depcor if dox > anaer else BRCON[1] * scrfac * depcor -# conc += releas -# return conc, releas - - -# #@jit(nopython=True) -# def decbal(TAMFG, PO4FG, decnit, decpo4, tam, no3, po4): -# ''' perform materials balance for transformation from organic to inorganic material by decay in reach water''' -# if TAMFG: -# tam += decnit # add nitrogen transformed to inorganic nitrogen by biomass decomposition -# else: -# no3 += decnit # add nitrogen transformed to inorganic nitrogen by biomass decomposition -# if PO4FG: # add phosphorus transformed to inorganic phosphorus by biomass decomposition to po4 state variable -# po4 += decpo4 -# return tam, no3, po4 - - -# #@jit(nopython=True) -# def sink (vol, avdepe, kset, conc, snkmat): -# ''' calculate quantity of material settling out of the control volume; determine the change in concentration as a result of sinking''' -# if kset > 0.0 and avdepe > 0.17: -# # calculate concentration change due to outgoing material; snkout is expressed in mass/liter/ivl; kset is expressed as ft/ivl and avdepe as feet -# snkout = conc * (kset / avdepe) if kset < avdepe else conc # calculate portion of material which settles out of the control volume during time step; snkout is expressed as mass/liter.ivl; conc is the concentration of material in the control volume -# conc -= snkout # calculate remaining concentration of material in the control volume -# snkmat = snkout * vol # find quantity of material that sinks out; units are mass.ft3/l.ivl in english system, and mass.m3/l.ivl in metric system -# else: -# snkout = 0.0 -# snkmat = 0.0 -# return conc, snkmat + return \ No newline at end of file diff --git a/HSP2/RQUAL_Class.py b/HSP2/RQUAL_Class.py new file mode 100644 index 00000000..0cade940 --- /dev/null +++ b/HSP2/RQUAL_Class.py @@ -0,0 +1,801 @@ +import numpy as np +from numpy import where, zeros, array +from math import log +import numba as nb +from numba.experimental import jitclass + +from HSP2.OXRX_Class import OXRX_Class +from HSP2.NUTRX_Class import NUTRX_Class +from HSP2.PLANK_Class import PLANK_Class +#from HSP2.PHCARB_Class import PHCARB_Class +from HSP2.utilities import make_numba_dict, initm + +spec = [ + ('OXRX', OXRX_Class.class_type.instance_type), + ('NUTRX', NUTRX_Class.class_type.instance_type), + ('PLANK', PLANK_Class.class_type.instance_type), + ('AFACT', nb.float64), + ('ALK', nb.float64[:]), + ('AVDEP', nb.float64[:]), + ('AVDEPE', nb.float64[:]), + ('AVVEL', nb.float64[:]), + ('AVVELE', nb.float64[:]), + ('BALCLA', nb.float64[:]), + ('BENAL1', nb.float64[:]), + ('BENRFG', nb.int32), + ('BOD', nb.float64[:]), + ('CO2', nb.float64[:]), + ('delt60', nb.float64), + ('delts', nb.float64), + ('DEPCOR', nb.float64[:]), + ('DOX', nb.float64[:]), + ('EOVOL', nb.float64[:,:]), + ('EROVOL', nb.float64[:]), + ('IBOD', nb.float64[:]), + ('ICO2', nb.float64[:]), + ('IDOX', nb.float64[:]), + ('INH4', nb.float64[:]), + ('INO2', nb.float64[:]), + ('INO3', nb.float64[:]), + ('IORC', nb.float64[:]), + ('IORN', nb.float64[:]), + ('IORP', nb.float64[:]), + ('IPHYT', nb.float64[:]), + ('IPO4', nb.float64[:]), + ('ISNH41', nb.float64[:]), + ('ISNH42', nb.float64[:]), + ('ISNH43', nb.float64[:]), + ('ISPO41', nb.float64[:]), + ('ISPO42', nb.float64[:]), + ('ISPO43', nb.float64[:]), + ('ITIC', nb.float64[:]), + ('IZOO', nb.float64[:]), + ('LKFG', nb.int32), + ('nexits', nb.int32), + ('NH3', nb.float64[:]), + ('NH4', nb.float64[:]), + ('NO2', nb.float64[:]), + ('NO3', nb.float64[:]), + ('NUADCN', nb.float64[:]), + ('NUADFG', nb.int32[:]), + ('NUADFX', nb.float64[:]), + ('NUTFG', nb.int32), + ('OBALCLA', nb.float64[:,:]), + ('OBENAL', nb.float64[:,:]), + ('OBOD', nb.float64[:,:]), + ('OCO2', nb.float64[:,:]), + ('ODOX', nb.float64[:,:]), + ('ONO2', nb.float64[:,:]), + ('ONO3', nb.float64[:,:]), + ('OPHYCLA', nb.float64[:,:]), + ('OPHYT', nb.float64[:,:]), + ('OPO4', nb.float64[:,:]), + ('OORC', nb.float64[:,:]), + ('OORN', nb.float64[:,:]), + ('OORP', nb.float64[:,:]), + ('ORC', nb.float64[:]), + ('ORN', nb.float64[:]), + ('ORP', nb.float64[:]), + ('OTAM', nb.float64[:,:]), + ('OTIC', nb.float64[:,:]), + ('OZOO', nb.float64[:,:]), + ('PH', nb.float64[:]), + ('PHFG', nb.int32), + ('PHYCLA', nb.float64[:]), + ('PHYTO', nb.float64[:]), + ('PKIF1', nb.float64[:]), + ('PKIF2', nb.float64[:]), + ('PKIF3', nb.float64[:]), + ('PKIF4', nb.float64[:]), + ('PKIF5', nb.float64[:]), + ('PLKFG', nb.int32), + ('PO4', nb.float64[:]), + ('POTBOD', nb.float64[:]), + ('PREC', nb.float64[:]), + ('RNH3', nb.float64[:]), + ('RNH4', nb.float64[:]), + ('RNO2', nb.float64[:]), + ('RNO3', nb.float64[:]), + ('RO', nb.float64[:]), + ('ROBALCLA', nb.float64[:]), + ('ROBENAL', nb.float64[:]), + ('ROBOD', nb.float64[:]), + ('ROCO2', nb.float64[:]), + ('RODOX', nb.float64[:]), + ('RONO2', nb.float64[:]), + ('RONO3', nb.float64[:]), + ('ROORC', nb.float64[:]), + ('ROORN', nb.float64[:]), + ('ROORP', nb.float64[:]), + ('ROPHYCLA', nb.float64[:]), + ('ROPHYT', nb.float64[:]), + ('ROPO4', nb.float64[:]), + ('ROSNH41', nb.float64[:]), + ('ROSNH42', nb.float64[:]), + ('ROSNH43', nb.float64[:]), + ('ROSPO41', nb.float64[:]), + ('ROSPO42', nb.float64[:]), + ('ROSPO43', nb.float64[:]), + ('ROTAM', nb.float64[:]), + ('ROTIC', nb.float64[:]), + ('ROTORC', nb.float64[:]), + ('ROTORN', nb.float64[:]), + ('ROTORP', nb.float64[:]), + ('ROZOO', nb.float64[:]), + ('RPO4', nb.float64[:]), + ('RTAM', nb.float64[:]), + ('SAREA', nb.float64[:]), + ('SATDO', nb.float64[:]), + ('SCRFAC', nb.float64[:]), + ('SEDFG', nb.int32), + ('simlen', nb.int32), + ('SOLRAD', nb.float64[:]), + ('SOVOL', nb.float64[:,:]), + ('SSED4', nb.float64[:]), + ('SROVOL', nb.float64[:]), + ('svol', nb.float64), + ('TAM', nb.float64[:]), + ('TBENAL1', nb.float64[:]), + ('TBENAL2', nb.float64[:]), + ('TIC', nb.float64[:]), + ('TN', nb.float64[:]), + ('TORC', nb.float64[:]), + ('TORN', nb.float64[:]), + ('TORP', nb.float64[:]), + ('TOTCO2', nb.float64[:]), + ('TP', nb.float64[:]), + ('TW', nb.float64[:]), + ('uunits', nb.int32), + ('vol', nb.float64), + ('VOL', nb.float64[:]), + ('WASH', nb.float64[:]), + ('WIND', nb.float64[:]), + ('ZOO', nb.float64[:]) +] + +@jitclass(spec) +class RQUAL_Class: + + #------------------------------------------------------------------- + # class initialization: + #------------------------------------------------------------------- + def __init__(self, siminfo, ui, ui_oxrx, ui_nutrx, ui_plank, ui_phcarb, ts): + + ''' Initialize instance variables for rqual wrapper ''' + + # simulation data: + delt60 = siminfo['delt'] / 60.0 # delt60 - simulation time interval in hours + simlen = int(siminfo['steps']) + + self.delt60 = delt60 + self.simlen = simlen + self.delts = siminfo['delt'] * 60 + self.uunits = int(siminfo['units']) + + # hydaulic results: + #(nexits, vol, VOL, self.SROVOL, self.EROVOL, self.SOVOL, self.EOVOL) = advectData + + self.AFACT = 43560.0 + if self.uunits == 2: + # si units conversion + self.AFACT = 1000000.0 + + nexits = int(ui['nexits']) + self.nexits = nexits + + self.VOL = ts['VOL'] * self.AFACT + self.SROVOL = ts['SROVOL'] + self.EROVOL = ts['EROVOL'] + self.SOVOL = zeros((self.simlen, nexits)) + self.EOVOL = zeros((simlen, nexits)) + for i in range(nexits): + self.SOVOL[:, i] = ts['SOVOL' + str(i + 1)] + self.EOVOL[:, i] = ts['EOVOL' + str(i + 1)] + + self.vol = ui['vol'] * self.AFACT + self.svol = self.vol + + # initialize flags: + self.BENRFG = int(ui['BENRFG']) # table-type benth-flag + + # table type ACTIVITY + self.NUTFG = int(ui['NUTFG']) + self.PLKFG = int(ui['PLKFG']) + self.PHFG = int(ui['PHFG']) + + self.SEDFG = int(ui['SEDFG']) + self.LKFG = int(ui['LKFG']) + + #TO-DO! - remove hardcode once PHCARB class is ready for testing + self.PHFG = 0 + + # get external time series + self.PREC = ts['PREC'] + self.SAREA = ts['SAREA'] #dbg * self.AFACT + + self.RO = ts['RO'] + self.AVDEP = ts['AVDEP'] + self.AVVEL = ts['AVVEL'] + self.TW = ts['TW'] + + if 'WIND' in ts: + self.WIND = ts['WIND'] * 1609.0 # miles/ivld to meters/ivld + else: + self.WIND = zeros(simlen) + + # initialize time series for physical variables: + self.AVDEPE = zeros(simlen) + self.AVVELE = zeros(simlen) + self.DEPCOR = zeros(simlen) + self.SCRFAC = zeros(simlen) + + for t in range(simlen): + # fix undefined temps if present + self.TW[t] = 20.0 if self.TW[t] < -100.0 else self.TW[t] + + # convert depth/velocity to english units, if necessary: + self.AVDEPE[t] = self.AVDEP[t] + self.AVVELE[t] = self.AVVEL[t] + + if self.uunits == 2: + self.AVDEPE[t] *= 3.28 + self.AVVELE[t] *= 3.28 + + # define conversion factor from mg/m2 to mg/l + if self.AVDEPE[t] > 0.0: + self.DEPCOR[t] = 3.28084e-3 / self.AVDEPE[t] + else: + self.DEPCOR[t] = -1.0e30 + + # calculate scouring factor + if self.BENRFG == 1: + scrvel = ui['SCRVEL'] # table-type scour-parms + scrmul = ui['SCRMUL'] # table-type scour-parms + + if self.AVVELE[t] > scrvel: + self.SCRFAC[t] = scrmul + else: + self.SCRFAC[t] = 1.0 + + #------------------------------------------------------- + # OXRX - initialize: + #------------------------------------------------------- + + if ('OXIF1' in ts): + self.IDOX = ts['OXIF1'] # optional, input flow + else: + self.IDOX = zeros(simlen) + + if ('OXIF2' in ts): + self.IBOD = ts['OXIF2'] # optional, input flow + else: + self.IBOD = zeros(simlen) + + # OXRX - instantiate: + self.OXRX = OXRX_Class(siminfo, self.nexits, self.vol, ui, ui_oxrx, ts) + + # OXRX - preallocate arrays for computed time series: + self.DOX = ts['DOX'] = zeros(simlen) # concentration, state variable + self.BOD = ts['BOD'] = zeros(simlen) # concentration, state variable + self.SATDO = ts['SATDO'] = zeros(simlen) # concentration, state variable + + self.RODOX = ts['OXCF11'] = zeros(simlen) # reach outflow of DOX + self.ROBOD = ts['OXCF12'] = zeros(simlen) # reach outflow of BOD + + if nexits > 1: + for i in range(nexits): + ts['OXCF2' + str(i + 1) + ' 1'] = zeros(simlen) # DOX outflow by exit + ts['OXCF2' + str(i + 1) + ' 2'] = zeros(simlen) # BOD outflow by exit + + #------------------------------------------------------- + # NUTRX - initialize: + #------------------------------------------------------- + if self.NUTFG == 1: + + # nutrient inflows: + self.INO3 = zeros(simlen); self.INO2 = zeros(simlen) + self.INH4 = zeros(simlen); self.IPO4 = zeros(simlen) + + if 'NUIF11' in ts: + self.INO3 = ts['NUIF11'] # optional, input + + if 'NUIF12' in ts: + self.INH4 = ts['NUIF12'] # optional, input + + if 'NUIF13' in ts: + self.INO2 = ts['NUIF13'] # optional, input + + if 'NUIF14' in ts: + self.IPO4 = ts['NUIF14'] # optional, input + + # sediment-adsorbed (NH4, PO4): + self.ISNH41 = zeros(simlen) + self.ISNH42 = zeros(simlen) + self.ISNH43 = zeros(simlen) + self.ISPO41 = zeros(simlen) + self.ISPO42 = zeros(simlen) + self.ISPO43 = zeros(simlen) + + if 'NUIF21 1' in ts: self.ISNH41 = ts['NUIF21 1'] + if 'NUIF22 1' in ts: self.ISNH42 = ts['NUIF22 1'] + if 'NUIF23 1' in ts: self.ISNH43 = ts['NUIF23 1'] + + if 'NUIF21 2' in ts: self.ISPO41 = ts['NUIF21 2'] + if 'NUIF22 2' in ts: self.ISPO42 = ts['NUIF22 2'] + if 'NUIF23 2' in ts: self.ISPO43 = ts['NUIF23 2'] + + # NUTRX - instantiate class: + self.NUTRX = NUTRX_Class(siminfo, self.nexits, self.vol, ui, ui_nutrx, ts, self.OXRX) + + # NUTRX - preallocate storage for computed time series + self.NO3 = ts['NO3'] = zeros(simlen) # concentration, state variable + self.TAM = ts['TAM'] = zeros(simlen) # concentration, state variable + self.NO2 = ts['NO2'] = zeros(simlen) # concentration, state variable + self.PO4 = ts['PO4'] = zeros(simlen) # concentration, state variable + self.NH4 = ts['NH4'] = zeros(simlen) # concentration, state variable + self.NH3 = ts['NH3'] = zeros(simlen) # concentration, state variable + + # inflows: + #self.NUIF11 = ts['NUIF1_NO3'] = zeros(simlen) # total outflow + #self.NUIF12 = ts['NUIF1_TAM'] = zeros(simlen) # total outflow + #self.NUIF13 = ts['NUIF1_NO2'] = zeros(simlen) # total outflow + #self.NUIF14 = ts['NUIF1_PO4'] = zeros(simlen) # total outflow + + # total outflows: + self.RONO3 = ts['NUCF11'] = zeros(simlen) # outflow + self.ROTAM = ts['NUCF12'] = zeros(simlen) # outflow + self.RONO2 = ts['NUCF13'] = zeros(simlen) # outflow + self.ROPO4 = ts['NUCF14'] = zeros(simlen) # outflow + + if self.NUTRX.ADNHFG > 0: + self.ROSNH41 = ts['NUCF21 1'] = zeros(simlen) # sand + self.ROSNH42 = ts['NUCF22 1'] = zeros(simlen) # silt + self.ROSNH43 = ts['NUCF23 1'] = zeros(simlen) # clay + + if self.NUTRX.ADPOFG > 0: + self.ROSPO41 = ts['NUCF21 2'] = zeros(simlen) # sand + self.ROSPO42 = ts['NUCF22 2'] = zeros(simlen) # silt + self.ROSPO43 = ts['NUCF23 2'] = zeros(simlen) # clay + + # exit outflows: + if nexits > 1: + for i in range(nexits): + ts['NUCF9' + str(i + 1) + ' 1'] = zeros(simlen) + ts['NUCF9' + str(i + 1) + ' 2'] = zeros(simlen) + ts['NUCF9' + str(i + 1) + ' 3'] = zeros(simlen) + ts['NUCF9' + str(i + 1) + ' 4'] = zeros(simlen) + + if self.NUTRX.ADNHFG > 0: + ts['OSNH4' + str(i + 1) + ' 1'] = zeros(simlen) # sand + ts['OSNH4' + str(i + 1) + ' 2'] = zeros(simlen) # silt + ts['OSNH4' + str(i + 1) + ' 3'] = zeros(simlen) # clay + + if self.NUTRX.ADPOFG > 0: + ts['OSPO4' + str(i + 1) + ' 1'] = zeros(simlen) # sand + ts['OSPO4' + str(i + 1) + ' 2'] = zeros(simlen) # silt + ts['OSPO4' + str(i + 1) + ' 3'] = zeros(simlen) # clay + + self.RNO3 = ts['RNO3'] = zeros(simlen) + self.RTAM = ts['RTAM'] = zeros(simlen) + self.RNO2 = ts['RNO2'] = zeros(simlen) + self.RPO4 = ts['RPO4'] = zeros(simlen) + self.RNH4 = ts['RNH4'] = zeros(simlen) + self.RNH3 = ts['RNH3'] = zeros(simlen) + + #------------------------------------------------------- + # PLANK - simulate biological components: + #------------------------------------------------------- + if self.PLKFG == 1: + + # get PLANK specific time series + self.IPHYT = zeros(simlen); self.IZOO = zeros(simlen) + self.IORN = zeros(simlen); self.IORP = zeros(simlen); self.IORC = zeros(simlen) + + if 'PKIF1' in ts: self.IPHYT = ts['PKIF1'] # optional input + if 'PKIF2' in ts: self.IZOO = ts['PKIF2'] # optional input + if 'PKIF3' in ts: self.IORN = ts['PKIF3'] # optional input + if 'PKIF4' in ts: self.IORP = ts['PKIF4'] # optional input + if 'PKIF5' in ts: self.IORC = ts['PKIF5'] # optional input + + self.WASH = zeros(simlen); self.SOLRAD = zeros(simlen) + if 'WASH' in ts: self.WASH = ts['WASH'] + if 'SOLRAD' in ts: self.SOLRAD = ts['SOLRAD'] + + # total suspended sediment conc: + self.SSED4 = zeros(simlen) + if 'SSED4' in ts: self.SSED4 = ts['SSED4'] + + # PLANK - instantiate class: + self.PLANK = PLANK_Class(siminfo, self.nexits, self.vol, ui, ui_plank, ts, self.OXRX, self.NUTRX) + + # PLANK - preallocate storage for computed time series + + self.PHYTO = ts['PHYTO'] = zeros(simlen) # concentration + self.ZOO = ts['ZOO'] = zeros(simlen) # concentration + self.BENAL1 = ts['BENAL1'] = zeros(simlen) # concentration + self.TBENAL1= ts['TBENAL1'] = zeros(simlen) # concentration + self.TBENAL2= ts['TBENAL2'] = zeros(simlen) # concentration + self.PHYCLA = ts['PHYCLA'] = zeros(simlen) # concentration + + self.ORN = ts['ORN'] = zeros(simlen) # state variable + self.ORP = ts['ORP'] = zeros(simlen) # state variable + self.ORC = ts['ORC'] = zeros(simlen) # state variable + self.TORN = ts['TORN'] = zeros(simlen) # state variable + self.TORP = ts['TORP'] = zeros(simlen) # state variable + self.TORC = ts['TORC'] = zeros(simlen) # state variable + self.POTBOD = ts['POTBOD'] = zeros(simlen) # state variable + self.TN = ts['TN'] = zeros(simlen) # state variable + self.TP = ts['TP'] = zeros(simlen) # state variable + + # inflows: + self.PKIF1 = ts['PKIF_PHYT'] = zeros(simlen) # total outflow + self.PKIF2 = ts['PKIF_ZOO'] = zeros(simlen) # total outflow + self.PKIF3 = ts['PKIF_ORN'] = zeros(simlen) # total outflow + self.PKIF4 = ts['PKIF_ORP'] = zeros(simlen) # total outflow + self.PKIF5 = ts['PKIF_ORC'] = zeros(simlen) # total outflow + + # outflows: + self.ROPHYT = ts['PKCF11'] = zeros(simlen) # total outflow + self.ROZOO = ts['PKCF12'] = zeros(simlen) # total outflow + self.ROORN = ts['PKCF13'] = zeros(simlen) # total outflow + self.ROORP = ts['PKCF14'] = zeros(simlen) # total outflow + self.ROORC = ts['PKCF15'] = zeros(simlen) # total outflow + + self.ROTORN = ts['ROTORN'] = zeros(simlen) # total outflow + self.ROTORP = ts['ROTORP'] = zeros(simlen) # total outflow + self.ROTORC = ts['ROTORC'] = zeros(simlen) # total outflow + #self.ROTN = ts['ROTN'] = zeros(simlen) # total outflow + #self.ROTP = ts['ROTP'] = zeros(simlen) # total outflow + + if nexits > 1: + for i in range(nexits): + ts['PKCF2' + str(i + 1) + ' 1'] = zeros(simlen) # OPHYT + ts['PKCF2' + str(i + 1) + ' 2'] = zeros(simlen) # OZOO + ts['PKCF2' + str(i + 1) + ' 3'] = zeros(simlen) # OORN + ts['PKCF2' + str(i + 1) + ' 4'] = zeros(simlen) # OORP + ts['PKCF2' + str(i + 1) + ' 5'] = zeros(simlen) # OORC + + #------------------------------------------------------- + # PHCARB - initialize: + #------------------------------------------------------- + if self.PHFG == 1: + + # get PHCARB() specific external time series + self.ALK = zeros(simlen) + self.ITIC = zeros(simlen) + self.ICO2 = zeros(simlen) + + if 'CON' in ts: self.ALK = ts['CON'] + if 'ITIC' in ts: self.ITIC = ts['ITIC'] + if 'ICO2' in ts: self.CO2 = ts['ICO2'] + + # preallocate output arrays for speed + self.PH = ts['PH'] = zeros(simlen) # state variable + self.TIC = ts['TIC'] = zeros(simlen) # state variable + self.CO2 = ts['CO2'] = zeros(simlen) # state variable + self.ROTIC = ts['ROTIC'] = zeros(simlen) # reach total outflow + self.ROCO2 = ts['ROCO2'] = zeros(simlen) # reach total outflow + self.OTIC = zeros((simlen, nexits)) # outflow by exit + self.OCO2 = zeros((simlen, nexits)) # outflow by exit + self.TOTCO2 = ts['TOTCO2'] = zeros(simlen) # ??? computed, but not returned??? + + for i in range(nexits): + ts['OTIC' + str(i + 1)] = zeros(simlen) + ts['OCO2' + str(i + 1)] = zeros(simlen) + + #PHCARB = PHCARB_Class(siminfo, self.nexits, self.vol, ui, ui_phcarb, ts) + + return + + def simulate(self, ts): + + for loop in range(self.simlen): + + #------------------------------------------------------- + # define inputs for current step: + #------------------------------------------------------- + ro = self.RO[loop] + avdepe = self.AVDEPE[loop] + avvele = self.AVVELE[loop] + scrfac = self.SCRFAC[loop] + + tw = self.TW[loop] + if self.uunits == 1: + tw = (tw - 32.0) * (5.0 / 9.0) + + wind = self.WIND[loop] + wind_r = wind + if self.LKFG == 0: wind_r = 0 + + depcor = self.DEPCOR[loop] + + svol = self.vol + self.svol = svol + + self.vol = self.VOL[loop] + advData = self.nexits, self.svol, self.vol, self.SROVOL[loop], self.EROVOL[loop], self.SOVOL[loop], self.EOVOL[loop] + + idox = self.IDOX[loop] + ibod = self.IBOD[loop] + + # define initial CO2 conentration value for NUTRX: #TO-DO! - needs implementation + co2 = -999.0 + #co2 = PHCARB.co2 + + # define initial pH concentration (for use in NUTRX): #TO-DO! - needs implementation + phval = 7.0 + + #if PHFG == 1: + # if PHFLAG == 1: + # phval = PHCARB.ph + + #------------------------------------------------------- + # OXRX - simulate do and bod balances: + #------------------------------------------------------- + self.OXRX.simulate(idox, ibod, wind_r, scrfac, avdepe, avvele, depcor, tw, advData) + + # initialize DO/BOD process quantities: + nitdox = 0.0 + denbod = 0.0 + phydox = 0.0 + zoodox = 0.0 + baldox = 0.0 + + #------------------------------------------------------- + # NUTRX - simulate primary nutrient (N/P) balances: + #------------------------------------------------------- + if self.NUTFG == 1: + + # sediment results: + depscr = zeros(5) + rosed = zeros(5) + osed = zeros((self.nexits,5)) + + if self.SEDFG == 1: + for j in range(1, 5): + rosed[j] = ts['ROSED' + str(j)][loop] + depscr[j] = ts['DEPSCR' + str(j)][loop] + + if self.nexits > 1: + for i in range(self.nexits): + osed[i,j] = ts['OSED' + str(i+1)][loop] + else: + osed[0,j] = rosed[j] + + # sediment-associated nutrient inflows: + isnh4 = zeros(5) + isnh4[1] = self.ISNH41[loop] + isnh4[2] = self.ISNH42[loop] + isnh4[3] = self.ISNH43[loop] + + for j in range(1,4): + isnh4[4] += isnh4[j] + + ispo4 = zeros(5) + ispo4[1] = self.ISPO41[loop] + ispo4[2] = self.ISPO42[loop] + ispo4[3] = self.ISPO43[loop] + + for j in range(1,4): + ispo4[4] += ispo4[j] + + # simulate nutrients: + self.OXRX = self.NUTRX.simulate(loop, tw, wind, phval, self.OXRX, + self.INO3[loop], self.INH4[loop], self.INO2[loop], self.IPO4[loop], isnh4, ispo4, + self.PREC[loop], self.SAREA[loop], scrfac, avdepe, depcor, depscr, rosed, osed, advData) + + + # update DO / BOD totals: + nitdox = self.NUTRX.nitdox + denbod = self.NUTRX.denbod + + #------------------------------------------------------- + # PLANK - simulate plankton components & associated reactions + #------------------------------------------------------- + if self.PLKFG == 1: + + co2 = 0.0 + #if self.PHFG == 1: co2 = PHCARB.co2 #TO-DO! + + (self.OXRX, self.NUTRX) \ + = self.PLANK.simulate(tw, phval, co2, self.SSED4[loop], self.OXRX, self.NUTRX, + self.IPHYT[loop], self.IZOO[loop], + self.IORN[loop], self.IORP[loop], self.IORC[loop], + self.WASH[loop], self.SOLRAD[loop], self.PREC[loop], + self.SAREA[loop], avdepe, avvele, depcor, ro, advData) + + phydox = self.PLANK.phydox + zoodox = self.PLANK.zoodox + baldox = self.PLANK.baldox + + #------------------------------------------------------- + # PHCARB - simulate: (TO-DO! - needs class implementation) + #------------------------------------------------------- + if self.PHFG == 1: + #self.PHCARB.simulate() + + # update pH and CO2 concentration for use in NUTRX/PLANK: + #if ui_nutrx['PHFLAG'] == 1: + # phval = PHCARB.ph + + #co2 = PHCARB.co2 + i = 0 + pass + + #------------------------------------------------------- + # NUTRX - update masses (TO-DO! - removed for now; fails on numba compile) + #------------------------------------------------------- + ''' + self.NUTRX.rno3 = self.NUTRX.no3 * self.vol + self.NUTRX.rtam = self.NUTRX.tam * self.vol + self.NUTRX.rno2 = self.NUTRX.no2 * self.vol + self.NUTRX.rpo4 = self.NUTRX.po4 * self.vol + self.NUTRX.rnh4 = self.NUTRX.nh4 * self.vol + self.NUTRX.rnh3 = self.NUTRX.nh3 * self.vol + + self.NUTRX.rrno3 = self.NUTRX.no3 * self.vol + self.NUTRX.rrtam = self.NUTRX.tam * self.vol + + if self.NUTRX.ADNHFG == 1: + self.NUTRX.rrtam += self.NUTRX.rsnh4[4] # add adsorbed suspended nh4 to dissolved + + self.NUTRX.rrno2 = self.NUTRX.no2 * self.vol + self.NUTRX.rrpo4 = self.NUTRX.po4 * self.vol + + if self.NUTRX.ADPOFG == 1: + self.NUTRX.rrpo4 += self.NUTRX.rspo4[4] # add adsorbed suspended po4 to dissolved + ''' + + #self.NUTRX.update_mass() + + #------------------------------------------------------- + # OXRX - finalize DO and calculate totals + #------------------------------------------------------- + + # check do level; if dox exceeds user specified level of supersaturation, then release excess do to the atmosphere + dox = self.OXRX.dox + doxs = dox + + if dox > (self.OXRX.supsat * self.OXRX.satdo): + dox = self.OXRX.supsat * self.OXRX.satdo + + self.OXRX.readox += (dox - doxs) * self.vol + self.OXRX.totdox = self.OXRX.readox + self.OXRX.boddox + self.OXRX.bendox \ + + nitdox + phydox + zoodox + baldox + + self.OXRX.dox = dox + self.OXRX.rdox = self.OXRX.dox * self.vol + self.OXRX.rbod = self.OXRX.bod * self.vol + + #self.OXRX.adjust_dox(nitdox, denbod, phydox, zoodox, baldox) + + # udate initial volume for next step: + #self.svol = self.vol + + #------------------------------------------------------- + # Store time series results (all WQ modules): + #------------------------------------------------------- + + # OXRX results: + self.DOX[loop] = self.OXRX.dox + self.BOD[loop] = self.OXRX.bod + self.SATDO[loop] = self.OXRX.satdo + + # outflows (convert to mass per interval (lb/ivld or kg/ivld)) + self.RODOX[loop] = self.OXRX.rodox * self.OXRX.conv + self.ROBOD[loop] = self.OXRX.robod * self.OXRX.conv + + if self.nexits > 1: + for i in range(self.nexits): + ts['OXCF2' + str(i + 1) + '1'][loop] = self.OXRX.odox[i] * self.OXRX.conv + ts['OXCF2' + str(i + 1) + '2'][loop] = self.OXRX.obod[i] * self.OXRX.conv + + # NUTRX results: + if self.NUTFG == 1: + self.NO3[loop] = self.NUTRX.no3 + self.TAM[loop] = self.NUTRX.tam + self.NO2[loop] = self.NUTRX.no2 + self.PO4[loop] = self.NUTRX.po4 + self.NH4[loop] = self.NUTRX.nh4 + self.NH3[loop] = self.NUTRX.nh3 + + # inflows (lb/ivld or kg/ivld): + + # outflows (convert to mass per interval (lb/ivld or kg/ivld)) + conv = self.NUTRX.conv + self.RONO3[loop] = self.NUTRX.rono3 * conv + self.ROTAM[loop] = self.NUTRX.rotam * conv + self.RONO2[loop] = self.NUTRX.rono2 * conv + self.ROPO4[loop] = self.NUTRX.ropo4 * conv + + if self.NUTRX.ADNHFG > 0: + self.ROSNH41[loop] = self.NUTRX.rosnh4[1] * conv + self.ROSNH42[loop] = self.NUTRX.rosnh4[2] * conv + self.ROSNH43[loop] = self.NUTRX.rosnh4[3] * conv + + if self.NUTRX.ADPOFG > 0: + self.ROSPO41[loop] = self.NUTRX.rospo4[1] * conv + self.ROSPO42[loop] = self.NUTRX.rospo4[2] * conv + self.ROSPO43[loop] = self.NUTRX.rospo4[3] * conv + + # exit outflows: + if self.nexits > 1: + for i in range(self.nexits): + ts['NUCF9' + str(i + 1) + ' 1'][loop] = self.NUTRX.ono3[i] * conv + ts['NUCF9' + str(i + 1) + ' 2'][loop] = self.NUTRX.otam[i] * conv + ts['NUCF9' + str(i + 1) + ' 3'][loop] = self.NUTRX.ono2[i] * conv + ts['NUCF9' + str(i + 1) + ' 4'][loop] = self.NUTRX.opo4[i] * conv + + if self.NUTRX.ADNHFG > 0: + ts['OSNH4' + str(i + 1) + ' 1'][loop] = self.NUTRX.osnh4[i,1] * conv # sand + ts['OSNH4' + str(i + 1) + ' 2'][loop] = self.NUTRX.osnh4[i,2] * conv # silt + ts['OSNH4' + str(i + 1) + ' 3'][loop] = self.NUTRX.osnh4[i,3] * conv # clay + + if self.NUTRX.ADPOFG > 0: + ts['OSPO4' + str(i + 1) + ' 1'][loop] = self.NUTRX.ospo4[i,1] * conv # sand + ts['OSPO4' + str(i + 1) + ' 2'][loop] = self.NUTRX.ospo4[i,2] * conv # silt + ts['OSPO4' + str(i + 1) + ' 3'][loop] = self.NUTRX.ospo4[i,3] * conv # clay + + # mass storages: + self.RNO3[loop] = self.NUTRX.no3 * self.vol + self.RTAM[loop] = self.NUTRX.tam * self.vol + self.RNO2[loop] = self.NUTRX.no2 * self.vol + self.RPO4[loop] = self.NUTRX.po4 * self.vol + self.RNH4[loop] = self.NUTRX.nh4 * self.vol + self.RNH3[loop] = self.NUTRX.nh3 * self.vol + + + # PLANK results: + if self.PLKFG == 1: + + self.PHYTO[loop] = self.PLANK.phyto + self.ZOO[loop] = self.PLANK.zoo + self.BENAL1[loop] = self.PLANK.benal[0] + self.TBENAL1[loop] = self.PLANK.tbenal[1] + self.TBENAL2[loop] = self.PLANK.tbenal[2] + self.PHYCLA[loop] = self.PLANK.phycla + + self.ORN[loop] = self.PLANK.orn + self.ORP[loop] = self.PLANK.orp + self.ORC[loop] = self.PLANK.orc + self.TORN[loop] = self.PLANK.torn + self.TORP[loop] = self.PLANK.torp + self.TORC[loop] = self.PLANK.torc + self.POTBOD[loop] = self.PLANK.potbod + self.TN[loop] = self.PLANK.tn + self.TP[loop] = self.PLANK.tp + + # inflows (lb/ivld or kg/ivld): + self.PKIF1[loop] = self.PLANK.iphyto + self.PKIF2[loop] = self.PLANK.izoo + self.PKIF3[loop] = self.PLANK.iorn + self.PKIF4[loop] = self.PLANK.iorp + self.PKIF5[loop] = self.PLANK.iorc + + # outflows (convert to mass per interval (lb/ivld or kg/ivld)) + conv = self.PLANK.conv + + self.ROPHYT[loop] = self.PLANK.rophyt * conv + self.ROZOO[loop] = self.PLANK.rozoo * conv + self.ROORN[loop] = self.PLANK.roorn * conv + self.ROORP[loop] = self.PLANK.roorp * conv + self.ROORC[loop] = self.PLANK.roorc * conv + + self.ROTORN[loop] = self.PLANK.rotorn * conv + self.ROTORP[loop] = self.PLANK.rotorp * conv + self.ROTORC[loop] = self.PLANK.rotorc * conv + + # exit outflows: + if self.nexits > 1: + for i in range(self.nexits): + ts['PKCF2' + str(i + 1) + ' 1'][loop] = self.PLANK.ophyt[i] * conv + ts['PKCF2' + str(i + 1) + ' 2'][loop] = self.PLANK.ozoo[i] * conv + ts['PKCF2' + str(i + 1) + ' 3'][loop] = self.PLANK.oorn[i] * conv + ts['PKCF2' + str(i + 1) + ' 4'][loop] = self.PLANK.oorp[i] * conv + ts['PKCF2' + str(i + 1) + ' 5'][loop] = self.PLANK.oorc[i] * conv + + # PHCARB results: + if self.PHFG == 1: + pass + + return \ No newline at end of file diff --git a/HSP2/RQUTIL.py b/HSP2/RQUTIL.py index 8a6aece1..1c58fd8b 100644 --- a/HSP2/RQUTIL.py +++ b/HSP2/RQUTIL.py @@ -1,7 +1,6 @@ +from numba import njit - - -#@jit(nopython=True) +@njit(cache=True) def benth (dox, anaer, BRCON, scrfac, depcor, conc): ''' simulate benthal release of constituent''' # calculate benthal release of constituent; release is a step function of aerobic/anaerobic conditions, and stream velocity; @@ -12,7 +11,7 @@ def benth (dox, anaer, BRCON, scrfac, depcor, conc): return conc, releas -#@jit(nopython=True) +@njit(cache=True) def decbal(TAMFG, PO4FG, decnit, decpo4, tam, no3, po4): ''' perform materials balance for transformation from organic to inorganic material by decay in reach water''' if TAMFG: @@ -24,9 +23,11 @@ def decbal(TAMFG, PO4FG, decnit, decpo4, tam, no3, po4): return tam, no3, po4 -#@jit(nopython=True) -def sink (vol, avdepe, kset, conc, snkmat): +@njit(cache=True) +def sink (vol, avdepe, kset, conc): ''' calculate quantity of material settling out of the control volume; determine the change in concentration as a result of sinking''' + snkmat = 0.0 + if kset > 0.0 and avdepe > 0.17: # calculate concentration change due to outgoing material; snkout is expressed in mass/liter/ivl; kset is expressed as ft/ivl and avdepe as feet snkout = conc * (kset / avdepe) if kset < avdepe else conc # calculate portion of material which settles out of the control volume during time step; snkout is expressed as mass/liter.ivl; conc is the concentration of material in the control volume diff --git a/HSP2/configuration.py b/HSP2/configuration.py index 63c835fb..158039be 100644 --- a/HSP2/configuration.py +++ b/HSP2/configuration.py @@ -26,6 +26,14 @@ from HSP2.SEDTRN import sedtrn, expand_SEDTRN_masslinks from HSP2.CONS import cons, expand_CONS_masslinks from HSP2.GQUAL import gqual, expand_GQUAL_masslinks +from HSP2.RQUAL import rqual +from HSP2.RQUAL import expand_OXRX_masslinks +from HSP2.RQUAL import expand_NUTRX_masslinks +from HSP2.RQUAL import expand_PLANK_masslinks + +#from HSP2.GENER import gener +from HSP2.COPY import Copy +from HSP2.GENER import Gener def noop (store, siminfo, ui, ts): ERRMSGS = [] @@ -34,13 +42,15 @@ def noop (store, siminfo, ui, ts): # Note: This is the ONLY place in HSP2 that defines activity execution order activities = { + 'COPY' : Copy, + 'GENER' : Gener, 'PERLND': {'ATEMP':atemp, 'SNOW':snow, 'PWATER':pwater, 'SEDMNT':sedmnt, 'PSTEMP':pstemp, 'PWTGAS':pwtgas, 'PQUAL':pqual, 'MSTLAY':noop, 'PEST':noop, 'NITR':noop, 'PHOS':noop, 'TRACER':noop}, 'IMPLND': {'ATEMP':atemp, 'SNOW':snow, 'IWATER':iwater, 'SOLIDS':solids, 'IWTGAS':iwtgas, 'IQUAL':iqual}, 'RCHRES': {'HYDR':hydr, 'ADCALC':adcalc, 'CONS':cons, 'HTRCH':htrch, - 'SEDTRN':sedtrn, 'GQUAL':gqual, 'OXRX':noop, 'NUTRX':noop, 'PLANK':noop, + 'SEDTRN':sedtrn, 'GQUAL':gqual, 'OXRX':noop, 'NUTRX':noop, 'PLANK':noop, 'RQUAL':rqual, 'PHCARB':noop}} def expand_masslinks(flags, uci, dat, recs): @@ -49,6 +59,11 @@ def expand_masslinks(flags, uci, dat, recs): recs = expand_CONS_masslinks(flags, uci, dat, recs) recs = expand_SEDTRN_masslinks(flags, uci, dat, recs) recs = expand_GQUAL_masslinks(flags, uci, dat, recs) + recs = expand_OXRX_masslinks(flags, uci, dat, recs) + recs = expand_NUTRX_masslinks(flags, uci, dat, recs) + recs = expand_PLANK_masslinks(flags, uci, dat, recs) + recs = expand_PHCARB_masslinks(flags, uci, dat, recs) + return recs # NOTE: the flowtype (Python set) at the top of utilities.py may need to be diff --git a/HSP2/main.py b/HSP2/main.py index 9163cb75..537bc294 100644 --- a/HSP2/main.py +++ b/HSP2/main.py @@ -11,12 +11,13 @@ from collections import defaultdict from datetime import datetime as dt import os -from HSP2.utilities import transform, versions +from HSP2.utilities import transform, versions, get_timeseries from HSP2.configuration import activities, noop, expand_masslinks +from typing import List def main(hdfname, saveall=False, jupyterlab=True): - '''Runs main HSP2 program. + """Runs main HSP2 program. Parameters ---------- @@ -25,126 +26,167 @@ def main(hdfname, saveall=False, jupyterlab=True): saveall: Boolean [optional] Default is False. Saves all calculated data ignoring SAVE tables. - ''' + """ if not os.path.exists(hdfname): - print(f'{hdfname} HDF5 File Not Found, QUITTING') - return + raise FileNotFoundError(f'{hdfname} HDF5 File Not Found') with HDFStore(hdfname, 'a') as store: msg = messages() msg(1, f'Processing started for file {hdfname}; saveall={saveall}') # read user control, parameters, states, and flags from HDF5 file - opseq, ddlinks, ddmasslinks, ddext_sources, uci, siminfo = get_uci(store) + opseq, ddlinks, ddmasslinks, ddext_sources, ddgener, uci, siminfo = get_uci(store) start, stop = siminfo['start'], siminfo['stop'] + copy_instances = {} + gener_instances = {} + # main processing loop msg(1, f'Simulation Start: {start}, Stop: {stop}') for _, operation, segment, delt in opseq.itertuples(): - msg(2, f'{operation} {segment} DELT(minutes): {delt}') - siminfo['delt'] = delt - siminfo['tindex'] = date_range(start, stop, freq=Minute(delt))[1:] - siminfo['steps'] = len(siminfo['tindex']) - - # now conditionally execute all activity modules for the op, segment - ts = get_timeseries(store,ddext_sources[(operation,segment)],siminfo) - flags = uci[(operation, 'GENERAL', segment)]['ACTIVITY'] - if operation == 'RCHRES': - get_flows(store, ts, flags, uci, segment, ddlinks, ddmasslinks, siminfo['steps'], msg) - - for activity, function in activities[operation].items(): - if function == noop or not flags[activity]: - continue - - msg(3, f'{activity}') - - 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'] - if operation == 'PERLND' and activity == 'PSTEMP': - # special exception here to make AIRTFG available - ui['PARAMETERS']['AIRTFG'] = flags['ATEMP'] - if operation == 'PERLND' and activity == 'PWTGAS': - # special exception here to make CSNOFG available - ui['PARAMETERS']['CSNOFG'] = uci[(operation, 'PWATER', segment)]['PARAMETERS']['CSNOFG'] + if operation == 'COPY': + copy_instances[segment] = activities[operation](store, siminfo, ddext_sources[(operation,segment)]) + elif operation == 'GENER': + gener_instances[segment] = activities[operation](segment, copy_instances, gener_instances, ddlinks, ddgener) + else: + msg(2, f'{operation} {segment} DELT(minutes): {delt}') + siminfo['delt'] = delt + siminfo['tindex'] = date_range(start, stop, freq=Minute(delt))[1:] + siminfo['steps'] = len(siminfo['tindex']) + + # now conditionally execute all activity modules for the op, segment + ts = get_timeseries(store,ddext_sources[(operation,segment)],siminfo) + ts = get_gener_timeseries(ts, gener_instances, ddlinks[segment]) + flags = uci[(operation, 'GENERAL', segment)]['ACTIVITY'] if operation == 'RCHRES': - if not 'PARAMETERS' in ui: - ui['PARAMETERS'] = {} - ui['PARAMETERS']['NEXITS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['NEXITS'] - if activity == 'ADCALC': - ui['PARAMETERS']['ADFG'] = flags['ADCALC'] - ui['PARAMETERS']['KS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['KS'] - ui['PARAMETERS']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - if activity == 'HTRCH': - ui['PARAMETERS']['ADFG'] = flags['ADCALC'] - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] - # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - if activity == 'CONS': - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] - if activity == 'SEDTRN': - ui['PARAMETERS']['ADFG'] = flags['ADCALC'] - 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'] - ui['PARAMETERS']['SEDFG'] = flags['SEDTRN'] - # ui['PARAMETERS']['REAMFG'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['REAMFG'] - ui['PARAMETERS']['HYDRFG'] = flags['HYDR'] - if flags['HYDR']: - ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] - ui['PARAMETERS']['AUX1FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX1FG'] - ui['PARAMETERS']['AUX2FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX2FG'] - ui['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN'] - ui['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH'] - if flags['OXRX']: - ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] - ui['PARAMETERS']['CFOREA'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['CFOREA'] - if flags['SEDTRN']: - ui['PARAMETERS']['SSED1'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED1'] - ui['PARAMETERS']['SSED2'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED2'] - ui['PARAMETERS']['SSED3'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED3'] - if flags['HTRCH']: - ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'HTRCH', segment)]['PARAMETERS']['CFSAEX'] - elif flags['PLANK']: - if 'CFSAEX' in uci[(operation, 'PLANK', segment)]['PARAMETERS']: - ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'PLANK', segment)]['PARAMETERS']['CFSAEX'] - - if activity in ['OXRX','NUTRX','PLANK']: - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] - if flags['HYDR']: - ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] - ui['FLAGS']['BENRFG'] = uci[(operation, 'RQUAL', segment)]['FLAGS']['BENRFG'] - - if activity == 'PLANK': + get_flows(store, ts, flags, uci, segment, ddlinks, ddmasslinks, siminfo['steps'], msg) + + for activity, function in activities[operation].items(): + if function == noop: #or not flags[activity]: + continue + + if (activity in flags) and (not flags[activity]): + continue + + msg(3, f'{activity}') + + 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'] + if operation == 'PERLND' and activity == 'PSTEMP': + # special exception here to make AIRTFG available + ui['PARAMETERS']['AIRTFG'] = flags['ATEMP'] + if operation == 'PERLND' and activity == 'PWTGAS': + # special exception here to make CSNOFG available + ui['PARAMETERS']['CSNOFG'] = uci[(operation, 'PWATER', segment)]['PARAMETERS']['CSNOFG'] + if operation == 'RCHRES': + if not 'PARAMETERS' in ui: + ui['PARAMETERS'] = {} + ui['PARAMETERS']['NEXITS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['NEXITS'] + if activity == 'ADCALC': + ui['PARAMETERS']['ADFG'] = flags['ADCALC'] + ui['PARAMETERS']['KS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['KS'] + ui['PARAMETERS']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] + ui['PARAMETERS']['ROS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['ROS'] + if activity == 'HTRCH': + ui['PARAMETERS']['ADFG'] = flags['ADCALC'] + ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] + # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] + if activity == 'CONS': + ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] + if activity == 'SEDTRN': + ui['PARAMETERS']['ADFG'] = flags['ADCALC'] + ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] + # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] ui['PARAMETERS']['HTFG'] = flags['HTRCH'] - ui['FLAGS']['BPCNTC'] = uci[(operation, 'NUTRX', segment)]['PARAMETERS']['BPCNTC'] - ui['FLAGS']['CVBO'] = uci[(operation, 'NUTRX', segment)]['PARAMETERS']['CVBO'] - ui['FLAGS']['CVBPC'] = uci[(operation, 'NUTRX', segment)]['PARAMETERS']['CVBPC'] - ui['FLAGS']['NH3FG'] = uci[(operation, 'NUTRX', segment)]['FLAGS']['NH3FG'] - ui['FLAGS']['PO4FG'] = uci[(operation, 'NUTRX', segment)]['FLAGS']['PO4FG'] - - if activity == 'RQUAL': - pass - - ############ calls activity function like snow() ############## - errors, errmessages = function(store, siminfo, ui, ts) - ############################################################### - - for errorcnt, errormsg in zip(errors, errmessages): - if errorcnt > 0: - msg(4, f'Error count {errorcnt}: {errormsg}') - if 'SAVE' in ui: - save_timeseries(store,ts,ui['SAVE'],siminfo,saveall,operation,segment,activity,jupyterlab) + 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'] + ui['PARAMETERS']['SEDFG'] = flags['SEDTRN'] + # ui['PARAMETERS']['REAMFG'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['REAMFG'] + ui['PARAMETERS']['HYDRFG'] = flags['HYDR'] + if flags['HYDR']: + ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] + ui['PARAMETERS']['AUX1FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX1FG'] + ui['PARAMETERS']['AUX2FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX2FG'] + ui['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN'] + ui['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH'] + if flags['OXRX']: + ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] + ui['PARAMETERS']['CFOREA'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['CFOREA'] + if flags['SEDTRN']: + ui['PARAMETERS']['SSED1'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED1'] + ui['PARAMETERS']['SSED2'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED2'] + ui['PARAMETERS']['SSED3'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED3'] + if flags['HTRCH']: + ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'HTRCH', segment)]['PARAMETERS']['CFSAEX'] + elif flags['PLANK']: + if 'CFSAEX' in uci[(operation, 'PLANK', segment)]['PARAMETERS']: + ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'PLANK', segment)]['PARAMETERS']['CFSAEX'] + + if activity == 'RQUAL': + # RQUAL inputs: + ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] + if flags['HYDR']: + ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] + + ui['FLAGS']['HTFG'] = flags['HTRCH'] + ui['FLAGS']['SEDFG'] = flags['SEDTRN'] + ui['FLAGS']['GQFG'] = flags['GQUAL'] + ui['FLAGS']['GQALFG4'] = uci[(operation, 'GQUAL', segment)]['GQUAL1']['QALFG4'] + ui['FLAGS']['OXFG'] = flags['OXFG'] + ui['FLAGS']['NUTFG'] = flags['NUTRX'] + ui['FLAGS']['PLKFG'] = flags['PLANK'] + ui['FLAGS']['PHFG'] = flags['PHCARB'] + + # OXRX module inputs: + ui_oxrx = uci[(operation, 'OXRX', segment)] + + if flags['HYDR']: + ui_oxrx['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN'] + ui_oxrx['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH'] + + if flags['HTRCH']: + ui_oxrx['PARAMETERS']['ELEV'] = uci[(operation, 'HTRCH', segment)]['PARAMETERS']['ELEV'] + + if flags['SEDTRN']: + ui['PARAMETERS']['SSED1'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED1'] + ui['PARAMETERS']['SSED2'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED2'] + ui['PARAMETERS']['SSED3'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED3'] + + # NUTRX, PLANK, PHCARB module inputs: + ui_nutrx = uci[(operation, 'NUTRX', segment)] + ui_plank = uci[(operation, 'PLANK', segment)] + ui_phcarb = uci[(operation, 'PHCARB', segment)] + + ############ calls activity function like snow() ############## + if operation not in ['COPY','GENER']: + if (activity != 'RQUAL'): + errors, errmessages = function(store, siminfo, ui, ts) + else: + errors, errmessages = function(store, siminfo, ui, ui_oxrx, ui_nutrx, ui_plank, ui_phcarb, ts) + ############################################################### + + for errorcnt, errormsg in zip(errors, errmessages): + if errorcnt > 0: + msg(4, f'Error count {errorcnt}: {errormsg}') + if 'SAVE' in ui: + save_timeseries(store,ts,ui['SAVE'],siminfo,saveall,operation,segment,activity,jupyterlab) + + if (activity == 'RQUAL'): + if 'SAVE' in ui_oxrx: save_timeseries(store,ts,ui_oxrx['SAVE'],siminfo,saveall,operation,segment,'OXRX',jupyterlab) + if 'SAVE' in ui_nutrx: save_timeseries(store,ts,ui_nutrx['SAVE'],siminfo,saveall,operation,segment,'NUTRX',jupyterlab) + if 'SAVE' in ui_plank: save_timeseries(store,ts,ui_plank['SAVE'],siminfo,saveall,operation,segment,'PLANK',jupyterlab) + #if 'SAVE' in ui_phcarb: save_timeseries(store,ts,ui_phcarb['SAVE'],siminfo,saveall,operation,segment,'PHCARB',jupyterlab) + msglist = msg(1, 'Done', final=True) df = DataFrame(msglist, columns=['logfile']) @@ -156,7 +198,6 @@ def main(hdfname, saveall=False, jupyterlab=True): print('\n\n', df) return - def messages(): '''Closure routine; msg() prints messages to screen and run log''' start = dt.now() @@ -180,6 +221,7 @@ def get_uci(store): ddlinks = defaultdict(list) ddmasslinks = defaultdict(list) ddext_sources = defaultdict(list) + ddgener =defaultdict(dict) siminfo = {} opseq = 0 @@ -196,8 +238,8 @@ def get_uci(store): if int(temp['Units']): siminfo['units'] = int(temp['Units']) elif module == 'LINKS': - for row in store[path].replace('na','').itertuples(): - ddlinks[row.TVOLNO].append(row) + for row in store[path].fillna('').itertuples(): + ddlinks[f'{row.TVOLNO}{row.TOPFST}'].append(row) elif module == 'MASS_LINKS': for row in store[path].replace('na','').itertuples(): ddmasslinks[row.MLNO].append(row) @@ -209,35 +251,11 @@ def get_uci(store): elif op in {'PERLND', 'IMPLND', 'RCHRES'}: for id, vdict in store[path].to_dict('index').items(): uci[(op, module, id)][s] = vdict - return opseq, ddlinks, ddmasslinks, ddext_sources, uci, siminfo - - -def get_timeseries(store, ext_sourcesdd, siminfo): - ''' makes timeseries for the current timestep and trucated to the sim interval''' - # explicit creation of Numba dictionary with signatures - ts = Dict.empty(key_type=types.unicode_type, value_type=types.float64[:]) - for row in ext_sourcesdd: - if row.SVOL == '*': - path = f'TIMESERIES/{row.SVOLNO}' - if path in store: - temp1 = store[path] - else: - print('Get Timeseries ERROR for', path) - continue - else: - temp1 = read_hdf(row.SVOL, path) - - if row.MFACTOR != 1.0: - temp1 *= row.MFACTOR - t = transform(temp1, row.TMEMN, row.TRAN, siminfo) - - tname = f'{row.TMEMN}{row.TMEMSB}' - if tname in ts: - ts[tname] += t - else: - ts[tname] = t - return ts - + elif op == 'GENER': + for row in store[path].itertuples(): + start, stop = row.OPNID.split() + for i in range(int(start), int(stop)): ddgener[module][f'G{i:03d}'] = row[2] + return opseq, ddlinks, ddmasslinks, ddext_sources, ddgener, uci, siminfo def save_timeseries(store, ts, savedict, siminfo, saveall, operation, segment, activity, jupyterlab=True): # save computed timeseries (at computation DELT) @@ -326,7 +344,9 @@ def get_flows(store, ts, flags, uci, segment, ddlinks, ddmasslinks, steps, msg): factor = afactr * mfactor # KLUDGE until remaining HSP2 modules are available. - if tmemn not in {'IVOL', 'ICON', 'IHEAT', 'ISED', 'ISED1', 'ISED2', 'ISED3', 'IDQAL', 'ISQAL1', 'ISQAL2', 'ISQAL3'}: + if tmemn not in {'IVOL', 'ICON', 'IHEAT', 'ISED', 'ISED1', 'ISED2', 'ISED3', + 'IDQAL', 'ISQAL1', 'ISQAL2', 'ISQAL3', + 'OXIF', 'NUIF1', 'NUIF2', 'PKIF'}: continue if (sgrpn == 'OFLOW' and smemn == 'OVOL') or (sgrpn == 'ROFLOW' and smemn == 'ROVOL'): sgrpn = 'HYDR' @@ -338,6 +358,13 @@ def get_flows(store, ts, flags, uci, segment, ddlinks, ddmasslinks, steps, msg): sgrpn = 'GQUAL' if (sgrpn == 'OFLOW' and smemn == 'OSQAL') or (sgrpn == 'ROFLOW' and smemn == 'ROSQAL'): sgrpn = 'GQUAL' + if (sgrpn == 'OFLOW' and smemn == 'OXCF2') or (sgrpn == 'ROFLOW' and smemn == 'OXCF1'): + sgrpn = 'OXRX' + if (sgrpn == 'OFLOW' and (smemn == 'NUCF9' or smemn == 'OSNH4' or smemn == 'OSPO4')) or (sgrpn == 'ROFLOW' and (smemn == 'NUCF1' or smemn == 'NUFCF2')): + sgrpn = 'NUTRX' + if (sgrpn == 'OFLOW' and smemn == 'PKCF2') or (sgrpn == 'ROFLOW' and smemn == 'PKCF1'): + sgrpn = 'PLANK' + if tmemn == 'ISED' or tmemn == 'ISQAL': tmemn = tmemn + tmemsb1 # need to add sand, silt, clay subscript @@ -401,6 +428,7 @@ def expand_timeseries_names(smemn, smemsb1, smemsb2, tmemn, tmemsb1, tmemsb2): else: smemn = 'CONS' + smemsb1 + '_ROCON' + # GQUAL: if tmemn == 'IDQAL': if tmemsb1 == '': tmemn = 'GQUAL1_IDQAL' @@ -420,8 +448,60 @@ def expand_timeseries_names(smemn, smemsb1, smemsb2, tmemn, tmemsb1, tmemsb2): if smemn == 'ROSQAL': smemn = 'GQUAL' + smemsb2 + '_ROSQAL' + smemsb1 # smemsb1 is ssc + # OXRX: + if smemn == 'OXCF1': + smemn = 'OXCF1' + smemsb1 + + if smemn == 'OXCF2': + smemn = 'OXCF2' + smemsb1 + ' ' + smemsb2 # smemsb1 is exit # + + if tmemn == 'OXIF': + tmemn = 'OXIF' + tmemsb1 + + # NUTRX - dissolved species: + if smemn == 'NUCF1': # total outflow + smemn = 'NUCF1' + smemsb1 + + if smemn == 'NUCF9': # exit-specific outflow + smemn = 'NUCF9' + smemsb1 + ' ' + smemsb2 # smemsb1 is exit # + + if tmemn == 'NUIF1': + tmemn = 'NUIF1' + tmemsb1 + + # NUTRX - particulate species: + if smemn == 'NUCF2': # total outflow + smemn = 'NUCF2' + smemsb1 + ' ' + smemsb2 # smemsb1 is sediment class + + if smemn == 'OSNH4' or smemn == 'OSPO4': # exit-specific outflow + smemn = smemn + smemsb1 + ' ' + smemsb2 # smemsb1 is exit #, smemsb2 is sed class + + if tmemn == 'NUIF2': + tmemn = 'NUIF2' + tmemsb1 + ' ' + tmemsb2 + + # PLANK: + if smemn == 'PKCF1': # total outflow + smemn = 'PKCF1' + smemsb1 # smemsb1 is species index + + if smemn == 'PKCF2': # exit-specific outflow + smemn = 'PKCF2' + smemsb1 + ' ' + smemsb2 # smemsb1 is exit #, smemsb2 is species index + + if tmemn == 'PKIF': + tmemn = 'PKIF' + tmemsb1 # smemsb1 is species index + return smemn, tmemn +def get_gener_timeseries(ts: Dict, gener_instances: Dict, ddlinks: List) -> Dict: + """ + Uses links tables to load necessary TimeSeries from Gener class instances to TS dictionary + """ + for link in ddlinks: + if link.SVOL == 'GENER': + gener = gener_instances[link.SVOLNO] + series = gener.get_ts() + ts[f'{link.TMEMN}{link.TMEMSB1}{link.TMEMSB2}'] = series + return ts + + ''' # This table defines the expansion to INFLOW, ROFLOW, OFLOW for RCHRES networks diff --git a/HSP2/utilities.py b/HSP2/utilities.py index 0c2b09d7..1790eb07 100644 --- a/HSP2/utilities.py +++ b/HSP2/utilities.py @@ -220,3 +220,29 @@ def versions(import_list=[]): data.extend([platform.platform(), platform.processor(), str(datetime.datetime.now())[0:19]]) return pandas.DataFrame(data, index=names, columns=['version']) + +def get_timeseries(store, ext_sourcesdd, siminfo): + """ makes timeseries for the current timestep and trucated to the sim interval""" + # explicit creation of Numba dictionary with signatures + ts = Dict.empty(key_type=types.unicode_type, value_type=types.float64[:]) + for row in ext_sourcesdd: + if row.SVOL == '*': + path = f'TIMESERIES/{row.SVOLNO}' + if path in store: + temp1 = store[path] + else: + print('Get Timeseries ERROR for', path) + continue + else: + temp1 = read_hdf(row.SVOL, path) + + if row.MFACTOR != 1.0: + temp1 *= row.MFACTOR + t = transform(temp1, row.TMEMN, row.TRAN, siminfo) + + tname = f'{row.TMEMN}{row.TMEMSB}' + if tname in ts: + ts[tname] += t + else: + ts[tname] = t + return ts diff --git a/HSP2tools/data/ParseTable.csv b/HSP2tools/data/ParseTable.csv index b48a170f..7835ac0d 100644 --- a/HSP2tools/data/ParseTable.csv +++ b/HSP2tools/data/ParseTable.csv @@ -2361,7 +2361,7 @@ 2302,NETWORK,na,GLOBAL,NETWORK,SVOLNO,C,6,11,na 2303,NETWORK,na,GLOBAL,NETWORK,SGRPN,C,11,18,na 2304,NETWORK,na,GLOBAL,NETWORK,SMEMN,C,18,24,na -2305,NETWORK,na,GLOBAL,NETWORK,SMEMSB1,C,25,26,na +2305,NETWORK,na,GLOBAL,NETWORK,SMEMSB1,C,24,26,na 2305,NETWORK,na,GLOBAL,NETWORK,SMEMSB2,C,27,28,na 2306,NETWORK,na,GLOBAL,NETWORK,MFACTOR,R,28,38,1 2307,NETWORK,na,GLOBAL,NETWORK,TRAN,C,38,43,na @@ -2391,3 +2391,7 @@ 2328,MASS-LINK,na,GLOBAL,MASS-LINK,TMEMN,C,65,71,na 2329,MASS-LINK,na,GLOBAL,MASS-LINK,TMEMSB1,C,72,73,na 2329,MASS-LINK,na,GLOBAL,MASS-LINK,TMEMSB2,C,74,75,na +2330,GENER,OPCODE,GENER,OPCODES,OPNID,C,0,10,0 +2331,GENER,OPCODE,GENER,OPCODES,OPCODE,I,11,16,0 +2332,GENER,PARM,GENER,INFO,OPNID,C,0,10,0 +2332,GENER,PARM,GENER,INFO,K,R,11,21,0.0E0 \ No newline at end of file diff --git a/HSP2tools/data/SaveTable.csv b/HSP2tools/data/SaveTable.csv index c4307e20..691fdf4a 100644 --- a/HSP2tools/data/SaveTable.csv +++ b/HSP2tools/data/SaveTable.csv @@ -369,64 +369,92 @@ RCHRES,GQUAL,SQDEC,0 RCHRES,GQUAL,TIQAL,0 RCHRES,GQUAL,TOSQAL,0 RCHRES,GQUAL,TROQAL,0 -RCHRES,OXRX,BOD,0 -RCHRES,OXRX,DOX,0 -RCHRES,OXRX,OXCF1,0 -RCHRES,OXRX,OXCF2,0 -RCHRES,OXRX,OXCF3,0 -RCHRES,OXRX,OXCF4,0 -RCHRES,OXRX,OXIF,0 +RCHRES,OXRX,BOD,1 +RCHRES,OXRX,DOX,1 RCHRES,OXRX,SATDO,0 -RCHRES,NUTRX,DNUST,0 -RCHRES,NUTRX,DNUST2,0 -RCHRES,NUTRX,NUADDR,0 -RCHRES,NUTRX,NUADEP,0 -RCHRES,NUTRX,NUADWT,0 -RCHRES,NUTRX,NUCF1,0 -RCHRES,NUTRX,NUCF2,0 -RCHRES,NUTRX,NUCF3,0 -RCHRES,NUTRX,NUCF4,0 -RCHRES,NUTRX,NUCF5,0 -RCHRES,NUTRX,NUCF6,0 -RCHRES,NUTRX,NUCF7,0 -RCHRES,NUTRX,NUCF8,0 -RCHRES,NUTRX,NUCF9,0 -RCHRES,NUTRX,NUIF1,0 -RCHRES,NUTRX,NUIF2,0 -RCHRES,NUTRX,NUST,0 -RCHRES,NUTRX,OSNH4,0 -RCHRES,NUTRX,OSPO4,0 -RCHRES,NUTRX,RSNH4,0 -RCHRES,NUTRX,RSPO4,0 -RCHRES,NUTRX,SNH4,0 -RCHRES,NUTRX,SPO4,0 -RCHRES,NUTRX,TNUCF1,0 -RCHRES,NUTRX,TNUCF2,0 -RCHRES,NUTRX,TNUIF,0 +RCHRES,OXRX,IDOX,0 +RCHRES,OXRX,IBOD,0 +RCHRES,OXRX,OXCF11,1 +RCHRES,OXRX,OXCF12,1 +RCHRES,OXRX,OXCF21 1,1 +RCHRES,OXRX,OXCF21 2,1 +RCHRES,OXRX,OXCF22 1,1 +RCHRES,OXRX,OXCF22 2,1 +RCHRES,OXRX,OXCF23 1,1 +RCHRES,OXRX,OXCF23 2,1 +RCHRES,OXRX,OXCF24 1,1 +RCHRES,OXRX,OXCF24 2,1 +RCHRES,OXRX,OXCF25 1,1 +RCHRES,OXRX,OXCF25 2,1 +RCHRES,OXRX,OXCF3_REAR,0 +RCHRES,OXRX,OXCF3_DEC,0 +RCHRES,OXRX,OXCF3_BENDO,0 +RCHRES,OXRX,OXCF3_NITR,0 +RCHRES,OXRX,OXCF3_PHYT,0 +RCHRES,OXRX,OXCF3_ZOO,0 +RCHRES,OXRX,OXCF3_BALG,0 +RCHRES,OXRX,OXCF3_TOTAL,0 +RCHRES,OXRX,OXCF4_DEC,0 +RCHRES,OXRX,OXCF4_BENR,0 +RCHRES,OXRX,OXCF4_SNK,0 +RCHRES,OXRX,OXCF4_PHYT,0 +RCHRES,OXRX,OXCF4_ZOO,0 +RCHRES,OXRX,OXCF4_BALG,0 +RCHRES,OXRX,OXCF4_TOTAL,0 +RCHRES,NUTRX,NO3,0 +RCHRES,NUTRX,TAM,0 +RCHRES,NUTRX,NO2,0 +RCHRES,NUTRX,PO4,0 +RCHRES,NUTRX,NH4,0 +RCHRES,NUTRX,NH3,0 +RCHRES,NUTRX,NUCF11,0 +RCHRES,NUTRX,NUCF12,0 +RCHRES,NUTRX,NUCF13,0 +RCHRES,NUTRX,NUCF14,0 +RCHRES,NUTRX,NUCF21 1,0 +RCHRES,NUTRX,NUCF22 1,0 +RCHRES,NUTRX,NUCF23 1,0 +RCHRES,NUTRX,NUCF21 2,0 +RCHRES,NUTRX,NUCF22 2,0 +RCHRES,NUTRX,NUCF23 2,0 +RCHRES,NUTRX,RNO3,0 +RCHRES,NUTRX,RTAM,0 +RCHRES,NUTRX,RNO2,0 +RCHRES,NUTRX,RPO4,0 +RCHRES,NUTRX,RNH4,0 +RCHRES,NUTRX,RNH3,0 RCHRES,PLANK,PHYTO,0 RCHRES,PLANK,ZOO,0 -RCHRES,PLANK,BENAL,0 -RCHRES,PLANK,TBENAL,0 +RCHRES,PLANK,BENAL1,0 +RCHRES,PLANK,BENAL2,0 +RCHRES,PLANK,BENAL3,0 +RCHRES,PLANK,BENAL4,0 +RCHRES,PLANK,TBENAL1,0 +RCHRES,PLANK,TBENAL2,0 RCHRES,PLANK,PHYCLA,0 RCHRES,PLANK,BALCLA,0 -RCHRES,PLANK,PKST3,0 -RCHRES,PLANK,PKST4,0 -RCHRES,PLANK,PKIF,0 -RCHRES,PLANK,TPKIF,0 -RCHRES,PLANK,PKCF1,0 -RCHRES,PLANK,TPKCF1,0 -RCHRES,PLANK,PKCF2,0 -RCHRES,PLANK,TPKCF2,0 -RCHRES,PLANK,PLADDR,0 -RCHRES,PLANK,PLADWT,0 -RCHRES,PLANK,PLADEP,0 -RCHRES,PLANK,PKCF5,0 -RCHRES,PLANK,PKCF6,0 -RCHRES,PLANK,PKCF7,0 -RCHRES,PLANK,TPKCF7,0 -RCHRES,PLANK,PKCF8,0 -RCHRES,PLANK,PKCF9,0 -RCHRES,PLANK,PKCF10,0 +RCHRES,PLANK,ORN,0 +RCHRES,PLANK,ORP,0 +RCHRES,PLANK,ORC,0 +RCHRES,PLANK,TORN,0 +RCHRES,PLANK,TORP,0 +RCHRES,PLANK,TORC,0 +RCHRES,PLANK,POTBOD,0 +RCHRES,PLANK,TN,0 +RCHRES,PLANK,TP,0 +RCHRES,PLANK,PKIF_PHYT,0 +RCHRES,PLANK,PKIF_ZOO,0 +RCHRES,PLANK,PKIF_ORN,0 +RCHRES,PLANK,PKIF_ORP,0 +RCHRES,PLANK,PKIF_ORC,0 +RCHRES,PLANK,PKCF11,0 +RCHRES,PLANK,PKCF12,0 +RCHRES,PLANK,PKCF13,0 +RCHRES,PLANK,PKCF14,0 +RCHRES,PLANK,PKCF15,0 +RCHRES,PLANK,ROTORN,0 +RCHRES,PLANK,ROTORP,0 +RCHRES,PLANK,ROTORC,0 RCHRES,PHCARB,PHCF1,0 RCHRES,PHCARB,PHCF2,0 RCHRES,PHCARB,PHCF3,0 diff --git a/HSP2tools/readUCI.py b/HSP2tools/readUCI.py index 32e51810..08538b1b 100644 --- a/HSP2tools/readUCI.py +++ b/HSP2tools/readUCI.py @@ -103,7 +103,7 @@ def fix_df(df, op, save, ddfaults, valid): ('RCHRES', 'BINARY-INFO')} -ops = {'PERLND','IMPLND','RCHRES'} +ops = {'PERLND', 'IMPLND', 'RCHRES', 'COPY', 'GENER'} conlike = {'CONS':'NCONS', 'PQUAL':'NQUAL', 'IQUAL':'NQUAL', 'GQUAL':'NQUAL'} def readUCI(uciname, hdfname): # create lookup dictionaries from 'ParseTable.csv' and 'rename.csv' @@ -149,6 +149,7 @@ def readUCI(uciname, hdfname): if line[0:9] == 'MASS-LINK': masslink(info, getlines(f)) if line[0:7] == 'FTABLES': ftables(info, getlines(f)) if line[0:3] == 'EXT': ext(info, getlines(f)) + if line[0:5] == 'GENER': gener(info, getlines(f)) if line[0:6] == 'PERLND': operation(info, getlines(f),'PERLND') if line[0:6] == 'IMPLND': operation(info, getlines(f),'IMPLND') if line[0:6] == 'RCHRES': operation(info, getlines(f),'RCHRES') @@ -612,3 +613,25 @@ def operation(info, llines, op): for name,value in savetable[op,activity].items(): df[name] = int(value) df.to_hdf(store, f'{op}/{activity}/SAVE', data_columns=True) + + +def copy(info, lines): + #placeholder - PRT - no sure I actually need to implement this method + pass + +def gener(info, lines): + store, parse, path, *_ = info + lst = [] + sub_blocks = ['OPCODE','PARM'] + current_block = '' + d = {} + for line in lines: + if line [2:5] == 'END': + df = DataFrame(lst, columns=d) + df.to_hdf(store, key=f'GENER/{current_block}', data_columns=True) + lst.clear() + elif any(s in line for s in sub_blocks): + current_block = [s for s in sub_blocks if s in line][0] + else: + d = parseD(line, parse['GENER',current_block]) + lst.append(d) diff --git a/tests/test10/HSP2results/test10.hdf b/tests/test10/HSP2results/test10.hdf new file mode 100644 index 00000000..c15440f6 Binary files /dev/null and b/tests/test10/HSP2results/test10.hdf differ