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

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