I have the following script:
import numpy as np
import matplotlib.pyplot as pl
time=np.linspace(0.,1825,1826)
#temp with peak at midpoint of year
T=23*np.sin((2*np.pi*time)/365)+1
#drainage basin area
A=8.094e10
#surface area of pond in mm
SA= 2.023e9
#hydraulic conductivity in mm/day/mm
K=2e-3
#GW level mm
Hgw=-2000
def precipPond(dayt):
'''
Input day with time notation for 365 days...
eg time[0] is dayt 1, time[1] is dayt 2, etc.
'''
#P varies from 0.2 mm/day to 1.8 mm/day with a peak at 120
P=np.sin(((2*np.pi*dayt)/365)+120)+0.8
#15 percent of all P flows into pond (mm)
infilP=(P*A*0.15)
return infilP
def evapoTrans(dayT):
'''
Input day with time notation for 365 days...
eg T[0] is temp at dayT 1, T[1] is temp at dayT 2, etc.
'''
#saturated vapor pressure for ET equation
et=0.611*(np.exp((17.3*dayT)/(dayT+237.3)))
#evapotranspiration
ET2=29.8*12*(et)/(dayT+273.2)
#return factor*ET2
return ET2
Hp=np.zeros(1826.)
Hp2=[]
precip=[]
evapoT=[]
def waterBudget(HpI,delt):
'''
Put in initial water level and change in time.
'''
#Save the initial water level in the first position
Hp[0]=HpI
steps=(len(time)-1)
factor=1
for i in range(steps):
#Get the mid point value, which will be used to estimate derivative
HpMid=Hp[i]+(delt/2)*(K*SA*((Hgw-Hp[i])/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*factor*SA))
#Find the derivative for dh/dt
Hp[i+1]=Hp[i]+delt*(K*SA*((Hgw-HpMid)/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*factor*SA))
if any(HpMid < 0):
factor=0
else:
factor=1
return Hp
ET=evapoTrans(T)
infilP=precipPond(time)
Hp=waterBudget(0,1)
I want the if statement in the loop (within the function waterBudget) to make ET zero when HpMid is below zero. I run this script and it makes factor 1 the whole time even when I've changed the script to make HpMid go below zero (**I didn't include that in here).