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

forrest@0	1	;********************************************************
forrest@0	2	; histogram normalized by rain and compute correleration
forrest@0	3	;********************************************************
forrest@0	4	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
forrest@0	5	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
forrest@0	6	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0	7
forrest@0	8	procedure pminmax(data:numeric,name:string)
forrest@0	9	begin
forrest@0	10	print ("min/max " + name + " = " + min(data) + "/" + max(data))
forrest@0	11	if(isatt(data,"units")) then
forrest@0	12	print (name + " units = " + data@units)
forrest@0	13	end if
forrest@0	14	end
forrest@0	15
forrest@0	16	; Main code.
forrest@0	17	begin
forrest@0	18
forrest@0	19	nclass = 20
forrest@0	20
forrest@0	21	plot_type = "ps"
forrest@0	22	plot_type_new = "png"
forrest@0	23
forrest@0	24	;************************************************
forrest@0	25	; read in data: model
forrest@0	26	;************************************************
forrest@0	27	;film = "b30.061n_1995-2004_MONS_climo_lnd.nc"
forrest@0	28	;model_name = "b30.061n"
forrest@0	29	;model_grid = "T31"
forrest@0	30
forrest@0	31	;film = "newcn05_ncep_1i_MONS_climo_lnd.nc"
forrest@0	32	;model_name = "newcn"
forrest@0	33	;model_grid = "1.9"
forrest@0	34
forrest@0	35	;film = "i01.06cn_1798-2004_MONS_climo.nc"
forrest@0	36	;model_name = "06cn"
forrest@0	37	;model_grid = "T42"
forrest@0	38
forrest@0	39	;film = "i01.06casa_1798-2004_MONS_climo.nc"
forrest@0	40	;model_name = "06casa"
forrest@0	41	;model_grid = "T42"
forrest@0	42
forrest@0	43	film = "i01.10cn_1948-2004_MONS_climo.nc"
forrest@0	44	model_name = "10cn"
forrest@0	45	model_grid = "T42"
forrest@0	46
forrest@0	47	;film = "i01.10casa_1948-2004_MONS_climo.nc"
forrest@0	48	;model_name = "10casa"
forrest@0	49	;model_grid = "T42"
forrest@0	50
forrest@0	51	html_name = "table.html." + model_name
forrest@0	52	html_new = html_name +".new"
forrest@0	53	system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
forrest@0	54	"mv -f "+html_new+" "+html_name)
forrest@0	55	;------------------------------------------------
forrest@0	56	dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0	57	fm = addfile(dirm+film,"r")
forrest@0	58
forrest@0	59	laimod = fm->TLAI
forrest@0	60
forrest@0	61	;************************************************
forrest@0	62	; read in data: observed
forrest@0	63	;************************************************
forrest@0	64
forrest@0	65	ob_name = "MODIS MOD 15A2 2000-2005"
forrest@0	66
forrest@0	67	diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
forrest@0	68	filo1 = "land_class_"+model_grid+".nc"
forrest@0	69	filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc"
forrest@0	70
forrest@0	71	fo1 = addfile(diro+filo1,"r")
forrest@0	72	fo2 = addfile(diro+filo2,"r")
forrest@0	73
forrest@0	74	classob = tofloat(fo1->LAND_CLASS)
forrest@0	75	laiob = fo2->LAI
forrest@0	76
forrest@0	77	;*******************************************************************
forrest@0	78	; for plotting table
forrest@0	79	;*******************************************************************
forrest@0	80	; table header name
forrest@0	81	table_header_name = "LAI"
forrest@0	82
forrest@0	83	; column (not including header column)
forrest@0	84
forrest@0	85	col_header1 = (/"Mean","Max","Phase","Growth"/)
forrest@0	86	col_header2 = (/"ob","model","M" \
forrest@0	87	,"ob","model","M" \
forrest@0	88	,"ob","model","M" \
forrest@0	89	,"ob","model","M" \
forrest@0	90	/)
forrest@0	91
forrest@0	92	ncol1 = dimsizes(col_header1)
forrest@0	93	ncol2 = dimsizes(col_header2)
forrest@0	94	ncol = ncol2
forrest@0	95
forrest@0	96	; row (not including header row)
forrest@0	97	row_header = (/"Water Bodies" \
forrest@0	98	,"Evergreen Needleleaf Forests" \
forrest@0	99	,"Evergreen Broadleaf Forests" \
forrest@0	100	,"Deciduous Needleleaf Forest" \
forrest@0	101	,"Deciduous Broadleaf Forests" \
forrest@0	102	,"Mixed Forests" \
forrest@0	103	,"Closed Bushlands" \
forrest@0	104	,"Open Bushlands" \
forrest@0	105	,"Woody Savannas (S. Hem.)" \
forrest@0	106	,"Savannas (S. Hem.)" \
forrest@0	107	,"Grasslands" \
forrest@0	108	,"Permanent Wetlands" \
forrest@0	109	,"Croplands" \
forrest@0	110	,"Urban and Built-Up" \
forrest@0	111	,"Cropland/Natural Vegetation Mosaic" \
forrest@0	112	,"Permanent Snow and Ice" \
forrest@0	113	,"Barren or Sparsely Vegetated" \
forrest@0	114	,"Unclassified" \
forrest@0	115	,"Woody Savannas (N. Hem.)" \
forrest@0	116	,"Savannas (N. Hem.)" \
forrest@0	117	,"All biome average" \
forrest@0	118	/)
forrest@0	119	nrow = dimsizes(row_header)
forrest@0	120
forrest@0	121	; arrays to be passed to table.
forrest@0	122	value = new ((/nrow, ncol/),string )
forrest@0	123
forrest@0	124	table_length = 0.995
forrest@0	125
forrest@0	126	; Table header
forrest@0	127	ncr1 = (/1,1/) ; 1 row, 1 column
forrest@0	128	xx1 = (/0.005,0.25/) ; Start and end X
forrest@0	129	yy1 = (/0.900,0.995/) ; Start and end Y
forrest@0	130	text1 = table_header_name
forrest@0	131	res1 = True
forrest@0	132	res1@txFontHeightF = 0.03
forrest@0	133	res1@gsFillColor = "CornFlowerBlue"
forrest@0	134
forrest@0	135	; Column header (equally space in x2)
forrest@0	136	ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
forrest@0	137	xx21 = (/xx1(1),0.995/) ; start from end of x1
forrest@0	138	yy21 = (/0.9475,0.995/) ; half of y1
forrest@0	139	text21 = col_header1
forrest@0	140	res21 = True
forrest@0	141	res21@txFontHeightF = 0.015
forrest@0	142	res21@gsFillColor = "Gray"
forrest@0	143
forrest@0	144	ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
forrest@0	145	xx22 = xx21 ; start from end of x1
forrest@0	146	yy22 = (/0.900,0.9475/) ; half of y1
forrest@0	147	text22 = col_header2
forrest@0	148	res22 = True
forrest@0	149	res22@txFontHeightF = 0.015
forrest@0	150	res22@gsFillColor = "Gray"
forrest@0	151
forrest@0	152	; Row header (equally space in y2)
forrest@0	153	ncr3 = (/nrow,1/) ; 20 rows, 1 columns
forrest@0	154	xx3 = xx1 ; same as x1
forrest@0	155	yy3 = (/1.0-table_length,0.900/) ; end at start of y1
forrest@0	156	text3 = row_header
forrest@0	157	res3 = True
forrest@0	158	res3@txFontHeightF = 0.01
forrest@0	159	res3@gsFillColor = "Gray"
forrest@0	160
forrest@0	161	; Main table body
forrest@0	162	ncr4 = (/nrow,ncol/) ; 5 rows, 5 columns
forrest@0	163	xx4 = xx21 ; Start and end x
forrest@0	164	yy4 = yy3 ; Start and end Y
forrest@0	165	text4 = new((/nrow,ncol/),string)
forrest@0	166
forrest@0	167	color_fill4 = new((/nrow,ncol/),string)
forrest@0	168	color_fill4 = "white"
forrest@0	169	color_fill4(nrow-1,:) = "CornFlowerBlue"
forrest@0	170
forrest@0	171	res4 = True ; Set up resource list
forrest@0	172	; res4@gsnDebug = True ; Useful to print NDC row,col values used.
forrest@0	173	res4@txFontHeightF = 0.015
forrest@0	174	res4@gsFillColor = color_fill4
forrest@0	175
forrest@0	176	delete (color_fill4)
forrest@0	177
forrest@0	178	;************************************************
forrest@0	179	; plot global land class: observed
forrest@0	180	;************************************************
forrest@0	181	;global res
forrest@0	182
forrest@0	183	resg = True ; Use plot options
forrest@0	184	resg@cnFillOn = True ; Turn on color fill
forrest@0	185	resg@gsnSpreadColors = True ; use full colormap
forrest@0	186	; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0	187	; resg@lbLabelAutoStride = True
forrest@0	188	resg@cnLinesOn = False ; Turn off contourn lines
forrest@0	189	resg@mpFillOn = False ; Turn off map fill
forrest@0	190
forrest@0	191	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	192	resg@cnMinLevelValF = 1. ; Min level
forrest@0	193	resg@cnMaxLevelValF = 19. ; Max level
forrest@0	194	resg@cnLevelSpacingF = 1. ; interval
forrest@0	195
forrest@0	196	;global contour ob
forrest@0	197	classob@_FillValue = 1.e+36
forrest@0	198	classob = where(classob.eq.0,classob@_FillValue,classob)
forrest@0	199
forrest@0	200	plot_name = "global_class_ob"
forrest@0	201	title = ob_name
forrest@0	202	resg@tiMainString = title
forrest@0	203
forrest@0	204	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	205	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	206
forrest@0	207	plot = gsn_csm_contour_map_ce(wks,classob,resg)
forrest@0	208	frame(wks)
forrest@0	209
forrest@0	210	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	211	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	212	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	213	"rm "+plot_name+"."+plot_type)
forrest@0	214
forrest@0	215	clear (wks)
forrest@0	216	;*******************************************************************
forrest@0	217	; Calculate "nice" bins for binning the data in equally spaced ranges
forrest@0	218	;********************************************************************
forrest@0	219	nclassn = nclass + 1
forrest@0	220	range = fspan(0,nclassn-1,nclassn)
forrest@0	221	; print (range)
forrest@0	222
forrest@0	223	; Use this range information to grab all the values in a
forrest@0	224	; particular range, and then take an average.
forrest@0	225
forrest@0	226	nr = dimsizes(range)
forrest@0	227	nx = nr-1
forrest@0	228	xvalues = new((/2,nx/),float)
forrest@0	229	xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
forrest@0	230	dx = xvalues(0,1) - xvalues(0,0) ; range width
forrest@0	231	dx4 = dx/4 ; 1/4 of the range
forrest@0	232	xvalues(1,:) = xvalues(0,:) - dx/5.
forrest@0	233	;-----------------------------------------------------------------
forrest@0	234	;(A) mean
forrest@0	235	;--------------------------------------------------------------------
forrest@0	236	; get data
forrest@0	237
forrest@0	238	laiob_mean = dim_avg_Wrap(laiob(lat\|:,lon\|:,time\|:))
forrest@0	239	laimod_mean = dim_avg_Wrap(laimod(lat\|:,lon\|:,time\|:))
forrest@0	240
forrest@0	241	DATA11_1D = ndtooned(classob)
forrest@0	242	DATA12_1D = ndtooned(laiob_mean)
forrest@0	243	DATA22_1D = ndtooned(laimod_mean)
forrest@0	244
forrest@0	245	yvalues = new((/2,nx/),float)
forrest@0	246	mn_yvalues = new((/2,nx/),float)
forrest@0	247	mx_yvalues = new((/2,nx/),float)
forrest@0	248
forrest@0	249	do nd=0,1
forrest@0	250
forrest@0	251	; See if we are doing model or observational data.
forrest@0	252
forrest@0	253	if(nd.eq.0) then
forrest@0	254	data_ob = DATA11_1D
forrest@0	255	data_mod = DATA12_1D
forrest@0	256	else
forrest@0	257	data_ob = DATA11_1D
forrest@0	258	data_mod = DATA22_1D
forrest@0	259	end if
forrest@0	260
forrest@0	261	; Loop through each range and check for values.
forrest@0	262
forrest@0	263	do i=0,nr-2
forrest@0	264	if (i.ne.(nr-2)) then
forrest@0	265	; print("")
forrest@0	266	; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0	267	idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))
forrest@0	268	else
forrest@0	269	; print("")
forrest@0	270	; print("In range ["+range(i)+",)")
forrest@0	271	idx = ind(data_ob.ge.range(i))
forrest@0	272	end if
forrest@0	273
forrest@0	274	; Calculate average, and get min and max.
forrest@0	275
forrest@0	276	if(.not.any(ismissing(idx))) then
forrest@0	277	yvalues(nd,i) = avg(data_mod(idx))
forrest@0	278	mn_yvalues(nd,i) = min(data_mod(idx))
forrest@0	279	mx_yvalues(nd,i) = max(data_mod(idx))
forrest@0	280	count = dimsizes(idx)
forrest@0	281	else
forrest@0	282	count = 0
forrest@0	283	yvalues(nd,i) = yvalues@_FillValue
forrest@0	284	mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0	285	mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0	286	end if
forrest@0	287
forrest@0	288	; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0	289	; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0	290
forrest@0	291	; Clean up for next time in loop.
forrest@0	292
forrest@0	293	delete(idx)
forrest@0	294	end do
forrest@0	295	delete(data_ob)
forrest@0	296	delete(data_mod)
forrest@0	297	end do
forrest@0	298	;-----------------------------------------------------------------
forrest@0	299	; compute correlation coef and M score
forrest@0	300
forrest@0	301	u = yvalues(0,:)
forrest@0	302	v = yvalues(1,:)
forrest@0	303
forrest@0	304	good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0	305	uu = u(good)
forrest@0	306	vv = v(good)
forrest@0	307
forrest@0	308	ccrMean = esccr(uu,vv,0)
forrest@0	309	; print (ccrMean)
forrest@0	310
forrest@0	311	; new eq
forrest@0	312	bias = sum(abs(vv-uu)/(vv+uu))
forrest@0	313	Mmean = (1.- (bias/dimsizes(uu)))*5.
forrest@0	314
forrest@0	315	M_lai_mean = sprintf("%.2f", Mmean)
forrest@0	316	system("sed s#M_lai_mean#"+M_lai_mean+"# "+html_name+" > "+html_new+";"+ \
forrest@0	317	"mv -f "+html_new+" "+html_name)
forrest@0	318	print (M_lai_mean)
forrest@0	319
forrest@0	320	do i=0,nrow-2
forrest@0	321	text4(i,0) = sprintf("%5.2f",u(i))
forrest@0	322	text4(i,1) = sprintf("%5.2f",v(i))
forrest@0	323	text4(i,2) = "-"
forrest@0	324	end do
forrest@0	325	text4(nrow-1,0) = sprintf("%5.2f",avg(u))
forrest@0	326	text4(nrow-1,1) = sprintf("%5.2f",avg(v))
forrest@0	327	text4(nrow-1,2) = sprintf("%5.2f",Mmean)
forrest@0	328
forrest@0	329	delete (u)
forrest@0	330	delete (v)
forrest@0	331	delete (uu)
forrest@0	332	delete (vv)
forrest@0	333	;--------------------------------------------------------------------
forrest@0	334	; histogram res
forrest@0	335
forrest@0	336	resm = True
forrest@0	337	resm@gsnMaximize = True
forrest@0	338	resm@gsnDraw = False
forrest@0	339	resm@gsnFrame = False
forrest@0	340	resm@xyMarkLineMode = "Markers"
forrest@0	341	resm@xyMarkerSizeF = 0.014
forrest@0	342	resm@xyMarker = 16
forrest@0	343	resm@xyMarkerColors = (/"Brown","Blue"/)
forrest@0	344	; resm@trYMinF = min(mn_yvalues) - 10.
forrest@0	345	; resm@trYMaxF = max(mx_yvalues) + 10.
forrest@0	346	resm@trYMinF = min(mn_yvalues) - 2
forrest@0	347	resm@trYMaxF = max(mx_yvalues) + 4
forrest@0	348
forrest@0	349	resm@tiYAxisString = "Mean LAI (Leaf Area Index)"
forrest@0	350	resm@tiXAxisString = "Land Cover Type"
forrest@0	351
forrest@0	352	max_bar = new((/2,nx/),graphic)
forrest@0	353	min_bar = new((/2,nx/),graphic)
forrest@0	354	max_cap = new((/2,nx/),graphic)
forrest@0	355	min_cap = new((/2,nx/),graphic)
forrest@0	356
forrest@0	357	lnres = True
forrest@0	358	line_colors = (/"brown","blue"/)
forrest@0	359	;------------------------------------------------------------------
forrest@0	360	; Start the graphics.
forrest@0	361
forrest@0	362	plot_name = "histogram_mean"
forrest@0	363	title = model_name + " vs Observed"
forrest@0	364	resm@tiMainString = title
forrest@0	365
forrest@0	366	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	367	;-----------------------------
forrest@0	368	; Add a boxed legend using the more simple method
forrest@0	369
forrest@0	370	resm@pmLegendDisplayMode = "Always"
forrest@0	371	; resm@pmLegendWidthF = 0.1
forrest@0	372	resm@pmLegendWidthF = 0.08
forrest@0	373	resm@pmLegendHeightF = 0.05
forrest@0	374	resm@pmLegendOrthogonalPosF = -1.17
forrest@0	375	; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0	376	; resm@pmLegendParallelPosF = 0.18
forrest@0	377	resm@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0	378
forrest@0	379	; resm@lgPerimOn = False
forrest@0	380	resm@lgLabelFontHeightF = 0.015
forrest@0	381	resm@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0	382	;-----------------------------
forrest@0	383	tRes = True
forrest@0	384	tRes@txFontHeightF = 0.025
forrest@0	385
forrest@0	386	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"
forrest@0	387
forrest@0	388	gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0	389
forrest@0	390	xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
forrest@0	391	;-------------------------------
forrest@0	392	;Attach the vertical bar and the horizontal cap line
forrest@0	393
forrest@0	394	do nd=0,1
forrest@0	395	lnres@gsLineColor = line_colors(nd)
forrest@0	396	do i=0,nx-1
forrest@0	397
forrest@0	398	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	399	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	400
forrest@0	401	; Attach the vertical bar, both above and below the marker.
forrest@0	402
forrest@0	403	x1 = xvalues(nd,i)
forrest@0	404	y1 = yvalues(nd,i)
forrest@0	405	y2 = mn_yvalues(nd,i)
forrest@0	406	min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	407
forrest@0	408	y2 = mx_yvalues(nd,i)
forrest@0	409	max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	410
forrest@0	411	; Attach the horizontal cap line, both above and below the marker.
forrest@0	412
forrest@0	413	x1 = xvalues(nd,i) - dx4
forrest@0	414	x2 = xvalues(nd,i) + dx4
forrest@0	415	y1 = mn_yvalues(nd,i)
forrest@0	416	min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	417
forrest@0	418	y1 = mx_yvalues(nd,i)
forrest@0	419	max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	420	end if
forrest@0	421	end do
forrest@0	422	end do
forrest@0	423
forrest@0	424	draw(xy)
forrest@0	425	frame(wks)
forrest@0	426
forrest@0	427	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	428	"rm "+plot_name+"."+plot_type)
forrest@0	429
forrest@0	430	clear (wks)
forrest@0	431
forrest@0	432	delete (DATA11_1D)
forrest@0	433	delete (DATA12_1D)
forrest@0	434	delete (DATA22_1D)
forrest@0	435	;delete (range)
forrest@0	436	;delete (xvalues)
forrest@0	437	delete (yvalues)
forrest@0	438	delete (mn_yvalues)
forrest@0	439	delete (mx_yvalues)
forrest@0	440	delete (good)
forrest@0	441	delete (max_bar)
forrest@0	442	delete (min_bar)
forrest@0	443	delete (max_cap)
forrest@0	444	delete (min_cap)
forrest@0	445	;-----------------------------------------------------------------
forrest@0	446	;global res
forrest@0	447
forrest@0	448	resg = True ; Use plot options
forrest@0	449	resg@cnFillOn = True ; Turn on color fill
forrest@0	450	resg@gsnSpreadColors = True ; use full colormap
forrest@0	451	; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0	452	; resg@lbLabelAutoStride = True
forrest@0	453	resg@cnLinesOn = False ; Turn off contourn lines
forrest@0	454	resg@mpFillOn = False ; Turn off map fill
forrest@0	455
forrest@0	456	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	457	resg@cnMinLevelValF = 0. ; Min level
forrest@0	458	resg@cnMaxLevelValF = 10. ; Max level
forrest@0	459	resg@cnLevelSpacingF = 1. ; interval
forrest@0	460
forrest@0	461	;global contour ob
forrest@0	462
forrest@0	463	delta = 0.00001
forrest@0	464	laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean)
forrest@0	465
forrest@0	466	plot_name = "global_mean_ob"
forrest@0	467	title = ob_name
forrest@0	468	resg@tiMainString = title
forrest@0	469
forrest@0	470	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	471	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	472
forrest@0	473	plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg)
forrest@0	474	frame(wks)
forrest@0	475
forrest@0	476	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	477	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	478	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	479	"rm "+plot_name+"."+plot_type)
forrest@0	480
forrest@0	481	clear (wks)
forrest@0	482	;------------------------------------------------------------------------
forrest@0	483	;global contour model
forrest@0	484
forrest@0	485	plot_name = "global_mean_model"
forrest@0	486	title = "Model " + model_name
forrest@0	487	resg@tiMainString = title
forrest@0	488
forrest@0	489	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	490	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	491
forrest@0	492	plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg)
forrest@0	493	frame(wks)
forrest@0	494
forrest@0	495	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	496	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	497	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	498	"rm "+plot_name+"."+plot_type)
forrest@0	499
forrest@0	500	clear (wks)
forrest@0	501	;-----------------------------------------------------------------
forrest@0	502	;(B) max
forrest@0	503	;--------------------------------------------------------------------
forrest@0	504	; get data
forrest@0	505
forrest@0	506	; observed
forrest@0	507	laiob_max = laiob(0,:,:)
forrest@0	508	s = laiob(:,0,0)
forrest@0	509	laiob_max@long_name = "Leaf Area Index Max"
forrest@0	510
forrest@0	511	dsizes_z = dimsizes(laiob)
forrest@0	512	nlat = dsizes_z(1)
forrest@0	513	nlon = dsizes_z(2)
forrest@0	514
forrest@0	515	do j = 0,nlat-1
forrest@0	516	do i = 0,nlon-1
forrest@0	517	s = laiob(:,j,i)
forrest@0	518	laiob_max(j,i) = max(s)
forrest@0	519	end do
forrest@0	520	end do
forrest@0	521
forrest@0	522	; print (min(y)+"/"+max(y))
forrest@0	523	delete (s)
forrest@0	524	delete (dsizes_z)
forrest@0	525	;-------------------------
forrest@0	526	; model
forrest@0	527	laimod_max = laimod(0,:,:)
forrest@0	528	s = laimod(:,0,0)
forrest@0	529	laimod_max@long_name = "Leaf Area Index Max"
forrest@0	530
forrest@0	531	dsizes_z = dimsizes(laimod)
forrest@0	532	nlat = dsizes_z(1)
forrest@0	533	nlon = dsizes_z(2)
forrest@0	534
forrest@0	535	do j = 0,nlat-1
forrest@0	536	do i = 0,nlon-1
forrest@0	537	s = laimod(:,j,i)
forrest@0	538	laimod_max(j,i) = max(s)
forrest@0	539	end do
forrest@0	540	end do
forrest@0	541
forrest@0	542	; print (min(laimod_max)+"/"+max(laimod_max))
forrest@0	543	delete (s)
forrest@0	544	delete (dsizes_z)
forrest@0	545	;------------------------
forrest@0	546	DATA11_1D = ndtooned(classob)
forrest@0	547	DATA12_1D = ndtooned(laiob_max)
forrest@0	548	DATA22_1D = ndtooned(laimod_max)
forrest@0	549
forrest@0	550	yvalues = new((/2,nx/),float)
forrest@0	551	mn_yvalues = new((/2,nx/),float)
forrest@0	552	mx_yvalues = new((/2,nx/),float)
forrest@0	553
forrest@0	554	do nd=0,1
forrest@0	555
forrest@0	556	; See if we are doing model or observational data.
forrest@0	557
forrest@0	558	if(nd.eq.0) then
forrest@0	559	data_ob = DATA11_1D
forrest@0	560	data_mod = DATA12_1D
forrest@0	561	else
forrest@0	562	data_ob = DATA11_1D
forrest@0	563	data_mod = DATA22_1D
forrest@0	564	end if
forrest@0	565
forrest@0	566	; Loop through each range and check for values.
forrest@0	567
forrest@0	568	do i=0,nr-2
forrest@0	569	if (i.ne.(nr-2)) then
forrest@0	570	; print("")
forrest@0	571	; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0	572	idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
forrest@0	573	else
forrest@0	574	; print("")
forrest@0	575	; print("In range ["+range(i)+",)")
forrest@0	576	idx = ind(range(i).le.data_ob)
forrest@0	577	end if
forrest@0	578
forrest@0	579	; Calculate average, and get min and max.
forrest@0	580
forrest@0	581	if(.not.any(ismissing(idx))) then
forrest@0	582	yvalues(nd,i) = avg(data_mod(idx))
forrest@0	583	mn_yvalues(nd,i) = min(data_mod(idx))
forrest@0	584	mx_yvalues(nd,i) = max(data_mod(idx))
forrest@0	585	count = dimsizes(idx)
forrest@0	586	else
forrest@0	587	count = 0
forrest@0	588	yvalues(nd,i) = yvalues@_FillValue
forrest@0	589	mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0	590	mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0	591	end if
forrest@0	592
forrest@0	593	; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0	594	; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0	595
forrest@0	596	; Clean up for next time in loop.
forrest@0	597
forrest@0	598	delete(idx)
forrest@0	599	end do
forrest@0	600	delete(data_ob)
forrest@0	601	delete(data_mod)
forrest@0	602	end do
forrest@0	603	;-----------------------------------------------------------------
forrest@0	604	; compute correlation coef and M score
forrest@0	605
forrest@0	606	u = yvalues(0,:)
forrest@0	607	v = yvalues(1,:)
forrest@0	608
forrest@0	609	good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0	610	uu = u(good)
forrest@0	611	vv = v(good)
forrest@0	612
forrest@0	613	ccrMax = esccr(uu,vv,0)
forrest@0	614	; print (ccrMax)
forrest@0	615
forrest@0	616	; new eq
forrest@0	617	bias = sum(abs(vv-uu)/(vv+uu))
forrest@0	618	Mmax = (1.- (bias/dimsizes(uu)))*5.
forrest@0	619
forrest@0	620	M_lai_max = sprintf("%.2f", Mmax)
forrest@0	621	system("sed s#M_lai_max#"+M_lai_max+"# "+html_name+" > "+html_new+";"+ \
forrest@0	622	"mv -f "+html_new+" "+html_name)
forrest@0	623	print (M_lai_max)
forrest@0	624
forrest@0	625	do i=0,nrow-2
forrest@0	626	text4(i,3) = sprintf("%5.2f",u(i))
forrest@0	627	text4(i,4) = sprintf("%5.2f",v(i))
forrest@0	628	text4(i,5) = "-"
forrest@0	629	end do
forrest@0	630	text4(nrow-1,3) = sprintf("%5.2f",avg(u))
forrest@0	631	text4(nrow-1,4) = sprintf("%5.2f",avg(v))
forrest@0	632	text4(nrow-1,5) = sprintf("%5.2f",Mmax)
forrest@0	633
forrest@0	634	delete (u)
forrest@0	635	delete (v)
forrest@0	636	delete (uu)
forrest@0	637	delete (vv)
forrest@0	638	;--------------------------------------------------------------------
forrest@0	639	; histogram res
forrest@0	640
forrest@0	641	resm = True
forrest@0	642	resm@gsnMaximize = True
forrest@0	643	resm@gsnDraw = False
forrest@0	644	resm@gsnFrame = False
forrest@0	645	resm@xyMarkLineMode = "Markers"
forrest@0	646	resm@xyMarkerSizeF = 0.014
forrest@0	647	resm@xyMarker = 16
forrest@0	648	resm@xyMarkerColors = (/"Brown","Blue"/)
forrest@0	649	; resm@trYMinF = min(mn_yvalues) - 10.
forrest@0	650	; resm@trYMaxF = max(mx_yvalues) + 10.
forrest@0	651	resm@trYMinF = min(mn_yvalues) - 2
forrest@0	652	resm@trYMaxF = max(mx_yvalues) + 4
forrest@0	653
forrest@0	654	resm@tiYAxisString = "Max LAI (Leaf Area Index)"
forrest@0	655	resm@tiXAxisString = "Land Cover Type"
forrest@0	656
forrest@0	657	max_bar = new((/2,nx/),graphic)
forrest@0	658	min_bar = new((/2,nx/),graphic)
forrest@0	659	max_cap = new((/2,nx/),graphic)
forrest@0	660	min_cap = new((/2,nx/),graphic)
forrest@0	661
forrest@0	662	lnres = True
forrest@0	663	line_colors = (/"brown","blue"/)
forrest@0	664	;------------------------------------------------------------------
forrest@0	665	; Start the graphics.
forrest@0	666
forrest@0	667	plot_name = "histogram_max"
forrest@0	668	title = model_name + " vs Observed"
forrest@0	669	resm@tiMainString = title
forrest@0	670
forrest@0	671	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	672	;-----------------------------
forrest@0	673	; Add a boxed legend using the more simple method
forrest@0	674
forrest@0	675	resm@pmLegendDisplayMode = "Always"
forrest@0	676	; resm@pmLegendWidthF = 0.1
forrest@0	677	resm@pmLegendWidthF = 0.08
forrest@0	678	resm@pmLegendHeightF = 0.05
forrest@0	679	resm@pmLegendOrthogonalPosF = -1.17
forrest@0	680	; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0	681	; resm@pmLegendParallelPosF = 0.18
forrest@0	682	resm@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0	683
forrest@0	684	; resm@lgPerimOn = False
forrest@0	685	resm@lgLabelFontHeightF = 0.015
forrest@0	686	resm@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0	687	;-----------------------------
forrest@0	688	tRes = True
forrest@0	689	tRes@txFontHeightF = 0.025
forrest@0	690
forrest@0	691	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
forrest@0	692
forrest@0	693	gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0	694
forrest@0	695	xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
forrest@0	696	;-------------------------------
forrest@0	697	;Attach the vertical bar and the horizontal cap line
forrest@0	698
forrest@0	699	do nd=0,1
forrest@0	700	lnres@gsLineColor = line_colors(nd)
forrest@0	701	do i=0,nx-1
forrest@0	702
forrest@0	703	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	704	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	705
forrest@0	706	; Attach the vertical bar, both above and below the marker.
forrest@0	707
forrest@0	708	x1 = xvalues(nd,i)
forrest@0	709	y1 = yvalues(nd,i)
forrest@0	710	y2 = mn_yvalues(nd,i)
forrest@0	711	min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	712
forrest@0	713	y2 = mx_yvalues(nd,i)
forrest@0	714	max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	715
forrest@0	716	; Attach the horizontal cap line, both above and below the marker.
forrest@0	717
forrest@0	718	x1 = xvalues(nd,i) - dx4
forrest@0	719	x2 = xvalues(nd,i) + dx4
forrest@0	720	y1 = mn_yvalues(nd,i)
forrest@0	721	min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	722
forrest@0	723	y1 = mx_yvalues(nd,i)
forrest@0	724	max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	725	end if
forrest@0	726	end do
forrest@0	727	end do
forrest@0	728
forrest@0	729	draw(xy)
forrest@0	730	frame(wks)
forrest@0	731
forrest@0	732	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	733	"rm "+plot_name+"."+plot_type)
forrest@0	734
forrest@0	735	clear (wks)
forrest@0	736
forrest@0	737	delete (DATA11_1D)
forrest@0	738	delete (DATA12_1D)
forrest@0	739	delete (DATA22_1D)
forrest@0	740	;delete (range)
forrest@0	741	;delete (xvalues)
forrest@0	742	delete (yvalues)
forrest@0	743	delete (mn_yvalues)
forrest@0	744	delete (mx_yvalues)
forrest@0	745	delete (good)
forrest@0	746	delete (max_bar)
forrest@0	747	delete (min_bar)
forrest@0	748	delete (max_cap)
forrest@0	749	delete (min_cap)
forrest@0	750	;-----------------------------------------------------------------
forrest@0	751	;global res
forrest@0	752
forrest@0	753	resg = True ; Use plot options
forrest@0	754	resg@cnFillOn = True ; Turn on color fill
forrest@0	755	resg@gsnSpreadColors = True ; use full colormap
forrest@0	756	; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0	757	; resg@lbLabelAutoStride = True
forrest@0	758	resg@cnLinesOn = False ; Turn off contourn lines
forrest@0	759	resg@mpFillOn = False ; Turn off map fill
forrest@0	760
forrest@0	761	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	762	resg@cnMinLevelValF = 0. ; Min level
forrest@0	763	resg@cnMaxLevelValF = 10. ; Max level
forrest@0	764	resg@cnLevelSpacingF = 1. ; interval
forrest@0	765
forrest@0	766	;global contour ob
forrest@0	767
forrest@0	768	delta = 0.00001
forrest@0	769	laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)
forrest@0	770
forrest@0	771	plot_name = "global_max_ob"
forrest@0	772	title = ob_name
forrest@0	773	resg@tiMainString = title
forrest@0	774
forrest@0	775	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	776	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	777
forrest@0	778	plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg)
forrest@0	779	frame(wks)
forrest@0	780
forrest@0	781	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	782	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	783	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	784	"rm "+plot_name+"."+plot_type)
forrest@0	785
forrest@0	786	clear (wks)
forrest@0	787	;------------------------------------------------------------------------
forrest@0	788	;global contour model
forrest@0	789
forrest@0	790	plot_name = "global_max_model"
forrest@0	791	title = "Model " + model_name
forrest@0	792	resg@tiMainString = title
forrest@0	793
forrest@0	794	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	795	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	796
forrest@0	797	plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg)
forrest@0	798	frame(wks)
forrest@0	799
forrest@0	800	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	801	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	802	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	803	"rm "+plot_name+"."+plot_type)
forrest@0	804
forrest@0	805	clear (wks)
forrest@0	806	;------------------------------------------------------------------------
forrest@0	807	;(C) phase
forrest@0	808	;--------------------------------------------------------------------
forrest@0	809	; get data
forrest@0	810
forrest@0	811	; observed
forrest@0	812	laiob_phase = laiob(0,:,:)
forrest@0	813	s = laiob(:,0,0)
forrest@0	814	laiob_phase@long_name = "Leaf Area Index Max Month"
forrest@0	815
forrest@0	816	dsizes_z = dimsizes(laiob)
forrest@0	817	nlat = dsizes_z(1)
forrest@0	818	nlon = dsizes_z(2)
forrest@0	819
forrest@0	820	do j = 0,nlat-1
forrest@0	821	do i = 0,nlon-1
forrest@0	822	s = laiob(:,j,i)
forrest@0	823	laiob_phase(j,i) = maxind(s) + 1
forrest@0	824	end do
forrest@0	825	end do
forrest@0	826
forrest@0	827	; print (min(laiob_phase)+"/"+max(laiob_phase))
forrest@0	828	delete (s)
forrest@0	829	delete (dsizes_z)
forrest@0	830	;-------------------------
forrest@0	831	; model
forrest@0	832	laimod_phase = laimod(0,:,:)
forrest@0	833	s = laimod(:,0,0)
forrest@0	834	laimod_phase@long_name = "Leaf Area Index Max Month"
forrest@0	835
forrest@0	836	dsizes_z = dimsizes(laimod)
forrest@0	837	nlat = dsizes_z(1)
forrest@0	838	nlon = dsizes_z(2)
forrest@0	839
forrest@0	840	do j = 0,nlat-1
forrest@0	841	do i = 0,nlon-1
forrest@0	842	s = laimod(:,j,i)
forrest@0	843	laimod_phase(j,i) = maxind(s) + 1
forrest@0	844	end do
forrest@0	845	end do
forrest@0	846
forrest@0	847	; print (min(laimod_phase)+"/"+max(laimod_phase))
forrest@0	848	delete (s)
forrest@0	849	delete (dsizes_z)
forrest@0	850	;------------------------
forrest@0	851	DATA11_1D = ndtooned(classob)
forrest@0	852	DATA12_1D = ndtooned(laiob_phase)
forrest@0	853	DATA22_1D = ndtooned(laimod_phase)
forrest@0	854
forrest@0	855	yvalues = new((/2,nx/),float)
forrest@0	856	mn_yvalues = new((/2,nx/),float)
forrest@0	857	mx_yvalues = new((/2,nx/),float)
forrest@0	858
forrest@0	859	do nd=0,1
forrest@0	860
forrest@0	861	; See if we are doing model or observational data.
forrest@0	862
forrest@0	863	if(nd.eq.0) then
forrest@0	864	data_ob = DATA11_1D
forrest@0	865	data_mod = DATA12_1D
forrest@0	866	else
forrest@0	867	data_ob = DATA11_1D
forrest@0	868	data_mod = DATA22_1D
forrest@0	869	end if
forrest@0	870
forrest@0	871	; Loop through each range and check for values.
forrest@0	872
forrest@0	873	do i=0,nr-2
forrest@0	874	if (i.ne.(nr-2)) then
forrest@0	875	; print("")
forrest@0	876	; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0	877	idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
forrest@0	878	else
forrest@0	879	; print("")
forrest@0	880	; print("In range ["+range(i)+",)")
forrest@0	881	idx = ind(range(i).le.data_ob)
forrest@0	882	end if
forrest@0	883
forrest@0	884	; Calculate average, and get min and max.
forrest@0	885
forrest@0	886	if(.not.any(ismissing(idx))) then
forrest@0	887	yvalues(nd,i) = avg(data_mod(idx))
forrest@0	888	mn_yvalues(nd,i) = min(data_mod(idx))
forrest@0	889	mx_yvalues(nd,i) = max(data_mod(idx))
forrest@0	890	count = dimsizes(idx)
forrest@0	891	else
forrest@0	892	count = 0
forrest@0	893	yvalues(nd,i) = yvalues@_FillValue
forrest@0	894	mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0	895	mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0	896	end if
forrest@0	897
forrest@0	898	; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0	899	; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0	900
forrest@0	901	; Clean up for next time in loop.
forrest@0	902
forrest@0	903	delete(idx)
forrest@0	904	end do
forrest@0	905	delete(data_ob)
forrest@0	906	delete(data_mod)
forrest@0	907	end do
forrest@0	908	;-----------------------------------------------------------------
forrest@0	909	; compute correlation coef and M score
forrest@0	910
forrest@0	911	u = yvalues(0,:)
forrest@0	912	v = yvalues(1,:)
forrest@0	913
forrest@0	914	good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0	915	uu = u(good)
forrest@0	916	vv = v(good)
forrest@0	917
forrest@0	918	ccrPhase = esccr(uu,vv,0)
forrest@0	919	; print (ccrPhase)
forrest@0	920
forrest@0	921	; old eq
forrest@0	922	; bias = abs(avg(vv)-avg(uu))
forrest@0	923	; new eq
forrest@0	924	bias = avg(abs(vv-uu))
forrest@0	925
forrest@0	926	bias = where((bias.gt. 6.),12.-bias,bias)
forrest@0	927	Mphase = ((6. - bias)/6.)*5.
forrest@0	928
forrest@0	929	M_lai_phase = sprintf("%.2f", Mphase)
forrest@0	930	system("sed s#M_lai_phase#"+M_lai_phase+"# "+html_name+" > "+html_new+";"+ \
forrest@0	931	"mv -f "+html_new+" "+html_name)
forrest@0	932	print (M_lai_phase)
forrest@0	933
forrest@0	934	do i=0,nrow-2
forrest@0	935	text4(i,6) = sprintf("%5.2f",u(i))
forrest@0	936	text4(i,7) = sprintf("%5.2f",v(i))
forrest@0	937	text4(i,8) = "-"
forrest@0	938	end do
forrest@0	939	text4(nrow-1,6) = sprintf("%5.2f",avg(u))
forrest@0	940	text4(nrow-1,7) = sprintf("%5.2f",avg(v))
forrest@0	941	text4(nrow-1,8) = sprintf("%5.2f",Mphase)
forrest@0	942
forrest@0	943	delete (u)
forrest@0	944	delete (v)
forrest@0	945	delete (uu)
forrest@0	946	delete (vv)
forrest@0	947	;--------------------------------------------------------------------
forrest@0	948	; histogram res
forrest@0	949
forrest@0	950	resm = True
forrest@0	951	resm@gsnMaximize = True
forrest@0	952	resm@gsnDraw = False
forrest@0	953	resm@gsnFrame = False
forrest@0	954	resm@xyMarkLineMode = "Markers"
forrest@0	955	resm@xyMarkerSizeF = 0.014
forrest@0	956	resm@xyMarker = 16
forrest@0	957	resm@xyMarkerColors = (/"Brown","Blue"/)
forrest@0	958	; resm@trYMinF = min(mn_yvalues) - 10.
forrest@0	959	; resm@trYMaxF = max(mx_yvalues) + 10.
forrest@0	960	resm@trYMinF = min(mn_yvalues) - 2
forrest@0	961	resm@trYMaxF = max(mx_yvalues) + 4
forrest@0	962
forrest@0	963	resm@tiYAxisString = "Max LAI (Leaf Area Index) Month"
forrest@0	964	resm@tiXAxisString = "Land Cover Type"
forrest@0	965
forrest@0	966	max_bar = new((/2,nx/),graphic)
forrest@0	967	min_bar = new((/2,nx/),graphic)
forrest@0	968	max_cap = new((/2,nx/),graphic)
forrest@0	969	min_cap = new((/2,nx/),graphic)
forrest@0	970
forrest@0	971	lnres = True
forrest@0	972	line_colors = (/"brown","blue"/)
forrest@0	973	;------------------------------------------------------------------
forrest@0	974	; Start the graphics.
forrest@0	975
forrest@0	976	plot_name = "histogram_phase"
forrest@0	977	title = model_name + " vs Observed"
forrest@0	978	resm@tiMainString = title
forrest@0	979
forrest@0	980	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	981	;-----------------------------
forrest@0	982	; Add a boxed legend using the more simple method
forrest@0	983
forrest@0	984	resm@pmLegendDisplayMode = "Always"
forrest@0	985	; resm@pmLegendWidthF = 0.1
forrest@0	986	resm@pmLegendWidthF = 0.08
forrest@0	987	resm@pmLegendHeightF = 0.05
forrest@0	988	resm@pmLegendOrthogonalPosF = -1.17
forrest@0	989	; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0	990	; resm@pmLegendParallelPosF = 0.18
forrest@0	991	resm@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0	992
forrest@0	993	; resm@lgPerimOn = False
forrest@0	994	resm@lgLabelFontHeightF = 0.015
forrest@0	995	resm@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0	996	;-----------------------------
forrest@0	997	tRes = True
forrest@0	998	tRes@txFontHeightF = 0.025
forrest@0	999
forrest@0	1000	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
forrest@0	1001
forrest@0	1002	gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0	1003
forrest@0	1004	xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
forrest@0	1005	;-------------------------------
forrest@0	1006	;Attach the vertical bar and the horizontal cap line
forrest@0	1007
forrest@0	1008	do nd=0,1
forrest@0	1009	lnres@gsLineColor = line_colors(nd)
forrest@0	1010	do i=0,nx-1
forrest@0	1011
forrest@0	1012	if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0	1013	.not.ismissing(mx_yvalues(nd,i))) then
forrest@0	1014
forrest@0	1015	; Attach the vertical bar, both above and below the marker.
forrest@0	1016
forrest@0	1017	x1 = xvalues(nd,i)
forrest@0	1018	y1 = yvalues(nd,i)
forrest@0	1019	y2 = mn_yvalues(nd,i)
forrest@0	1020	min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1021
forrest@0	1022	y2 = mx_yvalues(nd,i)
forrest@0	1023	max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0	1024
forrest@0	1025	; Attach the horizontal cap line, both above and below the marker.
forrest@0	1026
forrest@0	1027	x1 = xvalues(nd,i) - dx4
forrest@0	1028	x2 = xvalues(nd,i) + dx4
forrest@0	1029	y1 = mn_yvalues(nd,i)
forrest@0	1030	min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1031
forrest@0	1032	y1 = mx_yvalues(nd,i)
forrest@0	1033	max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0	1034	end if
forrest@0	1035	end do
forrest@0	1036	end do
forrest@0	1037
forrest@0	1038	draw(xy)
forrest@0	1039	frame(wks)
forrest@0	1040
forrest@0	1041	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1042	"rm "+plot_name+"."+plot_type)
forrest@0	1043
forrest@0	1044	clear (wks)
forrest@0	1045
forrest@0	1046	delete (DATA11_1D)
forrest@0	1047	delete (DATA12_1D)
forrest@0	1048	delete (DATA22_1D)
forrest@0	1049	;delete (range)
forrest@0	1050	;delete (xvalues)
forrest@0	1051	delete (yvalues)
forrest@0	1052	delete (mn_yvalues)
forrest@0	1053	delete (mx_yvalues)
forrest@0	1054	delete (good)
forrest@0	1055	delete (max_bar)
forrest@0	1056	delete (min_bar)
forrest@0	1057	delete (max_cap)
forrest@0	1058	delete (min_cap)
forrest@0	1059	;-----------------------------------------------------------------
forrest@0	1060	;global res
forrest@0	1061
forrest@0	1062	resg = True ; Use plot options
forrest@0	1063	resg@cnFillOn = True ; Turn on color fill
forrest@0	1064	resg@gsnSpreadColors = True ; use full colormap
forrest@0	1065	; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0	1066	; resg@lbLabelAutoStride = True
forrest@0	1067	resg@cnLinesOn = False ; Turn off contourn lines
forrest@0	1068	resg@mpFillOn = False ; Turn off map fill
forrest@0	1069
forrest@0	1070	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1071	resg@cnMinLevelValF = 1. ; Min level
forrest@0	1072	resg@cnMaxLevelValF = 12. ; Max level
forrest@0	1073	resg@cnLevelSpacingF = 1. ; interval
forrest@0	1074
forrest@0	1075	;global contour ob
forrest@0	1076
forrest@0	1077	delta = 0.00001
forrest@0	1078	laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)
forrest@0	1079
forrest@0	1080	plot_name = "global_phase_ob"
forrest@0	1081	title = ob_name
forrest@0	1082	resg@tiMainString = title
forrest@0	1083
forrest@0	1084	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1085	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1086
forrest@0	1087	plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
forrest@0	1088	frame(wks)
forrest@0	1089
forrest@0	1090	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1091	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1092	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	1093	"rm "+plot_name+"."+plot_type)
forrest@0	1094
forrest@0	1095	clear (wks)
forrest@0	1096	;------------------------------------------------------------------------
forrest@0	1097	;global contour model
forrest@0	1098
forrest@0	1099	plot_name = "global_phase_model"
forrest@0	1100	title = "Model " + model_name
forrest@0	1101	resg@tiMainString = title
forrest@0	1102
forrest@0	1103	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	1104	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1105
forrest@0	1106	delete (plot)
forrest@0	1107	plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
forrest@0	1108	frame(wks)
forrest@0	1109
forrest@0	1110	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1111	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1112	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	1113	"rm "+plot_name+"."+plot_type)
forrest@0	1114
forrest@0	1115	clear (wks)
forrest@0	1116	;-----------------------------------------------------------------
forrest@0	1117	;(D) grow day
forrest@0	1118	;--------------------------------------------------------------------
forrest@0	1119	; get data
forrest@0	1120
forrest@0	1121	day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
forrest@0	1122
forrest@0	1123	; observed
forrest@0	1124	laiob_grow = laiob(0,:,:)
forrest@0	1125	laiob_grow@long_name = "Days of Growing Season"
forrest@0	1126
forrest@0	1127	dsizes_z = dimsizes(laiob)
forrest@0	1128	ntime = dsizes_z(0)
forrest@0	1129	nlat = dsizes_z(1)
forrest@0	1130	nlon = dsizes_z(2)
forrest@0	1131
forrest@0	1132	do j = 0,nlat-1
forrest@0	1133	do i = 0,nlon-1
forrest@0	1134	nday = 0.
forrest@0	1135	do k = 0,ntime-1
forrest@0	1136	if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then
forrest@0	1137	nday = nday + day_of_data(k)
forrest@0	1138	end if
forrest@0	1139	end do
forrest@0	1140
forrest@0	1141	laiob_grow(j,i) = nday
forrest@0	1142	end do
forrest@0	1143	end do
forrest@0	1144
forrest@0	1145	; print (min(laiob_grow)+"/"+max(laiob_grow))
forrest@0	1146	;-------------------------
forrest@0	1147	; model
forrest@0	1148	laimod_grow = laimod(0,:,:)
forrest@0	1149	laimod_grow@long_name = "Days of Growing Season"
forrest@0	1150
forrest@0	1151	dsizes_z = dimsizes(laimod)
forrest@0	1152	ntime = dsizes_z(0)
forrest@0	1153	nlat = dsizes_z(1)
forrest@0	1154	nlon = dsizes_z(2)
forrest@0	1155
forrest@0	1156	do j = 0,nlat-1
forrest@0	1157	do i = 0,nlon-1
forrest@0	1158	nday = 0.
forrest@0	1159	do k = 0,ntime-1
forrest@0	1160	if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then
forrest@0	1161	nday = nday + day_of_data(k)
forrest@0	1162	end if
forrest@0	1163	end do
forrest@0	1164
forrest@0	1165	laimod_grow(j,i) = nday
forrest@0	1166	end do
forrest@0	1167	end do
forrest@0	1168
forrest@0	1169	; print (min(laimod_grow)+"/"+max(laimod_grow))
forrest@0	1170	;------------------------
forrest@0	1171	DATA11_1D = ndtooned(classob)
forrest@0	1172	DATA12_1D = ndtooned(laiob_grow)
forrest@0	1173	DATA22_1D = ndtooned(laimod_grow)
forrest@0	1174
forrest@0	1175	yvalues = new((/2,nx/),float)
forrest@0	1176	mn_yvalues = new((/2,nx/),float)
forrest@0	1177	mx_yvalues = new((/2,nx/),float)
forrest@0	1178
forrest@0	1179	do nd=0,1
forrest@0	1180
forrest@0	1181	; See if we are doing model or observational data.
forrest@0	1182
forrest@0	1183	if(nd.eq.0) then
forrest@0	1184	data_ob = DATA11_1D
forrest@0	1185	data_mod = DATA12_1D
forrest@0	1186	else
forrest@0	1187	data_ob = DATA11_1D
forrest@0	1188	data_mod = DATA22_1D
forrest@0	1189	end if
forrest@0	1190
forrest@0	1191	; Loop through each range and check for values.
forrest@0	1192
forrest@0	1193	do i=0,nr-2
forrest@0	1194	if (i.ne.(nr-2)) then
forrest@0	1195	; print("")
forrest@0	1196	; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0	1197	idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
forrest@0	1198	else
forrest@0	1199	; print("")
forrest@0	1200	; print("In range ["+range(i)+",)")
forrest@0	1201	idx = ind(range(i).le.data_ob)
forrest@0	1202	end if
forrest@0	1203
forrest@0	1204	; Calculate average, and get min and max.
forrest@0	1205
forrest@0	1206	if(.not.any(ismissing(idx))) then
forrest@0	1207	yvalues(nd,i) = avg(data_mod(idx))
forrest@0	1208	mn_yvalues(nd,i) = min(data_mod(idx))
forrest@0	1209	mx_yvalues(nd,i) = max(data_mod(idx))
forrest@0	1210	count = dimsizes(idx)
forrest@0	1211	else
forrest@0	1212	count = 0
forrest@0	1213	yvalues(nd,i) = yvalues@_FillValue
forrest@0	1214	mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0	1215	mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0	1216	end if
forrest@0	1217
forrest@0	1218	; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0	1219	; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0	1220
forrest@0	1221	; Clean up for next time in loop.
forrest@0	1222
forrest@0	1223	delete(idx)
forrest@0	1224	end do
forrest@0	1225	delete(data_ob)
forrest@0	1226	delete(data_mod)
forrest@0	1227	end do
forrest@0	1228	;-----------------------------------------------------------------
forrest@0	1229	; compute correlation coef and M score
forrest@0	1230
forrest@0	1231	u = yvalues(0,:)
forrest@0	1232	v = yvalues(1,:)
forrest@0	1233
forrest@0	1234	good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0	1235	uu = u(good)
forrest@0	1236	vv = v(good)
forrest@0	1237
forrest@0	1238	ccrGrow = esccr(uu,vv,0)
forrest@0	1239	; print (ccrGrow)
forrest@0	1240
forrest@0	1241	; new eq
forrest@0	1242	bias = sum(abs(vv-uu)/(vv+uu))
forrest@0	1243	Mgrow = (1.- (bias/dimsizes(uu)))*5.
forrest@0	1244
forrest@0	1245	M_lai_grow = sprintf("%.2f", Mgrow)
forrest@0	1246	system("sed s#M_lai_grow#"+M_lai_grow+"# "+html_name+" > "+html_new+";"+ \
forrest@0	1247	"mv -f "+html_new+" "+html_name)
forrest@0	1248	print (M_lai_grow)
forrest@0	1249
forrest@0	1250	do i=0,nrow-2
forrest@0	1251	text4(i,9) = sprintf("%5.2f",u(i))
forrest@0	1252	text4(i,10) = sprintf("%5.2f",v(i))
forrest@0	1253	text4(i,11) = "-"
forrest@0	1254	end do
forrest@0	1255	text4(nrow-1,9) = sprintf("%5.2f",avg(u))
forrest@0	1256	text4(nrow-1,10) = sprintf("%5.2f",avg(v))
forrest@0	1257	text4(nrow-1,11) = sprintf("%5.2f",Mgrow)
forrest@0	1258
forrest@0	1259	delete (u)
forrest@0	1260	delete (v)
forrest@0	1261	delete (uu)
forrest@0	1262	delete (vv)
forrest@0	1263	;--------------------------------------------------------------------
forrest@0	1264	; histogram res
forrest@0	1265
forrest@0	1266	resm = True
forrest@0	1267	resm@gsnMaximize = True
forrest@0	1268	resm@gsnDraw = False
forrest@0	1269	resm@gsnFrame = False
forrest@0	1270	resm@xyMarkLineMode = "Markers"
forrest@0	1271	resm@xyMarkerSizeF = 0.014
forrest@0	1272	resm@xyMarker = 16
forrest@0	1273	resm@xyMarkerColors = (/"Brown","Blue"/)
forrest@0	1274	; resm@trYMinF = min(mn_yvalues) - 10.
forrest@0	1275	; resm@trYMaxF = max(mx_yvalues) + 10.
forrest@0	1276	resm@trYMinF = min(mn_yvalues) - 2
forrest@0	1277	resm@trYMaxF = max(mx_yvalues) + 4
forrest@0	1278
forrest@0	1279	resm@tiYAxisString = "Days of Growing season"
forrest@0	1280	resm@tiXAxisString = "Land Cover Type"
forrest@0	1281
forrest@0	1282	max_bar = new((/2,nx/),graphic)
forrest@0	1283	min_bar = new((/2,nx/),graphic)
forrest@0	1284	max_cap = new((/2,nx/),graphic)
forrest@0	1285	min_cap = new((/2,nx/),graphic)
forrest@0	1286
forrest@0	1287	lnres = True
forrest@0	1288	line_colors = (/"brown","blue"/)
forrest@0	1289	;------------------------------------------------------------------
forrest@0	1290	; Start the graphics.
forrest@0	1291
forrest@0	1292	plot_name = "histogram_grow"
forrest@0	1293	title = model_name + " vs Observed"
forrest@0	1294	resm@tiMainString = title
forrest@0	1295
forrest@0	1296	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	1297	;-----------------------------
forrest@0	1298	; Add a boxed legend using the more simple method
forrest@0	1299
forrest@0	1300	resm@pmLegendDisplayMode = "Always"
forrest@0	1301	; resm@pmLegendWidthF = 0.1
forrest@0	1302	resm@pmLegendWidthF = 0.08
forrest@0	1303	resm@pmLegendHeightF = 0.05
forrest@0	1304	resm@pmLegendOrthogonalPosF = -1.17
forrest@0	1305	; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0	1306	; resm@pmLegendParallelPosF = 0.18
forrest@0	1307	resm@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0	1308
forrest@0	1309	; resm@lgPerimOn = False
forrest@0	1310	resm@lgLabelFontHeightF = 0.015
forrest@0	1311	resm@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0	1312	;-----------------------------
forrest@0	1313	tRes = True
forrest@0	1314	tRes@txFontHeightF = 0.025
forrest@0	1315
forrest@0	1316	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
forrest@0	1317
forrest@0	1318	gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0	1319
forrest@0	1320	xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
forrest@0	1321	;-------------------------------
forrest@0	1322	;Attach the vertical bar and the horizontal cap line
forrest@0	1323
forrest@0	1324	do nd=0,1
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	"rm "+plot_name+"."+plot_type)
forrest@0	1359
forrest@0	1360	clear (wks)
forrest@0	1361
forrest@0	1362	delete (DATA11_1D)
forrest@0	1363	delete (DATA12_1D)
forrest@0	1364	delete (DATA22_1D)
forrest@0	1365	;delete (range)
forrest@0	1366	;delete (xvalues)
forrest@0	1367	delete (yvalues)
forrest@0	1368	delete (mn_yvalues)
forrest@0	1369	delete (mx_yvalues)
forrest@0	1370	delete (good)
forrest@0	1371	delete (max_bar)
forrest@0	1372	delete (min_bar)
forrest@0	1373	delete (max_cap)
forrest@0	1374	delete (min_cap)
forrest@0	1375	;-----------------------------------------------------------------
forrest@0	1376	;global res
forrest@0	1377
forrest@0	1378	resg = True ; Use plot options
forrest@0	1379	resg@cnFillOn = True ; Turn on color fill
forrest@0	1380	resg@gsnSpreadColors = True ; use full colormap
forrest@0	1381	; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0	1382	; resg@lbLabelAutoStride = True
forrest@0	1383	resg@cnLinesOn = False ; Turn off contourn lines
forrest@0	1384	resg@mpFillOn = False ; Turn off map fill
forrest@0	1385
forrest@0	1386	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1387	resg@cnMinLevelValF = 60. ; Min level
forrest@0	1388	resg@cnMaxLevelValF = 360. ; Max level
forrest@0	1389	resg@cnLevelSpacingF = 20. ; interval
forrest@0	1390
forrest@0	1391	;global contour ob
forrest@0	1392
forrest@0	1393	laiob_grow@_FillValue = 1.e+36
forrest@0	1394	laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow)
forrest@0	1395
forrest@0	1396	plot_name = "global_grow_ob"
forrest@0	1397	title = ob_name
forrest@0	1398	resg@tiMainString = title
forrest@0	1399
forrest@0	1400	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1401	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1402
forrest@0	1403	plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
forrest@0	1404	frame(wks)
forrest@0	1405
forrest@0	1406	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1407	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1408	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	1409	"rm "+plot_name+"."+plot_type)
forrest@0	1410
forrest@0	1411	clear (wks)
forrest@0	1412	;------------------------------------------------------------------------
forrest@0	1413	;global contour model
forrest@0	1414
forrest@0	1415	laimod_grow@_FillValue = 1.e+36
forrest@0	1416	laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow)
forrest@0	1417
forrest@0	1418	plot_name = "global_grow_model"
forrest@0	1419	title = "Model " + model_name
forrest@0	1420	resg@tiMainString = title
forrest@0	1421
forrest@0	1422	wks = gsn_open_wks (plot_type,plot_name)
forrest@0	1423	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1424
forrest@0	1425	delete (plot)
forrest@0	1426	plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
forrest@0	1427	frame(wks)
forrest@0	1428
forrest@0	1429	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1430	"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1431	"rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0	1432	"rm "+plot_name+"."+plot_type)
forrest@0	1433
forrest@0	1434	clear (wks)
forrest@0	1435	;------------------------------------------------------------------------
forrest@0	1436	;global contour model vs ob
forrest@0	1437
forrest@0	1438	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1439	resg@cnMinLevelValF = 0. ; Min level
forrest@0	1440	resg@cnMaxLevelValF = 10. ; Max level
forrest@0	1441	resg@cnLevelSpacingF = 1. ; interval
forrest@0	1442
forrest@0	1443	plot_name = "global_mean_model_vs_ob"
forrest@0	1444
forrest@0	1445	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1446	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1447
forrest@0	1448	delete (plot)
forrest@0	1449	plot=new(3,graphic) ; create graphic array
forrest@0	1450
forrest@0	1451	resg@gsnFrame = False ; Do not draw plot
forrest@0	1452	resg@gsnDraw = False ; Do not advance frame
forrest@0	1453
forrest@0	1454	; plot correlation coef
forrest@0	1455
forrest@0	1456	gRes = True
forrest@0	1457	gRes@txFontHeightF = 0.02
forrest@0	1458	gRes@txAngleF = 90
forrest@0	1459
forrest@0	1460	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"
forrest@0	1461
forrest@0	1462	gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0	1463	;--------------------------------------------------------------------
forrest@0	1464
forrest@0	1465	;(a) ob
forrest@0	1466
forrest@0	1467	title = ob_name
forrest@0	1468	resg@tiMainString = title
forrest@0	1469
forrest@0	1470	plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg)
forrest@0	1471
forrest@0	1472	;(b) model
forrest@0	1473
forrest@0	1474	title = "Model "+ model_name
forrest@0	1475	resg@tiMainString = title
forrest@0	1476
forrest@0	1477	plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg)
forrest@0	1478
forrest@0	1479	;(c) model-ob
forrest@0	1480
forrest@0	1481	zz = laimod_mean
forrest@0	1482	zz = laimod_mean - laiob_mean
forrest@0	1483	title = "Model_"+model_name+" - Observed"
forrest@0	1484	resg@tiMainString = title
forrest@0	1485
forrest@0	1486	resg@cnMinLevelValF = -2. ; Min level
forrest@0	1487	resg@cnMaxLevelValF = 2. ; Max level
forrest@0	1488	resg@cnLevelSpacingF = 0.4 ; interval
forrest@0	1489
forrest@0	1490
forrest@0	1491	plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
forrest@0	1492
forrest@0	1493	pres = True ; panel plot mods desired
forrest@0	1494	; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0	1495	; indiv. plots in panel
forrest@0	1496	pres@gsnMaximize = True ; fill the page
forrest@0	1497
forrest@0	1498	gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
forrest@0	1499
forrest@0	1500	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1501	"rm "+plot_name+"."+plot_type)
forrest@0	1502
forrest@0	1503	clear (wks)
forrest@0	1504	delete (plot)
forrest@0	1505	;-----------------------------------------------------------------
forrest@0	1506	;global contour model vs ob
forrest@0	1507
forrest@0	1508	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1509	resg@cnMinLevelValF = 0. ; Min level
forrest@0	1510	resg@cnMaxLevelValF = 10. ; Max level
forrest@0	1511	resg@cnLevelSpacingF = 1. ; interval
forrest@0	1512
forrest@0	1513	plot_name = "global_max_model_vs_ob"
forrest@0	1514
forrest@0	1515	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1516	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1517
forrest@0	1518	; delete(plot)
forrest@0	1519	plot=new(3,graphic) ; create graphic array
forrest@0	1520
forrest@0	1521	resg@gsnFrame = False ; Do not draw plot
forrest@0	1522	resg@gsnDraw = False ; Do not advance frame
forrest@0	1523
forrest@0	1524	; plot correlation coef
forrest@0	1525
forrest@0	1526	gRes = True
forrest@0	1527	gRes@txFontHeightF = 0.02
forrest@0	1528	gRes@txAngleF = 90
forrest@0	1529
forrest@0	1530	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
forrest@0	1531
forrest@0	1532	gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0	1533	;--------------------------------------------------------------------
forrest@0	1534	;(a) ob
forrest@0	1535
forrest@0	1536	title = ob_name
forrest@0	1537	resg@tiMainString = title
forrest@0	1538
forrest@0	1539	plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg)
forrest@0	1540
forrest@0	1541	;(b) model
forrest@0	1542
forrest@0	1543	title = "Model "+ model_name
forrest@0	1544	resg@tiMainString = title
forrest@0	1545
forrest@0	1546	plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg)
forrest@0	1547
forrest@0	1548	;(c) model-ob
forrest@0	1549
forrest@0	1550	delete (zz)
forrest@0	1551	zz = laimod_max
forrest@0	1552	zz = laimod_max - laiob_max
forrest@0	1553	title = "Model_"+model_name+" - Observed"
forrest@0	1554	resg@tiMainString = title
forrest@0	1555
forrest@0	1556	resg@cnMinLevelValF = -6. ; Min level
forrest@0	1557	resg@cnMaxLevelValF = 6. ; Max level
forrest@0	1558	resg@cnLevelSpacingF = 1. ; interval
forrest@0	1559
forrest@0	1560
forrest@0	1561	plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
forrest@0	1562
forrest@0	1563	pres = True ; panel plot mods desired
forrest@0	1564	; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0	1565	; indiv. plots in panel
forrest@0	1566	pres@gsnMaximize = True ; fill the page
forrest@0	1567
forrest@0	1568	gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
forrest@0	1569
forrest@0	1570	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1571	"rm "+plot_name+"."+plot_type)
forrest@0	1572
forrest@0	1573	clear (wks)
forrest@0	1574	delete (plot)
forrest@0	1575	;-----------------------------------------------------------------
forrest@0	1576	;global contour model vs ob
forrest@0	1577
forrest@0	1578	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1579	resg@cnMinLevelValF = 1. ; Min level
forrest@0	1580	resg@cnMaxLevelValF = 12. ; Max level
forrest@0	1581	resg@cnLevelSpacingF = 1. ; interval
forrest@0	1582
forrest@0	1583	plot_name = "global_phase_model_vs_ob"
forrest@0	1584
forrest@0	1585	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1586	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1587
forrest@0	1588	; delete (plot)
forrest@0	1589	plot=new(3,graphic) ; create graphic array
forrest@0	1590
forrest@0	1591	resg@gsnFrame = False ; Do not draw plot
forrest@0	1592	resg@gsnDraw = False ; Do not advance frame
forrest@0	1593
forrest@0	1594	; plot correlation coef
forrest@0	1595
forrest@0	1596	gRes = True
forrest@0	1597	gRes@txFontHeightF = 0.02
forrest@0	1598	gRes@txAngleF = 90
forrest@0	1599
forrest@0	1600	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
forrest@0	1601
forrest@0	1602	gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0	1603	;--------------------------------------------------------------------
forrest@0	1604	;(a) ob
forrest@0	1605
forrest@0	1606	title = ob_name
forrest@0	1607	resg@tiMainString = title
forrest@0	1608
forrest@0	1609	plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
forrest@0	1610
forrest@0	1611	;(b) model
forrest@0	1612
forrest@0	1613	title = "Model "+ model_name
forrest@0	1614	resg@tiMainString = title
forrest@0	1615
forrest@0	1616	plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
forrest@0	1617
forrest@0	1618	;(c) model-ob
forrest@0	1619
forrest@0	1620	delete (zz)
forrest@0	1621	zz = laimod_phase
forrest@0	1622	zz = laimod_phase - laiob_phase
forrest@0	1623	title = "Model_"+model_name+" - Observed"
forrest@0	1624	resg@tiMainString = title
forrest@0	1625
forrest@0	1626	resg@cnMinLevelValF = -6. ; Min level
forrest@0	1627	resg@cnMaxLevelValF = 6. ; Max level
forrest@0	1628	resg@cnLevelSpacingF = 1. ; interval
forrest@0	1629
forrest@0	1630
forrest@0	1631	plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
forrest@0	1632
forrest@0	1633	pres = True ; panel plot mods desired
forrest@0	1634	; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0	1635	; indiv. plots in panel
forrest@0	1636	pres@gsnMaximize = True ; fill the page
forrest@0	1637
forrest@0	1638	gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
forrest@0	1639
forrest@0	1640	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1641	"rm "+plot_name+"."+plot_type)
forrest@0	1642
forrest@0	1643	clear (wks)
forrest@0	1644	delete (plot)
forrest@0	1645	;------------------------------------------------------------------
forrest@0	1646	;global contour model vs ob
forrest@0	1647
forrest@0	1648	resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0	1649	resg@cnMinLevelValF = 60. ; Min level
forrest@0	1650	resg@cnMaxLevelValF = 360. ; Max level
forrest@0	1651	resg@cnLevelSpacingF = 20. ; interval
forrest@0	1652
forrest@0	1653	plot_name = "global_grow_model_vs_ob"
forrest@0	1654
forrest@0	1655	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1656	gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0	1657
forrest@0	1658	; delete (plot)
forrest@0	1659	plot=new(3,graphic) ; create graphic array
forrest@0	1660
forrest@0	1661	resg@gsnFrame = False ; Do not draw plot
forrest@0	1662	resg@gsnDraw = False ; Do not advance frame
forrest@0	1663
forrest@0	1664	; plot correlation coef
forrest@0	1665
forrest@0	1666	gRes = True
forrest@0	1667	gRes@txFontHeightF = 0.02
forrest@0	1668	gRes@txAngleF = 90
forrest@0	1669
forrest@0	1670	correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
forrest@0	1671
forrest@0	1672	gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0	1673	;--------------------------------------------------------------------
forrest@0	1674	;(a) ob
forrest@0	1675
forrest@0	1676	title = ob_name
forrest@0	1677	resg@tiMainString = title
forrest@0	1678
forrest@0	1679	plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
forrest@0	1680
forrest@0	1681	;(b) model
forrest@0	1682
forrest@0	1683	title = "Model "+ model_name
forrest@0	1684	resg@tiMainString = title
forrest@0	1685
forrest@0	1686	plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
forrest@0	1687
forrest@0	1688	;(c) model-ob
forrest@0	1689
forrest@0	1690	delete (zz)
forrest@0	1691	zz = laimod_grow
forrest@0	1692	zz = laimod_grow - laiob_grow
forrest@0	1693	title = "Model_"+model_name+" - Observed"
forrest@0	1694	resg@tiMainString = title
forrest@0	1695
forrest@0	1696	resg@cnMinLevelValF = -100. ; Min level
forrest@0	1697	resg@cnMaxLevelValF = 100. ; Max level
forrest@0	1698	resg@cnLevelSpacingF = 20. ; interval
forrest@0	1699
forrest@0	1700	plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
forrest@0	1701
forrest@0	1702	pres = True ; panel plot mods desired
forrest@0	1703	; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0	1704	; indiv. plots in panel
forrest@0	1705	pres@gsnMaximize = True ; fill the page
forrest@0	1706
forrest@0	1707	gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
forrest@0	1708
forrest@0	1709	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1710	"rm "+plot_name+"."+plot_type)
forrest@0	1711
forrest@0	1712	clear (wks)
forrest@0	1713	delete (plot)
forrest@0	1714	;**************************************************
forrest@0	1715	; plot lai table
forrest@0	1716	;**************************************************
forrest@0	1717
forrest@0	1718	plot_name = "table_lai"
forrest@0	1719	wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0	1720
forrest@0	1721	gsn_table(wks,ncr1,xx1,yy1,text1,res1)
forrest@0	1722	gsn_table(wks,ncr21,xx21,yy21,text21,res21)
forrest@0	1723	gsn_table(wks,ncr22,xx22,yy22,text22,res22)
forrest@0	1724	gsn_table(wks,ncr3,xx3,yy3,text3,res3)
forrest@0	1725	gsn_table(wks,ncr4,xx4,yy4,text4,res4)
forrest@0	1726
forrest@0	1727	frame(wks)
forrest@0	1728
forrest@0	1729	system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0	1730	"rm "+plot_name+"."+plot_type)
forrest@0	1731	;-------------------------------------------------------------------
forrest@0	1732	temp_name = "lai." + model_name
forrest@0	1733	system("mkdir -p " + temp_name+";"+ \
forrest@0	1734	"cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
forrest@0	1735	"mv *.png " + temp_name +";"+ \
forrest@0	1736	"tar cf "+ temp_name +".tar " + temp_name)
forrest@0	1737	;-------------------------------------------------------------------
forrest@0	1738	end
forrest@0	1739

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