program daisy2
implicit none	!Tells the program not to automatically designate number type to a variable.

real wdaisy,wdgrowth,wdbeta,wdseeds,wddeath,wdtemp		!Variables whose values are real numbers
real bdaisy,bdgrowth,bdbeta,bdseeds,bddeath,bdtemp
real X,death_rate,dt,fert,daisyUT,daisyLT,alb_planet
real stefboltz,solconst,wdalbedo,bdalbedo,gralbedo,luminosity
real solflux,q,planet_temp,t,tmax
integer i,imax									!Variables whose values are integers
dimension wdaisy(10000),bdaisy(10000),X(10000)	!Variables whose values are arrays

! Program to calculate the percentage of "black" and "white" daisies covering a planet that
! has a surface comprised of only black and white daisies and bare ground, as well as the
! temperatures of the daisies and ground. Written by Kirsten Menking and Joel Dashnaw,
! Dept. of Geology and Geography, Vassar College, Poughkeepsie, NY  12604, July 21, 2003.
! Based on the following readings:
! Watson, A.J., and Lovelock, J.E., 1983, Biological homeostasis of the global environment:
!	the parable of Daisyworld, Tellus, v. 35B, p. 284-289.

! The variables used in this model are as follows:
! RESERVOIRS:
! wdaisy = The area of the planet covered by white daisies, measured as a fraction of the
!	total planetary area (unitless).
! bdaisy = The area of the planet covered by black daisies, measured as a fraction of the
!	total planetary area (unitless).

! FLUXES:
! wdgrowth = The white daisy growth flux is an increase in wdaisy over time.
! wddeath = The white daisy death flux is a decrease in wdaisy over time.
! bdgrowth = The black daisy growth flux is an increase in bdaisy over time.
! bddeath = The black daisy death flux is a decrease in bdaisy over time.

! CONVERTERS:
! wdbeta = WD_beta is the growth rate of white daisies, which is a function of WD_temp.
!	In this model, the daisies, like most Earth species, have very specific ecological tolerances.
!	Daisies cannot grow at temperatures colder than 5 degrees C or warmer than 40 degrees C.
!	Their growth is maximized at 22.5 degrees.
! wdtemp = Temperature of the white daisies population in degrees Celsius.
! bdbeta = BD_beta is the growth rate of black daisies, which is a function of BD_temp.
!	In this model, the daisies, like most Earth species, have very specific ecological tolerances.
!	Daisies cannot grow at temperatures colder than 5 degrees C or warmer than 40 degrees C.
!	Their growth is maximized at 22.5 degrees.
! bdtemp = Temperature of the black daisies population in degrees Celsius.
! X = The area of available fertile ground, not covered by either species, measured as a
!	fraction of the total planetary area (dimensionless).
! death_rate = The death rate of daisies (dimensionless).
! fert = The proportion of the planet's area which is fertile ground (dimensionless).
! daisyUT = The upper temperature limit of the daisies (40 degrees C).
! daisyLT = The lower temperature limit of the daisies (5 degrees C).
! alb_planet = The percentage of incoming solar radiation that is reflected back into space
!	(dimensionless).
! stefboltz = The Stefan Boltzmann constant {5.67e-8 (J/m^2*s*K^4)}.  The S-B constant relates
!	the energy contained in a black body to its temperature.
! solconst = The solar constant is the amount of light reaching the surface of the earth.  For
!	the solar constant, Watson and Lovelock (1983) use a value of 9.17*105 ergs/cm2s, 
!	or 917 J/(m^2*s).
! wdalbedo = The percentage of incoming solar radiation that is reflected back into space by the
!	white daisies (dimensionless).
! bdalbedo = The percentage of incoming solar radiation that is reflected back into space by the
!	black daisies (dimensionless).
! gralbedo = The percentage of incoming solar radiation that is reflected back into space by the
!	bare ground (dimensionless).
! luminosity = Luminosity takes into account an increase in the amount of solar energy given off
!	by the sun over time (dimensionless).  This equation begins with the sun at half its current
!	brightness and ramps up the insolation by 2% a year.
! solflux = Solar flux is the change in the amount of insolation reaching the earth over time.
! Q = Q is a positive constant that expresses the degree to which absorbed solar energy is
!	redistributed to the three types of surfaces on the planet (dimensionless).
! planet_temp = Temperature of the entire planet in degrees Celsius.

! OTHER VARIABLES: 
! i = Counting loop incremental
! imax = 
! t = Time
! tmax = The maximum time that the program runs to
! dt = Time step in years.

wdseeds=0.001		!Arbitrary value used to limit the amount of depletion of wdaisy reservoir
bdseeds=0.001		!Arbitrary value used to limit the amount of depletion of bdaisy reservoir
wdaisy(1)=0.001		!Initial value of wdaisy: unitless
bdaisy(1)=0.001		!Initial value of bdaisy: unitless
death_rate=0.3		!unitless, but essentially per time
fert=1.0			!unitless
daisyUT=40.0		!degC
daisyLT=5.0			!degC
dt=0.5				!timestep
	
solconst=917.0				!J/(m^2*s)
stefboltz=5.67e-8			!J/(m^2*s*degK^4)
wdalbedo=0.75				!unitless
bdalbedo=0.25				!unitless
gralbedo=0.5				!unitless
t=0.0						
tmax=200.0											
imax=int((tmax+dt)/dt)		

	do i=2,imax

	bdaisy(i)=bdaisy(i-1)
	wdaisy(i)=wdaisy(i-1)

	X(i)=fert-bdaisy(i)-wdaisy(i)	!unitless

	!********************planetary temperature************************
	alb_planet=(bdalbedo*bdaisy(i))+(wdalbedo*wdaisy(i))+((1.-bdaisy(i)-wdaisy(i))*gralbedo)

	luminosity=0.5+(0.02*t)			!unitless
	solflux=solconst*luminosity		!!J/(m^2*s)
	q=(0.2*solflux)/stefboltz		!unitless

	planet_temp=(((solflux*(1-alb_planet))/stefboltz)**0.25)-273.15

	!********************white daisy temp*****************************
	wdtemp=(((q*(alb_planet-wdalbedo))+((planet_temp+273.16)**4))**(0.25))-273.16

	!********************black daisy temp*****************************
	bdtemp=(((q*(alb_planet-bdalbedo))+((planet_temp+273.16)**4))**(0.25))-273.16

	!********************white daisy growth***************************
	if(wdtemp.gt.daisyUT.or.wdtemp.lt.daisyLT)then
	wdbeta=0.0
	else
	wdbeta=(1.0-0.003265*((22.5-wdtemp)**2))
	endif

	wdgrowth=wdaisy(i)*X(i)*wdbeta

	if(wdaisy(i).gt.wdseeds)then
	wddeath=wdaisy(i)*death_rate
	else
	wddeath=0.
	endif

! Quantity in each reservoir = quantity at previous timestep + (inflows - outflows)*dt	
	wdaisy(i)=wdaisy(i)+((wdgrowth-wddeath)*dt)

!	********************black daisy growth***************************
	if(bdtemp.gt.daisyUT.or.bdtemp.lt.daisyLT)then
	bdbeta=0.0
	else
	bdbeta=(1.0-0.003265*((22.5-bdtemp)**2))
	endif

	bdgrowth=bdaisy(i)*X(i)*bdbeta

	if(bdaisy(i).gt.bdseeds)then
	bddeath=bdaisy(i)*death_rate
	else
	bddeath=0.
	endif

! Quantity in each reservoir = quantity at previous timestep + (inflows - outflows)*dt
	bdaisy(i)=bdaisy(i)+((bdgrowth-bddeath)*dt)

	write(*,*)i,t,bdaisy(i),wdaisy(i)	!write statement
	t=t+dt
	enddo

end