program phosphorus

	! This model simulates the amount of Phosphorous in six major reservoirs of the closed system 
	! of the Earth.  Written by Kirsten Menking, Dept. of Geology and Geography, Vassar College,
	! Poughkeepsie, NY  12604 and based on a model shown on pg. 103 of Chameides, W.L., and 
	! Perdue, E.M., 1997, Biogeochemical cycles:  A computer-interactive study of Earth system 
	! science and global change, New York:  Oxford University Press, 224 p.

	! The variables used in the model are as follows:
	! RESERVOIRS:
	! landbio = the amount of Phosphorus (P) found in land plants and animals
	! docean = the amount of P found in the deep ocean
	! oceanbio = the amount of P found in marine plants and animals
	! seds = the amount of P found in marine sediments and sedimentary rocks
	! socean = the amount of P found in the surface ocean
	! soils = the amount of P found in soils

	! FLUXES:
	! decay = the transfer of P from land plants and animals to soils via decay of dead organic 
	! matter

	! mixing = the transfer of P from the surface ocean to the deep ocean

	! ocdecay = the loss of P from marine plants and animals to the deep ocean via decay 
	! of dead organic matter

	! ocprod = the transfer of P from the surface ocean into marine plants and animals via primary
	! production

	! sediment = the transfer of P from the deep ocean into oceanic sediments via phosphate 
	! mineralization

	! sinking = the transfer of P from marine plants and animals to the deep ocean via sinking
	! of organic particles in the water column

	! solwx = the transfer of P from soils to oceanic sediments via insolube weathering and
	! erosion processes

	! terrprod = the uptake of P from the soils by land plants

	! uplift = the transfer of P from deep ocean sediments to soils via uplift of sedimentary 
	! rocks and weathering

	! upwell = transfer of P from the deep to the surface ocean via upwelling

	! wxing = transfer of P from soils to oceanic sediments via weathering

	!OTHER VARIABLES
	! t = time in days
	! dt = time step in days
	! i = counting loop incrementer

	
	implicit none
	real docean,landbio,oceanbio,socean,soils
	real(kind=selected_real_kind(p=15)) seds !This tells the program that seds has to have at 
	!least 15 digits.  The reason this is needed is that the seds reservoir is enormous compared 
	!to the other reservoirs and not enough significant digits are retained to carry out the
	!calculations properly without this specification.
	real decay,mixing,ocdecay,ocprod,sediment,sinking,solwx,terrprod,uplift,upwell,wxing
	real dt,t
	integer i
	dimension soils(5000),landbio(5000),seds(5000),docean(5000)
	dimension socean(5000),oceanbio(5000)
	dimension uplift(5000),decay(5000),terrprod(5000),wxing(5000),solwx(5000)
	dimension sediment(5000),sinking(5000),mixing(5000),upwell(5000),ocprod(5000)
	dimension ocdecay(5000)

	open(unit=15,file='landres.txt',status='replace')	! Creates a file name for unit 15
	open(unit=16,file='oceres.txt',status='replace')	! Creates a file name for unit 16
	
	! Initial Values of Reservoirs*********************************************
	soils(1)=2.0e5		!Tg (Teragrams): 1Tg=1*10^12g
	landbio(1)=2600.	!Tg
	seds(1)=2.e9		!Tg
	docean(1)=1.e5		!Tg
	socean(1)=2800.		!Tg
	oceanbio(1)=44.		!Tg
	dt=0.03125			!timestep in days (1/32 of a day)
	t=1.				!day 1

	do i=2,3000
	
	! Equations for Flows*****************************************************
	! The equation for each flow is: reservoir * transfer coefficient of the flow. The transfer
	! coefficient of the flow is equal to (the flow out of the reservoir / contents of the 
	! reservoir when the system is at steady state). Unless a flow is altered, the value of the 
	! flow will be equal to the flow out of the reservoir at steady state, as the reservoir value 
	! is multiplied by the flow/steady state value of the reservoir.  
		! ex. reservoir A = 20.
		! flow from reservoir A to reservoir B = reservoir A * (flow out of reservoir/20.)
	
	decay(i)=(310./2600.)*landbio(i-1)		!Tg/yr
	mixing(i)=(18./2800.)*socean(i-1)		!Tg/yr
	ocdecay(i)=(940./44.)*oceanbio(i-1)		!Tg/yr
	ocprod(i)=(980./2800.)*socean(i-1)		!Tg/yr
	sediment(i)=(2./1.e5)*docean(i-1)		!Tg/yr
	sinking(i)=(40./44.)*oceanbio(i-1)		!Tg/yr
	solwx(i)=(2./2.e5)*soils(i-1)			!Tg/yr
	terrprod(i)=(310./2.e5)*soils(i-1)		!Tg/yr
	uplift(i)=(20./2.e9)*seds(i-1)			!Tg/yr
	upwell(i)=(56./1.e5)*docean(i-1)		!Tg/yr
	wxing(i)=(18./2.e5)*soils(i-1)			!Tg/yr
	


	! Equations to update the amounts of P in each reservoir ****************************
	! Amount of P in each reservoir = amount of P in previous timestep + (inflows - outflows)*dt
	
	docean(i)=docean(i-1)+((sinking(i-1)+mixing(i-1)-upwell(i-1)-sediment(i-1))*dt)
	landbio(i)=landbio(i-1)+((terrprod(i-1)-decay(i-1))*dt)
	oceanbio(i)=oceanbio(i-1)+((ocprod(i-1)-ocdecay(i-1)-sinking(i-1))*dt)
	seds(i)=seds(i-1)+((wxing(i-1)+sediment(i-1)-uplift(i-1))*dt)
	socean(i)=socean(i-1)+((solwx(i-1)+upwell(i-1)+ocdecay(i-1)-mixing(i-1)-ocprod(i-1))*dt)
	soils(i)=soils(i-1)+((uplift(i-1)+decay(i-1)-terrprod(i-1)-wxing(i-1)-solwx(i-1))*dt)

	write(15,*)t,soils(i),landbio(i),seds(i) ! write output file
	write(16,*)t,docean(i),socean(i),oceanbio(i) ! write output file

	t=t+dt ! update the total time elapsed for each incrementing of the counting loop

	enddo

	end