program Heat_Flow2
	implicit none

!	written by Kirsten Menking, Department of Geology and Geography, Vassar College.
!	This program calculates a temperature profile for permafrost regions based on the
!	assumption that heat flows only via conduction, and that conduction can be described
!	by Fourier's Law.  

	real::unit_area,air_tempK,layer_thickness,rock_density,rock_specific_heat
	real::t,dt,tmax,heat_capacity,sec_per_year,mean_heat_flow,k,geothermal_inflow
	real,dimension(50000)::layer1,layer2,layer3,layer4,layer5,layer6,layer7,layer8,layer9
	real,dimension(50000)::layer10,tz0,tz100,tz200,tz300,tz400,tz500,tz600,tz700,tz800,tz900
	real,dimension(50000)::to_air,q1,q2,q3,q4,q5,q6,q7,q8,q9
	integer::i,imax

	dt=1.0 ! time step in years

	unit_area=1.0 ! cross-sectional area of each permafrost layer, in square meters	
	air_tempK=268.15 ! air temperature in Kelvin, equivalent to -5 deg. C.
	layer_thickness=100.0 ! thickness of each permafrost layer, in meters
	rock_density=2700.0 ! density of each layer, in kg/m^3
	rock_specific_heat=800.0 ! specific heat of rock, in Joules/(kg*Kelvin)
	heat_capacity=layer_thickness*unit_area*rock_density*rock_specific_heat !volumetric heat
!	capacity of each layer of rock, in Joules/Kelvin
	sec_per_year=60.0*60.0*24.0*365.25 ! seconds per year
	mean_heat_flow=65.0e-3 ! Mean continental heat flow in W/m^2, use 101 for oceanic heat flow
	k=3.6*sec_per_year ! Thermal conductivity of average rock, in J/(m*K*yr)
!	average conductivity of salt = 6.1*sec_per_year, average conductivity of shale = 1.6*sec_per_year
	geothermal_inflow=mean_heat_flow*sec_per_year*unit_area ! inflow of heat into the bottom
!	layer of permafrost

	t=0.0 ! start time for model
	tmax=30000.0 ! end time for model
	imax=int((tmax+dt)/dt) ! maximum number of iterations

!Initial values of energy in each permafrost layer (Joules)***************************
	layer1(1)=0.
	layer2(1)=0.
	layer3(1)=0.
	layer4(1)=0.
	layer5(1)=0.
	layer6(1)=0.
	layer7(1)=0.
	layer8(1)=0.
	layer9(1)=0.
	layer10(1)=0.

!Alternative initial values of energy (Joules)*********************
!	layer1(1)=1155511.36	! these values are derived by running the model
!	layer2(1)=58140120105.89	! for 30,000 years with initial values of
!	layer3(1)=58364485003.41	! energy = 0.  They are the end values from
!	layer4(1)=58598012786.36	! that initial run and are used to start the
!	layer5(1)=58844971608.95	! second run of 30,000 years.  
!	layer6(1)=59109263261.61
!	layer7(1)=59394316740.35
!	layer8(1)=59702994712.71
!	layer9(1)=60037515431.60
!	layer10(1)=60399392227.20

	do i=2,imax

!Determination of layer temperatures (Kelvin)***************************
	tz0(i)=air_tempK
	tz100(i)=layer2(i-1)/heat_capacity
	tz200(i)=layer3(i-1)/heat_capacity
	tz300(i)=layer4(i-1)/heat_capacity
	tz400(i)=layer5(i-1)/heat_capacity
	tz500(i)=layer6(i-1)/heat_capacity
	tz600(i)=layer7(i-1)/heat_capacity
	tz700(i)=layer8(i-1)/heat_capacity
	tz800(i)=layer9(i-1)/heat_capacity
	tz900(i)=layer10(i-1)/heat_capacity


!Flows of energy from layer to layer (Joules/yr)****************************

		to_air(i)=layer1(i)
		q1(i)=k*(tz100(i)-tz0(i))*unit_area/layer_thickness
		q2(i)=k*(tz200(i)-tz100(i))*unit_area/layer_thickness
		q3(i)=k*(tz300(i)-tz200(i))*unit_area/layer_thickness
		q4(i)=k*(tz400(i)-tz300(i))*unit_area/layer_thickness
		q5(i)=k*(tz500(i)-tz400(i))*unit_area/layer_thickness
		q6(i)=k*(tz600(i)-tz500(i))*unit_area/layer_thickness
		q7(i)=k*(tz700(i)-tz600(i))*unit_area/layer_thickness
		q8(i)=k*(tz800(i)-tz700(i))*unit_area/layer_thickness
		q9(i)=k*(tz900(i)-tz800(i))*unit_area/layer_thickness
		
!Updating of energy held in each permafrost layer (Joules)********************

	layer10(i)=layer10(i-1)+(geothermal_inflow-q9(i))*dt	
	layer9(i)=layer9(i-1)+(q9(i)-q8(i))*dt
	layer8(i)=layer8(i-1)+(q8(i)-q7(i))*dt		
	layer7(i)=layer7(i-1)+(q7(i)-q6(i))*dt
	layer6(i)=layer6(i-1)+(q6(i)-q5(i))*dt		
	layer5(i)=layer5(i-1)+(q5(i)-q4(i))*dt		
	layer4(i)=layer4(i-1)+(q4(i)-q3(i))*dt		
	layer3(i)=layer3(i-1)+(q3(i)-q2(i))*dt		
	layer2(i)=layer2(i-1)+(q2(i)-q1(i))*dt		
	layer1(i)=layer1(i-1)+(q1(i)-to_air(i))*dt		

	write(15,*)i,tz0(i),tz100(i),tz200(i),tz300(i),tz400(i)
	write(16,*)tz500(i),tz600(i),tz700(i),tz800(i),tz900(i)
	enddo
	
end