program blackbody_plus_albedo

implicit none
real::t,tmax,dt,E_R,E_CS,E_SA,W_depth,W_density,W_SH,W_HC,SPY,SC,SB,Solar_to_Earth,Albedo
real,dimension(700)::EarthEnergy,TempK,TempC,Surfheat_to_Space
integer::i,imax

! Program to calculate the temperature of the surface of the Earth and one layer of atmosphere.
! Written by Joel Dashnaw and Kirsten Menking, Dept. of Geology and Geography, Vassar College,
! Poughkeepsie, NY  12604, July 7, 2003.  Based on the following readings:
! Few, A.A., 1996, System Behavior and System Modeling, Sausalito, CA: University Science 
!	Books, p. 27-37.
! Graedel, T.E., and Crutzen, P.J., 1993, Atmospheric Change:  An Earth System Perspective,
!	New York: W.H. Freeman and Company, 446 p.
! Harte, J., 1988, Consider a Spherical Cow, Sausalito, CA:  University Science Books, 
!	p. 69-72, p. 160-167.

! The variables used in this model are as follows:
! RESERVOIRS:
! EarthEnergy = This model is a "swamp model", which assumes that the whole of the Earth's 
!	surface is covered by a layer of water of a certain depth (in this case 1 meter).  This
!	is done to simplify the model, as the Earth's surface can then be considered to have uniform 
!	conditions.  Earth Energy reservoir is the amount of joules stored in the Earth's swamp
!	surface.  
		
! FLUXES:
! Solar_to_Earth = Incoming short-wave radiation emitted from the Sun to the surface of the
!	Earth over time.
! Surfheat_to_Space = Outgoing long-wave radiation emitted from the surface of the Earth to
!	Space over time.
	
! CONVERTERS:
! Albedo = The percentage of incoming solar radiation that is reflected back into space.
!	This fraction is multiplied by the solar constant and the cross section of
!	the Earth, and is a unitless number.  For the Earth, the albedo = 30%, however
!	this number can be changed to represent increases or decreases in cloud and
!	ice cover.
! E_R = Radius of the Earth
! E_CS = Cross section of the Earth
! E_SA = Surface Area of the Earth
! W_depth = Water depth.  The initial assumption of depth of the layer of water ("swamp") 
!	covering the Earth is 1.0 meter.  This may be changed for experimentation.
! W_density = Density of water
! W_SH = Specific Heat of water
! W_HC = Heat Capacity of water
! SPY = Seconds per year
! SC = The solar constant, which is the incoming solar flux multiplied by seconds per year.
! SB = The Stefan Boltzmann constant multiplied by seconds per year.
! TempK = The temperature of the Earth's surface in degrees Kelvin.
! TempC = The temperature of the Earth's surface in degrees Celsius.
	
! OTHER VARIABLES: 
! i = Counting loop incremental
! imax = 
! t = Time in years.
! tmax = 
! dt = Time step in years.

dt=(1./128.)
t=0.

!Initial values of reservoirs**************************************
EarthEnergy(1)=0.	!Joules
TempK(1)=0.			!degrees K: At time=1 initial value=0
TempC(1)=0.			!degrees K: At time=1 initial value=0

!Initial Values of Flows*******************************************
Surfheat_to_Space(1)=0.		!Joules/years: At time=1 initial value=0

!Equations and values of Program Variables************************
Albedo=0.3							!Unitless
E_R=6371.0e3						!meters
E_CS=3.1415*(E_R**2)				!meters^2
E_SA=4.0*3.1415*(E_R**2)			!meters^2
W_depth=1.0							!meters
W_density=1000.0					!Kilograms/meters^3
W_SH=4218.0							!Joules/Kilograms*degrees K
W_HC=W_depth*E_SA*W_density*W_SH	!Joules/degrees K
SPY=3.15576e7						!seconds/years
SC=1368.0*SPY						!Joules/meters^2*years
SB=5.67e-8*SPY						!Joules/meters^2*years*degrees K^4
Solar_to_Earth=SC*E_CS*(1-Albedo)	!Joules/years
tmax=5.								
imax=int((tmax+dt)/dt)				

	do i=2,imax
	
! Quantity in reservoir = quantity at previous timestep + (inflows - outflows)*dt

	EarthEnergy(i)=EarthEnergy(i-1)+(Solar_to_Earth-Surfheat_to_Space(i-1))*dt
	TempK(i)=EarthEnergy(i-1)/W_HC
	TempC(i)=TempK(i)-273.15
	Surfheat_to_Space(i)=E_SA*SB*(TempK(i)**4)
	
	t=t+dt
	write(*,*)i,t,TempK(i),TempC(i)		!Write output file

	enddo

end