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