diff -r 000000000000 -r 0c6405ab2ff4 ameriflux/22.plot.ncl
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ameriflux/22.plot.ncl Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,919 @@
+;************************************************************
+; merge annual
+; add table
+; add time series plot
+; sort by latitude in decending order (N->S)
+; find year_ob >= 4
+; add year_ob_i and year_ob_f
+; change long_name
+;************************************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+;************************************************************
+procedure set_line(lines:string,nline:integer,newlines:string)
+begin
+; add line to ascci/html file
+
+ nnewlines = dimsizes(newlines)
+ if(nline+nnewlines-1.ge.dimsizes(lines))
+ print("set_line: bad index, not setting anything.")
+ return
+ end if
+ lines(nline:nline+nnewlines-1) = newlines
+; print ("lines = " + lines(nline:nline+nnewlines-1))
+ nline = nline + nnewlines
+ return
+end
+;*************************************************************
+begin
+
+ plot_type = "ps"
+ plot_type_new = "png"
+
+;************************************************
+; model data info
+;************************************************
+ model = "cn"
+; model = "casa"
+
+ model_name = "i01.10" + model
+
+ ENERGY ="new"
+
+ nmonth = 12
+
+;************************************************
+; observed data info
+;************************************************
+
+ station = (/"ARM_Oklahoma" \
+ ,"ARM_Oklahoma_burn" \
+ ,"ARM_Oklahoma_control" \
+ ,"Atqasuk" \
+ ,"Audubon" \
+ ,"AustinCary" \
+ ,"Bartlett" \
+ ,"Bondville" \
+ ,"Brookings" \
+ ,"Donaldson" \
+ ,"Duke_Forest_Hardwoods" \
+ ,"Duke_Forest_Open_Field" \
+ ,"Duke_Forest_Pine" \
+ ,"Fermi_Ag" \
+ ,"Fermi_Prairie" \
+ ,"Flagstaff_Managed" \
+ ,"Flagstaff_Unmanaged" \
+ ,"Flagstaff_Wildfire" \
+ ,"FortPeck" \
+ ,"FreemanRanch_mesquite" \
+ ,"Goodwin_Creek" \
+ ,"HarvardForest" \
+ ,"HarvardForestHemlock" \
+ ,"HowlandForestMain" \
+ ,"HowlandForestWest" \
+ ,"Ivotuk" \
+ ,"KendallGrasslands" \
+ ,"KennedySpaceCenterPine" \
+ ,"KennedySpaceCenterScrub" \
+ ,"LittleProspect" \
+ ,"LostCreek" \
+ ,"Mead-irrigated" \
+ ,"Mead-irrigated-rotation" \
+ ,"Mead-rainfed" \
+ ,"Metolius_2nd_YoungPonderosaPine" \
+ ,"MetoliusEyerly" \
+ ,"MetoliusIntermediatePine" \
+ ,"MetoliusOldPonderosaPine" \
+ ,"MissouriOzark" \
+ ,"Mize" \
+ ,"MorganMonroe" \
+ ,"NiwotRidge" \
+ ,"NorthCarolina_cc" \
+ ,"NorthCarolina_lp" \
+ ,"ParkFalls" \
+ ,"Rayonier" \
+ ,"SantaRita" \
+ ,"SkyOaks_Old" \
+ ,"SkyOaks_PostFire" \
+ ,"SkyOaks_Young" \
+ ,"SylvaniaWilderness" \
+ ,"Toledo" \
+ ,"Tonzi" \
+ ,"UCI_1850" \
+ ,"UCI_1930" \
+ ,"UCI_1964" \
+ ,"UCI_1964wet" \
+ ,"UCI_1981" \
+ ,"UCI_1989" \
+ ,"UCI_1998" \
+ ,"UMBS" \
+ ,"Vaira" \
+ ,"WalkerBranch" \
+ ,"WillowCreek" \
+ ,"WindRiver" \
+ ,"Wisconsin_ihw" \
+ ,"Wisconsin_irp" \
+ ,"Wisconsin_mrp" \
+ ,"Wisconsin_myjp" \
+ ,"Wisconsin_pb" \
+ ,"Wisconsin_rpcc" \
+ ,"Wisconsin_yhw" \
+ ,"Wisconsin_yjp" \
+ ,"Wisconsin_yrp" \
+ /)
+
+ year_ob = (/"2003-2006" \
+ ,"2005-2006" \
+ ,"2005-2006" \
+ ,"1999-2006" \
+ ,"2002-2006" \
+ ,"2000-2005" \
+ ,"2004-2005" \
+ ,"1996-2006" \
+ ,"2004-2006" \
+ ,"1999-2004" \
+ ,"2003-2005" \
+ ,"2001-2005" \
+ ,"2001-2005" \
+ ,"2005-2006" \
+ ,"2004-2006" \
+ ,"2005-2006" \
+ ,"2005-2006" \
+ ,"2005-2006" \
+ ,"2000-2006" \
+ ,"2004-2006" \
+ ,"2002-2006" \
+ ,"1991-2004" \
+ ,"2004-2004" \
+ ,"1996-2004" \
+ ,"1999-2004" \
+ ,"2003-2006" \
+ ,"2004-2006" \
+ ,"2002-2002" \
+ ,"2000-2006" \
+ ,"2002-2005" \
+ ,"2001-2005" \
+ ,"2001-2005" \
+ ,"2001-2005" \
+ ,"2001-2005" \
+ ,"2004-2005" \
+ ,"2004-2005" \
+ ,"2003-2005" \
+ ,"1996-2000" \
+ ,"2004-2006" \
+ ,"1998-2004" \
+ ,"1999-2005" \
+ ,"1999-2003" \
+ ,"2005-2006" \
+ ,"2005-2006" \
+ ,"1996-2003" \
+ ,"1998-1998" \
+ ,"2004-2006" \
+ ,"1997-2006" \
+ ,"2004-2006" \
+ ,"1997-2006" \
+ ,"2002-2006" \
+ ,"2004-2005" \
+ ,"2001-2006" \
+ ,"2004-2005" \
+ ,"2001-2005" \
+ ,"2001-2005" \
+ ,"2002-2004" \
+ ,"2001-2005" \
+ ,"2001-2005" \
+ ,"2002-2005" \
+ ,"1999-2003" \
+ ,"2001-2006" \
+ ,"1995-1999" \
+ ,"1999-2005" \
+ ,"1999-2004" \
+ ,"2003-2003" \
+ ,"2003-2003" \
+ ,"2002-2005" \
+ ,"2004-2004" \
+ ,"2002-2002" \
+ ,"2005-2005" \
+ ,"2002-2002" \
+ ,"2004-2005" \
+ ,"2002-2002" \
+ /)
+
+ field = (/"NEE Flux" \
+ ,"Net Radiation" \
+ ,"Latent Heat" \
+ ,"Sensible Heat" \
+ ,"GPP Flux" \
+ ,"Respiration" \
+ /)
+
+ field_unit = (/"u mol/m2/s" \
+ ,"W/m2" \
+ ,"W/m2" \
+ ,"W/m2" \
+ ,"u mol/m2/s" \
+ ,"u mol/m2/s" \
+ /)
+
+ nstation = dimsizes(station)
+ nfield = dimsizes(field)
+
+;========================================================================
+; get observed info: number of year, first/last year, lat, lon
+; and annual data
+
+ dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
+
+ year_station = new ((/nstation/),integer) ; number of year
+ year_ob_i = new ((/nstation/),integer) ; first year
+ year_ob_f = new ((/nstation/),integer) ; last year
+ lat_ob = new ((/nstation/),float) ; latitude
+ lon_ob = new ((/nstation/),float) ; longitude
+
+ data_ob_ann = new ((/nfield, nmonth, nstation/),float)
+
+ do n = 0, nstation-1
+
+ diro = dir_root + station(n)+"/"
+ filo = year_ob(n)+"_L4_m.nc"
+ fo = addfile (diro+filo,"r")
+
+ lat_ob(n) = fo->lat
+ lon_ob(n) = fo->lon
+
+ year = fo->year
+
+ year_station(n) = dimsizes(year)
+ year_ob_i(n) = year(0)
+ year_ob_f(n) = year(year_station(n)-1)
+
+ delete (year)
+
+ data = fo->NEE_or_fMDS
+ data_ob_ann(0,:,n) = dim_avg(data(month|:,year|:))
+
+ data = fo->Rg_f
+ data_ob_ann(1,:,n) = dim_avg(data(month|:,year|:))
+
+ data = fo->LE_f
+ data_ob_ann(2,:,n) = dim_avg(data(month|:,year|:))
+
+ data = fo->H_f
+ data_ob_ann(3,:,n) = dim_avg(data(month|:,year|:))
+
+ data = fo->GPP_or_MDS
+ data_ob_ann(4,:,n) = dim_avg(data(month|:,year|:))
+
+ data = fo->Reco_or
+ data_ob_ann(5,:,n) = dim_avg(data(month|:,year|:))
+
+ delete (data)
+ delete (fo)
+ end do
+
+;print (year_ob_i)
+;print (year_ob_f)
+
+;--------------------------------------------------------------
+; find (# of year observed) >=4 and year_ob_i <= 2001
+
+ i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
+
+ station_long = station(i_long_ob)
+ lat_ob_long = lat_ob(i_long_ob)
+ lon_ob_long = lat_ob(i_long_ob)
+ year_ob_long = year_ob(i_long_ob)
+ year_ob_i_long = year_ob_i(i_long_ob)
+ year_ob_f_long = year_ob_f(i_long_ob)
+ year_station_long = year_station(i_long_ob)
+
+ nstation_long = dimsizes(station_long)
+
+; print (i_long_ob)
+; print (nstation_long)
+; print (station_long)
+; print (year_ob_i(i_long_ob))
+
+;=========================================================
+; get model data at observed lat-lon
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = model_name+"_ameriflux_1990-2004_monthly.nc"
+ fm = addfile (dirm+film,"r")
+
+ xm = fm->lon
+ ym = fm->lat
+ date = fm->date
+
+ date_dim = dimsizes(date)
+ nyear = date_dim(0)
+
+ data_mod = new ((/nfield,nyear,nmonth,nstation/),float)
+ data_mod_ann = new ((/nfield,nmonth,nstation/),float)
+ data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)
+
+; change to unit of observed (u mol/m2/s)
+; Model_units [=] gC/m2/s
+; 12. = molecular weight of C
+; u mol = 1e-6 mol
+
+ factor = 1.e6 /12.
+
+;------------------------------------------------------------
+; interpolate model data into observed station
+; note: model is 0-360E, 90S-90N
+
+; to be able to handle observation at (-89.98,-24.80)
+ ym(0) = -90.
+;------------------------------------------------------------
+
+if (ENERGY .eq. "old") then
+
+ data = fm->NEE
+ data_mod(0,:,:,:) = data(:,:,:) * factor
+
+; data = fm->LATENT
+ data1 = fm->FCEV
+ data2 = fm->FCTR
+ data3 = fm->FGEV
+ data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
+ delete (data1)
+ delete (data2)
+ delete (data3)
+
+; data = fm->SENSIBLE
+ data = fm->FSH
+ data_mod(3,:,:,:) = data(:,:,:)
+
+; data = fm->NETRAD
+ data1 = fm->FSA
+ data2 = fm->FIRA
+ data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:)
+ delete (data1)
+ delete (data2)
+
+else
+
+ data = fm->NEE
+ yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
+ data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
+ data_mod(0,:,:,:) = yy(:,:,:) * factor
+
+; printVarSummary(yy)
+
+ data = fm->NETRAD
+ yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
+ data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))
+ data_mod(1,:,:,:) = yy(:,:,:)
+
+ data = fm->LATENT
+ yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
+ data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))
+ data_mod(2,:,:,:) = yy(:,:,:)
+
+; data = fm->SENSIBLE
+ data = fm->FSH
+ yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
+ data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))
+ data_mod(3,:,:,:) = yy(:,:,:)
+
+end if
+
+ data = fm->GPP
+ yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
+ data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
+ data_mod(4,:,:,:) = yy(:,:,:) * factor
+
+ if (model .eq. "cn") then
+ data = fm->ER
+ else
+ data1 = fm->AR
+ data2 = fm->HR
+ data = data1 + data2
+
+ delete (data1)
+ delete (data2)
+ end if
+
+ yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
+ data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
+ data_mod(5,:,:,:) = yy(:,:,:) * factor
+
+ data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)
+
+ delete (data_mod)
+ delete (fm)
+ delete (data)
+ delete (yy)
+
+;*******************************************************************
+; for station line plot
+;*******************************************************************
+
+; for x-axis in xyplot
+ mon = ispan(1,12,1)
+ mon@long_name = "month"
+
+ res = True ; plot mods desired
+ res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
+ res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
+
+ res@tmXBFormat = "f" ; not to add trailing zeros
+
+;-------------------------------------------------------------------------
+; Add a boxed legend using the more simple method
+
+ res@pmLegendDisplayMode = "Always"
+; res@pmLegendWidthF = 0.1
+ res@pmLegendWidthF = 0.08
+ res@pmLegendHeightF = 0.06
+; res@pmLegendOrthogonalPosF = -1.17
+; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
+ res@pmLegendOrthogonalPosF = -0.30 ;(downward)
+
+; res@pmLegendParallelPosF = 0.18
+ res@pmLegendParallelPosF = 0.23 ;(rightward)
+
+; res@lgPerimOn = False
+ res@lgLabelFontHeightF = 0.015
+ res@xyExplicitLegendLabels = (/"observed",model_name/)
+;-------------------------------------------------------------------
+; for panel plot
+ res@gsnFrame = False ; Do not draw plot
+ res@gsnDraw = False ; Do not advance frame
+
+ pres = True ; panel plot mods desired
+ pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
+ ; indiv. plots in panel
+ pres@gsnMaximize = True ; fill the page
+;-------------------------------------------------------------------
+
+;==============================================================
+; get ob data at each site with long observation
+
+ dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
+
+ do n = 0,nstation_long-1
+; do n = 0,0
+
+ year_setback = 0
+
+ nyear = year_station_long(n)
+
+ if (year_ob_f_long(n).eq. 2006) then
+ year_setback = 2006 -2004
+ end if
+ if (year_ob_f_long(n).eq. 2005) then
+ year_setback = 2005 -2004
+ end if
+
+ ntime = (nyear - year_setback) * nmonth
+
+ data_ob = new ((/nfield, nyear, nmonth/),float)
+
+ diro = dir_root + station_long(n)+"/"
+ filo = year_ob_long(n)+"_L4_m.nc"
+ fo = addfile (diro+filo,"r")
+
+; print (diro)
+; print (filo)
+
+ data_ob(0,:,:) = fo->NEE_or_fMDS
+ data_ob(1,:,:) = fo->Rg_f
+ data_ob(2,:,:) = fo->LE_f
+ data_ob(3,:,:) = fo->H_f
+ data_ob(4,:,:) = fo->GPP_or_MDS
+ data_ob(5,:,:) = fo->Reco_or
+
+ delete (fo)
+
+ timeI = new((/ntime/),integer)
+ timeF = new((/ntime/),float)
+ timeI = ispan(1,ntime,1)
+ timeF = year_ob_i_long(n) + (timeI-1)/12.
+ timeF@long_name = "year"
+
+ plot_data = new((/2,ntime/),float)
+
+;----------------------------
+; for model_vs_ob
+
+ plot_name = station_long(n)+"_tseries_vs_ob"
+ title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"
+ res@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name)
+ plot=new(nfield,graphic) ; create graphic array
+
+ i_year_mod_i = year_ob_i_long(n) - 1990
+ i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
+
+ i_year_ob_f = nyear - year_setback - 1
+
+; print (nyear)
+; print (i_year_ob_f)
+; print (i_year_mod_i)
+; print (i_year_mod_f)
+
+ do i = 0,nfield-1
+ plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
+ plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))
+ plot_data@long_name = field(i)+" ("+field_unit(i)+")"
+ plot(i)=gsn_csm_xy(wks,timeF,plot_data,res) ; create plot
+ end do
+
+ gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+
+ clear (wks)
+ delete (plot)
+
+ delete (data_ob)
+ delete (timeI)
+ delete (timeF)
+ delete (plot_data)
+
+ end do
+
+;###################################################################
+; for the following tables,
+; sort by latitude in decending order (N->S)
+
+ isort = dim_pqsort(lat_ob_long,-1)
+
+ station_sort = station_long(isort)
+ year_ob_sort = year_ob_long(isort)
+ lat_ob_sort = lat_ob_long(isort)
+ lon_ob_sort = lon_ob_long(isort)
+
+; print(isort)
+; print(lat_ob_sort)
+
+;*******************************************************************
+; html table of site: observed
+;*******************************************************************
+ output_html = "tseries_vs_ob.html"
+
+ header = (/"" \
+ ,"
" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,"Timeseries at Site: "+model_name+" vs Observation
" \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site Name | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," Observed | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+; row of table
+
+ do n = 0,nstation_long-1
+
+ set_line(lines,nline,row_header)
+
+ txt0 = station_sort(n)
+ txt1 = sprintf("%5.2f", lat_ob_sort(n))
+ txt2 = sprintf("%5.2f", lon_ob_sort(n))
+ txt3 = year_ob_sort(n)
+
+ set_line(lines,nline,""+txt0+" | ")
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+ delete (idx)
+
+ delete (isort)
+ delete (station_sort)
+ delete (year_ob_sort)
+ delete (lat_ob_sort)
+ delete (lon_ob_sort)
+
+;************************************************************
+; compute annual cycle correlation coef and M score
+;************************************************************
+
+ score_max = 5.
+
+ ccr = new ((/nstation, nfield/),float)
+ M_score = new ((/nstation, nfield/),float)
+
+ do n=0,nstation-1
+ do m=0,nfield-1
+ ccr(n,m) = esccr(data_ob_ann(m,:,n),data_mod_ann(m,:,n),0)
+ bias = sum(abs(data_mod_ann(m,:,n)-data_ob_ann(m,:,n))/(abs(data_mod_ann(m,:,n))+abs(data_ob_ann(m,:,n))))
+ M_score(n,m) = (1. -(bias/nmonth)) * score_max
+ end do
+ end do
+
+ M_co2 = avg(M_score(:,0))
+ M_rad = avg(M_score(:,1))
+ M_lh = avg(M_score(:,2))
+ M_sh = avg(M_score(:,3))
+ M_gpp = avg(M_score(:,4))
+ M_er = avg(M_score(:,5))
+ M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er
+
+ M_energy_co2 = sprintf("%.2f", M_co2)
+ M_energy_rad = sprintf("%.2f", M_rad)
+ M_energy_lh = sprintf("%.2f", M_lh )
+ M_energy_sh = sprintf("%.2f", M_sh )
+ M_energy_gpp = sprintf("%.2f", M_gpp)
+ M_energy_er = sprintf("%.2f", M_er )
+ M_energy_all = sprintf("%.2f", M_all)
+
+;*******************************************************************
+; for station line plot
+;*******************************************************************
+
+; for x-axis in xyplot
+ mon = ispan(1,12,1)
+ mon@long_name = "month"
+
+;-------------------------------------------------------------------
+
+ plot_data = new((/2,nmonth/),float)
+ plot_data!0 = "case"
+ plot_data!1 = "month"
+
+ do n = 0,nstation-1
+;----------------------------
+; for observed
+
+ plot_name = station(n)+"_ob"
+ title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
+ res@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name)
+ plot=new(nfield,graphic) ; create graphic array
+
+ do i = 0,nfield-1
+ plot_data(0,:) = (/data_ob_ann(i,:,n)/)
+ plot_data@long_name = field(i)+" ("+field_unit(i)+")"
+ plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot
+ end do
+
+ gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+
+ clear (wks)
+ delete (plot)
+;----------------------------
+; for model_vs_ob
+
+ plot_name = station(n)+"_model_vs_ob"
+ title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
+ res@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name)
+ plot=new(nfield,graphic) ; create graphic array
+
+ do i = 0,nfield-1
+ plot_data(0,:) = (/data_ob_ann(i,:,n)/)
+ plot_data(1,:) = (/data_mod_ann(i,:,n)/)
+ plot_data@long_name = field(i)+" ("+field_unit(i)+")"
+ plot(i)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot
+ end do
+
+ gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+
+ clear (wks)
+ delete (plot)
+ end do
+
+;###################################################################
+; for the following tables,
+; sort by latitude in decending order (N->S)
+; sort by lat in decending order (N->S)
+
+ isort = dim_pqsort(lat_ob,-1)
+
+ station_sort = station(isort)
+ year_ob_sort = year_ob(isort)
+ lat_ob_sort = lat_ob(isort)
+ lon_ob_sort = lon_ob(isort)
+ M_score_sort = M_score(isort,:)
+
+; print(isort)
+; print(lat_ob_sort)
+;###################################################################
+;*******************************************************************
+; html table of site: observed
+;*******************************************************************
+ output_html = "line_ob.html"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,"Energy at Site: Observation
" \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site Name | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," Observed | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+; row of table
+
+ do n = 0,nstation-1
+ set_line(lines,nline,row_header)
+
+ txt0 = station_sort(n)
+ txt1 = sprintf("%5.2f", lat_ob_sort(n))
+ txt2 = sprintf("%5.2f", lon_ob_sort(n))
+ txt3 = year_ob_sort(n)
+
+ set_line(lines,nline,""+txt0+" | ")
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+ delete (idx)
+
+;*******************************************************************
+; score and line table : model vs observed
+;*******************************************************************
+ output_html = "score+line_vs_ob.html"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,"Energy at Site: Model "+model_name+"
" \
+ /)
+ footer = ""
+
+ delete (table_header)
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site Name | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," Observed | " \
+ ," CO2 Flux | " \
+ ," Net Radiation | " \
+ ," Latent Heat | " \
+ ," Sensible Heat | " \
+ ," GPP Glux | " \
+ ," Respiration | " \
+ ," Average | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+; row of table
+
+ do n = 0,nstation-1
+ set_line(lines,nline,row_header)
+
+ txt0 = station_sort(n)
+ txt1 = sprintf("%5.2f", lat_ob_sort(n))
+ txt2 = sprintf("%5.2f", lon_ob_sort(n))
+ txt3 = year_ob_sort(n)
+ txt4 = sprintf("%5.2f", M_score_sort(n,0))
+ txt5 = sprintf("%5.2f", M_score_sort(n,1))
+ txt6 = sprintf("%5.2f", M_score_sort(n,2))
+ txt7 = sprintf("%5.2f", M_score_sort(n,3))
+ txt8 = sprintf("%5.2f", M_score_sort(n,4))
+ txt9 = sprintf("%5.2f", M_score_sort(n,5))
+ txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))
+
+ set_line(lines,nline,""+txt0+" | ")
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+ set_line(lines,nline,""+txt7+" | ")
+ set_line(lines,nline,""+txt8+" | ")
+ set_line(lines,nline,""+txt9+" | ")
+ set_line(lines,nline,""+txt10+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+
+; last row, summary
+ set_line(lines,nline,row_header)
+
+ txt0 = "All_"+sprintf("%.0f", nstation)
+ txt1 = "-"
+ txt2 = "-"
+ txt3 = "-"
+ txt4 = M_energy_co2
+ txt5 = M_energy_rad
+ txt6 = M_energy_lh
+ txt7 = M_energy_sh
+ txt8 = M_energy_gpp
+ txt9 = M_energy_er
+ txt10 = M_energy_all
+
+ set_line(lines,nline,""+txt0+" | ")
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+ set_line(lines,nline,""+txt7+" | ")
+ set_line(lines,nline,""+txt8+" | ")
+ set_line(lines,nline,""+txt9+" | ")
+ set_line(lines,nline,""+txt10+" | ")
+
+ set_line(lines,nline,row_footer)
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+ delete (idx)
+
+;***************************************************************************
+; output plots
+;***************************************************************************
+ output_dir = model_name+"/ameriflux"
+
+ system("mv *.png *.html " + output_dir)
+;***************************************************************************
+end
+