PRO CH4_2BOX_MODERND, dev=dev

       ;IDL procedure to explore the effects of changing CH4 sources on
       ;concentration and isotopic gradients in a box model of the 
       ;atmosphere with 2 hemispheres.

CCG_OPENDEV, dev=dev, pen=pen, portrait=1

;PLOT OPTONS:
show=1          ;mixing ratios, isotopic ratios, differences
;show=2          ;isotopic forcing function

;variables:
;	Xn,Xs	Northern and southern hemisphere CH4 amounts  
;	I13n,I13s	Northern and southern hemisphere 13CH4 amounts
;	R13n,R13s	Isotopic ratios
;	Del13n,Del13s   Conventional delta values relative to PDB
;	I2n,I2s	Northern and southern hemisphere 13CH4 amounts
;	R2n,R2s	Isotopic ratios
;	Del2n,Del2s   Conventional delta values relative to PDB
;	Qs,Delqs,Rqs; Qn,Delqn,Rqn
;       Source, its del value, and its isotopic ratio 
;	exch    Interhemispheric mixing time in years
;	alpha13   Isotopic fractionation during destruction
;	alpha2   Isotopic fractionation during destruction

;Initial conditions and constants
step=0.20
time=0.D0
stop=510.

xran=[00,stop-10]         ;years to be plotted

nsteps=stop/step+1
steps=DINDGEN(nsteps)
PDB=0.01111       ;Isotopic standard defined as 13C/(12C+13C)
SMOW=0.00015573
exch=1.



;****************timee-varying parameters*********************************
Xn=DBLARR(nsteps)
Xs=DBLARR(nsteps)
I13n=DBLARR(nsteps)
I13s=DBLARR(nsteps)
I2n=DBLARR(nsteps)
I2s=DBLARR(nsteps)
Qn=DBLARR(nsteps)
Qs=DBLARR(nsteps)
R13qn=DBLARR(nsteps)
R13qs=DBLARR(nsteps)
R13n=DBLARR(nsteps)
R13s=DBLARR(nsteps)
Del13qn=DBLARR(nsteps)
Del13qs=DBLARR(nsteps)
Del13n=DBLARR(nsteps)
Del13s=DBLARR(nsteps)
R2qn=DBLARR(nsteps)
R2qs=DBLARR(nsteps)
R2n=DBLARR(nsteps)
R2s=DBLARR(nsteps)
Del2qn=DBLARR(nsteps)
Del2qs=DBLARR(nsteps)
Del2n=DBLARR(nsteps)
Del2s=DBLARR(nsteps)
t=DBLARR(nsteps)
bact=DBLARR(nsteps)
ff=DBLARR(nsteps)
bmb=DBLARR(nsteps)
Lfills=DBLARR(nsteps)
Qa1=DBLARR(nsteps)
Qnp=DBLARR(nsteps)
Qsp=DBLARR(nsteps)
Q1n=DBLARR(nsteps)
Q1s=DBLARR(nsteps)
Del13q1n=DBLARR(nsteps)
Del13q1s=DBLARR(nsteps)
R13q1n=DBLARR(nsteps)
R13q1s=DBLARR(nsteps)
Del2q1n=DBLARR(nsteps)
Del2q1s=DBLARR(nsteps)
R2q1n=DBLARR(nsteps)
R2q1s=DBLARR(nsteps)
;************************************************************************

       
;************source functions********
;hold first one hundred years constant so model can initialize smoothly


Bact(0:nsteps-1)=355.			;sources (Tg/yr)
;Bact=(228+12*cos((steps-100/step)/318))/1
;Bact=200+(205.-200.)/200.*steps
;Bact=239+1*exp((steps-100/step)/217)
;Bact(0:100/step)=240.
;Bact(200/step:500/step)=240.

FF(0:nsteps-1)=100.
;ff=(228+12*cos((steps-100/step)/318))/2
;FF=0+(100.-0.)/200.*steps*.2
;FF=-1+1*exp((steps-100/step)/217)
;FF(0:100/step)=0

BMB(0:nsteps-1)=56.
;BMB(501:nsteps-1)=10
;BMB=4.9+0.1*exp((steps-100/step)/267)
;BMB(0:500)=10
;BMB(0:200/step)=5

Lfills(0:nsteps-1)=35.

Qa1=(Bact+FF+BMB+Lfills)*1d12		;total global source in g(Ch4)/yr

Bact_N=0.716			;northern hemisphere fractions of sources
FF_N=0.95				;from Fung et al. 1991 based on zero deg equator
BMB_N=0.62
Lfills_N=0.96                    

bact_s=1-bact_n			;southern hemisphere fractions	
ff_s=1-ff_n
bmb_s=1-bmb_n
Lfills_s=1-lfills_N

Qnp=(Bact*Bact_N+FF*FF_N+BMB*BMB_N+Lfills*Lfills_n)/Qa1*1d12 	;n/s source distribution in g/yr
Qsp=1-Qnp

Qn1=Qa1*Qnp/(2.767e21/2)*1.E9	 ;sources mixing into each hemisphere
Qs1=Qa1*Qsp/(2.767e21/2)*1.E9
 
Qn=Qn1 	
Qs=Qs1
;********isotopic source functions

Del13qn=((-61*bact*bact_n)+(-40*ff*ff_n)+(-25*bmb*bmb_n)+(-50*lfills*lfills_n))/(qa1*qnp/1d12)
Del13qs=((-61*bact*bact_s)+(-40*ff*ff_s)+(-25*bmb*bmb_s)+(-50*lfills*lfills_s))/(qa1*qsp/1d12)

R13qn=(1.D0+Del13qn/1000.)*PDB
R13qs=(1.D0+Del13qs/1000.)*PDB


Del2qn=((-310*bact*bact_n)+(-175*ff*ff_n)+(-210*bmb*bmb_n)+(-293*lfills*lfills_n))/(qa1*qnp/1d12)
Del2qs=((-310*bact*bact_s)+(-175*ff*ff_s)+(-210*bmb*bmb_s)+(-293*lfills*lfills_s))/(qa1*qsp/1d12)

R2qn=(1.D0+Del2qn/1000.)*SMOW
R2qs=(1.D0+Del2qs/1000.)*SMOW
;**************sink functions and constants******************************

soil_n=7/2.767			;sink functions (note: soil sink is 0-th order)
soil_s=3/2.767			;ppb/yr
soil_n=0			;sink functions (note: soil sink is 0-th order)
soil_s=0
;oh=DBLARR(nsteps)+0.10215	;modern
oh=DBLARR(nsteps)+0.11236	;paleo	;Wang et al. ??
strat=DBLARR(nsteps)+0.00851	;weighted average of Cl, OH, and O(1D) after Hein et al.
;oh=0.11798-0.00562*cos((steps-100/step)*3.14/1000*10.745.)
;oh=.11049+((0.11049-0.12212)/200)*steps
oh(0:100/step)=0.11236
;oh(200/step:300/step)=0.12360

;************sink fractionation factors
a13_oh=0.9946			;Cantrell et al
a13_strat=0.988			;weighted average of OH, Cl, and O(1D) according to Hein 97
a13_soil=0.979			;King et al

;a2_oh=0.7450
a2_oh=0.7690
;a2_oh=0.8370
arr1=[0.745,0.769,0.837]			;in Snover et al, 2000
;a2_oh=MEAN(arr1)
a2_strat=0.840			;
a2_soil=0.919			;
;*************initial 13C and 12C
n=0
Xn(0)=770.
Xs(0)=720.
I13n(0)=0.01055*Xn(0)
I13s(0)=0.01055*Xs(0)
Del13n(n)=(I13n(n)/Xn(n)/PDB-1.)*1000.
Del13s(n)=(I13s(n)/Xs(n)/PDB-1.)*1000.
I2n(0)=0.00016*Xn(0)
I2s(0)=0.00016*Xs(0)
Del2n(n)=(I2n(n)/Xn(n)/SMOW-1.)*1000.
Del2s(n)=(I2s(n)/Xs(n)/SMOW-1.)*1000.

;********solve mass balance equations and iterate towards steady-state**********
  
WHILE time LT stop-0.0001 DO BEGIN
   
   n=n+1
   time=time+step
   t(n)=time

   Xn(n)=Xn(n-1)+step*(Qn(n)-soil_n-(oh(n)+strat(n))*Xn(n-1)-exch*(Xn(n-1)-Xs(n-1)))
   Xs(n)=Xs(n-1)+step*(Qs(n)-soil_s-(oh(n)+strat(n))*Xs(n-1)+exch*(Xn(n-1)-Xs(n-1)))
   I13n(n)=I13n(n-1)+step*(R13qn(n)*Qn(n)-a13_soil*soil_n*(I13n(n-1)/Xn(n-1))-(a13_oh*oh(n)+a13_strat*strat(n))*I13n(n-1)-exch*(I13n(n-1)-I13s(n-1)))
   I13s(n)=I13s(n-1)+step*(R13qs(n)*Qs(n)-a13_soil*soil_s*(I13s(n-1)/Xs(n-1))-(a13_oh*oh(n)+a13_strat*strat(n))*I13s(n-1)+exch*(I13n(n-1)-I13s(n-1)))
   Del13n(n)=(I13n(n)/Xn(n)/PDB-1.)*1000.
   Del13s(n)=(I13s(n)/Xs(n)/PDB-1.)*1000.
   I2n(n)=I2n(n-1)+step*(R2qn(n)*Qn(n)-a2_soil*soil_n*(I2n(n-1)/Xn(n-1))-(a2_oh*oh(n)+a2_strat*strat(n))*I2n(n-1)-exch*(I2n(n-1)-I2s(n-1)))
   I2s(n)=I2s(n-1)+step*(R2qs(n)*Qs(n)-a2_soil*soil_s*(I2s(n-1)/Xs(n-1))-(a2_oh*oh(n)+a2_strat*strat(n))*I2s(n-1)+exch*(I2n(n-1)-I2s(n-1)))
   Del2n(n)=(I2n(n)/Xn(n)/SMOW-1.)*1000.
   Del2s(n)=(I2s(n)/Xs(n)/SMOW-1.)*1000.
ENDWHILE

;PLOTS*********************************************************************
plottitle='Increase of BMB from 5 to 10 Tg/yr'
;plottitle='Exponential increase of BMB from 5 to 10 Tg/yrv'
;plottitle='LOWER Rq 4 PERMIL, Q UNCHANGED'
;plottitle='DOUBLE Q, LOWER Rq BY 4 PERMIL'

IF show EQ 1 THEN BEGIN

	PLOT,t,Xn, $
   		POSITION=[0.1,0.70,0.47,0.90], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		thick=1.5, $
   		XTICKNAME=replicate(' ',6), $
		;YRANGE=[400,800], $
   		YTICKS=4, $
		;YTICKNAME=[' ','500','1000','1500','2000'], $
		;YRANGE=[0,1200], $
		;YTICKS=6, $
		;YTICKNAME=[' ','200','400','600','800','1000','1200'], $
  		YCHARSIZE=0.9, $
   		YTITLE='ppb'

	OPLOT,t,Xs, $
   		COLOR=pen(2), $
   		THICK=1.5 ;, $
   		;LINESTYLE=2

	tarr=['NH','SH']
	larr=[0,2]
	;carr=[pen(1),pen(2)]
	carr=[pen(1),pen(2)]
	;CCG_LLEGEND,x=0.35,y=0.76,tarr=tarr,larr=larr,carr=carr
	;XYOUTS,/NORMAL,0.31,0.78,'mixing ratios:'
	CCG_LLEGEND,x=0.29,y=0.76,tarr=tarr,larr=larr,carr=carr
	XYOUTS,/NORMAL,0.25,0.78,'mixing ratios:'
	XYOUTS,/NORMAL,0.13,0.87,'A'

	;XYOUTS,/NORMAL,0.00,.94,plottitle,charsize=1.5
	XYOUTS,/NORMAL,0.1,.94,plottitle,charsize=1.5
	
	PLOT,t,Xn-Xs, $
   		/NOERASE, $
   		POSITION=[0.1,0.50,0.47,0.7], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		thick=1.5, $
   		XTICKNAME=replicate(' ',6), $
   		YRANGE=[30,60], $
   		YTICKS=4, $
   		YCHARSIZE=0.9, $
   		;YTICKNAME=[' ','50','100','150',' '], $
   		YTITLE='ppb'

	XYOUTS,/NORMAL,0.25,0.54,'difference NH-SH'
	XYOUTS,/NORMAL,0.13,0.67,'B'
  ;************13C**********************
	PLOT,t,Del13n, $
   		/NOERASE, $
   		POSITION=[0.1,0.30,0.47,0.5], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		thick=1.5, $
   		XTICKNAME=replicate(' ',6), $
		YRANGE=[MIN([Del13n,Del13s])-1,MAX([Del13n,Del13s])+1], $
   		;YRANGE=[-55,-54], $
   		;YTICKS=7, $
		;YTICKNAME=[' ','-49','-48','-47','-46','-45','-44',' '], $
   		;YMINOR=1, $
   		YCHARSIZE=0.9, $
   		;ystyle=1, $
   		YTITLE='permil (PDB)'
	OPLOT,t,Del13s, $
   		color=pen(2), $
   		thick=1.5

	tarr=['NH','SH']
	larr=[0,2]
	;carr=[pen(1),pen(2)]
	carr=[pen(1),pen(2)]
	CCG_LLEGEND,x=0.20,y=0.36,tarr=tarr,larr=larr,carr=carr
	XYOUTS,/NORMAL,0.15,0.38,'isotopic ratios:'
	;CCG_LLEGEND,x=0.29,y=0.36,tarr=tarr,larr=larr,carr=carr
	;XYOUTS,/NORMAL,0.25,0.38,'isotopic ratios:'
	XYOUTS,/NORMAL,0.13,0.47,'C'

	PLOT,t,Del13n-Del13s, $
   		/NOERASE, $
   		POSITION=[0.1,0.10,0.47,0.3], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		THICK=1.5, $
   		XTITLE='year', $
   		XCHARSIZE=1.0, $
   		YRANGE=[-1.,0.0], $
   		YTICKS=5, $
   		;YTICKNAME=['-1.0','-0.8','-0.6','-0.4','-0.2',' '], $
   		YMINOR=4, $
   		YCHARSIZE=0.9, $
   		YTITLE='permil'

	XYOUTS,/NORMAL,0.25,0.25,'difference NH-SH'
	XYOUTS,/NORMAL,0.13,0.27,'D'

	;************2H**********************
	PLOT,t,Del2n, $
   		/NOERASE, $
   		POSITION=[0.57,0.30,0.94,0.5], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		thick=1.5, $
   		XTICKNAME=replicate(' ',6), $
		;YRANGE=[-60,-45], $
   		YRANGE=[MIN([Del2n,Del2s])-1,MAX([Del2n,Del2s])+1], $
   		YTICKS=7, $
		;YTICKNAME=[' ','-49','-48','-47','-46','-45','-44',' '], $
   		YMINOR=1, $
   		YCHARSIZE=0.9, $
   		ystyle=1, $
   		YTITLE='permil (SMOW)'
	OPLOT,t,Del2s, $
   		color=pen(2), $
   		thick=1.5

	tarr=['NH','SH']
	larr=[0,2]
	;carr=[pen(1),pen(2)]
	carr=[pen(1),pen(2)]
	CCG_LLEGEND,x=0.20,y=0.36,tarr=tarr,larr=larr,carr=carr
	XYOUTS,/NORMAL,0.15,0.38,'isotopic ratios:'
	;CCG_LLEGEND,x=0.29,y=0.36,tarr=tarr,larr=larr,carr=carr
	;XYOUTS,/NORMAL,0.25,0.38,'isotopic ratios:'
	XYOUTS,/NORMAL,0.13,0.47,'C'

	PLOT,t,Del2n-Del2s, $
   		/NOERASE, $
   		POSITION=[0.57,0.10,0.94,0.3], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		THICK=1.5, $
   		XTITLE='year', $
   		XCHARSIZE=1.0, $
   		YRANGE=[-10.,10.0], $
   		YTICKS=5, $
   		;YTICKNAME=['-1.0','-0.8','-0.6','-0.4','-0.2',' '], $
   		YMINOR=4, $
   		YCHARSIZE=0.9, $
   		YTITLE='permil'

	XYOUTS,/NORMAL,0.25,0.25,'difference NH-SH'
	XYOUTS,/NORMAL,0.13,0.27,'D'
	
ENDIF ELSE BEGIN
	
	;***********************************************************************
	;PLOT ISOTOPIC FORCING FUNCTIONS

   	;Tease apart factors making up the first quasi-steady state term
	PLOT,t,Delq-Deln, $ 
   		POSITION=[0.1,0.70,0.47,0.9], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		XTICKNAME=replicate(' ',6), $
   		YTICKNAME=[' ','-5','-4','-3','-2','-1','0'], $
   		YMINOR=2, $
   		YCHARSIZE=0.9, $
   		YTITLE='permil'
   
	XYOUTS,/NORMAL,0.13,0.87,'A'

	PLOT,t,Q/(2*Xn)/exch, $ 
   		/NOERASE, $
   		POSITION=[0.1,0.50,0.47,0.7], $
   		YRANGE=[0,1], $
   		YTICKS=10, $
   		YMINOR=1, $
   		YTITLE='fraction', $
   		YCHARSIZE=0.9, $
   		YTICKNAME=[' ',' ','0.2',' ','0.4',' ','0.6',' ','0.8',' ',' '], $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		XTICKNAME=replicate(' ',6)
   
	XYOUTS,/NORMAL,0.13,0.67,'B'

	PLOT,t,Q*(Delq-Deln)/(2*Xn)/exch, $ 
		/NOERASE, $
		POSITION=[0.1,0.30,0.47,0.5], $
		YRANGE=[-1.0,0], $
		YTICKS=10, $
		YMINOR=1, $
		YTICKNAME=[' ',' ','-0.8',' ','-0.6',' ','-0.4',' ','-0.2',' ',' '], $
		YCHARSIZE=0.9, $
		YTITLE='permil', $
		XRANGE=xran, $
		;XSTYLE=1, $
		XTICKNAME=replicate(' ',6)
   
	XYOUTS,/NORMAL,0.13,0.47,'C'
 
	PLOT,t,Deln-Dels-Q*(Delq-Deln)/(2*Xn)/exch, $ 
   		/NOERASE, $
   		POSITION=[0.1,0.10,0.47,0.3], $
   		YRANGE=[-0.2,0.8], $
   		YTICKS=10, $
   		YMINOR=1, $
   		YTICKNAME=['-0.2',' ','0',' ','0.2',' ','0.4',' ','0.6',' ',' '], $
   		YCHARSIZE=0.9, $
   		YTITLE='permil', $
   		XRANGE=xran, $
		;XSTYLE=1, $
   		XCHARSIZE=1.0, $
   		XTITLE='year'
   
	XYOUTS,/NORMAL,0.13,0.27,'D'
 
ENDELSE

;**************************WRITE DATA TO FILE************************

;CCG_FWRITE, file='~john/idl/ch4model/data/trial',nc=13,t,Xn,Xs,deln,dels, $
;	bact,ff,bmb,lfills,oh,strat,qa1,qnp

;*************************************

CCG_CLOSEDEV, dev=dev
stop
END