1 ;************************************************************
5 ; sort by latitude in decending order (N->S)
7 ; add year_ob_i and year_ob_f
9 ;************************************************************
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
12 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
13 ;************************************************************
14 procedure set_line(lines:string,nline:integer,newlines:string)
16 ; add line to ascci/html file
18 nnewlines = dimsizes(newlines)
19 if(nline+nnewlines-1.ge.dimsizes(lines))
20 print("set_line: bad index, not setting anything.")
23 lines(nline:nline+nnewlines-1) = newlines
24 ; print ("lines = " + lines(nline:nline+nnewlines-1))
25 nline = nline + nnewlines
28 ;*************************************************************
34 ;************************************************
36 ;************************************************
40 model_name = "i01.10" + model
46 ;************************************************
48 ;************************************************
50 station = (/"ARM_Oklahoma" \
51 ,"ARM_Oklahoma_burn" \
52 ,"ARM_Oklahoma_control" \
60 ,"Duke_Forest_Hardwoods" \
61 ,"Duke_Forest_Open_Field" \
65 ,"Flagstaff_Managed" \
66 ,"Flagstaff_Unmanaged" \
67 ,"Flagstaff_Wildfire" \
69 ,"FreemanRanch_mesquite" \
72 ,"HarvardForestHemlock" \
73 ,"HowlandForestMain" \
74 ,"HowlandForestWest" \
76 ,"KendallGrasslands" \
77 ,"KennedySpaceCenterPine" \
78 ,"KennedySpaceCenterScrub" \
82 ,"Mead-irrigated-rotation" \
84 ,"Metolius_2nd_YoungPonderosaPine" \
86 ,"MetoliusIntermediatePine" \
87 ,"MetoliusOldPonderosaPine" \
100 ,"SylvaniaWilderness" \
126 year_ob = (/"2003-2006" \
202 field = (/"NEE Flux" \
210 field_unit = (/"u mol/m2/s" \
218 nstation = dimsizes(station)
219 nfield = dimsizes(field)
221 ;========================================================================
222 ; get observed info: number of year, first/last year, lat, lon
225 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
227 year_station = new ((/nstation/),integer) ; number of year
228 year_ob_i = new ((/nstation/),integer) ; first year
229 year_ob_f = new ((/nstation/),integer) ; last year
230 lat_ob = new ((/nstation/),float) ; latitude
231 lon_ob = new ((/nstation/),float) ; longitude
233 data_ob_ann = new ((/nfield, nmonth, nstation/),float)
237 diro = dir_root + station(n)+"/"
238 filo = year_ob(n)+"_L4_m.nc"
239 fo = addfile (diro+filo,"r")
246 year_station(n) = dimsizes(year)
247 year_ob_i(n) = year(0)
248 year_ob_f(n) = year(year_station(n)-1)
252 data = fo->NEE_or_fMDS
253 data_ob_ann(0,:,n) = dim_avg(data(month|:,year|:))
256 data_ob_ann(1,:,n) = dim_avg(data(month|:,year|:))
259 data_ob_ann(2,:,n) = dim_avg(data(month|:,year|:))
262 data_ob_ann(3,:,n) = dim_avg(data(month|:,year|:))
264 data = fo->GPP_or_MDS
265 data_ob_ann(4,:,n) = dim_avg(data(month|:,year|:))
268 data_ob_ann(5,:,n) = dim_avg(data(month|:,year|:))
277 ;--------------------------------------------------------------
278 ; find (# of year observed) >=4 and year_ob_i <= 2001
280 i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
282 station_long = station(i_long_ob)
283 lat_ob_long = lat_ob(i_long_ob)
284 lon_ob_long = lat_ob(i_long_ob)
285 year_ob_long = year_ob(i_long_ob)
286 year_ob_i_long = year_ob_i(i_long_ob)
287 year_ob_f_long = year_ob_f(i_long_ob)
288 year_station_long = year_station(i_long_ob)
290 nstation_long = dimsizes(station_long)
293 ; print (nstation_long)
294 ; print (station_long)
295 ; print (year_ob_i(i_long_ob))
297 ;=========================================================
298 ; get model data at observed lat-lon
300 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
301 film = model_name+"_ameriflux_1990-2004_monthly.nc"
302 fm = addfile (dirm+film,"r")
308 date_dim = dimsizes(date)
311 data_mod = new ((/nfield,nyear,nmonth,nstation/),float)
312 data_mod_ann = new ((/nfield,nmonth,nstation/),float)
313 data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)
315 ; change to unit of observed (u mol/m2/s)
316 ; Model_units [=] gC/m2/s
317 ; 12. = molecular weight of C
322 ;------------------------------------------------------------
323 ; interpolate model data into observed station
324 ; note: model is 0-360E, 90S-90N
326 ; to be able to handle observation at (-89.98,-24.80)
328 ;------------------------------------------------------------
330 if (ENERGY .eq. "old") then
333 data_mod(0,:,:,:) = data(:,:,:) * factor
339 data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
344 ; data = fm->SENSIBLE
346 data_mod(3,:,:,:) = data(:,:,:)
351 data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:)
358 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
359 data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
360 data_mod(0,:,:,:) = yy(:,:,:) * factor
362 ; printVarSummary(yy)
365 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
366 data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))
367 data_mod(1,:,:,:) = yy(:,:,:)
370 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
371 data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))
372 data_mod(2,:,:,:) = yy(:,:,:)
374 ; data = fm->SENSIBLE
376 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
377 data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))
378 data_mod(3,:,:,:) = yy(:,:,:)
383 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
384 data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
385 data_mod(4,:,:,:) = yy(:,:,:) * factor
387 if (model .eq. "cn") then
398 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
399 data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
400 data_mod(5,:,:,:) = yy(:,:,:) * factor
402 data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)
409 ;*******************************************************************
410 ; for station line plot
411 ;*******************************************************************
413 ; for x-axis in xyplot
415 mon@long_name = "month"
417 res = True ; plot mods desired
418 res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
419 res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
421 res@tmXBFormat = "f" ; not to add trailing zeros
423 ;-------------------------------------------------------------------------
424 ; Add a boxed legend using the more simple method
426 res@pmLegendDisplayMode = "Always"
427 ; res@pmLegendWidthF = 0.1
428 res@pmLegendWidthF = 0.08
429 res@pmLegendHeightF = 0.06
430 ; res@pmLegendOrthogonalPosF = -1.17
431 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
432 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
434 ; res@pmLegendParallelPosF = 0.18
435 res@pmLegendParallelPosF = 0.23 ;(rightward)
437 ; res@lgPerimOn = False
438 res@lgLabelFontHeightF = 0.015
439 res@xyExplicitLegendLabels = (/"observed",model_name/)
440 ;-------------------------------------------------------------------
442 res@gsnFrame = False ; Do not draw plot
443 res@gsnDraw = False ; Do not advance frame
445 pres = True ; panel plot mods desired
446 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
447 ; indiv. plots in panel
448 pres@gsnMaximize = True ; fill the page
449 ;-------------------------------------------------------------------
451 ;==============================================================
452 ; get ob data at each site with long observation
454 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
456 do n = 0,nstation_long-1
461 nyear = year_station_long(n)
463 if (year_ob_f_long(n).eq. 2006) then
464 year_setback = 2006 -2004
466 if (year_ob_f_long(n).eq. 2005) then
467 year_setback = 2005 -2004
470 ntime = (nyear - year_setback) * nmonth
472 data_ob = new ((/nfield, nyear, nmonth/),float)
474 diro = dir_root + station_long(n)+"/"
475 filo = year_ob_long(n)+"_L4_m.nc"
476 fo = addfile (diro+filo,"r")
481 data_ob(0,:,:) = fo->NEE_or_fMDS
482 data_ob(1,:,:) = fo->Rg_f
483 data_ob(2,:,:) = fo->LE_f
484 data_ob(3,:,:) = fo->H_f
485 data_ob(4,:,:) = fo->GPP_or_MDS
486 data_ob(5,:,:) = fo->Reco_or
490 timeI = new((/ntime/),integer)
491 timeF = new((/ntime/),float)
492 timeI = ispan(1,ntime,1)
493 timeF = year_ob_i_long(n) + (timeI-1)/12.
494 timeF@long_name = "year"
496 plot_data = new((/2,ntime/),float)
498 ;----------------------------
501 plot_name = station_long(n)+"_tseries_vs_ob"
502 title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"
503 res@tiMainString = title
505 wks = gsn_open_wks (plot_type,plot_name)
506 plot=new(nfield,graphic) ; create graphic array
508 i_year_mod_i = year_ob_i_long(n) - 1990
509 i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
511 i_year_ob_f = nyear - year_setback - 1
514 ; print (i_year_ob_f)
515 ; print (i_year_mod_i)
516 ; print (i_year_mod_f)
519 plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
520 plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))
521 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
522 plot(i)=gsn_csm_xy(wks,timeF,plot_data,res) ; create plot
525 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
527 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
528 "rm "+plot_name+"."+plot_type)
540 ;###################################################################
541 ; for the following tables,
542 ; sort by latitude in decending order (N->S)
544 isort = dim_pqsort(lat_ob_long,-1)
546 station_sort = station_long(isort)
547 year_ob_sort = year_ob_long(isort)
548 lat_ob_sort = lat_ob_long(isort)
549 lon_ob_sort = lon_ob_long(isort)
554 ;*******************************************************************
555 ; html table of site: observed
556 ;*******************************************************************
557 output_html = "tseries_vs_ob.html"
559 header = (/"<HTML>" \
561 ,"<TITLE>CLAMP metrics</TITLE>" \
563 ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \
568 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
570 ," <th bgcolor=DDDDDD >Site Name</th>" \
571 ," <th bgcolor=DDDDDD >Latitude</th>" \
572 ," <th bgcolor=DDDDDD >Longitude</th>" \
573 ," <th bgcolor=DDDDDD >Observed</th>" \
576 table_footer = "</table>"
580 lines = new(50000,string)
583 set_line(lines,nline,header)
584 set_line(lines,nline,table_header)
585 ;-----------------------------------------------
588 do n = 0,nstation_long-1
590 set_line(lines,nline,row_header)
592 txt0 = station_sort(n)
593 txt1 = sprintf("%5.2f", lat_ob_sort(n))
594 txt2 = sprintf("%5.2f", lon_ob_sort(n))
595 txt3 = year_ob_sort(n)
597 set_line(lines,nline,"<th><a href="+txt0+"_tseries_vs_ob.png>"+txt0+"</a></th>")
598 set_line(lines,nline,"<th>"+txt1+"</th>")
599 set_line(lines,nline,"<th>"+txt2+"</th>")
600 set_line(lines,nline,"<th>"+txt3+"</th>")
602 set_line(lines,nline,row_footer)
604 ;-----------------------------------------------
605 set_line(lines,nline,table_footer)
606 set_line(lines,nline,footer)
608 ; Now write to an HTML file.
609 idx = ind(.not.ismissing(lines))
610 if(.not.any(ismissing(idx))) then
611 asciiwrite(output_html,lines(idx))
618 delete (station_sort)
619 delete (year_ob_sort)
623 ;************************************************************
624 ; compute annual cycle correlation coef and M score
625 ;************************************************************
629 ccr = new ((/nstation, nfield/),float)
630 M_score = new ((/nstation, nfield/),float)
634 ccr(n,m) = esccr(data_ob_ann(m,:,n),data_mod_ann(m,:,n),0)
635 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))))
636 M_score(n,m) = (1. -(bias/nmonth)) * score_max
640 M_co2 = avg(M_score(:,0))
641 M_rad = avg(M_score(:,1))
642 M_lh = avg(M_score(:,2))
643 M_sh = avg(M_score(:,3))
644 M_gpp = avg(M_score(:,4))
645 M_er = avg(M_score(:,5))
646 M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er
648 M_energy_co2 = sprintf("%.2f", M_co2)
649 M_energy_rad = sprintf("%.2f", M_rad)
650 M_energy_lh = sprintf("%.2f", M_lh )
651 M_energy_sh = sprintf("%.2f", M_sh )
652 M_energy_gpp = sprintf("%.2f", M_gpp)
653 M_energy_er = sprintf("%.2f", M_er )
654 M_energy_all = sprintf("%.2f", M_all)
656 ;*******************************************************************
657 ; for station line plot
658 ;*******************************************************************
660 ; for x-axis in xyplot
662 mon@long_name = "month"
664 ;-------------------------------------------------------------------
666 plot_data = new((/2,nmonth/),float)
668 plot_data!1 = "month"
671 ;----------------------------
674 plot_name = station(n)+"_ob"
675 title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
676 res@tiMainString = title
678 wks = gsn_open_wks (plot_type,plot_name)
679 plot=new(nfield,graphic) ; create graphic array
682 plot_data(0,:) = (/data_ob_ann(i,:,n)/)
683 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
684 plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot
687 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
689 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
690 "rm "+plot_name+"."+plot_type)
694 ;----------------------------
697 plot_name = station(n)+"_model_vs_ob"
698 title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
699 res@tiMainString = title
701 wks = gsn_open_wks (plot_type,plot_name)
702 plot=new(nfield,graphic) ; create graphic array
705 plot_data(0,:) = (/data_ob_ann(i,:,n)/)
706 plot_data(1,:) = (/data_mod_ann(i,:,n)/)
707 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
708 plot(i)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot
711 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
713 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
714 "rm "+plot_name+"."+plot_type)
720 ;###################################################################
721 ; for the following tables,
722 ; sort by latitude in decending order (N->S)
723 ; sort by lat in decending order (N->S)
725 isort = dim_pqsort(lat_ob,-1)
727 station_sort = station(isort)
728 year_ob_sort = year_ob(isort)
729 lat_ob_sort = lat_ob(isort)
730 lon_ob_sort = lon_ob(isort)
731 M_score_sort = M_score(isort,:)
735 ;###################################################################
736 ;*******************************************************************
737 ; html table of site: observed
738 ;*******************************************************************
739 output_html = "line_ob.html"
741 header = (/"<HTML>" \
743 ,"<TITLE>CLAMP metrics</TITLE>" \
745 ,"<H1>Energy at Site: Observation</H1>" \
750 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
752 ," <th bgcolor=DDDDDD >Site Name</th>" \
753 ," <th bgcolor=DDDDDD >Latitude</th>" \
754 ," <th bgcolor=DDDDDD >Longitude</th>" \
755 ," <th bgcolor=DDDDDD >Observed</th>" \
758 table_footer = "</table>"
762 lines = new(50000,string)
765 set_line(lines,nline,header)
766 set_line(lines,nline,table_header)
767 ;-----------------------------------------------
771 set_line(lines,nline,row_header)
773 txt0 = station_sort(n)
774 txt1 = sprintf("%5.2f", lat_ob_sort(n))
775 txt2 = sprintf("%5.2f", lon_ob_sort(n))
776 txt3 = year_ob_sort(n)
778 set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")
779 set_line(lines,nline,"<th>"+txt1+"</th>")
780 set_line(lines,nline,"<th>"+txt2+"</th>")
781 set_line(lines,nline,"<th>"+txt3+"</th>")
783 set_line(lines,nline,row_footer)
785 ;-----------------------------------------------
786 set_line(lines,nline,table_footer)
787 set_line(lines,nline,footer)
789 ; Now write to an HTML file.
790 idx = ind(.not.ismissing(lines))
791 if(.not.any(ismissing(idx))) then
792 asciiwrite(output_html,lines(idx))
798 ;*******************************************************************
799 ; score and line table : model vs observed
800 ;*******************************************************************
801 output_html = "score+line_vs_ob.html"
803 header = (/"<HTML>" \
805 ,"<TITLE>CLAMP metrics</TITLE>" \
807 ,"<H1>Energy at Site: Model "+model_name+"</H1>" \
811 delete (table_header)
813 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
815 ," <th bgcolor=DDDDDD >Site Name</th>" \
816 ," <th bgcolor=DDDDDD >Latitude</th>" \
817 ," <th bgcolor=DDDDDD >Longitude</th>" \
818 ," <th bgcolor=DDDDDD >Observed</th>" \
819 ," <th bgcolor=DDDDDD >CO2 Flux</th>" \
820 ," <th bgcolor=DDDDDD >Net Radiation</th>" \
821 ," <th bgcolor=DDDDDD >Latent Heat</th>" \
822 ," <th bgcolor=DDDDDD >Sensible Heat</th>" \
823 ," <th bgcolor=DDDDDD >GPP Glux</th>" \
824 ," <th bgcolor=DDDDDD >Respiration</th>" \
825 ," <th bgcolor=DDDDDD >Average</th>" \
828 table_footer = "</table>"
832 lines = new(50000,string)
835 set_line(lines,nline,header)
836 set_line(lines,nline,table_header)
837 ;-----------------------------------------------
841 set_line(lines,nline,row_header)
843 txt0 = station_sort(n)
844 txt1 = sprintf("%5.2f", lat_ob_sort(n))
845 txt2 = sprintf("%5.2f", lon_ob_sort(n))
846 txt3 = year_ob_sort(n)
847 txt4 = sprintf("%5.2f", M_score_sort(n,0))
848 txt5 = sprintf("%5.2f", M_score_sort(n,1))
849 txt6 = sprintf("%5.2f", M_score_sort(n,2))
850 txt7 = sprintf("%5.2f", M_score_sort(n,3))
851 txt8 = sprintf("%5.2f", M_score_sort(n,4))
852 txt9 = sprintf("%5.2f", M_score_sort(n,5))
853 txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))
855 set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")
856 set_line(lines,nline,"<th>"+txt1+"</th>")
857 set_line(lines,nline,"<th>"+txt2+"</th>")
858 set_line(lines,nline,"<th>"+txt3+"</th>")
859 set_line(lines,nline,"<th>"+txt4+"</th>")
860 set_line(lines,nline,"<th>"+txt5+"</th>")
861 set_line(lines,nline,"<th>"+txt6+"</th>")
862 set_line(lines,nline,"<th>"+txt7+"</th>")
863 set_line(lines,nline,"<th>"+txt8+"</th>")
864 set_line(lines,nline,"<th>"+txt9+"</th>")
865 set_line(lines,nline,"<th>"+txt10+"</th>")
867 set_line(lines,nline,row_footer)
871 set_line(lines,nline,row_header)
873 txt0 = "All_"+sprintf("%.0f", nstation)
885 set_line(lines,nline,"<th>"+txt0+"</th>")
886 set_line(lines,nline,"<th>"+txt1+"</th>")
887 set_line(lines,nline,"<th>"+txt2+"</th>")
888 set_line(lines,nline,"<th>"+txt3+"</th>")
889 set_line(lines,nline,"<th>"+txt4+"</th>")
890 set_line(lines,nline,"<th>"+txt5+"</th>")
891 set_line(lines,nline,"<th>"+txt6+"</th>")
892 set_line(lines,nline,"<th>"+txt7+"</th>")
893 set_line(lines,nline,"<th>"+txt8+"</th>")
894 set_line(lines,nline,"<th>"+txt9+"</th>")
895 set_line(lines,nline,"<th>"+txt10+"</th>")
897 set_line(lines,nline,row_footer)
898 ;-----------------------------------------------
899 set_line(lines,nline,table_footer)
900 set_line(lines,nline,footer)
902 ; Now write to an HTML file.
903 idx = ind(.not.ismissing(lines))
904 if(.not.any(ismissing(idx))) then
905 asciiwrite(output_html,lines(idx))
911 ;***************************************************************************
913 ;***************************************************************************
914 output_dir = model_name+"/ameriflux"
916 system("mv *.png *.html " + output_dir)
917 ;***************************************************************************