forrest@0: ;************************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;************************************************************
forrest@0: procedure set_line(lines:string,nline:integer,newlines:string)
forrest@0: begin
forrest@0: ; add line to ascci/html file
forrest@0:
forrest@0: nnewlines = dimsizes(newlines)
forrest@0: if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0: print("set_line: bad index, not setting anything.")
forrest@0: return
forrest@0: end if
forrest@0: lines(nline:nline+nnewlines-1) = newlines
forrest@0: ; print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0: nline = nline + nnewlines
forrest@0: return
forrest@0: end
forrest@0: ;*************************************************************
forrest@0: begin
forrest@0:
forrest@0: plot_type = "ps"
forrest@0: plot_type_new = "png"
forrest@0:
forrest@0: ;------------------------------------------------------
forrest@0: ; edit table.html of current model for movel1_vs_model2
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: html_name2 = compare+"/table.html"
forrest@0: html_new2 = html_name2 +".new"
forrest@0: end if
forrest@0:
forrest@0: ;------------------------------------------------------
forrest@0: ; edit table.html for current model
forrest@0:
forrest@0: html_name = model_name+"/table.html"
forrest@0: html_new = html_name +".new"
forrest@0:
forrest@0: ;------------------------------------------------------
forrest@0:
forrest@0: nmonth = 12
forrest@0:
forrest@0: ; for nee, gpp, and ar
forrest@0: ; observed unit is gC/m2/day
forrest@0: ; model unit is gC/m2/s
forrest@0: ; to change to observed unit,
forrest@0:
forrest@0: factor_flux = 86400.
forrest@0:
forrest@0: ; for incident solar radiation,
forrest@0: ; observed Rg_f unit is MJ/m2/day
forrest@0: ; model (FSDS) unit is W/m2
forrest@0: ; to change to model unit,
forrest@0:
forrest@0: factor_rad = 1.e6/86400.
forrest@0:
forrest@0: ;************************************************
forrest@0: ; observed data info
forrest@0: ;************************************************
forrest@0:
forrest@0: station = (/"ARM_Oklahoma" \
forrest@0: ,"ARM_Oklahoma_burn" \
forrest@0: ,"ARM_Oklahoma_control" \
forrest@0: ,"Atqasuk" \
forrest@0: ,"Audubon" \
forrest@0: ,"AustinCary" \
forrest@0: ,"Bartlett" \
forrest@0: ,"Bondville" \
forrest@0: ,"Brookings" \
forrest@0: ,"Donaldson" \
forrest@0: ,"Duke_Forest_Hardwoods" \
forrest@0: ,"Duke_Forest_Open_Field" \
forrest@0: ,"Duke_Forest_Pine" \
forrest@0: ,"Fermi_Ag" \
forrest@0: ,"Fermi_Prairie" \
forrest@0: ,"Flagstaff_Managed" \
forrest@0: ,"Flagstaff_Unmanaged" \
forrest@0: ,"Flagstaff_Wildfire" \
forrest@0: ,"FortPeck" \
forrest@0: ,"FreemanRanch_mesquite" \
forrest@0: ,"Goodwin_Creek" \
forrest@0: ,"HarvardForest" \
forrest@0: ,"HarvardForestHemlock" \
forrest@0: ,"HowlandForestMain" \
forrest@0: ,"HowlandForestWest" \
forrest@0: ,"Ivotuk" \
forrest@0: ,"KendallGrasslands" \
forrest@0: ,"KennedySpaceCenterPine" \
forrest@0: ,"KennedySpaceCenterScrub" \
forrest@0: ,"LittleProspect" \
forrest@0: ,"LostCreek" \
forrest@0: ,"Mead-irrigated" \
forrest@0: ,"Mead-irrigated-rotation" \
forrest@0: ,"Mead-rainfed" \
forrest@0: ,"Metolius_2nd_YoungPonderosaPine" \
forrest@0: ,"MetoliusEyerly" \
forrest@0: ,"MetoliusIntermediatePine" \
forrest@0: ,"MetoliusOldPonderosaPine" \
forrest@0: ,"MissouriOzark" \
forrest@0: ,"Mize" \
forrest@0: ,"MorganMonroe" \
forrest@0: ,"NiwotRidge" \
forrest@0: ,"NorthCarolina_cc" \
forrest@0: ,"NorthCarolina_lp" \
forrest@0: ,"ParkFalls" \
forrest@0: ,"Rayonier" \
forrest@0: ,"SantaRita" \
forrest@0: ,"SkyOaks_Old" \
forrest@0: ,"SkyOaks_PostFire" \
forrest@0: ,"SkyOaks_Young" \
forrest@0: ,"SylvaniaWilderness" \
forrest@0: ,"Toledo" \
forrest@0: ,"Tonzi" \
forrest@0: ,"UCI_1850" \
forrest@0: ,"UCI_1930" \
forrest@0: ,"UCI_1964" \
forrest@0: ,"UCI_1964wet" \
forrest@0: ,"UCI_1981" \
forrest@0: ,"UCI_1989" \
forrest@0: ,"UCI_1998" \
forrest@0: ,"UMBS" \
forrest@0: ,"Vaira" \
forrest@0: ,"WalkerBranch" \
forrest@0: ,"WillowCreek" \
forrest@0: ,"WindRiver" \
forrest@0: ,"Wisconsin_ihw" \
forrest@0: ,"Wisconsin_irp" \
forrest@0: ,"Wisconsin_mrp" \
forrest@0: ,"Wisconsin_myjp" \
forrest@0: ,"Wisconsin_pb" \
forrest@0: ,"Wisconsin_rpcc" \
forrest@0: ,"Wisconsin_yhw" \
forrest@0: ,"Wisconsin_yjp" \
forrest@0: ,"Wisconsin_yrp" \
forrest@0: /)
forrest@0:
forrest@0: field = (/"NEE Flux" \
forrest@0: ,"Shortwave Incoming" \
forrest@0: ,"Latent Heat" \
forrest@0: ,"Sensible Heat" \
forrest@0: ,"GPP Flux" \
forrest@0: ,"Respiration" \
forrest@0: /)
forrest@0:
forrest@0: field_unit = (/"gC/m2/day" \
forrest@0: ,"W/m2" \
forrest@0: ,"W/m2" \
forrest@0: ,"W/m2" \
forrest@0: ,"gC/m2/day" \
forrest@0: ,"gC/m2/day" \
forrest@0: /)
forrest@0:
forrest@0: nstation = dimsizes(station)
forrest@0: nfield = dimsizes(field)
forrest@0:
forrest@0: ;========================================================================
forrest@0: ; get observed info: number of year, first/last year, lat, lon
forrest@0: ; and annual data
forrest@0:
forrest@0: dir_root = diro + "ameriflux/"
forrest@0:
forrest@0: year_station = new ((/nstation/),integer) ; number of year
forrest@0: year_ob = new ((/nstation/),string) ; observed year
forrest@0: year_ob_i = new ((/nstation/),integer) ; first year
forrest@0: year_ob_f = new ((/nstation/),integer) ; last year
forrest@0: lat_ob = new ((/nstation/),float) ; latitude
forrest@0: lon_ob = new ((/nstation/),float) ; longitude
forrest@0:
forrest@0: data_ob_ann = new ((/nfield, nmonth, nstation/),float)
forrest@0:
forrest@0: do n = 0, nstation-1
forrest@0:
forrest@0: dir_f = dir_root + station(n)+"/"
forrest@0: fil_f = "timeseries_L4_m.nc"
forrest@0: fo = addfile (dir_f+fil_f,"r")
forrest@0:
forrest@0: lat_ob(n) = fo->lat
forrest@0: lon_ob(n) = fo->lon
forrest@0:
forrest@0: year = fo->year
forrest@0:
forrest@0: year_station(n) = dimsizes(year)
forrest@0: year_ob_i(n) = year(0)
forrest@0: year_ob_f(n) = year(year_station(n)-1)
forrest@0: year_ob(n) = year_ob_i(n) + "-" + year_ob_f(n)
forrest@0:
forrest@0: delete (year)
forrest@0:
forrest@0: data = fo->NEE_or_fMDS
forrest@0: data_ob_ann(0,:,n) = dim_avg(data(month|:,year|:))
forrest@0:
forrest@0: data = fo->Rg_f
forrest@0: data_ob_ann(1,:,n) = dim_avg(data(month|:,year|:)) * factor_rad
forrest@0:
forrest@0: data = fo->LE_f
forrest@0: data_ob_ann(2,:,n) = dim_avg(data(month|:,year|:))
forrest@0:
forrest@0: data = fo->H_f
forrest@0: data_ob_ann(3,:,n) = dim_avg(data(month|:,year|:))
forrest@0:
forrest@0: data = fo->GPP_or_MDS
forrest@0: data_ob_ann(4,:,n) = dim_avg(data(month|:,year|:))
forrest@0:
forrest@0: data = fo->Reco_or
forrest@0: data_ob_ann(5,:,n) = dim_avg(data(month|:,year|:))
forrest@0:
forrest@0: delete (data)
forrest@0: delete (fo)
forrest@0: end do
forrest@0:
forrest@0: ;--------------------------------------------------------------
forrest@0: ; find (# of year observed) >=4 and year_ob_i <= 2001
forrest@0:
forrest@0: i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
forrest@0:
forrest@0: station_long = station(i_long_ob)
forrest@0: lat_ob_long = lat_ob(i_long_ob)
forrest@0: lon_ob_long = lat_ob(i_long_ob)
forrest@0: year_ob_long = year_ob(i_long_ob)
forrest@0: year_ob_i_long = year_ob_i(i_long_ob)
forrest@0: year_ob_f_long = year_ob_f(i_long_ob)
forrest@0: year_station_long = year_station(i_long_ob)
forrest@0:
forrest@0: nstation_long = dimsizes(station_long)
forrest@0:
forrest@0: ;=========================================================
forrest@0: ; get model data at observed lat-lon
forrest@0:
forrest@0: fm = addfile (dirm+film8,"r")
forrest@0:
forrest@0: xm = fm->lon
forrest@0: ym = fm->lat
forrest@0: date = fm->date
forrest@0:
forrest@0: date_dim = dimsizes(date)
forrest@0: nyear = date_dim(0)
forrest@0:
forrest@0: data_mod = new ((/nfield,nyear,nmonth,nstation/),float)
forrest@0: data_mod_ann = new ((/nfield,nmonth,nstation/),float)
forrest@0: data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)
forrest@0:
forrest@0: ; change to unit of observed (u mol/m2/s)
forrest@0: ; Model_units [=] gC/m2/s
forrest@0: ; 12. = molecular weight of C
forrest@0: ; u mol = 1e-6 mol
forrest@0:
forrest@0: factor = 1.e6 /12.
forrest@0:
forrest@0: ;------------------------------------------------------------
forrest@0: ; interpolate model data into observed station
forrest@0: ; note: model is 0-360E, 90S-90N
forrest@0:
forrest@0: ; to be able to handle observation at (-89.98,-24.80)
forrest@0: ym(0) = -90.
forrest@0: ;------------------------------------------------------------
forrest@0:
forrest@0: if (ENERGY .eq. "old") then
forrest@0:
forrest@0: data = fm->NEE
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0: data_mod(0,:,:,:) = yy(:,:,:) * factor_flux
forrest@0:
forrest@0: ;;data = fm->NETRAD
forrest@0: ; data = fm->FSA
forrest@0: ; data1 = fm->FIRA
forrest@0: ; data = data - data1
forrest@0: ; delete (data1)
forrest@0:
forrest@0: data = fm->FSDS
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0: data_mod(1,:,:,:) = yy(:,:,:)
forrest@0:
forrest@0:
forrest@0: ; data = fm->LATENT
forrest@0: data = fm->FCEV
forrest@0: data1 = fm->FCTR
forrest@0: data2 = fm->FGEV
forrest@0: data = data + data1 + data2
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0: data_mod(2,:,:,:) = yy(:,:,:)
forrest@0: delete (data1)
forrest@0: delete (data2)
forrest@0:
forrest@0: ; data = fm->SENSIBLE
forrest@0: data = fm->FSH
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0: data_mod(3,:,:,:) = yy(:,:,:)
forrest@0:
forrest@0: else
forrest@0:
forrest@0: data = fm->NEE
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0: data_mod(0,:,:,:) = yy(:,:,:) * factor_flux
forrest@0:
forrest@0: ; data = fm->NETRAD
forrest@0: data = fm->FSDS
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0: data_mod(1,:,:,:) = yy(:,:,:)
forrest@0:
forrest@0: data = fm->LATENT
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0: data_mod(2,:,:,:) = yy(:,:,:)
forrest@0:
forrest@0: ; data = fm->SENSIBLE
forrest@0: data = fm->FSH
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0: data_mod(3,:,:,:) = yy(:,:,:)
forrest@0:
forrest@0: end if
forrest@0:
forrest@0: data = fm->GPP
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0: data_mod(4,:,:,:) = yy(:,:,:) * factor_flux
forrest@0:
forrest@0: data = fm->ER
forrest@0: yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0: data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0: data_mod(5,:,:,:) = yy(:,:,:) * factor_flux
forrest@0:
forrest@0: data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)
forrest@0:
forrest@0: delete (data_mod)
forrest@0: delete (fm)
forrest@0: delete (data)
forrest@0: delete (yy)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; for station line plot
forrest@0: ;*******************************************************************
forrest@0:
forrest@0: ; for x-axis in xyplot
forrest@0: mon = ispan(1,12,1)
forrest@0: mon@long_name = "month"
forrest@0:
forrest@0: res = True ; plot mods desired
forrest@0: res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
forrest@0: res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
forrest@0:
forrest@0: res@tmXBFormat = "f" ; not to add trailing zeros
forrest@0:
forrest@0: ;-------------------------------------------------------------------------
forrest@0: ; Add a boxed legend using the more simple method
forrest@0:
forrest@0: res@pmLegendDisplayMode = "Always"
forrest@0: ; res@pmLegendWidthF = 0.1
forrest@0: res@pmLegendWidthF = 0.08
forrest@0: res@pmLegendHeightF = 0.06
forrest@0: ; res@pmLegendOrthogonalPosF = -1.17
forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0: res@pmLegendOrthogonalPosF = -0.30 ;(downward)
forrest@0:
forrest@0: ; res@pmLegendParallelPosF = 0.18
forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward)
forrest@0:
forrest@0: ; res@lgPerimOn = False
forrest@0: res@lgLabelFontHeightF = 0.015
forrest@0: res@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0: ;-------------------------------------------------------------------
forrest@0: ; for panel plot
forrest@0: res@gsnFrame = False ; Do not draw plot
forrest@0: res@gsnDraw = False ; Do not advance frame
forrest@0:
forrest@0: pres = True ; panel plot mods desired
forrest@0: pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0: ; indiv. plots in panel
forrest@0: pres@gsnMaximize = True ; fill the page
forrest@0: ;-------------------------------------------------------------------
forrest@0:
forrest@0: ; change longitude from 0-360 to -180-180
forrest@0: lon_ob = where(lon_ob .gt. 180.,lon_ob-360., lon_ob)
forrest@0: lon_ob_long = where(lon_ob_long .gt. 180.,lon_ob_long-360., lon_ob_long)
forrest@0:
forrest@0: ;==============================================================
forrest@0: ; get ob data at each site with long observation
forrest@0:
forrest@0: do n = 0,nstation_long-1
forrest@0:
forrest@0: ;##################################################################
forrest@0: ; hardwired: model up to year 2004
forrest@0: ; observed up to year 2006
forrest@0:
forrest@0: year_setback = 0
forrest@0:
forrest@0: nyear = year_station_long(n)
forrest@0:
forrest@0: if (year_ob_f_long(n).eq. 2006) then
forrest@0: year_setback = 2006 -2004
forrest@0: end if
forrest@0: if (year_ob_f_long(n).eq. 2005) then
forrest@0: year_setback = 2005 -2004
forrest@0: end if
forrest@0: ;##################################################################
forrest@0:
forrest@0: ntime = (nyear - year_setback) * nmonth
forrest@0:
forrest@0: data_ob = new ((/nfield, nyear, nmonth/),float)
forrest@0:
forrest@0: dir_f = dir_root + station_long(n)+"/"
forrest@0: fil_f = "timeseries_L4_m.nc"
forrest@0: fo = addfile (dir_f+fil_f,"r")
forrest@0:
forrest@0: data_ob(0,:,:) = fo->NEE_or_fMDS
forrest@0: data_ob(1,:,:) = fo->Rg_f
forrest@0: data_ob(2,:,:) = fo->LE_f
forrest@0: data_ob(3,:,:) = fo->H_f
forrest@0: data_ob(4,:,:) = fo->GPP_or_MDS
forrest@0: data_ob(5,:,:) = fo->Reco_or
forrest@0:
forrest@0: data_ob(1,:,:) = data_ob(1,:,:) * factor_rad
forrest@0:
forrest@0: delete (fo)
forrest@0:
forrest@0: timeI = new((/ntime/),integer)
forrest@0: timeF = new((/ntime/),float)
forrest@0: timeI = ispan(1,ntime,1)
forrest@0: timeF = year_ob_i_long(n) + (timeI-1)/12.
forrest@0: timeF@long_name = "year"
forrest@0:
forrest@0: plot_data = new((/2,ntime/),float)
forrest@0:
forrest@0: ;----------------------------
forrest@0: ; for model_vs_ob
forrest@0:
forrest@0: plot_name = station_long(n)+"_tseries_vs_ob"
forrest@0: title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"
forrest@0: res@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name)
forrest@0: plot=new(nfield,graphic) ; create graphic array
forrest@0:
forrest@0: i_year_mod_i = year_ob_i_long(n) - 1990
forrest@0: i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
forrest@0:
forrest@0: i_year_ob_f = nyear - year_setback - 1
forrest@0:
forrest@0: ; print (nyear)
forrest@0: ; print (i_year_ob_f)
forrest@0: ; print (i_year_mod_i)
forrest@0: ; print (i_year_mod_f)
forrest@0:
forrest@0: do i = 0,nfield-1
forrest@0: plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
forrest@0: plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))
forrest@0: plot_data@long_name = field(i)+" ("+field_unit(i)+")"
forrest@0: plot(i)=gsn_csm_xy(wks,timeF,plot_data,res) ; create plot
forrest@0: end do
forrest@0:
forrest@0: gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
forrest@0:
forrest@0: delete (wks)
forrest@0: delete (plot)
forrest@0:
forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@0:
forrest@0: delete (data_ob)
forrest@0: delete (timeI)
forrest@0: delete (timeF)
forrest@0: delete (plot_data)
forrest@0:
forrest@0: end do
forrest@0:
forrest@0: ;###################################################################
forrest@0: ; for the following tables,
forrest@0: ; sort by latitude in decending order (N->S)
forrest@0:
forrest@0: isort = dim_pqsort(lat_ob_long,-1)
forrest@0:
forrest@0: station_sort = station_long(isort)
forrest@0: year_ob_sort = year_ob_long(isort)
forrest@0: lat_ob_sort = lat_ob_long(isort)
forrest@0: lon_ob_sort = lon_ob_long(isort)
forrest@0:
forrest@0: ; print(isort)
forrest@0: ; print(lat_ob_sort)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; html table of site: observed
forrest@0: ;*******************************************************************
forrest@0: output_html = "tseries_vs_ob.html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"
" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"Timeseries at Site: "+model_name+" vs Observation
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," Site Name | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@0: ," Observed | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ; row of table
forrest@0:
forrest@0: do n = 0,nstation_long-1
forrest@0:
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt0 = station_sort(n)
forrest@0: txt1 = sprintf("%5.2f", lat_ob_sort(n))
forrest@0: txt2 = sprintf("%5.2f", lon_ob_sort(n))
forrest@0: txt3 = year_ob_sort(n)
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0: delete (idx)
forrest@0:
forrest@0: delete (isort)
forrest@0: delete (station_sort)
forrest@0: delete (year_ob_sort)
forrest@0: delete (lat_ob_sort)
forrest@0: delete (lon_ob_sort)
forrest@0:
forrest@0: ;************************************************************
forrest@0: ; compute annual cycle correlation coef and M score
forrest@0: ;************************************************************
forrest@0:
forrest@0: score_max = 1.
forrest@0:
forrest@0: ccr = new ((/nstation, nfield/),float)
forrest@0: M_score = new ((/nstation, nfield/),float)
forrest@0:
forrest@0: do n=0,nstation-1
forrest@0: do m=0,nfield-1
forrest@0: ccr(n,m) = esccr(data_ob_ann(m,:,n),data_mod_ann(m,:,n),0)
forrest@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))))
forrest@0: M_score(n,m) = (1. -(bias/nmonth)) * score_max
forrest@0: end do
forrest@0: end do
forrest@0:
forrest@0: M_nee = avg(M_score(:,0))
forrest@0: M_rad = avg(M_score(:,1))
forrest@0: M_lh = avg(M_score(:,2))
forrest@0: M_sh = avg(M_score(:,3))
forrest@0: M_gpp = avg(M_score(:,4))
forrest@0: M_er = avg(M_score(:,5))
forrest@0: M_all = M_nee+ M_rad +M_lh + M_sh + M_gpp + M_er
forrest@0:
forrest@0: M_ameriflux_nee = sprintf("%.2f", M_nee)
forrest@0: M_ameriflux_rad = sprintf("%.2f", M_rad)
forrest@0: M_ameriflux_lh = sprintf("%.2f", M_lh )
forrest@0: M_ameriflux_sh = sprintf("%.2f", M_sh )
forrest@0: M_ameriflux_gpp = sprintf("%.2f", M_gpp)
forrest@0: M_ameriflux_er = sprintf("%.2f", M_er )
forrest@0: M_ameriflux_all = sprintf("%.2f", M_all)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; for station line plot
forrest@0: ;*******************************************************************
forrest@0:
forrest@0: ; for x-axis in xyplot
forrest@0: mon = ispan(1,12,1)
forrest@0: mon@long_name = "month"
forrest@0:
forrest@0: ;-------------------------------------------------------------------
forrest@0:
forrest@0: plot_data = new((/2,nmonth/),float)
forrest@0: plot_data!0 = "case"
forrest@0: plot_data!1 = "month"
forrest@0:
forrest@0: do n = 0,nstation-1
forrest@0: ;----------------------------
forrest@0: ; for observed
forrest@0:
forrest@0: plot_name = station(n)+"_ob"
forrest@0: title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
forrest@0: res@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name)
forrest@0: plot=new(nfield,graphic) ; create graphic array
forrest@0:
forrest@0: do i = 0,nfield-1
forrest@0: plot_data(0,:) = (/data_ob_ann(i,:,n)/)
forrest@0: plot_data@long_name = field(i)+" ("+field_unit(i)+")"
forrest@0: plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot
forrest@0: end do
forrest@0:
forrest@0: gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
forrest@0:
forrest@0: delete (wks)
forrest@0: delete (plot)
forrest@0:
forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@0:
forrest@0: ;----------------------------
forrest@0: ; for model_vs_ob
forrest@0:
forrest@0: plot_name = station(n)+"_model_vs_ob"
forrest@0: title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
forrest@0: res@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name)
forrest@0: plot=new(nfield,graphic) ; create graphic array
forrest@0:
forrest@0: do i = 0,nfield-1
forrest@0: plot_data(0,:) = (/data_ob_ann(i,:,n)/)
forrest@0: plot_data(1,:) = (/data_mod_ann(i,:,n)/)
forrest@0: plot_data@long_name = field(i)+" ("+field_unit(i)+")"
forrest@0: plot(i)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot
forrest@0: end do
forrest@0:
forrest@0: gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
forrest@0:
forrest@0: delete (wks)
forrest@0: delete (plot)
forrest@0:
forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@0: end do
forrest@0:
forrest@0: ;###################################################################
forrest@0: ; for the following tables,
forrest@0: ; sort by latitude in decending order (N->S)
forrest@0: ; sort by lat in decending order (N->S)
forrest@0:
forrest@0: isort = dim_pqsort(lat_ob,-1)
forrest@0:
forrest@0: station_sort = station(isort)
forrest@0: year_ob_sort = year_ob(isort)
forrest@0: lat_ob_sort = lat_ob(isort)
forrest@0: lon_ob_sort = lon_ob(isort)
forrest@0: M_score_sort = M_score(isort,:)
forrest@0:
forrest@0: ; print(isort)
forrest@0: ; print(lat_ob_sort)
forrest@0: ;###################################################################
forrest@0: ;*******************************************************************
forrest@0: ; html table of site: observed
forrest@0: ;*******************************************************************
forrest@0: output_html = "line_ob.html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"Energy at Site: Observation
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," Site Name | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@0: ," Observed | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ; row of table
forrest@0:
forrest@0: do n = 0,nstation-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt0 = station_sort(n)
forrest@0: txt1 = sprintf("%5.2f", lat_ob_sort(n))
forrest@0: txt2 = sprintf("%5.2f", lon_ob_sort(n))
forrest@0: txt3 = year_ob_sort(n)
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0: delete (idx)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; score and line table : model vs observed
forrest@0: ;*******************************************************************
forrest@0: output_html = "score+line_vs_ob.html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"Energy at Site: Model "+model_name+"
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: delete (table_header)
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," Site Name | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@0: ," Observed | " \
forrest@0: ," NEE Flux | " \
forrest@0: ," Shortwave Incoming | " \
forrest@0: ," Latent Heat | " \
forrest@0: ," Sensible Heat | " \
forrest@0: ," GPP Flux | " \
forrest@0: ," Respiration | " \
forrest@0: ," Average | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ; row of table
forrest@0:
forrest@0: do n = 0,nstation-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt0 = station_sort(n)
forrest@0: txt1 = sprintf("%5.2f", lat_ob_sort(n))
forrest@0: txt2 = sprintf("%5.2f", lon_ob_sort(n))
forrest@0: txt3 = year_ob_sort(n)
forrest@0: txt4 = sprintf("%5.2f", M_score_sort(n,0))
forrest@0: txt5 = sprintf("%5.2f", M_score_sort(n,1))
forrest@0: txt6 = sprintf("%5.2f", M_score_sort(n,2))
forrest@0: txt7 = sprintf("%5.2f", M_score_sort(n,3))
forrest@0: txt8 = sprintf("%5.2f", M_score_sort(n,4))
forrest@0: txt9 = sprintf("%5.2f", M_score_sort(n,5))
forrest@0: txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0: set_line(lines,nline,""+txt7+" | ")
forrest@0: set_line(lines,nline,""+txt8+" | ")
forrest@0: set_line(lines,nline,""+txt9+" | ")
forrest@0: set_line(lines,nline,""+txt10+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0:
forrest@0: ; last row, summary
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt0 = "All_"+sprintf("%.0f", nstation)
forrest@0: txt1 = "-"
forrest@0: txt2 = "-"
forrest@0: txt3 = "-"
forrest@0: txt4 = M_ameriflux_nee
forrest@0: txt5 = M_ameriflux_rad
forrest@0: txt6 = M_ameriflux_lh
forrest@0: txt7 = M_ameriflux_sh
forrest@0: txt8 = M_ameriflux_gpp
forrest@0: txt9 = M_ameriflux_er
forrest@0: txt10 = M_ameriflux_all
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0: set_line(lines,nline,""+txt7+" | ")
forrest@0: set_line(lines,nline,""+txt8+" | ")
forrest@0: set_line(lines,nline,""+txt9+" | ")
forrest@0: set_line(lines,nline,""+txt10+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0: delete (idx)
forrest@0:
forrest@0: ;**************************************************************************************
forrest@0: ; update score
forrest@0: ;**************************************************************************************
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_ameriflux_nee/s/M_ameriflux_nee/"+M_ameriflux_nee+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: "sed -e '1,/M_ameriflux_rad/s/M_ameriflux_rad/"+M_ameriflux_rad+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: "sed -e '1,/M_ameriflux_lh/s/M_ameriflux_lh/"+M_ameriflux_lh+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: "sed -e '1,/M_ameriflux_sh/s/M_ameriflux_sh/"+M_ameriflux_sh+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: "sed -e '1,/M_ameriflux_gpp/s/M_ameriflux_gpp/"+M_ameriflux_gpp+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: "sed -e '1,/M_ameriflux_er/s/M_ameriflux_er/"+M_ameriflux_er+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_ameriflux_nee#"+M_ameriflux_nee+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name+";"+ \
forrest@0: "sed s#M_ameriflux_rad#"+M_ameriflux_rad+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name+";"+ \
forrest@0: "sed s#M_ameriflux_lh#"+M_ameriflux_lh+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name+";"+ \
forrest@0: "sed s#M_ameriflux_sh#"+M_ameriflux_sh+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name+";"+ \
forrest@0: "sed s#M_ameriflux_gpp#"+M_ameriflux_gpp+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name+";"+ \
forrest@0: "sed s#M_ameriflux_er#"+M_ameriflux_er+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ; add total score and write to file
forrest@0: ;***************************************************************************
forrest@0: M_total = M_ameriflux_all
forrest@0:
forrest@0: asciiwrite("M_save.ameriflux", M_total)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ; output plot and html
forrest@0: ;***************************************************************************
forrest@0: output_dir = model_name+"/ameriflux"
forrest@0:
forrest@0: system("mv *.png *.html " + output_dir)
forrest@0: ;***************************************************************************
forrest@0:
forrest@0: end
forrest@0: