diff -r 000000000000 -r 0c6405ab2ff4 ameriflux/25.plot.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ameriflux/25.plot.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,863 @@ +;************************************************************ +; change data name to timeseries_L4_m.nc +; move data +; 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 + +; for nee, gpp, and ar +; observed unit is gC/m2/day +; model unit is gC/m2/s +; to change to observed unit, + + factor_flux = 86400. + +; for incident solar radiation, +; observed Rg_f unit is MJ/m2/day +; model (FSDS) unit is W/m2 +; to change to model unit, + + factor_rad = 1.e6/86400. + +; model unit is gC/m2/s +; change to unit of observed (u mol/m2/s) +; 12. = molecular weight of C +; u mol = 1e-6 mol +; factor_flux = 1.e6 /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" \ + /) + + field = (/"NEE Flux" \ + ,"Shortwave Incoming" \ + ,"Latent Heat" \ + ,"Sensible Heat" \ + ,"GPP Flux" \ + ,"Respiration" \ + /) + + field_unit = (/"gC/m2/day" \ + ,"W/m2" \ + ,"W/m2" \ + ,"W/m2" \ + ,"gC/m2/day" \ + ,"gC/m2/day" \ + /) + + 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_data/observed/ameriflux/" + + year_station = new ((/nstation/),integer) ; number of year + year_ob = new ((/nstation/),string) ; observed 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 = "timeseries_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) + year_ob(n) = year_ob_i(n) + "-" + year_ob_f(n) + + 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|:))*factor_rad + + 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) + +;------------------------------------------------------------ +; 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_flux + data_mod(0,:,:,:) = yy(:,:,:) * factor_flux + +; printVarSummary(yy) + +; data = fm->NETRAD + data = fm->FSDS + 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_flux + data_mod(4,:,:,:) = yy(:,:,:) * factor_flux + + 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_flux + data_mod(5,:,:,:) = yy(:,:,:) * factor_flux + + 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_data/observed/ameriflux/" + + 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 = "timeseries_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 + + data_ob(1,:,:) = data_ob(1,:,:) * factor_rad + + 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 = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Site NameLatitudeLongitudeObserved
" + 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 = 1. + + 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 = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Site NameLatitudeLongitudeObserved
" + 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 = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Site NameLatitudeLongitudeObservedCO2 FluxNet RadiationLatent HeatSensible HeatGPP GluxRespirationAverage
" + 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 +