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

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