1 ;************************************************************
4 ; sort by latitude in decending order (N->S)
6 ; add year_ob_i and year_ob_f
8 ;************************************************************
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
12 ;************************************************************
13 procedure set_line(lines:string,nline:integer,newlines:string)
15 ; add line to ascci/html file
17 nnewlines = dimsizes(newlines)
18 if(nline+nnewlines-1.ge.dimsizes(lines))
19 print("set_line: bad index, not setting anything.")
22 lines(nline:nline+nnewlines-1) = newlines
23 ; print ("lines = " + lines(nline:nline+nnewlines-1))
24 nline = nline + nnewlines
27 ;*************************************************************
33 ;************************************************
35 ;************************************************
39 model_name = "i01.10" + model
41 ; change to unit of observed (u mol/m2/s)
42 ; Model_units [=] gC/m2/s
43 ; 12. = molecular weight of C
49 ;************************************************
51 ;************************************************
53 station = (/"ARM_Oklahoma" \
54 ,"ARM_Oklahoma_burn" \
55 ,"ARM_Oklahoma_control" \
63 ,"Duke_Forest_Hardwoods" \
64 ,"Duke_Forest_Open_Field" \
68 ,"Flagstaff_Managed" \
69 ,"Flagstaff_Unmanaged" \
70 ,"Flagstaff_Wildfire" \
72 ,"FreemanRanch_mesquite" \
75 ,"HarvardForestHemlock" \
76 ,"HowlandForestMain" \
77 ,"HowlandForestWest" \
79 ,"KendallGrasslands" \
80 ,"KennedySpaceCenterPine" \
81 ,"KennedySpaceCenterScrub" \
85 ,"Mead-irrigated-rotation" \
87 ,"Metolius_2nd_YoungPonderosaPine" \
89 ,"MetoliusIntermediatePine" \
90 ,"MetoliusOldPonderosaPine" \
101 ,"SkyOaks_PostFire" \
103 ,"SylvaniaWilderness" \
129 year_ob = (/"2003-2006" \
357 field = (/"NEE Flux" \
365 field_unit = (/"u mol/m2/s" \
373 nfield = dimsizes(field)
375 ;========================================================================
376 ; find (# of year observed) >=4 and year_ob_i <= 2001
378 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
380 nstation = dimsizes(station)
382 year_station = new ((/nstation/),integer)
386 diro = dir_root + station(n)+"/"
387 filo = year_ob_i(n)+"-"+year_ob_f(n)+"_L4_m.nc"
388 fo = addfile (diro+filo,"r")
391 year_station(n) = dimsizes(year)
398 i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
400 station_long = station(i_long_ob)
401 year_ob_long = year_ob(i_long_ob)
402 year_ob_i_long = year_ob_i(i_long_ob)
403 year_ob_f_long = year_ob_f(i_long_ob)
404 year_station_long = year_station(i_long_ob)
406 nstation_long = dimsizes(station_long)
407 ; print (nstation_long)
409 ; print (year_ob_i(i_long_ob))
410 ; print (station_long)
412 ;========================================================================
413 ; get long observed lat and lon
415 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
417 lat_ob_long = new ((/nstation_long/),float)
418 lon_ob_long = new ((/nstation_long/),float)
420 do n = 0,nstation_long-1
422 diro = dir_root + station_long(n)+"/"
423 filo = year_ob_i_long(n)+"-"+year_ob_f_long(n)+"_L4_m.nc"
424 fo = addfile (diro+filo,"r")
429 lon_ob_long(n) = fo->lon
430 lat_ob_long(n) = fo->lat
436 ;=========================================================
437 ; get model data at observed lat-lon
439 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
440 film = model_name+"_ameriflux_1990-2004_monthly.nc"
441 fm = addfile (dirm+film,"r")
447 date_dim = dimsizes(date)
452 data_mod = new ((/nfield,nyear,nmonth,nstation_long/),float)
454 ;************************************************************
455 ; interpolate model data into observed station
456 ; note: model is 0-360E, 90S-90N
457 ;************************************************************
459 ; to be able to handle observation at (-89.98,-24.80)
462 if (ENERGY .eq. "old") then
465 data_mod(0,:,:,:) = data(:,:,:) * factor
471 data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
476 ; data = fm->SENSIBLE
478 data_mod(3,:,:,:) = data(:,:,:)
483 data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:)
490 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob_long,lat_ob_long,0)
491 data_mod(0,:,:,:) = yy(:,:,:) * factor
494 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob_long,lat_ob_long,0)
495 data_mod(1,:,:,:) = yy(:,:,:)
498 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob_long,lat_ob_long,0)
499 data_mod(2,:,:,:) = yy(:,:,:)
501 ; data = fm->SENSIBLE
503 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob_long,lat_ob_long,0)
504 data_mod(3,:,:,:) = yy(:,:,:)
509 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob_long,lat_ob_long,0)
510 data_mod(4,:,:,:) = yy(:,:,:) * factor
512 if (model .eq. "cn") then
523 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob_long,lat_ob_long,0)
524 data_mod(5,:,:,:) = yy(:,:,:) * factor
529 ;*******************************************************************
530 ; for station line plot
531 ;*******************************************************************
533 ; for x-axis in xyplot
535 mon@long_name = "month"
537 res = True ; plot mods desired
538 res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
539 res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
541 res@tmXBFormat = "f" ; not to add trailing zeros
543 ;-------------------------------------------------------------------------
544 ; Add a boxed legend using the more simple method
546 res@pmLegendDisplayMode = "Always"
547 ; res@pmLegendWidthF = 0.1
548 res@pmLegendWidthF = 0.08
549 res@pmLegendHeightF = 0.06
550 ; res@pmLegendOrthogonalPosF = -1.17
551 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
552 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
554 ; res@pmLegendParallelPosF = 0.18
555 res@pmLegendParallelPosF = 0.23 ;(rightward)
557 ; res@lgPerimOn = False
558 res@lgLabelFontHeightF = 0.015
559 res@xyExplicitLegendLabels = (/"observed",model_name/)
560 ;-------------------------------------------------------------------
562 res@gsnFrame = False ; Do not draw plot
563 res@gsnDraw = False ; Do not advance frame
565 pres = True ; panel plot mods desired
566 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
567 ; indiv. plots in panel
568 pres@gsnMaximize = True ; fill the page
569 ;-------------------------------------------------------------------
571 ;==============================================================
572 ; get ob data at each site
574 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
576 do n = 0,nstation_long-1
581 nyear = year_station_long(n)
583 if (year_ob_f_long(n).eq. 2006) then
584 year_setback = 2006 -2004
586 if (year_ob_f_long(n).eq. 2005) then
587 year_setback = 2005 -2004
590 ntime = (nyear - year_setback) * nmonth
593 data_ob = new ((/nfield, nyear, nmonth/),float)
595 diro = dir_root + station_long(n)+"/"
596 filo = year_ob_i_long(n)+"-"+year_ob_f_long(n)+"_L4_m.nc"
597 fo = addfile (diro+filo,"r")
602 data_ob(0,:,:) = fo->NEE_or_fMDS
603 data_ob(1,:,:) = fo->Rg_f
604 data_ob(2,:,:) = fo->LE_f
605 data_ob(3,:,:) = fo->H_f
606 data_ob(4,:,:) = fo->GPP_or_MDS
607 data_ob(5,:,:) = fo->Reco_or
611 timeI = new((/ntime/),integer)
612 timeF = new((/ntime/),float)
613 timeI = ispan(1,ntime,1)
614 timeF = year_ob_i_long(n) + (timeI-1)/12.
615 timeF@long_name = "year"
617 plot_data = new((/2,ntime/),float)
619 ;----------------------------
622 plot_name = station_long(n)+"_tseries_vs_ob"
623 title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"
624 res@tiMainString = title
626 wks = gsn_open_wks (plot_type,plot_name)
627 plot=new(nfield,graphic) ; create graphic array
629 i_year_mod_i = year_ob_i_long(n) - 1990
630 i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
632 i_year_ob_f = nyear - year_setback - 1
635 ; print (i_year_ob_f)
636 ; print (i_year_mod_i)
637 ; print (i_year_mod_f)
640 plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
641 plot_data(1,:) = ndtooned(data_mod(i,i_year_mod_i:i_year_mod_f,:,n))
642 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
643 plot(i)=gsn_csm_xy(wks,timeF,plot_data,res) ; create plot
646 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
648 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
649 "rm "+plot_name+"."+plot_type)
661 ;###################################################################
662 ; for the following tables,
663 ; sort by latitude in decending order (N->S)
665 isort = dim_pqsort(lat_ob_long,-1)
667 station_sort = station_long(isort)
668 year_ob_sort = year_ob_long(isort)
669 lat_ob_sort = lat_ob_long(isort)
670 lon_ob_sort = lon_ob_long(isort)
675 ;*******************************************************************
676 ; html table of site: observed
677 ;*******************************************************************
678 output_html = "tseries_vs_ob.html"
680 header = (/"<HTML>" \
682 ,"<TITLE>CLAMP metrics</TITLE>" \
684 ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \
689 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
691 ," <th bgcolor=DDDDDD >Site Name</th>" \
692 ," <th bgcolor=DDDDDD >Latitude</th>" \
693 ," <th bgcolor=DDDDDD >Longitude</th>" \
694 ," <th bgcolor=DDDDDD >Observed</th>" \
697 table_footer = "</table>"
701 lines = new(50000,string)
704 set_line(lines,nline,header)
705 set_line(lines,nline,table_header)
706 ;-----------------------------------------------
709 do n = 0,nstation_long-1
711 set_line(lines,nline,row_header)
713 txt0 = station_sort(n)
714 txt1 = sprintf("%5.2f", lat_ob_sort(n))
715 txt2 = sprintf("%5.2f", lon_ob_sort(n))
716 txt3 = year_ob_sort(n)
718 set_line(lines,nline,"<th><a href="+txt0+"_tseries_vs_ob.png>"+txt0+"</a></th>")
719 set_line(lines,nline,"<th>"+txt1+"</th>")
720 set_line(lines,nline,"<th>"+txt2+"</th>")
721 set_line(lines,nline,"<th>"+txt3+"</th>")
723 set_line(lines,nline,row_footer)
725 ;-----------------------------------------------
726 set_line(lines,nline,table_footer)
727 set_line(lines,nline,footer)
729 ; Now write to an HTML file.
730 idx = ind(.not.ismissing(lines))
731 if(.not.any(ismissing(idx))) then
732 asciiwrite(output_html,lines(idx))
738 ;***************************************************************************
740 ;***************************************************************************
741 output_dir = model_name+"/ameriflux"
743 system("mv *.png *.html " + output_dir)
744 ;***************************************************************************