program two_layer_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,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.
! 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) !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))*dt
TempK(i)=EarthEnergy(i-1)/W_HC
TempC(i)=TempK(i)-273.15
Surf_to_Layer1(i)=E_SA*SB*(TempK(i)**4)
Layer1(i)=Layer1(i-1)+(Surf_to_Layer1(i-1)+Layer2_to_Layer1(i-1)-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)-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