program complex_atmosphere

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::A_density,A_SH,L1_R,L1_depth,L1_SA,L1_HC,L2_R,L2_depth,L2_SA,L2_HC
real::Loss_to_Space,Convection,LH_to_Layer1,LH_to_Layer2,Solar_to_Layer1,Solar_to_Layer2
real,dimension(700)::EarthEnergy,Layer1,Layer2,TempK,TempC,Layer1_TempK,Layer1_TempC
real,dimension(700)::Layer2_TempK,Layer2_TempC,Surf_to_Layer1,Layer1_to_Surf
real,dimension(700)::Layer1_to_Layer2,Layer2_to_Layer1,Layer2_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.
! Layer1 = The amount of joules stored in atmosphere layer one of an atmosphere comprised of
!	two layers.
! Layer2 = The amount of joules stored in atmosphere layer two of an atmosphere comprised of
!	two layers.
		
! FLUXES:
! Solar_to_Earth = Incoming short-wave radiation emitted from the Sun to the surface of the
!	Earth over time.
! Surf_to_Layer1 = Long-wave radiation emitted from the surface of the Earth to Layer1 over time.
! Layer1_to_Surf = Long-wave radiation emitted from Layer1 to the surface of the Earth over time.
! Layer1_to_Layer2 = Long-wave radiation emitted from Layer1 to Layer2 over time.
! Layer2_to_Layer1 = Long-wave radiation emitted from Layer2 to Layer1 over time.
! Layer2_to_Space = Long-wave radiation emitted from Layer2 to Space over time.

! Loss_to_Space = This flow represents 20 W/m^2 of emitted infrared heat that is radiated 
!	directly to space without being absorbed into the atmosphere.
! Convection = This flow takes into account heat flow in the form of convection from the 
!	Earth's surface to the atmosphere.  This is about 17 W/m^2.
! LH_to_Layer1 = This flow takes into account 80 W/m^2 of heat added to the atmosphere from
!	the vaporization of water at the surface.  This quantity is divided equally between the
!	two layers of the atmosphere as the amount of water vapor in the atmosphere is distributed 
!	evenly between the two.  Thus, the flow is as follows:
!	Latent_heat_to_Layer_1 = 40*Seconds_per_Year*Earth_Surface_Area
!	For more information, refer to the reading by Harte, J., 1988.
! LH_to_Layer2 = This flow takes into account 80 W/m^2 of heat added to the atmosphere from
!	the vaporization of water at the surface.  This quantity is divided equally between the
!	two layers of the atmosphere as the amount of water vapor in the atmosphere is distributed
!	evenly between the two.  Thus, the flow is as follows:
!	Latent_heat_to_Layer_2 = 40*Seconds_per_Year*Earth_Surface_Area
!	For more information, refer to the reading by Harte, J., 1988.
! Solar_to_Layer1 = Incoming solar energy absorbed directly into atmosphere layer 1.
!	About 30% of the 86 W/m^2 absorbed into the atmosphere is absorbed into the lower atmosphere
!	layer (Layer 1).  See the following reading:
!	Harte, J., 1988, Consider a Spherical Cow, Sausalito, CA:  University Science Books, p. 69-72, 
!	p. 160-167.
! Solar_to_Layer2 = Incoming solar energy absorbed directly into atmosphere layer 2.
!	About 70% of the 86 W/m^2 absorbed into the atmosphere is absorbed into the upperr atmosphere
!	layer (Layer 2).  See the following reading:
!	Harte, J., 1988, Consider a Spherical Cow, Sausalito, CA:  University Science Books, p. 69-72, 
!	p. 160-167.
	
! 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
! A_density = Density of Air
! A_SH = Specific Heat of Air
! L1_depth = Depth of atmosphere layer one in the two layer atmosphere.  Initially this value 
!	is 5 kilometers, however this number can be changed for experimentation.
! L1_R = Radius of the Earth plus the depth of atmosphere layer one
! L1_SA = Surface Area of atmosphere layer one in the two layer atmosphere
! L1_HC = Heat Capacity of atmosphere layer one in the two layer atmosphere
! L2_depth = Depth of atmosphere layer two in the two layer atmosphere.  Initially this value
!	is 5 kilometers, however this number can be changed for experimentation.
! L2_R = Radius of the Earth plus the depth of atmosphere layer one and atmosphere layer two
! L2_SA = Surface area of atmosphere layer two in the two layer atmosphere
! L2_HC = Heat Capacity of atmosphere layer two in the two layer atmosphere
! SPY = Seconds per year
! SC = Solar Constant, which is the incoming solar flux multiplied by seconds per year
! SB = Stefan Boltzmann constant multiplied by seconds per year
! TempK = Temperature of the Earth's surface in degrees Kelvin
! TempC = 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
Layer1(1)=0.		!Joules
Layer2(1)=0.		!Joules

TempK(1)=0.			!degrees K: At time=1 initial value=0
TempC(1)=0.			!degrees C: At time=1 initial value=0
Layer1_TempK(1)=0.	!degrees K: At time=1 initial value=0
Layer1_TempC(1)=0.	!degrees C: At time=1 initial value=0
Layer2_TempK(1)=0.	!degrees K: At time=1 initial value=0
Layer2_TempC(1)=0.	!degrees C: At time=1 initial value=0

!Initial Values of Flows*******************************************
Surf_to_Layer1(1)=0.		!Joules/years: At time=1 initial value=0
Layer1_to_Surf(1)=0.		!Joules/years: At time=1 initial value=0
Layer1_to_Layer2(1)=0.		!Joules/years: At time=1 initial value=0
Layer2_to_Layer1(1)=0.		!Joules/years: At time=1 initial value=0
Layer2_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
A_density=1.3							!Kilograms/meters^3
A_SH=750.0								!Joules/Kilograms*degrees K
L1_depth=5.0e3							!meters
L1_R=6371.0e3+L1_depth					!meters
L1_SA=4.0*3.1415*(L1_R**2)				!meters^2
L1_HC=L1_depth*L1_SA*A_density*A_SH		!Joules/degrees K
L2_depth=5.0e3							!meters
L2_R=6371.0e3+L1_depth+L2_depth			!meters
L2_SA=4.0*3.1415*(L2_R**2)				!meters^2
L2_HC=L2_depth*L2_SA*A_density*A_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))-(86.0*SPY*E_SA)		!Joules/years
Loss_to_Space=20.0*SPY*E_SA				!Joules/years
Convection=17.0*SPY*E_SA				!Joules/years
LH_to_Layer1=40.0*SPY*E_SA				!Joules/years
LH_to_Layer2=40.0*SPY*E_SA				!Joules/years
Solar_to_Layer1=0.3*86.0*SPY*E_SA		!Joules/years
Solar_to_Layer2=0.7*86.0*SPY*E_SA		!Joules/years
tmax=5.									
imax=int((tmax+dt)/dt)					

	do i=2,imax
	
! Quantity in each reservoir = quantity at previous timestep + (inflows - outflows)*dt	
	EarthEnergy(i)=EarthEnergy(i-1)+(Solar_to_Earth+Layer1_to_Surf(i-1)-Surf_to_Layer1(i-1)-Loss_to_Space-LH_to_Layer1-LH_to_Layer2-Convection)*dt
		TempK(i)=EarthEnergy(i-1)/W_HC
		TempC(i)=TempK(i)-273.15
		Surf_to_Layer1(i)=(E_SA*SB*(TempK(i)**4))-(20.0*SPY*E_SA)

	Layer1(i)=Layer1(i-1)+(Surf_to_Layer1(i-1)+Layer2_to_Layer1(i-1)+LH_to_Layer1+Solar_to_Layer1+Convection-Layer1_to_Surf(i-1)-Layer1_to_Layer2(i-1))*dt
		Layer1_TempK(i)=Layer1(i-1)/L1_HC
		Layer1_TempC(i)=Layer1_TempK(i)-273.15
		Layer1_to_Layer2(i)=L1_SA*SB*(Layer1_TempK(i)**4)
		Layer1_to_Surf(i)=L1_SA*SB*(Layer1_TempK(i)**4)
	
	Layer2(i)=Layer2(i-1)+(Layer1_to_Layer2(i-1)+LH_to_Layer2+Solar_to_Layer2-Layer2_to_Layer1(i-1)-Layer2_to_Space(i-1))*dt
		Layer2_TempK(i)=Layer2(i-1)/L2_HC
		Layer2_TempC(i)=Layer2_TempK(i)-273.15
		Layer2_to_Space(i)=L2_SA*SB*(Layer2_TempK(i)**4)
		layer2_to_Layer1(i)=L2_SA*SB*(Layer2_TempK(i)**4)

	t=t+dt
	
	write(*,*)i,t,TempK(i),TempC(i)					!Write output file
	write(*,*)Layer1_TempK(i),Layer1_TempC(i)		!Write output file
	write(*,*)Layer2_TempK(i),Layer2_TempC(i)		!Write output file
	
	enddo

end