program lake_chain
	implicit none

	real ovol,odepth,oarea,oerate,o_spill,oevap
	real cvol,cdepth,carea,cerate,c_spill,cevap
	real svol,sdepth,sarea,serate,s_spill,sevap
	real pvol,pdepth,parea,perate,p_spill,pevap
	real mvol,mdepth,marea,merate,mevap
	real dt,runoff, time
	integer t
	dimension ovol(40000),cvol(40000),svol(40000),pvol(40000),mvol(40000)

!	****************initial conditions***********************

	ovol(1)=0.	!owens lake volume (m^3)
	odepth=0.	!owens lake depth (m)
	oarea=0.	!owens lake surface area (m^2)
	oerate=1.34 !owens evaporation rate (m/yr)
	o_spill=0.	!overflow to china lake (m^3)

	cvol(1)=0.		!china lake volume (m^3)
	cdepth=0.	!china lake depth (m)
	carea=0.	!china lake surface area (m^2)
	cerate=1.41	!china evaporation rate (m/yr)
	c_spill=0.	!overflow to searles lake (m^3)

	svol(1)=0.	!searles lake depth (m)
	sarea=0.	!searles lake area (m^2)
	serate=1.65 !searles evaporation rate (m/yr)
	s_spill=0.	!overflow to panamint (m^3)

	pvol(1)=0.		!panamint lake volume (m^3)
	pdepth=0.	!panamint lake depth (m)
	parea=0.	!panamint lake area (m^2)
	perate=1.65 !panamint evaporation rate (m/yr)
	p_spill=0.	!overflow to lake manly (m^3)

	mvol(1)=0.		!lake manly volume (m^3)
	mdepth=0.	!lake manly depth (m)
	marea=0.	!lake manly area (m^2)
	merate=1.97 !manly evaporation rate (m/yr)

	runoff=1.*3.98e8	!owens river flow = 3.98e8 m^3/yr

	dt=0.5	! model timestep = 0.5 years
	time=0.	!starting time

!	*****************end initial conditions*********************

	do t=2,600
	time=time+dt

!	**************************start owens***************************

	oevap=oarea*oerate						!calculate the evaporative outflow

	ovol(t)=ovol(t-1)+((runoff-oevap)*dt)	!update the volume

	if(ovol(t).gt.25.0e9)then				!determine whether overflow	is necessary		
	o_spill=ovol(t)-25.0e9					
	ovol(t)=25.0e9							!if the lake is overflowing, set the
	endif									!volume to be the maximum volume
	
	if(ovol(t).le.0)then					!make sure the area and depth are >=0.
	oarea=0.
	odepth=0.
	else
	oarea=75991.*(ovol(t)**0.37636)
	odepth=3.5837+(4.009e-9*ovol(t))-(9.5021e-20*(ovol(t)**2))+(1.2319e-30*(ovol(t)**3))
	endif

!	**************************start china***************************

	cevap=carea*cerate

	cvol(t)=cvol(t-1)+o_spill-(cevap*dt)	!note that o_spill is a volume.  Therefore, it
											!is not necessary to multiply it by dt.
	if(cvol(t).gt.1.2e9)then				
	c_spill=cvol(t)-1.2e9
	cvol(t)=1.2e9
	endif

	if(cvol(t).le.0)then
	carea=0.
	cdepth=0.
	else
	carea=1.238e8+(0.047279*cvol(t))-(1.4788e-12*(cvol(t)**2))
	cdepth=0.98637+(5.9886e-9*cvol(t))-(3.7316e-19*(cvol(t)**2))+(1.2192e-29*(cvol(t)**3))
	endif

!	********************start searles********************************


	sevap=sarea*serate

	svol(t)=svol(t-1)+c_spill-(sevap*dt)	!note that c_spill is a volume.  Therefore, it
											!is not necessary to multiply it by dt.
	if(svol(t).gt.76.0e9)then				
	s_spill=svol(t)-76.0e9
	svol(t)=76.0e9
	endif

	if(svol(t).le.0)then
	sarea=0.
	sdepth=0.
	else
	sarea=1.691e8+(0.008888*svol(t))-(3.4502e-14*(svol(t)**2))
	sdepth=6.3665+(4.5202e-9*svol(t))-(4.2614e-20*(svol(t)**2))+(2.1146e-31*(svol(t)**3))
	endif

!	********************start panamint******************************

	pevap=parea*perate

	pvol(t)=pvol(t-1)+s_spill-(pevap*dt)	!note that s_spill is a volume.  Therefore, it
											!is not necessary to multiply it by dt.
	if(pvol(t).gt.117.0e9)then
	p_spill=pvol(t)-117.0e9
	pvol(t)=117.0e9
	endif

	if(pvol(t).le.0)then
	parea=0.
	pdepth=0.
	else
	parea=1.1632e8+(0.0085949*pvol(t))-(2.4725e-14*(pvol(t)**2))
	pdepth=11.93+(5.4291e-9*pvol(t))-(4.6059e-20*(pvol(t)**2))+(1.7191e-31*(pvol(t)**3))
	endif

!	********************start lake manly****************************

	mevap=marea*merate

	mvol(t)=mvol(t-1)+p_spill-(mevap*dt)	!note that p_spill is a volume.  Therefore, it
											!is not necessary to multiply it by dt.
	if(mvol(t).le.0)then
	marea=0.
	mdepth=0.
	else
	marea=2.1382e8+(0.023849*mvol(t))-(1.9685e-13*(mvol(t)**2))+(6.5806e-25*(mvol(t)**3))
	mdepth=14.782+(2.1584e-9*mvol(t))-(1.3675e-20*(mvol(t)**2))+(4.2019e-32*(mvol(t)**3))
	endif

	write(17,'(f5.1,5(1x,f5.1))')time,odepth,cdepth,sdepth,pdepth,mdepth	!write data

	enddo

	end