Index: wflow-py/wflow/reservoir_Sw.py =================================================================== diff -u -rad5e049f1419c25bc30c20d2d0e1bcc56dcb4381 -r12ea40dc08628f654753679e0972e87b7bb12f7a --- wflow-py/wflow/reservoir_Sw.py (.../reservoir_Sw.py) (revision ad5e049f1419c25bc30c20d2d0e1bcc56dcb4381) +++ wflow-py/wflow/reservoir_Sw.py (.../reservoir_Sw.py) (revision 12ea40dc08628f654753679e0972e87b7bb12f7a) @@ -19,21 +19,23 @@ from copy import copy as copylist try: - from wflow.wf_DynamicFramework import * + from wflow.wf_DynamicFramework import * except ImportError: - from wf_DynamicFramework import * + from wf_DynamicFramework import * import scipy import JarvisCoefficients + def selectSwR(i): if i == 1: - name = 'snow' + name = "snow" if i == 2: - name = 'snowHour' + name = "snowHour" return name + def snow_no_reservoir(self, k): """ This function is used when no snow zone reservoir is used and only @@ -49,234 +51,315 @@ """ try: - JarvisCoefficients.calcEpSnow(self,k) + JarvisCoefficients.calcEpSnow(self, k) except: - JarvisCoefficients.calcEpSnowHour(self,k) + JarvisCoefficients.calcEpSnowHour(self, k) self.PotEvaporation = self.EpHour - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Qw_[k] = max(self.PrecipitationSnow, 0) self.Ew_[k] = 0. self.Sw[k] = 0. - self.wbSw_[k] = self.Precipitation - self.Ew_[k] - self.Qw_[k] - self.Sw[k] + self.Sw_t[k] - -def snow(self,k): + self.wbSw_[k] = ( + self.Precipitation - self.Ew_[k] - self.Qw_[k] - self.Sw[k] + self.Sw_t[k] + ) + + +def snow(self, k): """ - snow melt based on degree day factor and - - Code for ini-file: 1 """ - JarvisCoefficients.calcEpSnow(self,k) + JarvisCoefficients.calcEpSnow(self, k) self.PotEvaporation = self.EpHour - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - - self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0) + + self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]), 0) self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw - -def snow_rain(self,k): + + +def snow_rain(self, k): """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature - - Code for ini-file: 6 """ - - JarvisCoefficients.calcEpSnow(self,k) - #self.PotEvaporation = self.EpHour - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + + JarvisCoefficients.calcEpSnow(self, k) + # self.PotEvaporation = self.EpHour + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - + self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k]) - self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0) + self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]), 0) self.Qw1 = max(min(self.Fm2 * (self.Temperature - self.Tm[k]), self.Sw[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw -def snow_rain_hourlyEp(self,k): + +def snow_rain_hourlyEp(self, k): """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature - - Code for ini-file: 6 """ - - JarvisCoefficients.calcEpSnowHour(self,k) - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + + JarvisCoefficients.calcEpSnowHour(self, k) + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - + self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k]) - self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0) + self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]), 0) self.Qw1 = max(min(self.Fm2 * (self.Temperature - self.Tm[k]), self.Sw[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw -def snow_rain_Tsurf(self,k): + +def snow_rain_Tsurf(self, k): """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature - - Code for ini-file: 3 """ - - JarvisCoefficients.calcEpSnow(self,k) - #self.PotEvaporation = self.EpHour - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + + JarvisCoefficients.calcEpSnow(self, k) + # self.PotEvaporation = self.EpHour + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - + self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k]) - self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0) + self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]), 0) self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw -def snow_rain_TsurfAir(self,k): + +def snow_rain_TsurfAir(self, k): """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature - - Code for ini-file: 4 """ - JarvisCoefficients.calcEpSnow(self,k) - #self.PotEvaporation = self.EpHour - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + JarvisCoefficients.calcEpSnow(self, k) + # self.PotEvaporation = self.EpHour + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow self.Temp = (self.TempSurf + self.Temperature) / 2 - + self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k]) - self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]),0) + self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]), 0) self.Qw1 = max(min(self.Fm2 * (self.Temp - self.Tm[k]), self.Sw[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw -def snow_rain_Tsurf_noEw(self,k): + +def snow_rain_Tsurf_noEw(self, k): """ - snow melt based on degree day factor and minimum surface temperature - meltfactor increases with temperature - - Code for ini-file: 5 """ - JarvisCoefficients.calcEpSnow(self,k) - #self.PotEvaporation = self.EpHour - self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0),0) - + JarvisCoefficients.calcEpSnow(self, k) + # self.PotEvaporation = self.EpHour + self.PotEvaporation = cover(ifthenelse(self.EpHour > 0, self.EpHour, 0), 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - + self.Fm2 = max(self.Fm[k] * self.Precipitation, self.Fm[k]) self.Ew1 = 0 self.Qw1 = max(min(self.Fm2 * (self.TempSurf - self.Tm[k]), self.Sw[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw -def snowHour(self,k): + +def snowHour(self, k): """ - snow melt based on degree day factor and - for hourly input data - Code for ini-file: 2 """ -# JarvisCoefficients.calcEpSnowHour(self,k) -# self.PotEvaporation = self.EpHour -# self.PotEvaporation = ifthenelse(self.EpHour > 0, self.EpHour, 0) - + # JarvisCoefficients.calcEpSnowHour(self,k) + # self.PotEvaporation = self.EpHour + # self.PotEvaporation = ifthenelse(self.EpHour > 0, self.EpHour, 0) + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - - self.Ew1 = max(min(self.PotEvaporation,self.Sw[k]),0) -# self.Ew1 = 0 + + self.Ew1 = max(min(self.PotEvaporation, self.Sw[k]), 0) + # self.Ew1 = 0 self.Qw1 = max(self.Fm[k] * (self.Temperature - self.Tm[k]), 0) - + self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew1 - self.Qw1 self.Sw_diff = ifthenelse(self.Sw[k] < 0, self.Sw[k], 0) - self.Ew = self.Ew1 + (self.Ew1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff - self.Qw = self.Qw1 + (self.Qw1/ifthenelse(self.Ew1 + self.Qw1 > 0 , self.Ew1 + self.Qw1 , 1)) * self.Sw_diff + self.Ew = ( + self.Ew1 + + (self.Ew1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) + self.Qw = ( + self.Qw1 + + (self.Qw1 / ifthenelse(self.Ew1 + self.Qw1 > 0, self.Ew1 + self.Qw1, 1)) + * self.Sw_diff + ) self.Sw[k] = self.Sw_t[k] + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0 , self.Sw[k]) - self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) + self.Sw[k] = ifthenelse(self.Sw[k] < 0, 0, self.Sw[k]) + self.Sw_diff2 = ifthen(self.Sw[k] < 0, self.Sw[k]) -# if any(pcr2numpy(self.Sw[k],nan) > 0): -# pdb.set_trace() - self.wbSw_[k] = self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] - + # if any(pcr2numpy(self.Sw[k],nan) > 0): + # pdb.set_trace() + self.wbSw_[k] = ( + self.PrecipitationSnow - self.Ew - self.Qw - self.Sw[k] + self.Sw_t[k] + ) + self.Ew_[k] = self.Ew self.Qw_[k] = self.Qw -# if any(pcr2numpy(self.Qw,nan) > 0): + + +# if any(pcr2numpy(self.Qw,nan) > 0): # pdb.set_trace()