clamp: npp/99.all.ncl.1@4be95183fbcd (annotated)

forrest@0	1	;*****************************************************
forrest@0	2	; combine scatter, histogram, global and zonal plots
forrest@0	3	; compute all correlation coef and M score
forrest@0	4	; add 1-to-1 line in scatter plots
forrest@0	5	;*****************************************************
forrest@0	6	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
forrest@0	7	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
forrest@0	8	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0	9	;*****************************************************
forrest@0	10	procedure get_bin(RAIN1_1D[]:numeric, NPP1_1D[]:numeric \
forrest@0	11	,RAIN2_1D[]:numeric, NPP2_1D[]:numeric \
forrest@0	12	,xvalues[][]:numeric, yvalues[][]:numeric \
forrest@0	13	,mn_yvalues[][]:numeric, mx_yvalues[][]:numeric \
forrest@0	14	,dx4[1]:numeric \
forrest@0	15	)
forrest@0	16	begin
forrest@0	17	; Calculaee "nice" bins for binning the data in equally spaced ranges.
forrest@0	18	; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
forrest@0	19	; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
forrest@0	20
forrest@0	21	nbins = 15 ; Number of bins to use.
forrest@0	22
forrest@0	23	nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
forrest@0	24	nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
forrest@0	25	range = fspan(nicevals(0),nicevals(1),nvals)
forrest@0	26
forrest@0	27	; Use this range information to grab all the values in a
forrest@0	28	; particular range, and then take an average.
forrest@0	29
forrest@0	30	nr = dimsizes(range)
forrest@0	31	nx = nr-1
forrest@0	32	; print (nx)
forrest@0	33
forrest@0	34	; xvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0	35	xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
forrest@0	36	dx = xvalues(0,1) - xvalues(0,0) ; range width
forrest@0	37	dx4 = dx/4 ; 1/4 of the range
forrest@0	38	xvalues(1,:) = xvalues(0,:) - dx/5.
forrest@0	39	; yvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0	40	; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0	41	; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0	42
forrest@0	43	do nd=0,1
forrest@0	44
forrest@0	45	; See if we are doing model or observational data.
forrest@0	46
forrest@0	47	if(nd.eq.0) then
forrest@0	48	data = RAIN1_1D
forrest@0	49	npp_data = NPP1_1D
forrest@0	50	else
forrest@0	51	data = RAIN2_1D
forrest@0	52	npp_data = NPP2_1D
forrest@0	53	end if
forrest@0	54
forrest@0	55	; Loop through each range and check for values.
forrest@0	56
forrest@0	57	do i=0,nr-2
forrest@0	58	if (i.ne.(nr-2)) then
forrest@0	59	; print("")
forrest@0	60	; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0	61	idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
forrest@0	62	else
forrest@0	63	; print("")
forrest@0	64	; print("In range ["+range(i)+",)")
forrest@0	65	idx = ind(range(i).le.data)
forrest@0	66	end if
forrest@0	67
forrest@0	68	; Calculate average, and get min and max.
forrest@0	69
forrest@0	70	if(.not.any(ismissing(idx))) then
forrest@0	71	yvalues(nd,i) = avg(npp_data(idx))
forrest@0	72	mn_yvalues(nd,i) = min(npp_data(idx))
forrest@0	73	mx_yvalues(nd,i) = max(npp_data(idx))
forrest@0	74	count = dimsizes(idx)
forrest@0	75	else
forrest@0	76	count = 0
forrest@0	77	yvalues(nd,i) = yvalues@_FillValue
forrest@0	78	mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0	79	mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0	80	end if
forrest@0	81
forrest@0	82	; Print out information.
forrest@0	83
forrest@0	84	; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0	85	; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0	86
forrest@0	87
forrest@0	88	; Clean up for next time in loop.
forrest@0	89
forrest@0	90	delete(idx)
forrest@0	91	end do
forrest@0	92	delete(data)
forrest@0	93	delete(npp_data)
forrest@0	94	end do
forrest@0	95	end
forrest@0	96	;****************************************************************************
forrest@0	97	; Main code.
forrest@0	98	;****************************************************************************
forrest@0	99	begin
forrest@0	100
forrest@0	101	plot_type = "ps"
forrest@0	102	plot_type_new = "png"
forrest@0	103
forrest@0	104	;************************************************
forrest@0	105	; read model data
forrest@0	106	;************************************************
forrest@0	107
forrest@0	108	;film = "i01.06cn_1798-2004_ANN_climo.nc"
forrest@0	109	;model_name = "06cn"
forrest@0	110	;model_grid = "T42"
forrest@0	111
forrest@0	112	;film = "i01.06casa_1798-2004_ANN_climo.nc"
forrest@0	113	;model_name = "06casa"
forrest@0	114	;model_grid = "T42"
forrest@0	115
forrest@0	116	film = "i01.10casa_1948-2004_ANN_climo.nc"
forrest@0	117	model_name = "10casa"
forrest@0	118	model_grid = "T42"
forrest@0	119
forrest@0	120	;film = "i01.10cn_1948-2004_ANN_climo.nc"
forrest@0	121	;model_name = "10cn"
forrest@0	122	;model_grid = "T42"
forrest@0	123	;--------------------------------------------------
forrest@0	124	dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0	125	fm = addfile (dirm+film,"r")
forrest@0	126
forrest@0	127	nppmod0 = fm->NPP
forrest@0	128	rainmod0 = fm->RAIN
forrest@0	129	xm = fm->lon
forrest@0	130	ym = fm->lat
forrest@0	131
forrest@0	132	;************************************************
forrest@0	133	; read ob data
forrest@0	134	;************************************************
forrest@0	135
forrest@0	136	;(1) data at 81 sites
forrest@0	137	dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
forrest@0	138	fil81 = "data.81.nc"
forrest@0	139	f81 = addfile (dir81+fil81,"r")
forrest@0	140
forrest@0	141	id81 = f81->SITE_ID
forrest@0	142	npp81 = f81->TNPP_C
forrest@0	143	rain81 = tofloat(f81->PREC_ANN)
forrest@0	144	x81 = f81->LONG_DD
forrest@0	145	y81 = f81->LAT_DD
forrest@0	146
forrest@0	147	id81@long_name = "SITE_ID"
forrest@0	148
forrest@0	149	;change longitude from (-180,180) to (0,360)
forrest@0	150	;for model data interpolation
forrest@0	151
forrest@0	152	do i= 0,dimsizes(x81)-1
forrest@0	153	if (x81(i) .lt. 0.) then
forrest@0	154	x81(i) = x81(i)+ 360.
forrest@0	155	end if
forrest@0	156	end do
forrest@0	157	;print (x81)
forrest@0	158	;-------------------------------------
forrest@0	159	;(2) data at 933 sites
forrest@0	160	dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
forrest@0	161	fil933 = "data.933.nc"
forrest@0	162	f933 = addfile (dir933+fil933,"r")
forrest@0	163
forrest@0	164	id933 = f933->SITE_ID
forrest@0	165	npp933 = f933->TNPP_C
forrest@0	166	rain933 = f933->PREC
forrest@0	167	x933 = f933->LONG_DD
forrest@0	168	y933 = f933->LAT_DD
forrest@0	169
forrest@0	170	id933@long_name = "SITE_ID"
forrest@0	171
forrest@0	172	;change longitude from (-180,180) to (0,360)
forrest@0	173	;for model data interpolation
forrest@0	174
forrest@0	175	do i= 0,dimsizes(x933)-1
forrest@0	176	if (x933(i) .lt. 0.) then
forrest@0	177	x933(i) = x933(i)+ 360.
forrest@0	178	end if
forrest@0	179	end do
forrest@0	180	;print (x933)
forrest@0	181	;----------------------------------------
forrest@0	182	;(3) global data, interpolated into model grid
forrest@0	183	dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
forrest@0	184	filglobe = "Npp_"+model_grid+"_mean.nc"
forrest@0	185	fglobe = addfile (dirglobe+filglobe,"r")
forrest@0	186
forrest@0	187	nppglobe0 = fglobe->NPP
forrest@0	188	nppglobe = nppglobe0
forrest@0	189
forrest@0	190	;***********************************************************************
forrest@0	191	; interpolate model data into ob sites
forrest@0	192	;***********************************************************************
forrest@0	193
forrest@0	194	nppmod = nppmod0(0,:,:)
forrest@0	195	rainmod = rainmod0(0,:,:)
forrest@0	196	delete (nppmod0)
forrest@0	197	delete (rainmod0)
forrest@0	198
forrest@0	199	nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
forrest@0	200
forrest@0	201	nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
forrest@0	202
forrest@0	203	rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
forrest@0	204
forrest@0	205	rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
forrest@0	206
forrest@0	207	;**********************************************************
forrest@0	208	; unified units
forrest@0	209	;**********************************************************
forrest@0	210	; Units for these variables are:
forrest@0	211	;
forrest@0	212	; rain81 : mm/year
forrest@0	213	; rainmod : mm/s
forrest@0	214	; npp81 : g C/m^2/year
forrest@0	215	; nppmod81: g C/m^2/s
forrest@0	216	; nppglobe: g C/m^2/year
forrest@0	217	;
forrest@0	218	; We want to convert these to "m/year" and "g C/m^2/year".
forrest@0	219
forrest@0	220	nsec_per_year = 606024*365
forrest@0	221
forrest@0	222	rain81 = rain81 / 1000.
forrest@0	223	rainmod81 = (rainmod81/ 1000.) * nsec_per_year
forrest@0	224	nppmod81 = nppmod81 * nsec_per_year
forrest@0	225
forrest@0	226	rain933 = rain933 / 1000.
forrest@0	227	rainmod933 = (rainmod933/ 1000.) * nsec_per_year
forrest@0	228	nppmod933 = nppmod933 * nsec_per_year
forrest@0	229
forrest@0	230	nppmod = nppmod * nsec_per_year
forrest@0	231
forrest@0	232	npp81@units = "gC/m^2/yr"
forrest@0	233	nppmod81@units = "gC/m^2/yr"
forrest@0	234	npp933@units = "gC/m^2/yr"
forrest@0	235	nppmod933@units = "gC/m^2/yr"
forrest@0	236	nppmod@units = "gC/m^2/yr"
forrest@0	237	nppglobe@units = "gC/m^2/yr"
forrest@0	238	rain81@units = "m/yr"
forrest@0	239	rainmod81@units = "m/yr"
forrest@0	240	rain933@units = "m/yr"
forrest@0	241	rainmod933@units = "m/yr"
forrest@0	242
forrest@0	243	npp81@long_name = "NPP (gC/m2/year)"
forrest@0	244	npp933@long_name = "NPP (gC/m2/year)"
forrest@0	245	nppmod81@long_name = "NPP (gC/m2/year)"
forrest@0	246	nppmod933@long_name = "NPP (gC/m2/year)"
forrest@0	247	nppmod@long_name = "NPP (gC/m2/year)"
forrest@0	248	nppglobe@long_name = "NPP (gC/m2/year)"
forrest@0	249	rain81@long_name = "PREC (m/year)"
forrest@0	250	rain933@long_name = "PREC (m/year)"
forrest@0	251	rainmod81@long_name = "PREC (m/year)"
forrest@0	252	rainmod933@long_name = "PREC (m/year)"
forrest@0	253
forrest@0	254	;*******************************************************************
forrest@0	255	;(A)-1 table of site81 (1) (the first 41 sites)
forrest@0	256	;*******************************************************************
forrest@0	257
forrest@0	258	table_length = 0.95
forrest@0	259
forrest@0	260	; table header name
forrest@0	261	table_header_name = "Site ID"
forrest@0	262
forrest@0	263	; column (not including header column)
forrest@0	264	col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
forrest@0	265	ncol = dimsizes(col_header)
forrest@0	266
forrest@0	267	; row (not including header row)
forrest@0	268	nrow = 41
forrest@0	269	row_header = new ((/nrow/),string )
forrest@0	270	row_header(0:nrow-1) = id81(0:nrow-1)
forrest@0	271
forrest@0	272	; Table header
forrest@0	273	ncr1 = (/1,1/) ; 1 row, 1 column
forrest@0	274	x1 = (/0.005,0.15/) ; Start and end X
forrest@0	275	y1 = (/0.900,0.95/) ; Start and end Y
forrest@0	276	text1 = table_header_name
forrest@0	277	res1 = True
forrest@0	278	res1@txFontHeightF = 0.015
forrest@0	279	res1@gsFillColor = "CornFlowerBlue"
forrest@0	280
forrest@0	281	; Column header (equally space in x2)
forrest@0	282	ncr2 = (/1,ncol/) ; 1 rows, ncol columns
forrest@0	283	x2 = (/x1(1),0.995/) ; start from end of x1
forrest@0	284	y2 = y1 ; same as y1
forrest@0	285	text2 = col_header
forrest@0	286	res2 = True
forrest@0	287	res2@txFontHeightF = 0.015
forrest@0	288	res2@gsFillColor = "Gray"
forrest@0	289
forrest@0	290	; Row header (equally space in y2)
forrest@0	291	ncr3 = (/nrow,1/) ; nrow rows, 1 columns
forrest@0	292	x3 = x1 ; same as x1
forrest@0	293	y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
forrest@0	294	text3 = row_header
forrest@0	295	res3 = True
forrest@0	296	res3@txFontHeightF = 0.010
forrest@0	297	res3@gsFillColor = "Gray"
forrest@0	298
forrest@0	299	; Main table body
forrest@0	300	ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
forrest@0	301	x4 = x2 ; Start and end x
forrest@0	302	y4 = y3 ; Start and end Y
forrest@0	303	text4 = new((/nrow,ncol/),string)
forrest@0	304
forrest@0	305	color_fill4 = new((/nrow,ncol/),string)
forrest@0	306	color_fill4 = "white"
forrest@0	307	; color_fill4(:,ncol-1) = "grey"
forrest@0	308	; color_fill4(nrow-1,:) = "green"
forrest@0	309
forrest@0	310	res4 = True ; Set up resource list
forrest@0	311	; res4@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	312	res4@txFontHeightF = 0.015
forrest@0	313	res4@gsFillColor = color_fill4
forrest@0	314
forrest@0	315	delete (color_fill4)
forrest@0	316	;-------------------------------------------------------------------
forrest@0	317	; for table value
forrest@0	318
forrest@0	319	do n = 0,nrow-1
forrest@0	320	text4(n,0) = sprintf("%5.2f", y81(n))
forrest@0	321	text4(n,1) = sprintf("%5.2f", x81(n))
forrest@0	322	text4(n,2) = sprintf("%5.2f", npp81(n))
forrest@0	323	text4(n,3) = sprintf("%5.2f", rain81(n))
forrest@0	324	end do
forrest@0	325	;-------------------------------------------------------------------
forrest@0	326
forrest@0	327	plot_name = "table_site81_ob1"
forrest@0	328
forrest@0	329	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	330	;------------------------------------------
forrest@0	331	; for table title
forrest@0	332
forrest@0	333	gRes = True
forrest@0	334	gRes@txFontHeightF = 0.02
forrest@0	335	; gRes@txAngleF = 90
forrest@0	336
forrest@0	337	title_text = "Observation at 81 Sites (1)"
forrest@0	338
forrest@0	339	gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
forrest@0	340	;------------------------------------------
forrest@0	341
forrest@0	342	gsn_table(wks,ncr1,x1,y1,text1,res1)
forrest@0	343	gsn_table(wks,ncr2,x2,y2,text2,res2)
forrest@0	344	gsn_table(wks,ncr3,x3,y3,text3,res3)
forrest@0	345	gsn_table(wks,ncr4,x4,y4,text4,res4)
forrest@0	346
forrest@0	347	frame(wks)
forrest@0	348	clear (wks)
forrest@0	349
forrest@0	350	delete (row_header)
forrest@0	351	delete (res3)
forrest@0	352	delete (ncr3)
forrest@0	353	delete (text3)
forrest@0	354	delete (res4)
forrest@0	355	delete (ncr4)
forrest@0	356	delete (text4)
forrest@0	357
forrest@0	358	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	359	system("rm "+plot_name+"."+plot_type)
forrest@0	360	system("rm "+plot_name+"-1."+plot_type_new)
forrest@0	361	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	362
forrest@0	363	;*******************************************************************
forrest@0	364	;(A)-2 table of site81 (2) (the last 40 sites)
forrest@0	365	;*******************************************************************
forrest@0	366
forrest@0	367	table_length = 0.95
forrest@0	368
forrest@0	369	; table header name
forrest@0	370	table_header_name = "Site ID"
forrest@0	371
forrest@0	372	; column (not including header column)
forrest@0	373	col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
forrest@0	374	ncol = dimsizes(col_header)
forrest@0	375
forrest@0	376	; row (not including header row)
forrest@0	377	nrow = 40
forrest@0	378	row_header = new ((/nrow/),string )
forrest@0	379	row_header(0:nrow-1) = id81(41:41+nrow-1)
forrest@0	380
forrest@0	381	; Table header
forrest@0	382	ncr1 = (/1,1/) ; 1 row, 1 column
forrest@0	383	x1 = (/0.005,0.15/) ; Start and end X
forrest@0	384	y1 = (/0.900,0.95/) ; Start and end Y
forrest@0	385	text1 = table_header_name
forrest@0	386	res1 = True
forrest@0	387	res1@txFontHeightF = 0.015
forrest@0	388	res1@gsFillColor = "CornFlowerBlue"
forrest@0	389
forrest@0	390	; Column header (equally space in x2)
forrest@0	391	ncr2 = (/1,ncol/) ; 1 rows, ncol columns
forrest@0	392	x2 = (/x1(1),0.995/) ; start from end of x1
forrest@0	393	y2 = y1 ; same as y1
forrest@0	394	text2 = col_header
forrest@0	395	res2 = True
forrest@0	396	res2@txFontHeightF = 0.015
forrest@0	397	res2@gsFillColor = "Gray"
forrest@0	398
forrest@0	399	; Row header (equally space in y2)
forrest@0	400	ncr3 = (/nrow,1/) ; nrow rows, 1 columns
forrest@0	401	x3 = x1 ; same as x1
forrest@0	402	y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
forrest@0	403	text3 = row_header
forrest@0	404	res3 = True
forrest@0	405	res3@txFontHeightF = 0.010
forrest@0	406	res3@gsFillColor = "Gray"
forrest@0	407
forrest@0	408	; Main table body
forrest@0	409	ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
forrest@0	410	x4 = x2 ; Start and end x
forrest@0	411	y4 = y3 ; Start and end Y
forrest@0	412	text4 = new((/nrow,ncol/),string)
forrest@0	413
forrest@0	414	color_fill4 = new((/nrow,ncol/),string)
forrest@0	415	color_fill4 = "white"
forrest@0	416	; color_fill4(:,ncol-1) = "grey"
forrest@0	417	; color_fill4(nrow-1,:) = "green"
forrest@0	418
forrest@0	419	res4 = True ; Set up resource list
forrest@0	420	; res4@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	421	res4@txFontHeightF = 0.015
forrest@0	422	res4@gsFillColor = color_fill4
forrest@0	423
forrest@0	424	delete (color_fill4)
forrest@0	425	;-------------------------------------------------------------------
forrest@0	426	; for table value
forrest@0	427
forrest@0	428	do n = 0,nrow-1
forrest@0	429	text4(n,0) = sprintf("%5.2f", y81(n+41))
forrest@0	430	text4(n,1) = sprintf("%5.2f", x81(n+41))
forrest@0	431	text4(n,2) = sprintf("%5.2f", npp81(n+41))
forrest@0	432	text4(n,3) = sprintf("%5.2f", rain81(n+41))
forrest@0	433	end do
forrest@0	434	;-------------------------------------------------------------------
forrest@0	435
forrest@0	436	plot_name = "table_site81_ob2"
forrest@0	437
forrest@0	438	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	439	;------------------------------------------
forrest@0	440	; for table title
forrest@0	441
forrest@0	442	gRes = True
forrest@0	443	gRes@txFontHeightF = 0.02
forrest@0	444	; gRes@txAngleF = 90
forrest@0	445
forrest@0	446	title_text = "Observation at 81 Sites (2)"
forrest@0	447
forrest@0	448	gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
forrest@0	449	;------------------------------------------
forrest@0	450
forrest@0	451	gsn_table(wks,ncr1,x1,y1,text1,res1)
forrest@0	452	gsn_table(wks,ncr2,x2,y2,text2,res2)
forrest@0	453	gsn_table(wks,ncr3,x3,y3,text3,res3)
forrest@0	454	gsn_table(wks,ncr4,x4,y4,text4,res4)
forrest@0	455
forrest@0	456	frame(wks)
forrest@0	457	clear (wks)
forrest@0	458
forrest@0	459	delete (row_header)
forrest@0	460	delete (res3)
forrest@0	461	delete (ncr3)
forrest@0	462	delete (text3)
forrest@0	463	delete (res4)
forrest@0	464	delete (ncr4)
forrest@0	465	delete (text4)
forrest@0	466
forrest@0	467	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	468	system("rm "+plot_name+"."+plot_type)
forrest@0	469	system("rm "+plot_name+"-1."+plot_type_new)
forrest@0	470	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	471
forrest@0	472	;**************************************************************************
forrest@0	473	;(A)-3 scatter plot, ob 933
forrest@0	474	;**************************************************************************
forrest@0	475
forrest@0	476	plot_name = "scatter_ob_933"
forrest@0	477	title = "Observed 933 sites"
forrest@0	478
forrest@0	479	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	480	res = True ; plot mods desired
forrest@0	481	res@tiMainString = title ; add title
forrest@0	482	res@xyMarkLineModes = "Markers" ; choose which have markers
forrest@0	483	res@xyMarkers = 16 ; choose type of marker
forrest@0	484	res@xyMarkerColor = "red" ; Marker color
forrest@0	485	res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
forrest@0	486	res@tmLabelAutoStride = True ; nice tick mark labels
forrest@0	487
forrest@0	488	plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
forrest@0	489	frame(wks)
forrest@0	490	clear (wks)
forrest@0	491
forrest@0	492	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	493	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	494	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	495	system("rm "+plot_name+"."+plot_type)
forrest@0	496
forrest@0	497	;**************************************************************************
forrest@0	498	;(A)-4 scatter plot, model 933
forrest@0	499	;**************************************************************************
forrest@0	500
forrest@0	501	plot_name = "scatter_model_933"
forrest@0	502	title = "Model "+ model_name +" 933 sites"
forrest@0	503
forrest@0	504	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	505	res = True ; plot mods desired
forrest@0	506	res@tiMainString = title ; add title
forrest@0	507	res@xyMarkLineModes = "Markers" ; choose which have markers
forrest@0	508	res@xyMarkers = 16 ; choose type of marker
forrest@0	509	res@xyMarkerColor = "red" ; Marker color
forrest@0	510	res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
forrest@0	511	res@tmLabelAutoStride = True ; nice tick mark labels
forrest@0	512
forrest@0	513	plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
forrest@0	514	frame(wks)
forrest@0	515	clear (wks)
forrest@0	516
forrest@0	517	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	518	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	519	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	520	system("rm "+plot_name+"."+plot_type)
forrest@0	521
forrest@0	522	;***************************************************************************
forrest@0	523	;(A)-5 scatter plot, model vs ob 81
forrest@0	524	;***************************************************************************
forrest@0	525
forrest@0	526	plot_name = "scatter_model_vs_ob_81"
forrest@0	527	title = model_name +" vs ob 81 sites"
forrest@0	528
forrest@0	529	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	530	res = True ; plot mods desired
forrest@0	531	res@tiMainString = title ; add title
forrest@0	532	res@xyMarkLineModes = "Markers" ; choose which have markers
forrest@0	533	res@xyMarkers = 16 ; choose type of marker
forrest@0	534	res@xyMarkerColor = "red" ; Marker color
forrest@0	535	res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
forrest@0	536	res@tmLabelAutoStride = True ; nice tick mark labels
forrest@0	537
forrest@0	538	res@gsnDraw = False
forrest@0	539	res@gsnFrame = False ; don't advance frame yet
forrest@0	540	;-------------------------------
forrest@0	541	;compute correlation coef. and M
forrest@0	542	ccr81 = esccr(nppmod81,npp81,0)
forrest@0	543	;print (ccr81)
forrest@0	544
forrest@0	545	;new eq
forrest@0	546	bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
forrest@0	547	M81 = (1. - (bias/dimsizes(y81)))*5.
forrest@0	548	print (M81)
forrest@0	549
forrest@0	550	tRes = True
forrest@0	551	tRes@txFontHeightF = 0.025
forrest@0	552
forrest@0	553	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
forrest@0	554
forrest@0	555	gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
forrest@0	556	;--------------------------------
forrest@0	557	plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
forrest@0	558	;-------------------------------
forrest@0	559	; add polyline
forrest@0	560
forrest@0	561	dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
forrest@0	562	;-------------------------------
forrest@0	563	draw (plot)
forrest@0	564	frame(wks)
forrest@0	565	clear (wks)
forrest@0	566	delete (dum)
forrest@0	567
forrest@0	568	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	569	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	570	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	571	system("rm "+plot_name+"."+plot_type)
forrest@0	572
forrest@0	573	;***************************************************************************
forrest@0	574	;(A)-6 scatter plot, model vs ob 933
forrest@0	575	;***************************************************************************
forrest@0	576
forrest@0	577	plot_name = "scatter_model_vs_ob_933"
forrest@0	578	title = model_name +" vs ob 933 sites"
forrest@0	579
forrest@0	580	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	581	res = True ; plot mods desired
forrest@0	582	res@tiMainString = title ; add title
forrest@0	583	res@xyMarkLineModes = "Markers" ; choose which have markers
forrest@0	584	res@xyMarkers = 16 ; choose type of marker
forrest@0	585	res@xyMarkerColor = "red" ; Marker color
forrest@0	586	res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
forrest@0	587	res@tmLabelAutoStride = True ; nice tick mark labels
forrest@0	588
forrest@0	589	res@gsnDraw = False
forrest@0	590	res@gsnFrame = False ; don't advance frame yet
forrest@0	591	;-------------------------------
forrest@0	592	;compute correlation coef. and M
forrest@0	593	ccr933 = esccr(nppmod933,npp933,0)
forrest@0	594	;print (ccr933)
forrest@0	595
forrest@0	596	;new eq
forrest@0	597	bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
forrest@0	598	M933 = (1. - (bias/dimsizes(y933)))*5.
forrest@0	599	print (M933)
forrest@0	600
forrest@0	601	tRes = True
forrest@0	602	tRes@txFontHeightF = 0.025
forrest@0	603
forrest@0	604	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
forrest@0	605
forrest@0	606	gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
forrest@0	607	;--------------------------------
forrest@0	608	plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
forrest@0	609	;-------------------------------
forrest@0	610	; add polyline
forrest@0	611
forrest@0	612	dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
forrest@0	613	;-------------------------------
forrest@0	614	draw (plot)
forrest@0	615	frame(wks)
forrest@0	616	clear (wks)
forrest@0	617	delete (dum)
forrest@0	618
forrest@0	619	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	620	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	621	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	622	system("rm "+plot_name+"."+plot_type)
forrest@0	623
forrest@0	624	;*******************************************************************
forrest@0	625	;(A)-7 table of site81 (1) (the first 41 sites), model vs ob
forrest@0	626	;*******************************************************************
forrest@0	627
forrest@0	628	table_length = 0.95
forrest@0	629
forrest@0	630	; table header name
forrest@0	631	table_header_name = "Site ID"
forrest@0	632
forrest@0	633	; row (not including header row)
forrest@0	634	nrow = 41
forrest@0	635	row_header = new ((/nrow/),string )
forrest@0	636	row_header(0:nrow-1) = id81(0:nrow-1)
forrest@0	637
forrest@0	638	; Table header
forrest@0	639	ncr1 = (/1,1/) ; 1 row, 1 column
forrest@0	640	x1 = (/0.005,0.15/) ; Start and end X
forrest@0	641	y1 = (/0.85 ,0.95/) ; Start and end Y
forrest@0	642	text1 = table_header_name
forrest@0	643	res1 = True
forrest@0	644	res1@txFontHeightF = 0.015
forrest@0	645	res1@gsFillColor = "CornFlowerBlue"
forrest@0	646
forrest@0	647	; Column header1 (equally space in x2)
forrest@0	648
forrest@0	649	delete (col_header)
forrest@0	650	delete (ncr2)
forrest@0	651	delete (text2)
forrest@0	652
forrest@0	653	col_header = (/"Latitude","Longitude"/)
forrest@0	654	ncol = dimsizes(col_header)
forrest@0	655
forrest@0	656	ncr2 = (/1,ncol/) ; 1 rows, ncol columns
forrest@0	657	x2 = (/x1(1),0.35/) ; start from end of x1
forrest@0	658	y2 = y1 ; same as y1
forrest@0	659	text2 = col_header
forrest@0	660	res2 = True
forrest@0	661	res2@txFontHeightF = 0.015
forrest@0	662	res2@gsFillColor = "Gray"
forrest@0	663
forrest@0	664	; Column header2 (equally space in x2)
forrest@0	665
forrest@0	666	col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
forrest@0	667	col_header2 = (/"ob",model_name \
forrest@0	668	,"ob",model_name \
forrest@0	669	/)
forrest@0	670
forrest@0	671	ncol1 = dimsizes(col_header1)
forrest@0	672	ncol2 = dimsizes(col_header2)
forrest@0	673
forrest@0	674	ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
forrest@0	675	x21 = (/x2(1),0.995/) ; start from end of x2
forrest@0	676	yhalf = (y1(0)+y1(1))*0.5
forrest@0	677	y21 = (/yhalf,y1(1)/) ; top half of y1
forrest@0	678	text21 = col_header1
forrest@0	679	res21 = True
forrest@0	680	res21@txFontHeightF = 0.015
forrest@0	681	res21@gsFillColor = "Gray"
forrest@0	682
forrest@0	683	ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
forrest@0	684	x22 = x21 ; start from end of x1
forrest@0	685	y22 = (/y1(0),yhalf/) ; bottom half of y1
forrest@0	686	text22 = col_header2
forrest@0	687	res22 = True
forrest@0	688	res22@txFontHeightF = 0.015
forrest@0	689	res22@gsFillColor = "Gray"
forrest@0	690
forrest@0	691	; Row header (equally space in y2)
forrest@0	692	ncr3 = (/nrow,1/) ; nrow rows, 1 columns
forrest@0	693	x3 = x1 ; same as x1
forrest@0	694	y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
forrest@0	695	text3 = row_header
forrest@0	696	res3 = True
forrest@0	697	res3@txFontHeightF = 0.010
forrest@0	698	res3@gsFillColor = "Gray"
forrest@0	699
forrest@0	700	; Main table body-1
forrest@0	701	ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
forrest@0	702	x4 = x2 ; Start and end x
forrest@0	703	y4 = y3 ; Start and end Y
forrest@0	704	text4 = new((/nrow,ncol/),string)
forrest@0	705
forrest@0	706	color_fill4 = new((/nrow,ncol/),string)
forrest@0	707	color_fill4 = "white"
forrest@0	708	; color_fill4(:,ncol-1) = "grey"
forrest@0	709	; color_fill4(nrow-1,:) = "green"
forrest@0	710
forrest@0	711	res4 = True ; Set up resource list
forrest@0	712	; res4@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	713	res4@txFontHeightF = 0.015
forrest@0	714	res4@gsFillColor = color_fill4
forrest@0	715
forrest@0	716	delete (color_fill4)
forrest@0	717
forrest@0	718	; Main table body-2
forrest@0	719	ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
forrest@0	720	x42 = x21 ; Start and end x
forrest@0	721	y42 = y3 ; Start and end Y
forrest@0	722	text42 = new((/nrow,ncol2/),string)
forrest@0	723
forrest@0	724	color_fill42 = new((/nrow,ncol2/),string)
forrest@0	725	color_fill42 = "white"
forrest@0	726	; color_fill42(:,ncol-1) = "grey"
forrest@0	727	; color_fill42(nrow-1,:) = "green"
forrest@0	728
forrest@0	729	res42 = True ; Set up resource list
forrest@0	730	; res42@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	731	res42@txFontHeightF = 0.015
forrest@0	732	res42@gsFillColor = color_fill42
forrest@0	733
forrest@0	734	delete (color_fill42)
forrest@0	735	;-------------------------------------------------------------------
forrest@0	736	; for table value
forrest@0	737
forrest@0	738	do n = 0,nrow-1
forrest@0	739	text4(n,0) = sprintf("%5.2f", y81(n))
forrest@0	740	text4(n,1) = sprintf("%5.2f", x81(n))
forrest@0	741	end do
forrest@0	742
forrest@0	743	do n = 0,nrow-1
forrest@0	744	text42(n,0) = sprintf("%5.2f", npp81(n))
forrest@0	745	text42(n,1) = sprintf("%5.2f", nppmod81(n))
forrest@0	746	text42(n,2) = sprintf("%5.2f", rain81(n))
forrest@0	747	text42(n,3) = sprintf("%5.2f", rainmod81(n))
forrest@0	748	end do
forrest@0	749	;---------------------------------------------------------------------------
forrest@0	750
forrest@0	751	plot_name = "table_site81_model_vs_ob1"
forrest@0	752
forrest@0	753	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	754	;------------------------------------------
forrest@0	755	; for table title
forrest@0	756
forrest@0	757	gRes = True
forrest@0	758	gRes@txFontHeightF = 0.02
forrest@0	759	; gRes@txAngleF = 90
forrest@0	760
forrest@0	761	title_text = "Model vs Observation at 81 Sites (1)"
forrest@0	762
forrest@0	763	gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
forrest@0	764	;------------------------------------------
forrest@0	765
forrest@0	766	gsn_table(wks,ncr1,x1,y1,text1,res1)
forrest@0	767	gsn_table(wks,ncr2,x2,y2,text2,res2)
forrest@0	768	gsn_table(wks,ncr21,x21,y21,text21,res21)
forrest@0	769	gsn_table(wks,ncr22,x22,y22,text22,res22)
forrest@0	770	gsn_table(wks,ncr3,x3,y3,text3,res3)
forrest@0	771	gsn_table(wks,ncr4,x4,y4,text4,res4)
forrest@0	772	gsn_table(wks,ncr42,x42,y42,text42,res42)
forrest@0	773
forrest@0	774	frame(wks)
forrest@0	775	clear (wks)
forrest@0	776
forrest@0	777	delete (col_header)
forrest@0	778	delete (row_header)
forrest@0	779	delete (res1)
forrest@0	780	delete (ncr1)
forrest@0	781	delete (text1)
forrest@0	782	delete (res2)
forrest@0	783	delete (ncr2)
forrest@0	784	delete (text2)
forrest@0	785	delete (res21)
forrest@0	786	delete (ncr21)
forrest@0	787	delete (text21)
forrest@0	788	delete (res22)
forrest@0	789	delete (ncr22)
forrest@0	790	delete (text22)
forrest@0	791	delete (res3)
forrest@0	792	delete (ncr3)
forrest@0	793	delete (text3)
forrest@0	794	delete (res4)
forrest@0	795	delete (ncr4)
forrest@0	796	delete (text4)
forrest@0	797	delete (res42)
forrest@0	798	delete (ncr42)
forrest@0	799	delete (text42)
forrest@0	800
forrest@0	801	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	802	system("rm "+plot_name+"."+plot_type)
forrest@0	803	system("rm "+plot_name+"-1."+plot_type_new)
forrest@0	804	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	805
forrest@0	806	;*******************************************************************
forrest@0	807	;(A)-8 table of site81 (2) (the last 40 sites), model vs ob
forrest@0	808	;*******************************************************************
forrest@0	809
forrest@0	810	table_length = 0.95
forrest@0	811
forrest@0	812	; table header name
forrest@0	813	table_header_name = "Site ID"
forrest@0	814
forrest@0	815	; row (not including header row)
forrest@0	816	nrow = 40
forrest@0	817	row_header = new ((/nrow/),string )
forrest@0	818	row_header(0:nrow-1) = id81(41:41+nrow-1)
forrest@0	819
forrest@0	820	; Table header
forrest@0	821	ncr1 = (/1,1/) ; 1 row, 1 column
forrest@0	822	x1 = (/0.005,0.15/) ; Start and end X
forrest@0	823	y1 = (/0.85 ,0.95/) ; Start and end Y
forrest@0	824	text1 = table_header_name
forrest@0	825	res1 = True
forrest@0	826	res1@txFontHeightF = 0.015
forrest@0	827	res1@gsFillColor = "CornFlowerBlue"
forrest@0	828
forrest@0	829	; Column header1 (equally space in x2)
forrest@0	830
forrest@0	831	col_header = (/"Latitude","Longitude"/)
forrest@0	832	ncol = dimsizes(col_header)
forrest@0	833
forrest@0	834	ncr2 = (/1,ncol/) ; 1 rows, ncol columns
forrest@0	835	x2 = (/x1(1),0.35/) ; start from end of x1
forrest@0	836	y2 = y1 ; same as y1
forrest@0	837	text2 = col_header
forrest@0	838	res2 = True
forrest@0	839	res2@txFontHeightF = 0.015
forrest@0	840	res2@gsFillColor = "Gray"
forrest@0	841
forrest@0	842	; Column header2 (equally space in x2)
forrest@0	843
forrest@0	844	; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
forrest@0	845	col_header2 = (/"ob",model_name \
forrest@0	846	,"ob",model_name \
forrest@0	847	/)
forrest@0	848
forrest@0	849	ncol1 = dimsizes(col_header1)
forrest@0	850	ncol2 = dimsizes(col_header2)
forrest@0	851
forrest@0	852	ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
forrest@0	853	x21 = (/x2(1),0.995/) ; start from end of x2
forrest@0	854	yhalf = (y1(0)+y1(1))*0.5
forrest@0	855	y21 = (/yhalf,y1(1)/) ; top half of y1
forrest@0	856	text21 = col_header1
forrest@0	857	res21 = True
forrest@0	858	res21@txFontHeightF = 0.015
forrest@0	859	res21@gsFillColor = "Gray"
forrest@0	860
forrest@0	861	ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
forrest@0	862	x22 = x21 ; start from end of x1
forrest@0	863	y22 = (/y1(0),yhalf/) ; bottom half of y1
forrest@0	864	text22 = col_header2
forrest@0	865	res22 = True
forrest@0	866	res22@txFontHeightF = 0.015
forrest@0	867	res22@gsFillColor = "Gray"
forrest@0	868
forrest@0	869	; Row header (equally space in y2)
forrest@0	870	ncr3 = (/nrow,1/) ; nrow rows, 1 columns
forrest@0	871	x3 = x1 ; same as x1
forrest@0	872	y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
forrest@0	873	text3 = row_header
forrest@0	874	res3 = True
forrest@0	875	res3@txFontHeightF = 0.010
forrest@0	876	res3@gsFillColor = "Gray"
forrest@0	877
forrest@0	878	; Main table body-1
forrest@0	879	ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
forrest@0	880	x4 = x2 ; Start and end x
forrest@0	881	y4 = y3 ; Start and end Y
forrest@0	882	text4 = new((/nrow,ncol/),string)
forrest@0	883
forrest@0	884	color_fill4 = new((/nrow,ncol/),string)
forrest@0	885	color_fill4 = "white"
forrest@0	886	; color_fill4(:,ncol-1) = "grey"
forrest@0	887	; color_fill4(nrow-1,:) = "green"
forrest@0	888
forrest@0	889	res4 = True ; Set up resource list
forrest@0	890	; res4@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	891	res4@txFontHeightF = 0.015
forrest@0	892	res4@gsFillColor = color_fill4
forrest@0	893
forrest@0	894	delete (color_fill4)
forrest@0	895
forrest@0	896	; Main table body-2
forrest@0	897	ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
forrest@0	898	x42 = x21 ; Start and end x
forrest@0	899	y42 = y3 ; Start and end Y
forrest@0	900	text42 = new((/nrow,ncol2/),string)
forrest@0	901
forrest@0	902	color_fill42 = new((/nrow,ncol2/),string)
forrest@0	903	color_fill42 = "white"
forrest@0	904	; color_fill42(:,ncol-1) = "grey"
forrest@0	905	; color_fill42(nrow-1,:) = "green"
forrest@0	906
forrest@0	907	res42 = True ; Set up resource list
forrest@0	908	; res42@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	909	res42@txFontHeightF = 0.015
forrest@0	910	res42@gsFillColor = color_fill42
forrest@0	911
forrest@0	912	delete (color_fill42)
forrest@0	913	;-------------------------------------------------------------------
forrest@0	914	; for table value
forrest@0	915
forrest@0	916	do n = 0,nrow-1
forrest@0	917	text4(n,0) = sprintf("%5.2f", y81(n+41))
forrest@0	918	text4(n,1) = sprintf("%5.2f", x81(n+41))
forrest@0	919	end do
forrest@0	920
forrest@0	921	do n = 0,nrow-1
forrest@0	922	text42(n,0) = sprintf("%5.2f", npp81(n+41))
forrest@0	923	text42(n,1) = sprintf("%5.2f", nppmod81(n+41))
forrest@0	924	text42(n,2) = sprintf("%5.2f", rain81(n+41))
forrest@0	925	text42(n,3) = sprintf("%5.2f", rainmod81(n+41))
forrest@0	926	end do
forrest@0	927	;-------------------------------------------------------------------
forrest@0	928
forrest@0	929	plot_name = "table_site81_model_vs_ob2"
forrest@0	930
forrest@0	931	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	932	;------------------------------------------
forrest@0	933	; for table title
forrest@0	934
forrest@0	935	gRes = True
forrest@0	936	gRes@txFontHeightF = 0.02
forrest@0	937	; gRes@txAngleF = 90
forrest@0	938
forrest@0	939	title_text = "Model vs Observation at 81 Sites (2)"
forrest@0	940
forrest@0	941	gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
forrest@0	942	;------------------------------------------
forrest@0	943
forrest@0	944	gsn_table(wks,ncr1,x1,y1,text1,res1)
forrest@0	945	gsn_table(wks,ncr2,x2,y2,text2,res2)
forrest@0	946	gsn_table(wks,ncr21,x21,y21,text21,res21)
forrest@0	947	gsn_table(wks,ncr22,x22,y22,text22,res22)
forrest@0	948	gsn_table(wks,ncr3,x3,y3,text3,res3)
forrest@0	949	gsn_table(wks,ncr4,x4,y4,text4,res4)
forrest@0	950	gsn_table(wks,ncr42,x42,y42,text42,res42)
forrest@0	951
forrest@0	952	frame(wks)
forrest@0	953	clear (wks)
forrest@0	954
forrest@0	955	delete (col_header)
forrest@0	956	delete (row_header)
forrest@0	957	delete (res1)
forrest@0	958	delete (ncr1)
forrest@0	959	delete (text1)
forrest@0	960	delete (res2)
forrest@0	961	delete (ncr2)
forrest@0	962	delete (text2)
forrest@0	963	delete (res21)
forrest@0	964	delete (ncr21)
forrest@0	965	delete (text21)
forrest@0	966	delete (res22)
forrest@0	967	delete (ncr22)
forrest@0	968	delete (text22)
forrest@0	969	delete (res3)
forrest@0	970	delete (ncr3)
forrest@0	971	delete (text3)
forrest@0	972	delete (res4)
forrest@0	973	delete (ncr4)
forrest@0	974	delete (text4)
forrest@0	975	delete (res42)
forrest@0	976	delete (ncr42)
forrest@0	977	delete (text42)
forrest@0	978
forrest@0	979	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	980	system("rm "+plot_name+"."+plot_type)
forrest@0	981	system("rm "+plot_name+"-1."+plot_type_new)
forrest@0	982	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	983
forrest@0	984	;**************************************************************************
forrest@0	985	;(B) histogram 81
forrest@0	986	;**************************************************************************
forrest@0	987	;get data
forrest@0	988	RAIN1_1D = ndtooned(rain81)
forrest@0	989	RAIN2_1D = ndtooned(rainmod81)
forrest@0	990	NPP1_1D = ndtooned(npp81)
forrest@0	991	NPP2_1D = ndtooned(nppmod81)
forrest@0	992
forrest@0	993	; number of bin
forrest@0	994	nx = 8
forrest@0	995
forrest@0	996	xvalues = new((/2,nx/),float)
forrest@0	997	yvalues = new((/2,nx/),float)
forrest@0	998	mn_yvalues = new((/2,nx/),float)
forrest@0	999	mx_yvalues = new((/2,nx/),float)
forrest@0	1000	dx4 = new((/1/),float)
forrest@0	1001
forrest@0	1002	get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
forrest@0	1003	,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
forrest@0	1004
forrest@0	1005	;----------------------------------------
forrest@0	1006	;compute correlation coeff and M score
forrest@0	1007
forrest@0	1008	u = yvalues(0,:)
forrest@0	1009	v = yvalues(1,:)
forrest@0	1010
forrest@0	1011	good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0	1012	uu = u(good)
forrest@0	1013	vv = v(good)
forrest@0	1014
forrest@0	1015	ccr81h = esccr(uu,vv,0)
forrest@0	1016	;print (ccr81h)
forrest@0	1017
forrest@0	1018	;new eq
forrest@0	1019	bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
forrest@0	1020	M81h = (1.- (bias/dimsizes(uu)))*5.
forrest@0	1021	print (M81h)
forrest@0	1022
forrest@0	1023	delete (u)
forrest@0	1024	delete (v)
forrest@0	1025	delete (uu)
forrest@0	1026	delete (vv)
forrest@0	1027	;----------------------------------------------------------------------
forrest@0	1028	; histogram res
forrest@0	1029	resh = True
forrest@0	1030	resh@gsnMaximize = True
forrest@0	1031	resh@gsnDraw = False
forrest@0	1032	resh@gsnFrame = False
forrest@0	1033	resh@xyMarkLineMode = "Markers"
forrest@0	1034	resh@xyMarkerSizeF = 0.014
forrest@0	1035	resh@xyMarker = 16
forrest@0	1036	resh@xyMarkerColors = (/"Brown","Blue"/)
forrest@0	1037	resh@trYMinF = min(mn_yvalues) - 10.
forrest@0	1038	resh@trYMaxF = max(mx_yvalues) + 10.
forrest@0	1039
forrest@0	1040	resh@tiYAxisString = "NPP (g C/m2/year)"
forrest@0	1041	resh@tiXAxisString = "Precipitation (m/year)"
forrest@0	1042
forrest@0	1043	max_bar = new((/2,nx/),graphic)
forrest@0	1044	min_bar = new((/2,nx/),graphic)
forrest@0	1045	max_cap = new((/2,nx/),graphic)
forrest@0	1046	min_cap = new((/2,nx/),graphic)
forrest@0	1047
forrest@0	1048	lnres = True
forrest@0	1049	line_colors = (/"brown","blue"/)
forrest@0	1050
forrest@0	1051	;**************************************************************************
forrest@0	1052	;(B)-1 histogram plot, ob 81 site
forrest@0	1053	;**************************************************************************
forrest@0	1054
forrest@0	1055	plot_name = "histogram_ob_81"
forrest@0	1056	title = "Observed 81 site"
forrest@0	1057	resh@tiMainString = title
forrest@0	1058
forrest@0	1059	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	1060
forrest@0	1061	xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
forrest@0	1062
forrest@0	1063	;-------------------------------
forrest@0	1064	;Attach the vertical bar and the horizontal cap line
forrest@0	1065
forrest@0	1066	delete (x1)
forrest@0	1067	delete (x2)
forrest@0	1068	delete (y1)
forrest@0	1069	delete (y2)
forrest@0	1070
forrest@0	1071	do nd=0,0
forrest@0	1072	lnres@gsLineColor = line_colors(nd)
forrest@0	1073	do i=0,nx-1
forrest@0	1074
forrest@0	1075	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	1076	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	1077	;
forrest@0	1078	; Attach the vertical bar, both above and below the marker.
forrest@0	1079	;
forrest@0	1080	x1 = xvalues(nd,i)
forrest@0	1081	y1 = yvalues(nd,i)
forrest@0	1082	y2 = mn_yvalues(nd,i)
forrest@0	1083	plx = (/x1,x1/)
forrest@0	1084	ply = (/y1,y2/)
forrest@0	1085	min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0	1086
forrest@0	1087	y2 = mx_yvalues(nd,i)
forrest@0	1088	plx = (/x1,x1/)
forrest@0	1089	ply = (/y1,y2/)
forrest@0	1090	max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0	1091	;
forrest@0	1092	; Attach the horizontal cap line, both above and below the marker.
forrest@0	1093	;
forrest@0	1094	x1 = xvalues(nd,i) - dx4
forrest@0	1095	x2 = xvalues(nd,i) + dx4
forrest@0	1096	y1 = mn_yvalues(nd,i)
forrest@0	1097	plx = (/x1,x2/)
forrest@0	1098	ply = (/y1,y1/)
forrest@0	1099	min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0	1100
forrest@0	1101	y1 = mx_yvalues(nd,i)
forrest@0	1102	plx = (/x1,x2/)
forrest@0	1103	ply = (/y1,y1/)
forrest@0	1104	max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0	1105	end if
forrest@0	1106	end do
forrest@0	1107	end do
forrest@0	1108
forrest@0	1109	draw(xy)
forrest@0	1110	frame(wks)
forrest@0	1111	clear (wks)
forrest@0	1112
forrest@0	1113	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1114	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1115	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1116	system("rm "+plot_name+"."+plot_type)
forrest@0	1117
forrest@0	1118	;****************************************************************************
forrest@0	1119	;(B)-2 histogram plot, model vs ob 81 site
forrest@0	1120	;****************************************************************************
forrest@0	1121
forrest@0	1122	plot_name = "histogram_model_vs_ob_81"
forrest@0	1123	title = model_name+ " vs Observed 81 site"
forrest@0	1124	resh@tiMainString = title
forrest@0	1125
forrest@0	1126	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1127
forrest@0	1128	;-----------------------------
forrest@0	1129	; Add a boxed legend using the more simple method, which won't have
forrest@0	1130	; vertical lines going through the markers.
forrest@0	1131
forrest@0	1132	resh@pmLegendDisplayMode = "Always"
forrest@0	1133	; resh@pmLegendWidthF = 0.1
forrest@0	1134	resh@pmLegendWidthF = 0.08
forrest@0	1135	resh@pmLegendHeightF = 0.05
forrest@0	1136	resh@pmLegendOrthogonalPosF = -1.17
forrest@0	1137	; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0	1138	; resh@pmLegendParallelPosF = 0.18
forrest@0	1139	resh@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0	1140
forrest@0	1141	; resh@lgPerimOn = False
forrest@0	1142	resh@lgLabelFontHeightF = 0.015
forrest@0	1143	resh@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0	1144	;-----------------------------
forrest@0	1145	tRes = True
forrest@0	1146	tRes@txFontHeightF = 0.025
forrest@0	1147
forrest@0	1148	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
forrest@0	1149
forrest@0	1150	gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0	1151
forrest@0	1152	xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
forrest@0	1153	;-------------------------------
forrest@0	1154	;Attach the vertical bar and the horizontal cap line
forrest@0	1155
forrest@0	1156	do nd=0,1
forrest@0	1157	lnres@gsLineColor = line_colors(nd)
forrest@0	1158	do i=0,nx-1
forrest@0	1159
forrest@0	1160	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	1161	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	1162	;
forrest@0	1163	; Attach the vertical bar, both above and below the marker.
forrest@0	1164	;
forrest@0	1165	x1 = xvalues(nd,i)
forrest@0	1166	y1 = yvalues(nd,i)
forrest@0	1167	y2 = mn_yvalues(nd,i)
forrest@0	1168	min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1169
forrest@0	1170	y2 = mx_yvalues(nd,i)
forrest@0	1171	max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1172	;
forrest@0	1173	; Attach the horizontal cap line, both above and below the marker.
forrest@0	1174	;
forrest@0	1175	x1 = xvalues(nd,i) - dx4
forrest@0	1176	x2 = xvalues(nd,i) + dx4
forrest@0	1177	y1 = mn_yvalues(nd,i)
forrest@0	1178	min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1179
forrest@0	1180	y1 = mx_yvalues(nd,i)
forrest@0	1181	max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1182	end if
forrest@0	1183	end do
forrest@0	1184	end do
forrest@0	1185
forrest@0	1186	draw(xy)
forrest@0	1187	frame(wks)
forrest@0	1188	clear(wks)
forrest@0	1189
forrest@0	1190	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1191	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1192	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1193	system("rm "+plot_name+"."+plot_type)
forrest@0	1194
forrest@0	1195	delete (RAIN1_1D)
forrest@0	1196	delete (RAIN2_1D)
forrest@0	1197	delete (NPP1_1D)
forrest@0	1198	delete (NPP2_1D)
forrest@0	1199	;delete (range)
forrest@0	1200	delete (xvalues)
forrest@0	1201	delete (yvalues)
forrest@0	1202	delete (mn_yvalues)
forrest@0	1203	delete (mx_yvalues)
forrest@0	1204	delete (good)
forrest@0	1205	delete (max_bar)
forrest@0	1206	delete (min_bar)
forrest@0	1207	delete (max_cap)
forrest@0	1208	delete (min_cap)
forrest@0	1209
forrest@0	1210	;**************************************************************************
forrest@0	1211	;(B) histogram 933
forrest@0	1212	;**************************************************************************
forrest@0	1213
forrest@0	1214	; get data
forrest@0	1215	RAIN1_1D = ndtooned(rain933)
forrest@0	1216	RAIN2_1D = ndtooned(rainmod933)
forrest@0	1217	NPP1_1D = ndtooned(npp933)
forrest@0	1218	NPP2_1D = ndtooned(nppmod933)
forrest@0	1219
forrest@0	1220	; number of bin
forrest@0	1221	nx = 10
forrest@0	1222
forrest@0	1223	xvalues = new((/2,nx/),float)
forrest@0	1224	yvalues = new((/2,nx/),float)
forrest@0	1225	mn_yvalues = new((/2,nx/),float)
forrest@0	1226	mx_yvalues = new((/2,nx/),float)
forrest@0	1227	dx4 = new((/1/),float)
forrest@0	1228
forrest@0	1229	get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
forrest@0	1230	,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
forrest@0	1231
forrest@0	1232	;----------------------------------------
forrest@0	1233	;compute correlation coeff and M score
forrest@0	1234
forrest@0	1235	u = yvalues(0,:)
forrest@0	1236	v = yvalues(1,:)
forrest@0	1237
forrest@0	1238	good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0	1239	uu = u(good)
forrest@0	1240	vv = v(good)
forrest@0	1241
forrest@0	1242	ccr933h = esccr(uu,vv,0)
forrest@0	1243	;print (ccr933h)
forrest@0	1244
forrest@0	1245	;new eq
forrest@0	1246	bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
forrest@0	1247	M933h = (1.- (bias/dimsizes(uu)))*5.
forrest@0	1248	print (M933h)
forrest@0	1249
forrest@0	1250	delete (u)
forrest@0	1251	delete (v)
forrest@0	1252	delete (uu)
forrest@0	1253	delete (vv)
forrest@0	1254	;----------------------------------------------------------------------
forrest@0	1255	; histogram res
forrest@0	1256	delete (resh)
forrest@0	1257	resh = True
forrest@0	1258	resh@gsnMaximize = True
forrest@0	1259	resh@gsnDraw = False
forrest@0	1260	resh@gsnFrame = False
forrest@0	1261	resh@xyMarkLineMode = "Markers"
forrest@0	1262	resh@xyMarkerSizeF = 0.014
forrest@0	1263	resh@xyMarker = 16
forrest@0	1264	resh@xyMarkerColors = (/"Brown","Blue"/)
forrest@0	1265	resh@trYMinF = min(mn_yvalues) - 10.
forrest@0	1266	resh@trYMaxF = max(mx_yvalues) + 10.
forrest@0	1267
forrest@0	1268	resh@tiYAxisString = "NPP (g C/m2/year)"
forrest@0	1269	resh@tiXAxisString = "Precipitation (m/year)"
forrest@0	1270
forrest@0	1271	max_bar = new((/2,nx/),graphic)
forrest@0	1272	min_bar = new((/2,nx/),graphic)
forrest@0	1273	max_cap = new((/2,nx/),graphic)
forrest@0	1274	min_cap = new((/2,nx/),graphic)
forrest@0	1275
forrest@0	1276	lnres = True
forrest@0	1277	line_colors = (/"brown","blue"/)
forrest@0	1278
forrest@0	1279	;**************************************************************************
forrest@0	1280	;(B)-3 histogram plot, ob 933 site
forrest@0	1281	;**************************************************************************
forrest@0	1282
forrest@0	1283	plot_name = "histogram_ob_933"
forrest@0	1284	title = "Observed 933 site"
forrest@0	1285	resh@tiMainString = title
forrest@0	1286
forrest@0	1287	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	1288
forrest@0	1289	xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
forrest@0	1290
forrest@0	1291	;-------------------------------
forrest@0	1292	;Attach the vertical bar and the horizontal cap line
forrest@0	1293
forrest@0	1294	do nd=0,0
forrest@0	1295	lnres@gsLineColor = line_colors(nd)
forrest@0	1296	do i=0,nx-1
forrest@0	1297
forrest@0	1298	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	1299	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	1300
forrest@0	1301	; Attach the vertical bar, both above and below the marker.
forrest@0	1302
forrest@0	1303	x1 = xvalues(nd,i)
forrest@0	1304	y1 = yvalues(nd,i)
forrest@0	1305	y2 = mn_yvalues(nd,i)
forrest@0	1306	min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1307
forrest@0	1308	y2 = mx_yvalues(nd,i)
forrest@0	1309	max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1310
forrest@0	1311	; Attach the horizontal cap line, both above and below the marker.
forrest@0	1312
forrest@0	1313	x1 = xvalues(nd,i) - dx4
forrest@0	1314	x2 = xvalues(nd,i) + dx4
forrest@0	1315	y1 = mn_yvalues(nd,i)
forrest@0	1316	min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1317
forrest@0	1318	y1 = mx_yvalues(nd,i)
forrest@0	1319	max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1320	end if
forrest@0	1321	end do
forrest@0	1322	end do
forrest@0	1323
forrest@0	1324	draw(xy)
forrest@0	1325	frame(wks)
forrest@0	1326	delete (xy)
forrest@0	1327	clear (wks)
forrest@0	1328
forrest@0	1329	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1330	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1331	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1332	system("rm "+plot_name+"."+plot_type)
forrest@0	1333
forrest@0	1334	;****************************************************************************
forrest@0	1335	;(B)-4 histogram plot, model vs ob 933 site
forrest@0	1336	;****************************************************************************
forrest@0	1337
forrest@0	1338	plot_name = "histogram_model_vs_ob_933"
forrest@0	1339	title = model_name+ " vs Observed 933 site"
forrest@0	1340	resh@tiMainString = title
forrest@0	1341
forrest@0	1342	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1343
forrest@0	1344	;-----------------------------
forrest@0	1345	; Add a boxed legend using the more simple method, which won't have
forrest@0	1346	; vertical lines going through the markers.
forrest@0	1347
forrest@0	1348	resh@pmLegendDisplayMode = "Always"
forrest@0	1349	; resh@pmLegendWidthF = 0.1
forrest@0	1350	resh@pmLegendWidthF = 0.08
forrest@0	1351	resh@pmLegendHeightF = 0.05
forrest@0	1352	resh@pmLegendOrthogonalPosF = -1.17
forrest@0	1353	; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0	1354	; resh@pmLegendParallelPosF = 0.18
forrest@0	1355	resh@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0	1356
forrest@0	1357	; resh@lgPerimOn = False
forrest@0	1358	resh@lgLabelFontHeightF = 0.015
forrest@0	1359	resh@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0	1360	;-----------------------------
forrest@0	1361	tRes = True
forrest@0	1362	tRes@txFontHeightF = 0.025
forrest@0	1363
forrest@0	1364	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
forrest@0	1365
forrest@0	1366	gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0	1367
forrest@0	1368	xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
forrest@0	1369	;-------------------------------
forrest@0	1370	;Attach the vertical bar and the horizontal cap line
forrest@0	1371
forrest@0	1372	do nd=0,1
forrest@0	1373	lnres@gsLineColor = line_colors(nd)
forrest@0	1374	do i=0,nx-1
forrest@0	1375
forrest@0	1376	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	1377	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	1378	;
forrest@0	1379	; Attach the vertical bar, both above and below the marker.
forrest@0	1380	;
forrest@0	1381	x1 = xvalues(nd,i)
forrest@0	1382	y1 = yvalues(nd,i)
forrest@0	1383	y2 = mn_yvalues(nd,i)
forrest@0	1384	min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1385
forrest@0	1386	y2 = mx_yvalues(nd,i)
forrest@0	1387	max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1388	;
forrest@0	1389	; Attach the horizontal cap line, both above and below the marker.
forrest@0	1390	;
forrest@0	1391	x1 = xvalues(nd,i) - dx4
forrest@0	1392	x2 = xvalues(nd,i) + dx4
forrest@0	1393	y1 = mn_yvalues(nd,i)
forrest@0	1394	min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1395
forrest@0	1396	y1 = mx_yvalues(nd,i)
forrest@0	1397	max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1398	end if
forrest@0	1399	end do
forrest@0	1400	end do
forrest@0	1401
forrest@0	1402	draw(xy)
forrest@0	1403	frame(wks)
forrest@0	1404	delete(xy)
forrest@0	1405	clear(wks)
forrest@0	1406
forrest@0	1407	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1408	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1409	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1410	system("rm "+plot_name+"."+plot_type)
forrest@0	1411
forrest@0	1412	;***************************************************************************
forrest@0	1413	;(C) global contour
forrest@0	1414	;***************************************************************************
forrest@0	1415
forrest@0	1416	;res
forrest@0	1417	resg = True ; Use plot options
forrest@0	1418	resg@cnFillOn = True ; Turn on color fill
forrest@0	1419	resg@gsnSpreadColors = True ; use full colormap
forrest@0	1420	; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0	1421	; resg@lbLabelAutoStride = True
forrest@0	1422	resg@cnLinesOn = False ; Turn off contourn lines
forrest@0	1423	resg@mpFillOn = False ; Turn off map fill
forrest@0	1424
forrest@0	1425	resg@gsnSpreadColors = True ; use full colormap
forrest@0	1426	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1427	resg@cnMinLevelValF = 0. ; Min level
forrest@0	1428	resg@cnMaxLevelValF = 2200. ; Max level
forrest@0	1429	resg@cnLevelSpacingF = 200. ; interval
forrest@0	1430
forrest@0	1431	;****************************************************************************
forrest@0	1432	;(C)-1 global contour plot, ob
forrest@0	1433	;****************************************************************************
forrest@0	1434
forrest@0	1435	delta = 0.00000000001
forrest@0	1436	nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
forrest@0	1437
forrest@0	1438	plot_name = "global_ob"
forrest@0	1439	title = "Observed MODIS MOD 17"
forrest@0	1440	resg@tiMainString = title
forrest@0	1441
forrest@0	1442	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1443	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1444
forrest@0	1445	plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
forrest@0	1446
forrest@0	1447	frame(wks)
forrest@0	1448	clear (wks)
forrest@0	1449
forrest@0	1450	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1451	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1452	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1453	system("rm "+plot_name+"."+plot_type)
forrest@0	1454
forrest@0	1455	;****************************************************************************
forrest@0	1456	;(C)-2 global contour plot, model
forrest@0	1457	;****************************************************************************
forrest@0	1458
forrest@0	1459	plot_name = "global_model"
forrest@0	1460	title = "Model "+ model_name
forrest@0	1461	resg@tiMainString = title
forrest@0	1462
forrest@0	1463	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1464	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1465
forrest@0	1466	plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
forrest@0	1467
forrest@0	1468	frame(wks)
forrest@0	1469	clear (wks)
forrest@0	1470
forrest@0	1471	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1472	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1473	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1474	system("rm "+plot_name+"."+plot_type)
forrest@0	1475
forrest@0	1476	;****************************************************************************
forrest@0	1477	;(C)-3 global contour plot, model vs ob
forrest@0	1478	;****************************************************************************
forrest@0	1479
forrest@0	1480	plot_name = "global_model_vs_ob"
forrest@0	1481
forrest@0	1482	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1483	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1484
forrest@0	1485	delete (plot)
forrest@0	1486	plot=new(3,graphic) ; create graphic array
forrest@0	1487
forrest@0	1488	resg@gsnFrame = False ; Do not draw plot
forrest@0	1489	resg@gsnDraw = False ; Do not advance frame
forrest@0	1490
forrest@0	1491	; compute correlation coef and M score
forrest@0	1492
forrest@0	1493	uu1 = ndtooned(nppmod)
forrest@0	1494	vv1 = ndtooned(nppglobe)
forrest@0	1495
forrest@0	1496	delete (good)
forrest@0	1497	good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
forrest@0	1498
forrest@0	1499	ug = uu1(good)
forrest@0	1500	vg = vv1(good)
forrest@0	1501
forrest@0	1502	ccrG = esccr(ug,vg,0)
forrest@0	1503	; print (ccrG)
forrest@0	1504
forrest@0	1505	MG = (ccrGccrG) 5.0
forrest@0	1506	print (MG)
forrest@0	1507
forrest@0	1508	; plot correlation coef
forrest@0	1509
forrest@0	1510	gRes = True
forrest@0	1511	gRes@txFontHeightF = 0.02
forrest@0	1512	gRes@txAngleF = 90
forrest@0	1513
forrest@0	1514	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
forrest@0	1515
forrest@0	1516	gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0	1517	;--------------------------------------------------------------------
forrest@0	1518	;(a) ob
forrest@0	1519
forrest@0	1520	title = "Observed MODIS MOD 17"
forrest@0	1521	resg@tiMainString = title
forrest@0	1522
forrest@0	1523	plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
forrest@0	1524
forrest@0	1525	;(b) model
forrest@0	1526
forrest@0	1527	title = "Model "+ model_name
forrest@0	1528	resg@tiMainString = title
forrest@0	1529
forrest@0	1530	plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
forrest@0	1531
forrest@0	1532	;(c) model-ob
forrest@0	1533
forrest@0	1534	zz = nppmod
forrest@0	1535	zz = nppmod - nppglobe
forrest@0	1536	title = "Model_"+model_name+" - Observed"
forrest@0	1537
forrest@0	1538	resg@cnMinLevelValF = -500 ; Min level
forrest@0	1539	resg@cnMaxLevelValF = 500. ; Max level
forrest@0	1540	resg@cnLevelSpacingF = 50. ; interval
forrest@0	1541	resg@tiMainString = title
forrest@0	1542
forrest@0	1543	plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
forrest@0	1544
forrest@0	1545	pres = True ; panel plot mods desired
forrest@0	1546	; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0	1547	; indiv. plots in panel
forrest@0	1548	pres@gsnMaximize = True ; fill the page
forrest@0	1549
forrest@0	1550	gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
forrest@0	1551
forrest@0	1552	frame (wks)
forrest@0	1553	clear (wks)
forrest@0	1554	delete (plot)
forrest@0	1555
forrest@0	1556	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1557	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1558	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1559	system("rm "+plot_name+"."+plot_type)
forrest@0	1560
forrest@0	1561	;***************************************************************************
forrest@0	1562	;(D)-1 zonal line plot, ob
forrest@0	1563	;***************************************************************************
forrest@0	1564
forrest@0	1565	vv = zonalAve(nppglobe)
forrest@0	1566	vv@long_name = nppglobe@long_name
forrest@0	1567
forrest@0	1568	plot_name = "zonal_ob"
forrest@0	1569	title = "Observed MODIS MOD 17"
forrest@0	1570
forrest@0	1571	resz = True
forrest@0	1572	resz@tiMainString = title
forrest@0	1573
forrest@0	1574	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1575
forrest@0	1576	plot = gsn_csm_xy (wks,ym,vv,resz)
forrest@0	1577
forrest@0	1578	frame(wks)
forrest@0	1579	clear (wks)
forrest@0	1580
forrest@0	1581	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1582	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1583	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1584	system("rm "+plot_name+"."+plot_type)
forrest@0	1585
forrest@0	1586	;****************************************************************************
forrest@0	1587	;(D)-2 zonal line plot, model vs ob
forrest@0	1588	;****************************************************************************
forrest@0	1589
forrest@0	1590	s = new ((/2,dimsizes(ym)/), float)
forrest@0	1591
forrest@0	1592	s(0,:) = zonalAve(nppglobe)
forrest@0	1593	s(1,:) = zonalAve(nppmod)
forrest@0	1594
forrest@0	1595	s@long_name = nppglobe@long_name
forrest@0	1596	;-------------------------------------------
forrest@0	1597	; compute correlation coef and M score
forrest@0	1598
forrest@0	1599	ccrZ = esccr(s(1,:), s(0,:),0)
forrest@0	1600	; print (ccrZ)
forrest@0	1601
forrest@0	1602	MZ = (ccrZccrZ) 5.0
forrest@0	1603	print (MZ)
forrest@0	1604	;-------------------------------------------
forrest@0	1605	plot_name = "zonal_model_vs_ob"
forrest@0	1606	title = "Zonal Average"
forrest@0	1607	resz@tiMainString = title
forrest@0	1608
forrest@0	1609	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	1610
forrest@0	1611	; resz@vpHeightF = 0.4 ; change aspect ratio of plot
forrest@0	1612	; resz@vpWidthF = 0.7
forrest@0	1613
forrest@0	1614	resz@xyMonoLineColor = "False" ; want colored lines
forrest@0	1615	resz@xyLineColors = (/"black","red"/) ; colors chosen
forrest@0	1616	; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
forrest@0	1617	resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
forrest@0	1618	resz@xyDashPatterns = (/0.,0./) ; make all lines solid
forrest@0	1619
forrest@0	1620	resz@tiYAxisString = s@long_name ; add a axis title
forrest@0	1621	resz@txFontHeightF = 0.0195 ; change title font heights
forrest@0	1622
forrest@0	1623	; Legent
forrest@0	1624	resz@pmLegendDisplayMode = "Always" ; turn on legend
forrest@0	1625	resz@pmLegendSide = "Top" ; Change location of
forrest@0	1626	; resz@pmLegendParallelPosF = .45 ; move units right
forrest@0	1627	resz@pmLegendParallelPosF = .82 ; move units right
forrest@0	1628	resz@pmLegendOrthogonalPosF = -0.4 ; move units down
forrest@0	1629
forrest@0	1630	resz@pmLegendWidthF = 0.10 ; Change width and
forrest@0	1631	resz@pmLegendHeightF = 0.10 ; height of legend.
forrest@0	1632	resz@lgLabelFontHeightF = .02 ; change font height
forrest@0	1633	; resz@lgTitleOn = True ; turn on legend title
forrest@0	1634	; resz@lgTitleString = "Example" ; create legend title
forrest@0	1635	; resz@lgTitleFontHeightF = .025 ; font of legend title
forrest@0	1636	resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
forrest@0	1637	;--------------------------------------------------------------------
forrest@0	1638	zRes = True
forrest@0	1639	zRes@txFontHeightF = 0.025
forrest@0	1640
forrest@0	1641	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
forrest@0	1642
forrest@0	1643	gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
forrest@0	1644	;--------------------------------------------------------------------
forrest@0	1645
forrest@0	1646	plot = gsn_csm_xy (wks,ym,s,resz)
forrest@0	1647
forrest@0	1648	frame(wks)
forrest@0	1649	clear (wks)
forrest@0	1650
forrest@0	1651	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0	1652	system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0	1653	system("rm "+plot_name+"-*."+plot_type_new)
forrest@0	1654	system("rm "+plot_name+"."+plot_type)
forrest@0	1655
forrest@0	1656	;***************************************************************************
forrest@0	1657	; tar up output plots
forrest@0	1658	;***************************************************************************
forrest@0	1659
forrest@0	1660	temp_name = "temp." + model_name
forrest@0	1661	system("mkdir -p " + temp_name)
forrest@0	1662	system("mv *.png " + temp_name)
forrest@0	1663	system("tar cf "+ temp_name +".tar " + temp_name)
forrest@0	1664
forrest@0	1665	;***************************************************************************
forrest@0	1666	end
forrest@0	1667

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--