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
43 ;************************************************
45 ;************************************************
47 station = (/"ARM_Oklahoma" \
48 ,"ARM_Oklahoma_burn" \
49 ,"ARM_Oklahoma_control" \
57 ,"Duke_Forest_Hardwoods" \
58 ,"Duke_Forest_Open_Field" \
62 ,"Flagstaff_Managed" \
63 ,"Flagstaff_Unmanaged" \
64 ,"Flagstaff_Wildfire" \
66 ,"FreemanRanch_mesquite" \
69 ,"HarvardForestHemlock" \
70 ,"HowlandForestMain" \
71 ,"HowlandForestWest" \
73 ,"KendallGrasslands" \
74 ,"KennedySpaceCenterPine" \
75 ,"KennedySpaceCenterScrub" \
79 ,"Mead-irrigated-rotation" \
81 ,"Metolius_2nd_YoungPonderosaPine" \
83 ,"MetoliusIntermediatePine" \
84 ,"MetoliusOldPonderosaPine" \
97 ,"SylvaniaWilderness" \
123 year_ob = (/"2003-2006" \
199 field = (/"NEE Flux" \
207 field_unit = (/"u mol/m2/s" \
215 nstation = dimsizes(station)
216 nfield = dimsizes(field)
218 ;========================================================================
219 ; get observed info: number of year, first/last year, lat, lon
221 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
223 year_station = new ((/nstation/),integer) ; number of year
224 year_ob_i = new ((/nstation/),integer) ; first year
225 year_ob_f = new ((/nstation/),integer) ; last year
226 lat_ob = new ((/nstation/),float) ; latitude
227 lon_ob = new ((/nstation/),float) ; longitude
231 diro = dir_root + station(n)+"/"
232 filo = year_ob(n)+"_L4_m.nc"
233 fo = addfile (diro+filo,"r")
240 year_station(n) = dimsizes(year)
241 year_ob_i(n) = year(0)
242 year_ob_f(n) = year(year_station(n)-1)
252 ;--------------------------------------------------------------
253 ; find (# of year observed) >=4 and year_ob_i <= 2001
255 i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
257 station_long = station(i_long_ob)
258 lat_ob_long = lat_ob(i_long_ob)
259 lon_ob_long = lat_ob(i_long_ob)
260 year_ob_long = year_ob(i_long_ob)
261 year_ob_i_long = year_ob_i(i_long_ob)
262 year_ob_f_long = year_ob_f(i_long_ob)
263 year_station_long = year_station(i_long_ob)
265 nstation_long = dimsizes(station_long)
268 ; print (nstation_long)
269 ; print (station_long)
270 ; print (year_ob_i(i_long_ob))
272 ;=========================================================
273 ; get model data at observed lat-lon
275 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
276 film = model_name+"_ameriflux_1990-2004_monthly.nc"
277 fm = addfile (dirm+film,"r")
283 date_dim = dimsizes(date)
288 data_mod = new ((/nfield,nyear,nmonth,nstation/),float)
289 data_mod_ann = new ((/nfield,nmonth,nstation/),float)
290 data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)
292 ; change to unit of observed (u mol/m2/s)
293 ; Model_units [=] gC/m2/s
294 ; 12. = molecular weight of C
299 ;------------------------------------------------------------
300 ; interpolate model data into observed station
301 ; note: model is 0-360E, 90S-90N
303 ; to be able to handle observation at (-89.98,-24.80)
305 ;------------------------------------------------------------
307 if (ENERGY .eq. "old") then
310 data_mod(0,:,:,:) = data(:,:,:) * factor
316 data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
321 ; data = fm->SENSIBLE
323 data_mod(3,:,:,:) = data(:,:,:)
328 data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:)
335 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
336 data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
337 data_mod(0,:,:,:) = yy(:,:,:) * factor
339 ; printVarSummary(yy)
342 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
343 data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))
344 data_mod(1,:,:,:) = yy(:,:,:)
347 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
348 data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))
349 data_mod(2,:,:,:) = yy(:,:,:)
351 ; data = fm->SENSIBLE
353 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
354 data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))
355 data_mod(3,:,:,:) = yy(:,:,:)
360 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
361 data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
362 data_mod(4,:,:,:) = yy(:,:,:) * factor
364 if (model .eq. "cn") then
375 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
376 data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
377 data_mod(5,:,:,:) = yy(:,:,:) * factor
379 data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)
386 ;*******************************************************************
387 ; for station line plot
388 ;*******************************************************************
390 ; for x-axis in xyplot
392 mon@long_name = "month"
394 res = True ; plot mods desired
395 res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
396 res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
398 res@tmXBFormat = "f" ; not to add trailing zeros
400 ;-------------------------------------------------------------------------
401 ; Add a boxed legend using the more simple method
403 res@pmLegendDisplayMode = "Always"
404 ; res@pmLegendWidthF = 0.1
405 res@pmLegendWidthF = 0.08
406 res@pmLegendHeightF = 0.06
407 ; res@pmLegendOrthogonalPosF = -1.17
408 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
409 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
411 ; res@pmLegendParallelPosF = 0.18
412 res@pmLegendParallelPosF = 0.23 ;(rightward)
414 ; res@lgPerimOn = False
415 res@lgLabelFontHeightF = 0.015
416 res@xyExplicitLegendLabels = (/"observed",model_name/)
417 ;-------------------------------------------------------------------
419 res@gsnFrame = False ; Do not draw plot
420 res@gsnDraw = False ; Do not advance frame
422 pres = True ; panel plot mods desired
423 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
424 ; indiv. plots in panel
425 pres@gsnMaximize = True ; fill the page
426 ;-------------------------------------------------------------------
428 ;==============================================================
429 ; get ob data at each site with long observation
431 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
433 do n = 0,nstation_long-1
438 nyear = year_station_long(n)
440 if (year_ob_f_long(n).eq. 2006) then
441 year_setback = 2006 -2004
443 if (year_ob_f_long(n).eq. 2005) then
444 year_setback = 2005 -2004
447 ntime = (nyear - year_setback) * nmonth
450 data_ob = new ((/nfield, nyear, nmonth/),float)
452 diro = dir_root + station_long(n)+"/"
453 filo = year_ob_long(n)+"_L4_m.nc"
454 fo = addfile (diro+filo,"r")
459 data_ob(0,:,:) = fo->NEE_or_fMDS
460 data_ob(1,:,:) = fo->Rg_f
461 data_ob(2,:,:) = fo->LE_f
462 data_ob(3,:,:) = fo->H_f
463 data_ob(4,:,:) = fo->GPP_or_MDS
464 data_ob(5,:,:) = fo->Reco_or
468 timeI = new((/ntime/),integer)
469 timeF = new((/ntime/),float)
470 timeI = ispan(1,ntime,1)
471 timeF = year_ob_i_long(n) + (timeI-1)/12.
472 timeF@long_name = "year"
474 plot_data = new((/2,ntime/),float)
476 ;----------------------------
479 plot_name = station_long(n)+"_tseries_vs_ob"
480 title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"
481 res@tiMainString = title
483 wks = gsn_open_wks (plot_type,plot_name)
484 plot=new(nfield,graphic) ; create graphic array
486 i_year_mod_i = year_ob_i_long(n) - 1990
487 i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
489 i_year_ob_f = nyear - year_setback - 1
492 ; print (i_year_ob_f)
493 ; print (i_year_mod_i)
494 ; print (i_year_mod_f)
497 plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
498 plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))
499 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
500 plot(i)=gsn_csm_xy(wks,timeF,plot_data,res) ; create plot
503 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
505 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
506 "rm "+plot_name+"."+plot_type)
518 ;###################################################################
519 ; for the following tables,
520 ; sort by latitude in decending order (N->S)
522 isort = dim_pqsort(lat_ob_long,-1)
524 station_sort = station_long(isort)
525 year_ob_sort = year_ob_long(isort)
526 lat_ob_sort = lat_ob_long(isort)
527 lon_ob_sort = lon_ob_long(isort)
532 ;*******************************************************************
533 ; html table of site: observed
534 ;*******************************************************************
535 output_html = "tseries_vs_ob.html"
537 header = (/"<HTML>" \
539 ,"<TITLE>CLAMP metrics</TITLE>" \
541 ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \
546 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
548 ," <th bgcolor=DDDDDD >Site Name</th>" \
549 ," <th bgcolor=DDDDDD >Latitude</th>" \
550 ," <th bgcolor=DDDDDD >Longitude</th>" \
551 ," <th bgcolor=DDDDDD >Observed</th>" \
554 table_footer = "</table>"
558 lines = new(50000,string)
561 set_line(lines,nline,header)
562 set_line(lines,nline,table_header)
563 ;-----------------------------------------------
566 do n = 0,nstation_long-1
568 set_line(lines,nline,row_header)
570 txt0 = station_sort(n)
571 txt1 = sprintf("%5.2f", lat_ob_sort(n))
572 txt2 = sprintf("%5.2f", lon_ob_sort(n))
573 txt3 = year_ob_sort(n)
575 set_line(lines,nline,"<th><a href="+txt0+"_tseries_vs_ob.png>"+txt0+"</a></th>")
576 set_line(lines,nline,"<th>"+txt1+"</th>")
577 set_line(lines,nline,"<th>"+txt2+"</th>")
578 set_line(lines,nline,"<th>"+txt3+"</th>")
580 set_line(lines,nline,row_footer)
582 ;-----------------------------------------------
583 set_line(lines,nline,table_footer)
584 set_line(lines,nline,footer)
586 ; Now write to an HTML file.
587 idx = ind(.not.ismissing(lines))
588 if(.not.any(ismissing(idx))) then
589 asciiwrite(output_html,lines(idx))
595 ;***************************************************************************
597 ;***************************************************************************
598 output_dir = model_name+"/ameriflux"
600 system("mv *.png *.html " + output_dir)
601 ;***************************************************************************