diff -r 000000000000 -r 0c6405ab2ff4 lai/99.all.ncl.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lai/99.all.ncl.1 Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,1763 @@ +;******************************************************** +; histogram normalized by rain and compute correleration +;******************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" + +procedure pminmax(data:numeric,name:string) +begin + print ("min/max " + name + " = " + min(data) + "/" + max(data)) + if(isatt(data,"units")) then + print (name + " units = " + data@units) + end if +end + +; Main code. +begin + + nclass = 20 + + plot_type = "ps" + plot_type_new = "png" + +;************************************************ +; read in data: model +;************************************************ +;film = "b30.061n_1995-2004_MONS_climo_lnd.nc" +;model_name = "b30.061n" +;model_grid = "T31" + +;film = "newcn05_ncep_1i_MONS_climo_lnd.nc" +;model_name = "newcn" +;model_grid = "1.9" + +;film = "i01.06cn_1798-2004_MONS_climo.nc" +;model_name = "06cn" +;model_grid = "T42" + +;film = "i01.06casa_1798-2004_MONS_climo.nc" +;model_name = "06casa" +;model_grid = "T42" + + film = "i01.10cn_1948-2004_MONS_climo.nc" + model_name = "10cn" + model_grid = "T42" + +;film = "i01.10casa_1948-2004_MONS_climo.nc" +;model_name = "10casa" +;model_grid = "T42" + + system("sed s#model_name#"+model_name+"# table.html > table.html.new") + system("mv -f table.html.new table.html") +;------------------------------------------------ + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + fm = addfile(dirm+film,"r") + + laimod = fm->TLAI + +;************************************************ +; read in data: observed +;************************************************ + + ob_name = "MODIS MOD 15A2 2000-2005" + + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" + filo1 = "land_class_"+model_grid+".nc" + filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc" + + fo1 = addfile(diro+filo1,"r") + fo2 = addfile(diro+filo2,"r") + + classob = tofloat(fo1->LAND_CLASS) + laiob = fo2->LAI + +;******************************************************************* +; for plotting table +;******************************************************************* +; table header name + table_header_name = "LAI" + +; column (not including header column) + + col_header1 = (/"Mean","Max","Phase","Growth"/) + col_header2 = (/"ob","model","M" \ + ,"ob","model","M" \ + ,"ob","model","M" \ + ,"ob","model","M" \ + /) + + ncol1 = dimsizes(col_header1) + ncol2 = dimsizes(col_header2) + ncol = ncol2 + +; row (not including header row) + row_header = (/"Water Bodies" \ + ,"Evergreen Needleleaf Forests" \ + ,"Evergreen Broadleaf Forests" \ + ,"Deciduous Needleleaf Forest" \ + ,"Deciduous Broadleaf Forests" \ + ,"Mixed Forests" \ + ,"Closed Bushlands" \ + ,"Open Bushlands" \ + ,"Woody Savannas (S. Hem.)" \ + ,"Savannas (S. Hem.)" \ + ,"Grasslands" \ + ,"Permanent Wetlands" \ + ,"Croplands" \ + ,"Urban and Built-Up" \ + ,"Cropland/Natural Vegetation Mosaic" \ + ,"Permanent Snow and Ice" \ + ,"Barren or Sparsely Vegetated" \ + ,"Unclassified" \ + ,"Woody Savannas (N. Hem.)" \ + ,"Savannas (N. Hem.)" \ + ,"All biome average" \ + /) + nrow = dimsizes(row_header) + +; arrays to be passed to table. + value = new ((/nrow, ncol/),string ) + + table_length = 0.995 + +; Table header + ncr1 = (/1,1/) ; 1 row, 1 column + xx1 = (/0.005,0.25/) ; Start and end X + yy1 = (/0.900,0.995/) ; Start and end Y + text1 = table_header_name + res1 = True + res1@txFontHeightF = 0.03 + res1@gsFillColor = "CornFlowerBlue" + +; Column header (equally space in x2) + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns + xx21 = (/xx1(1),0.995/) ; start from end of x1 + yy21 = (/0.9475,0.995/) ; half of y1 + text21 = col_header1 + res21 = True + res21@txFontHeightF = 0.015 + res21@gsFillColor = "Gray" + + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns + xx22 = xx21 ; start from end of x1 + yy22 = (/0.900,0.9475/) ; half of y1 + text22 = col_header2 + res22 = True + res22@txFontHeightF = 0.015 + res22@gsFillColor = "Gray" + +; Row header (equally space in y2) + ncr3 = (/nrow,1/) ; 20 rows, 1 columns + xx3 = xx1 ; same as x1 + yy3 = (/1.0-table_length,0.900/) ; end at start of y1 + text3 = row_header + res3 = True + res3@txFontHeightF = 0.01 + res3@gsFillColor = "Gray" + +; Main table body + ncr4 = (/nrow,ncol/) ; 5 rows, 5 columns + xx4 = xx21 ; Start and end x + yy4 = yy3 ; Start and end Y + text4 = new((/nrow,ncol/),string) + + color_fill4 = new((/nrow,ncol/),string) + color_fill4 = "white" + color_fill4(nrow-1,:) = "CornFlowerBlue" + + res4 = True ; Set up resource list +; res4@gsnDebug = True ; Useful to print NDC row,col values used. + res4@txFontHeightF = 0.015 + res4@gsFillColor = color_fill4 + + delete (color_fill4) + +;************************************************ +; plot global land class: observed +;************************************************ +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 1. ; Min level + resg@cnMaxLevelValF = 19. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +;global contour ob + classob@_FillValue = 1.e+36 + classob = where(classob.eq.0,classob@_FillValue,classob) + + plot_name = "global_class_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,classob,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;******************************************************************* +; Calculate "nice" bins for binning the data in equally spaced ranges +;******************************************************************** + nclassn = nclass + 1 + range = fspan(0,nclassn-1,nclassn) +; print (range) + +; Use this range information to grab all the values in a +; particular range, and then take an average. + + nr = dimsizes(range) + nx = nr-1 + xvalues = new((/2,nx/),float) + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. + dx = xvalues(0,1) - xvalues(0,0) ; range width + dx4 = dx/4 ; 1/4 of the range + xvalues(1,:) = xvalues(0,:) - dx/5. +;----------------------------------------------------------------- +;(A) mean +;-------------------------------------------------------------------- +; get data + + laiob_mean = dim_avg_Wrap(laiob(lat|:,lon|:,time|:)) + laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:)) + + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_mean) + DATA22_1D = ndtooned(laimod_mean) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(data_ob.ge.range(i)) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrMean = esccr(uu,vv,0) +; print (ccrMean) + +; new eq + bias = sum(abs(vv-uu)/(vv+uu)) + Mmean = (1.- (bias/dimsizes(uu)))*5. + + M_lai_mean = sprintf("%.2f", Mmean) + system("sed s#M_lai_mean#"+M_lai_mean+"# table.html > table.html.new") + system("mv -f table.html.new table.html") + print (M_lai_mean) + + do i=0,nrow-2 + text4(i,0) = sprintf("%5.2f",u(i)) + text4(i,1) = sprintf("%5.2f",v(i)) + text4(i,2) = "-" + end do + text4(nrow-1,0) = sprintf("%5.2f",avg(u)) + text4(nrow-1,1) = sprintf("%5.2f",avg(v)) + text4(nrow-1,2) = sprintf("%5.2f",Mmean) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Mean LAI (Leaf Area Index)" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_mean" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 10. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +;global contour ob + + delta = 0.00001 + laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean) + + plot_name = "global_mean_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + plot_name = "global_mean_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;----------------------------------------------------------------- +;(B) max +;-------------------------------------------------------------------- +; get data + +; observed + laiob_max = laiob(0,:,:) + s = laiob(:,0,0) + laiob_max@long_name = "Leaf Area Index Max" + + dsizes_z = dimsizes(laiob) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laiob(:,j,i) + laiob_max(j,i) = max(s) + end do + end do + +; print (min(y)+"/"+max(y)) + delete (s) + delete (dsizes_z) +;------------------------- +; model + laimod_max = laimod(0,:,:) + s = laimod(:,0,0) + laimod_max@long_name = "Leaf Area Index Max" + + dsizes_z = dimsizes(laimod) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laimod(:,j,i) + laimod_max(j,i) = max(s) + end do + end do + +; print (min(laimod_max)+"/"+max(laimod_max)) + delete (s) + delete (dsizes_z) +;------------------------ + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_max) + DATA22_1D = ndtooned(laimod_max) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data_ob) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrMax = esccr(uu,vv,0) +; print (ccrMax) + +; new eq + bias = sum(abs(vv-uu)/(vv+uu)) + Mmax = (1.- (bias/dimsizes(uu)))*5. + + M_lai_max = sprintf("%.2f", Mmax) + system("sed s#M_lai_max#"+M_lai_max+"# table.html > table.html.new") + system("mv -f table.html.new table.html") + print (M_lai_max) + + do i=0,nrow-2 + text4(i,3) = sprintf("%5.2f",u(i)) + text4(i,4) = sprintf("%5.2f",v(i)) + text4(i,5) = "-" + end do + text4(nrow-1,3) = sprintf("%5.2f",avg(u)) + text4(nrow-1,4) = sprintf("%5.2f",avg(v)) + text4(nrow-1,5) = sprintf("%5.2f",Mmax) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Max LAI (Leaf Area Index)" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_max" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 10. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +;global contour ob + + delta = 0.00001 + laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max) + + plot_name = "global_max_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + plot_name = "global_max_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;(C) phase +;-------------------------------------------------------------------- +; get data + +; observed + laiob_phase = laiob(0,:,:) + s = laiob(:,0,0) + laiob_phase@long_name = "Leaf Area Index Max Month" + + dsizes_z = dimsizes(laiob) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laiob(:,j,i) + laiob_phase(j,i) = maxind(s) + 1 + end do + end do + +; print (min(laiob_phase)+"/"+max(laiob_phase)) + delete (s) + delete (dsizes_z) +;------------------------- +; model + laimod_phase = laimod(0,:,:) + s = laimod(:,0,0) + laimod_phase@long_name = "Leaf Area Index Max Month" + + dsizes_z = dimsizes(laimod) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laimod(:,j,i) + laimod_phase(j,i) = maxind(s) + 1 + end do + end do + +; print (min(laimod_phase)+"/"+max(laimod_phase)) + delete (s) + delete (dsizes_z) +;------------------------ + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_phase) + DATA22_1D = ndtooned(laimod_phase) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data_ob) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrPhase = esccr(uu,vv,0) +; print (ccrPhase) + +; old eq +; bias = abs(avg(vv)-avg(uu)) +; new eq + bias = avg(abs(vv-uu)) + + bias = where((bias.gt. 6.),12.-bias,bias) + Mphase = ((6. - bias)/6.)*5. + + M_lai_phase = sprintf("%.2f", Mphase) + system("sed s#M_lai_phase#"+M_lai_phase+"# table.html > table.html.new") + system("mv -f table.html.new table.html") + print (M_lai_phase) + + do i=0,nrow-2 + text4(i,6) = sprintf("%5.2f",u(i)) + text4(i,7) = sprintf("%5.2f",v(i)) + text4(i,8) = "-" + end do + text4(nrow-1,6) = sprintf("%5.2f",avg(u)) + text4(nrow-1,7) = sprintf("%5.2f",avg(v)) + text4(nrow-1,8) = sprintf("%5.2f",Mphase) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Max LAI (Leaf Area Index) Month" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_phase" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 1. ; Min level + resg@cnMaxLevelValF = 12. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +;global contour ob + + delta = 0.00001 + laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase) + + plot_name = "global_phase_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + plot_name = "global_phase_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;----------------------------------------------------------------- +;(D) grow day +;-------------------------------------------------------------------- +; get data + + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) + +; observed + laiob_grow = laiob(0,:,:) + laiob_grow@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(laiob) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + + laiob_grow(j,i) = nday + end do + end do + +; print (min(laiob_grow)+"/"+max(laiob_grow)) +;------------------------- +; model + laimod_grow = laimod(0,:,:) + laimod_grow@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(laimod) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + + laimod_grow(j,i) = nday + end do + end do + +; print (min(laimod_grow)+"/"+max(laimod_grow)) +;------------------------ + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_grow) + DATA22_1D = ndtooned(laimod_grow) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data_ob) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrGrow = esccr(uu,vv,0) +; print (ccrGrow) + +; new eq + bias = sum(abs(vv-uu)/(vv+uu)) + Mgrow = (1.- (bias/dimsizes(uu)))*5. + + M_lai_grow = sprintf("%.2f", Mgrow) + system("sed s#M_lai_grow#"+M_lai_grow+"# table.html > table.html.new") + system("mv -f table.html.new table.html") + print (M_lai_grow) + + do i=0,nrow-2 + text4(i,9) = sprintf("%5.2f",u(i)) + text4(i,10) = sprintf("%5.2f",v(i)) + text4(i,11) = "-" + end do + text4(nrow-1,9) = sprintf("%5.2f",avg(u)) + text4(nrow-1,10) = sprintf("%5.2f",avg(v)) + text4(nrow-1,11) = sprintf("%5.2f",Mgrow) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Days of Growing season" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_grow" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 60. ; Min level + resg@cnMaxLevelValF = 360. ; Max level + resg@cnLevelSpacingF = 20. ; interval + +;global contour ob + + laiob_grow@_FillValue = 1.e+36 + laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow) + + plot_name = "global_grow_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + laimod_grow@_FillValue = 1.e+36 + laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow) + + plot_name = "global_grow_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model vs ob + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 10. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + plot_name = "global_mean_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- + +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg) + +;(c) model-ob + + zz = laimod_mean + zz = laimod_mean - laiob_mean + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -2. ; Min level + resg@cnMaxLevelValF = 2. ; Max level + resg@cnLevelSpacingF = 0.4 ; interval + + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + + pres = True ; panel plot mods desired + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;----------------------------------------------------------------- +;global contour model vs ob + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 10. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + plot_name = "global_max_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + +; delete(plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg) + +;(c) model-ob + + delete (zz) + zz = laimod_max + zz = laimod_max - laiob_max + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -6. ; Min level + resg@cnMaxLevelValF = 6. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + + pres = True ; panel plot mods desired + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;----------------------------------------------------------------- +;global contour model vs ob + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 1. ; Min level + resg@cnMaxLevelValF = 12. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + plot_name = "global_phase_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + +; delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg) + +;(c) model-ob + + delete (zz) + zz = laimod_phase + zz = laimod_phase - laiob_phase + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -6. ; Min level + resg@cnMaxLevelValF = 6. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + +; pres = True ; panel plot mods desired +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel +; pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;------------------------------------------------------------------ +;global contour model vs ob + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 60. ; Min level + resg@cnMaxLevelValF = 360. ; Max level + resg@cnLevelSpacingF = 20. ; interval + + plot_name = "global_grow_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + +; delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg) + +;(c) model-ob + + delete (zz) + zz = laimod_grow + zz = laimod_grow - laiob_grow + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -100. ; Min level + resg@cnMaxLevelValF = 100. ; Max level + resg@cnLevelSpacingF = 20. ; interval + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + + pres = True ; panel plot mods desired + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;************************************************** +; plot lai table +;************************************************** + + plot_name = "table_lai" + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + + gsn_table(wks,ncr1,xx1,yy1,text1,res1) + gsn_table(wks,ncr21,xx21,yy21,text21,res21) + gsn_table(wks,ncr22,xx22,yy22,text22,res22) + gsn_table(wks,ncr3,xx3,yy3,text3,res3) + gsn_table(wks,ncr4,xx4,yy4,text4,res4) + + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) +;------------------------------------------------------------------- + temp_name = "lai." + model_name + system("mkdir -p " + temp_name) + system("cp table.html " + temp_name) + system("mv *.png " + temp_name) + system("tar cf "+ temp_name +".tar " + temp_name) +;------------------------------------------------------------------- +end +