##EMULATE NOAA RADIATIVE FORCING GRAPHIC FROM GHG CONCENTRATIONS


##RADIATIVE FORCING CALCULATION FUNCTION

  #the function below is under construction for various methods
   #only the NOAA method is used in the present example
####

 force=function(x,method="hansen88",scenario="B",base=c(278,700,270,0,0,0)) {
   x=unlist(x)	 
   name0=c("CO2","CH4","N2O","CFC11","CFC12","otg")
    force=rep(NA,7);
   names(base)=name0

  if(method=="hansen88") { 	
	#g=function (M,N) (0.394*M^0.66+0.16*M*exp(-1.6*M) )/(1+0.169*M^0.62) +
        #          1.556*log(1 + N^0.77* (109.8+3.5*N)/(100+0.14*N^2) ) -0.14*log(1+.636*(M*N)^.75 + 0.007*M*(M*N)^1.52) 
			#this formula won't work 
	g=function (M,N) 0.47*log(1 + 2.01E-5 *(M*N)^0.75 + 5.31E-15*M*(M*N)^1.52)
	h=function(x) {h=  log(1+ 1.2* x +0.005*x^2+1.4E-6*x^3);h}
	force[1] = 3.35*(h(x["CO2"] )-h(base["CO2"])) #3.35 from IPCC TAR simplified expressions
	force[2] = 3.35* ( 0.01074*(x["CH4"]^.5 - base["CH4"]^.5) - (g(x["CH4"],N=base["N2O"]) - g(base["CH4"],base["N2O"])))
	force[3] =  3.35*( 0.0418*(x["N2O"]^.5 - base["N2O"]^.5) - (g(M=base["CH4"],x["N2O"]) - g(base["CH4"],base["N2O"])))
	force[4] = 3.35* 0.066*(x["CFC11"]-base["CFC11"])  #0.2211
	force[5] =3.35* 0.084*(x["CFC12"] -base["CFC12"])  # 0.2814
	}

		#TAR has a donstant of 3.35 in its formula for CO2; IPCC 1990 page 52 provides figure of 3.35 (Lacis, pers comm to 
	#convert forcing as a temperature change to forcing as a change in net flux at the tropopause after allowing for stratopsheric temperatyere change.
	#These expessions should be considered as global mean forcings; they implicitly include the radiative effects of global mean cloud cover". 
	#this applies to all terms I guess

  if(method=="ipcc90") { 	
	g=function (M,N) 0.47*log(1 + 2.01E-5 *(M*N)^0.75 + 5.31E-15*M*(M*N)^1.52)
	force[1] =6.3*log(x["CO2"]/base["CO2"])
 	force[2] = 0.036*(x["CH4"]^.5 - base["CH4"]^.5) - (g(x["CH4"],N=base["N2O"]) - g(base["CH4"],base["N2O"]))
	force[3] =  0.14*(x["N2O"]^.5 - base["N2O"]^.5) - (g(M=base["CH4"],x["N2O"]) - g(base["CH4"],base["N2O"]))
	force[4] =0.22*(x["CFC11"]-base["CFC11"])
	force[5] =0.28*(x["CFC12"] -base["CFC12"])
	force[6] = 0.25*(x["otg"]-base["otg"])
	#strat=0.011*(x["CH4"]^.5-base["CH4"]^,5)
	#ozone=0.02*(x["ozone"]-base["ozone"])
	}

   if(method=="noaa"){
  	g=function (M,N) 0.47*log(1 + 2.01E-5 *(M*N)^0.75 + 5.31E-15*M*(M*N)^1.52)
	force[1] =5.35*log(x["CO2"]/base["CO2"])
 	force[2] = 0.036*(x["CH4"]^.5 - base["CH4"]^.5) - (g(x["CH4"],N=base["N2O"]) - g(base["CH4"],base["N2O"])) # in ppb
	force[3] =  0.12*(x["N2O"]^.5 - base["N2O"]^.5) - (g(M=base["CH4"],x["N2O"]) - g(base["CH4"],base["N2O"])) #in ppb
	force[4] =0.25*(x["CFC11"]-base["CFC11"])
	force[5] =0.32*(x["CFC12"] -base["CFC12"])
	force[6] = 0.25*(x["otg"]-base["otg"])
	}

 if(method=="wigley"){
   	#g=function (M,N) 0.47*log(1 + 2.01E-5 *(M*N)^0.75 + 5.31E-15*M*(M*N)^1.52)
	force[1] =6.333*log(x["CO2"]/base["CO2"])
 	force[2] = 0.0398*(x["CH4"]^.5 - base["CH4"]^.5) 
	force[3] =  0.105*(x["N2O"]^.5 - base["N2O"]^.5) 
	force[4] =0.27*(x["CFC11"]-base["CFC11"])
	force[5] =0.31*(x["CFC12"] -base["CFC12"])
	force[6] = 0.25*(x["otg"]-base["otg"])
	}

 if(method=="hansen98"){
	base=c(290.8,837,278.2,0,0,0);names(base)=name0  #1880 values
	g=function (M,N) 0.5*log(1 + 2.0E-5 *(M*N)^0.75 )
	force[1] =5.04*( log(x["CO2"] +0.0005*x["CO2"]^2)- log ( base["CO2"] +0.0005*base["CO2"]^2))
	force[2] = 0.04*(x["CH4"]^.5 - base["CH4"]^.5) - (g(x["CH4"],N=base["N2O"]) - g(base["CH4"],base["N2O"]))
	force[3] = 0.14*(x["N2O"]^.5 - base["N2O"]^.5) - (g(M=base["CH4"],x["N2O"]) - g(base["CH4"],base["N2O"]))
	force[4] =0.25*(x["CFC11"]-base["CFC11"])
	force[5] =0.30*(x["CFC12"] -base["CFC12"])
	force[6] = 0
	}

	force[7]=sum(force,na.rm=T)
	force=t(as.matrix(force))
	force
	}
################



##LOAD GHG CONCENTRATION INFO (GISS VERSION)	
	giss=read.table("http://www.climateaudit.org/data/hansen/giss_ghg.2007.dat",sep="\t",header=TRUE)
 		#this table was previously collated from three sources at NASA GISS:
    		#1850-2000 from http://data.giss.nasa.gov/modelforce/ghgases/GHGs.1850-2000.txt
		#2001-2004 from http://data.giss.nasa.gov/modelforce/ghgases/GHGs.Obsv.2001-2004.txt
		#2005-2006 from individual files
	names(giss)[7]="otg"

##CALCULATE RF USING NOAA FORMULAE
   #see webpage http://www.esrl.noaa.gov/gmd/aggi/
   #base is set at 1750 levels in this comparison
 
	  test=ts(array(NA,dim=c(nrow(giss),7)),start=1850)
	  for(j in 1:nrow(giss)){ test[j,]=force(giss[j,2:7],method="noaa") }
	    test[,7]=apply(test[,1:6],1,sum,na.rm=T)
	    dimnames(test)[[2]]=c(name0,"total")

	X=test[,c(6,4,5,3:1,7)] ;Y=X;year=c(time(X));N=length(year)
	X[is.na(X)]=0
	for (j in 1:nrow(Y)) Y[j,]=cumsum(X[j,]) 
 
##PLOT GRAPHIC
	ylim0=c(0,3)	#for (i in 1:ncol(X)){
	  plot(year,Y[,6],type="l",xlim=c(1978,2006),ylim=ylim0,ylab="RF (wm-2)",xaxs="i",axes=FALSE)
	  axis(side=2,las=1,font=2,cex=.7);box()
	  axis(side=1,font=2)
	 i=6; polygon(xy.coords(c(year[1],year,year[N]),c(0,Y[,i],0)),col="cyan")
	 i=5; polygon(xy.coords(c(year[1],year,year[N]),c(0,Y[,i],0)),col="palegreen")
	 i=4; polygon(xy.coords(c(year[1],year,year[N]),c(0,Y[,i],0)),col="magenta")
	 i=3; polygon(xy.coords(c(year[1],year,year[N]),c(0,Y[,i],0)),col="yellow")
	 i=2; polygon(xy.coords(c(year[1],year,year[N]),c(0,Y[,i],0)),col="blue")
 	i=1; polygon(xy.coords(c(year[1],year,year[N]),c(0,Y[,i],0)),col="red")
	legend(1978,3.1,fill=c("cyan","palegreen","magenta","yellow","blue","red"),
		legend=c("CO2","CH4","N2O","CFC12","CFC11","OTG"),cex=.8)
	title(main="NOAA Radiative Forcing",cex=.8)