clamp: all/11.ameriflux.ncl@0c6405ab2ff4

     1 ;************************************************************

     2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

     3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     5 ;************************************************************

     6 procedure set_line(lines:string,nline:integer,newlines:string)

     7 begin

     8 ; add line to ascci/html file

    10   nnewlines = dimsizes(newlines)

    11   if(nline+nnewlines-1.ge.dimsizes(lines))

    12     print("set_line: bad index, not setting anything.")

    13     return

    14   end if

    15   lines(nline:nline+nnewlines-1) = newlines

    16 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    17   nline = nline + nnewlines

    18   return

    19 end

    20 ;*************************************************************

    21 begin

    23   plot_type     = "ps"

    24   plot_type_new = "png"

    26 ;------------------------------------------------------

    27 ; edit table.html of current model for movel1_vs_model2

    29  if (isvar("compare")) then

    30     html_name2 = compare+"/table.html"

    31     html_new2  = html_name2 +".new"

    32  end if

    34 ;------------------------------------------------------

    35 ; edit table.html for current model

    37  html_name = model_name+"/table.html"

    38  html_new  = html_name +".new"

    40 ;------------------------------------------------------

    42   nmonth = 12

    44 ; for nee, gpp, and ar

    45 ; observed unit is gC/m2/day

    46 ; model    unit is gC/m2/s

    47 ; to change to observed unit,

    49   factor_flux = 86400.

    51 ; for incident solar radiation,

    52 ; observed Rg_f unit is MJ/m2/day

    53 ; model (FSDS)  unit is  W/m2

    54 ; to change to model unit,

    56   factor_rad = 1.e6/86400.

    58 ;************************************************

    59 ; observed data info

    60 ;************************************************

    62  station = (/"ARM_Oklahoma" \

    63             ,"ARM_Oklahoma_burn" \

    64             ,"ARM_Oklahoma_control" \

    65             ,"Atqasuk" \

    66             ,"Audubon" \

    67             ,"AustinCary" \

    68             ,"Bartlett" \

    69             ,"Bondville" \

    70             ,"Brookings" \

    71             ,"Donaldson" \

    72             ,"Duke_Forest_Hardwoods" \

    73             ,"Duke_Forest_Open_Field" \

    74             ,"Duke_Forest_Pine" \

    75             ,"Fermi_Ag" \

    76             ,"Fermi_Prairie" \

    77             ,"Flagstaff_Managed" \

    78             ,"Flagstaff_Unmanaged" \

    79             ,"Flagstaff_Wildfire" \

    80             ,"FortPeck" \

    81             ,"FreemanRanch_mesquite" \

    82             ,"Goodwin_Creek" \

    83             ,"HarvardForest" \

    84             ,"HarvardForestHemlock" \

    85             ,"HowlandForestMain" \

    86             ,"HowlandForestWest" \

    87             ,"Ivotuk" \

    88             ,"KendallGrasslands" \

    89             ,"KennedySpaceCenterPine" \

    90             ,"KennedySpaceCenterScrub" \

    91             ,"LittleProspect" \

    92             ,"LostCreek" \

    93             ,"Mead-irrigated" \

    94             ,"Mead-irrigated-rotation" \

    95             ,"Mead-rainfed" \

    96             ,"Metolius_2nd_YoungPonderosaPine" \

    97             ,"MetoliusEyerly" \

    98             ,"MetoliusIntermediatePine" \

    99             ,"MetoliusOldPonderosaPine" \

   100             ,"MissouriOzark" \

   101             ,"Mize" \

   102             ,"MorganMonroe" \

   103             ,"NiwotRidge" \

   104             ,"NorthCarolina_cc" \

   105             ,"NorthCarolina_lp" \

   106             ,"ParkFalls" \

   107             ,"Rayonier" \

   108             ,"SantaRita" \

   109             ,"SkyOaks_Old" \

   110             ,"SkyOaks_PostFire" \

   111             ,"SkyOaks_Young" \

   112             ,"SylvaniaWilderness" \

   113             ,"Toledo" \

   114             ,"Tonzi" \

   115             ,"UCI_1850" \

   116             ,"UCI_1930" \

   117             ,"UCI_1964" \

   118             ,"UCI_1964wet" \

   119             ,"UCI_1981" \

   120             ,"UCI_1989" \

   121             ,"UCI_1998" \

   122             ,"UMBS" \

   123             ,"Vaira" \

   124             ,"WalkerBranch" \

   125             ,"WillowCreek" \

   126             ,"WindRiver" \

   127             ,"Wisconsin_ihw" \

   128             ,"Wisconsin_irp" \

   129             ,"Wisconsin_mrp" \

   130             ,"Wisconsin_myjp" \

   131             ,"Wisconsin_pb" \

   132             ,"Wisconsin_rpcc" \

   133             ,"Wisconsin_yhw" \

   134             ,"Wisconsin_yjp" \

   135             ,"Wisconsin_yrp" \

   136             /)

   138  field   = (/"NEE Flux" \

   139             ,"Shortwave Incoming" \

   140             ,"Latent Heat" \

   141             ,"Sensible Heat" \

   142             ,"GPP Flux" \

   143             ,"Respiration" \

   144             /)

   146  field_unit = (/"gC/m2/day" \

   147                ,"W/m2" \

   148                ,"W/m2" \

   149                ,"W/m2" \

   150                ,"gC/m2/day" \

   151                ,"gC/m2/day" \

   152                /)

   154  nstation = dimsizes(station)

   155  nfield   = dimsizes(field)

   157 ;========================================================================

   158 ; get observed info: number of year, first/last year, lat, lon

   159 ; and annual data

   161  dir_root = diro + "ameriflux/"

   163  year_station = new ((/nstation/),integer) ; number of year

   164  year_ob      = new ((/nstation/),string)  ; observed year

   165  year_ob_i    = new ((/nstation/),integer) ; first year

   166  year_ob_f    = new ((/nstation/),integer) ; last year

   167  lat_ob       = new ((/nstation/),float)   ; latitude

   168  lon_ob       = new ((/nstation/),float)   ; longitude

   170  data_ob_ann  = new ((/nfield, nmonth, nstation/),float)

   172  do n = 0, nstation-1

   174     dir_f = dir_root + station(n)+"/"

   175     fil_f = "timeseries_L4_m.nc"

   176     fo   = addfile (dir_f+fil_f,"r")

   178     lat_ob(n) = fo->lat

   179     lon_ob(n) = fo->lon

   181     year      = fo->year

   183     year_station(n) = dimsizes(year)

   184     year_ob_i(n)    = year(0)

   185     year_ob_f(n)    = year(year_station(n)-1)

   186     year_ob(n)      = year_ob_i(n) + "-" + year_ob_f(n)

   188     delete (year)

   190     data = fo->NEE_or_fMDS

   191     data_ob_ann(0,:,n) = dim_avg(data(month|:,year|:))

   193     data = fo->Rg_f

   194     data_ob_ann(1,:,n) = dim_avg(data(month|:,year|:)) * factor_rad

   196     data = fo->LE_f

   197     data_ob_ann(2,:,n) = dim_avg(data(month|:,year|:))

   199     data = fo->H_f

   200     data_ob_ann(3,:,n) = dim_avg(data(month|:,year|:))

   202     data = fo->GPP_or_MDS

   203     data_ob_ann(4,:,n) = dim_avg(data(month|:,year|:))

   205     data = fo->Reco_or

   206     data_ob_ann(5,:,n) = dim_avg(data(month|:,year|:))

   208     delete (data)

   209     delete (fo)

   210  end do

   212 ;--------------------------------------------------------------

   213 ; find (# of year observed) >=4 and year_ob_i <= 2001

   215  i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)

   217  station_long      = station(i_long_ob)

   218  lat_ob_long       = lat_ob(i_long_ob)

   219  lon_ob_long       = lat_ob(i_long_ob)

   220  year_ob_long      = year_ob(i_long_ob)

   221  year_ob_i_long    = year_ob_i(i_long_ob)

   222  year_ob_f_long    = year_ob_f(i_long_ob)

   223  year_station_long = year_station(i_long_ob)

   225  nstation_long     = dimsizes(station_long)

   227 ;=========================================================

   228 ;   get model data at observed lat-lon

   230     fm   = addfile (dirm+film8,"r")

   232     xm   = fm->lon

   233     ym   = fm->lat

   234     date = fm->date

   236     date_dim = dimsizes(date)

   237     nyear  = date_dim(0)

   239     data_mod      = new ((/nfield,nyear,nmonth,nstation/),float)

   240     data_mod_ann  = new ((/nfield,nmonth,nstation/),float)

   241     data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)

   243 ;   change to unit of observed (u mol/m2/s)

   244 ;   Model_units [=] gC/m2/s

   245 ;   12. = molecular weight of C

   246 ;   u mol = 1e-6 mol

   248     factor = 1.e6 /12.

   250 ;------------------------------------------------------------

   251 ;   interpolate model data into observed station

   252 ;   note: model is 0-360E, 90S-90N

   254 ;   to be able to handle observation at (-89.98,-24.80)

   255     ym(0) = -90.

   256 ;------------------------------------------------------------

   258 if (ENERGY .eq. "old") then

   260   data = fm->NEE

   261   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   262   data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   263   data_mod(0,:,:,:) = yy(:,:,:) * factor_flux

   265 ;;data  = fm->NETRAD

   266 ; data  = fm->FSA

   267 ; data1 = fm->FIRA

   268 ; data  = data - data1

   269 ; delete (data1)

   271   data  = fm->FSDS

   272   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   273   data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   274   data_mod(1,:,:,:) = yy(:,:,:)

   277 ; data  = fm->LATENT

   278   data  = fm->FCEV

   279   data1 = fm->FCTR

   280   data2 = fm->FGEV

   281   data  = data + data1 + data2

   282   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   283   data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   284   data_mod(2,:,:,:) = yy(:,:,:)

   285   delete (data1)

   286   delete (data2)

   288 ; data = fm->SENSIBLE

   289   data  = fm->FSH

   290   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   291   data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   292   data_mod(3,:,:,:) = yy(:,:,:)

   294 else

   296   data = fm->NEE

   297   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   298   data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   299   data_mod(0,:,:,:) = yy(:,:,:) * factor_flux

   301 ; data = fm->NETRAD

   302   data = fm->FSDS

   303   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   304   data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   305   data_mod(1,:,:,:) = yy(:,:,:)

   307   data = fm->LATENT

   308   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   309   data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   310   data_mod(2,:,:,:) = yy(:,:,:)

   312 ; data = fm->SENSIBLE

   313   data = fm->FSH

   314   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   315   data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   316   data_mod(3,:,:,:) = yy(:,:,:)

   318 end if

   320   data = fm->GPP

   321   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   322   data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   323   data_mod(4,:,:,:) = yy(:,:,:) * factor_flux

   325   data = fm->ER

   326   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   327   data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   328   data_mod(5,:,:,:) = yy(:,:,:) * factor_flux

   330   data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)

   332   delete (data_mod)

   333   delete (fm)

   334   delete (data)

   335   delete (yy)

   337 ;*******************************************************************

   338 ; for station line plot

   339 ;*******************************************************************

   341 ; for x-axis in xyplot

   342   mon = ispan(1,12,1)

   343   mon@long_name = "month"

   345   res                   = True               ; plot mods desired

   346   res@xyLineThicknesses = (/2.0,2.0/)        ; make 2nd lines thicker

   347   res@xyLineColors      = (/"blue","red"/)   ; line color (ob,model)

   349   res@tmXBFormat  = "f"                      ; not to add trailing zeros

   351 ;-------------------------------------------------------------------------

   352 ; Add a boxed legend using the more simple method

   354   res@pmLegendDisplayMode    = "Always"

   355 ; res@pmLegendWidthF         = 0.1

   356   res@pmLegendWidthF         = 0.08

   357   res@pmLegendHeightF        = 0.06

   358 ; res@pmLegendOrthogonalPosF = -1.17

   359 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   360   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   362 ; res@pmLegendParallelPosF   =  0.18

   363   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   365 ; res@lgPerimOn             = False

   366   res@lgLabelFontHeightF     = 0.015

   367   res@xyExplicitLegendLabels = (/"observed",model_name/)

   368 ;-------------------------------------------------------------------

   369 ; for panel plot

   370   res@gsnFrame     = False                   ; Do not draw plot

   371   res@gsnDraw      = False                   ; Do not advance frame

   373   pres                            = True     ; panel plot mods desired

   374   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   375                                              ; indiv. plots in panel

   376   pres@gsnMaximize                = True     ; fill the page

   377 ;-------------------------------------------------------------------

   379 ; change longitude from 0-360 to -180-180

   380   lon_ob      = where(lon_ob      .gt. 180.,lon_ob-360.,      lon_ob)

   381   lon_ob_long = where(lon_ob_long .gt. 180.,lon_ob_long-360., lon_ob_long)

   383 ;==============================================================

   384 ; get ob data at each site with long observation

   386   do n = 0,nstation_long-1

   388 ;##################################################################

   389 ; hardwired: model    up to year 2004

   390 ;            observed up to year 2006

   392     year_setback = 0

   394     nyear = year_station_long(n)

   396     if (year_ob_f_long(n).eq. 2006) then

   397        year_setback = 2006 -2004

   398     end if

   399     if (year_ob_f_long(n).eq. 2005) then

   400        year_setback = 2005 -2004

   401     end if

   402 ;##################################################################

   404     ntime = (nyear - year_setback) * nmonth

   406     data_ob   = new ((/nfield, nyear, nmonth/),float)

   408     dir_f = dir_root + station_long(n)+"/"

   409     fil_f = "timeseries_L4_m.nc"

   410     fo    = addfile (dir_f+fil_f,"r")

   412     data_ob(0,:,:) = fo->NEE_or_fMDS

   413     data_ob(1,:,:) = fo->Rg_f

   414     data_ob(2,:,:) = fo->LE_f

   415     data_ob(3,:,:) = fo->H_f

   416     data_ob(4,:,:) = fo->GPP_or_MDS

   417     data_ob(5,:,:) = fo->Reco_or

   419     data_ob(1,:,:) = data_ob(1,:,:) * factor_rad

   421     delete (fo)

   423     timeI = new((/ntime/),integer)

   424     timeF = new((/ntime/),float)

   425     timeI = ispan(1,ntime,1)

   426     timeF = year_ob_i_long(n) + (timeI-1)/12.

   427     timeF@long_name = "year"

   429     plot_data = new((/2,ntime/),float)

   431 ;----------------------------

   432 ; for model_vs_ob

   434     plot_name = station_long(n)+"_tseries_vs_ob"

   435     title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"

   436     res@tiMainString = title

   438     wks = gsn_open_wks (plot_type,plot_name)

   439     plot=new(nfield,graphic)                         ; create graphic array

   441     i_year_mod_i = year_ob_i_long(n) - 1990

   442     i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback

   444     i_year_ob_f =  nyear - year_setback - 1

   446 ;   print (nyear)

   447 ;   print (i_year_ob_f)

   448 ;   print (i_year_mod_i)

   449 ;   print (i_year_mod_f)

   451     do i = 0,nfield-1

   452        plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))

   453        plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))

   454        plot_data@long_name = field(i)+" ("+field_unit(i)+")"

   455        plot(i)=gsn_csm_xy(wks,timeF,plot_data,res)            ; create plot

   456     end do

   458     gsn_panel(wks,plot,(/3,2/),pres)                    ; create panel plot

   460     delete (wks)

   461     delete (plot)

   463     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   464            "rm "+plot_name+"."+plot_type)

   466     delete (data_ob)

   467     delete (timeI)

   468     delete (timeF)

   469     delete (plot_data)

   471  end do

   473 ;###################################################################

   474 ; for the following tables,

   475 ; sort by latitude in decending order (N->S)

   477   isort = dim_pqsort(lat_ob_long,-1)

   479   station_sort = station_long(isort)

   480   year_ob_sort = year_ob_long(isort)

   481   lat_ob_sort  = lat_ob_long(isort)

   482   lon_ob_sort  = lon_ob_long(isort)

   484 ; print(isort)

   485 ; print(lat_ob_sort)

   487 ;*******************************************************************

   488 ; html table of site: observed

   489 ;*******************************************************************

   490   output_html = "tseries_vs_ob.html"

   492   header = (/"<HTML>" \

   493             ,"<HEAD>" \

   494             ,"<TITLE>CLAMP metrics</TITLE>" \

   495             ,"</HEAD>" \

   496             ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \

   497             /)

   498   footer = "</HTML>"

   500   table_header = (/ \

   501         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \

   502        ,"<tr>" \

   503        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   504        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   505        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   506        ,"   <th bgcolor=DDDDDD >Observed</th>" \

   507        ,"</tr>" \

   508        /)

   509   table_footer = "</table>"

   510   row_header = "<tr>"

   511   row_footer = "</tr>"

   513   lines = new(50000,string)

   514   nline = 0

   516   set_line(lines,nline,header)

   517   set_line(lines,nline,table_header)

   518 ;-----------------------------------------------

   519 ; row of table

   521   do n = 0,nstation_long-1

   523      set_line(lines,nline,row_header)

   525      txt0 = station_sort(n)

   526      txt1 = sprintf("%5.2f", lat_ob_sort(n))

   527      txt2 = sprintf("%5.2f", lon_ob_sort(n))

   528      txt3 = year_ob_sort(n)

   530      set_line(lines,nline,"<th><a href="+txt0+"_tseries_vs_ob.png>"+txt0+"</a></th>")

   531      set_line(lines,nline,"<th>"+txt1+"</th>")

   532      set_line(lines,nline,"<th>"+txt2+"</th>")

   533      set_line(lines,nline,"<th>"+txt3+"</th>")

   535      set_line(lines,nline,row_footer)

   536   end do

   537 ;-----------------------------------------------

   538   set_line(lines,nline,table_footer)

   539   set_line(lines,nline,footer)

   541 ; Now write to an HTML file.

   542   idx = ind(.not.ismissing(lines))

   543   if(.not.any(ismissing(idx))) then

   544     asciiwrite(output_html,lines(idx))

   545   else

   546    print ("error?")

   547   end if

   548   delete (idx)

   550   delete (isort)

   551   delete (station_sort)

   552   delete (year_ob_sort)

   553   delete (lat_ob_sort)

   554   delete (lon_ob_sort)

   556 ;************************************************************

   557 ; compute annual cycle correlation coef and M score

   558 ;************************************************************

   560  score_max = 1.

   562  ccr     = new ((/nstation, nfield/),float)

   563  M_score = new ((/nstation, nfield/),float)

   565  do n=0,nstation-1

   566  do m=0,nfield-1

   567     ccr(n,m) = esccr(data_ob_ann(m,:,n),data_mod_ann(m,:,n),0)

   568     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))))

   569     M_score(n,m) = (1. -(bias/nmonth)) * score_max

   570  end do

   571  end do

   573  M_nee = avg(M_score(:,0))

   574  M_rad = avg(M_score(:,1))

   575  M_lh  = avg(M_score(:,2))

   576  M_sh  = avg(M_score(:,3))

   577  M_gpp = avg(M_score(:,4))

   578  M_er  = avg(M_score(:,5))

   579  M_all = M_nee+ M_rad +M_lh + M_sh + M_gpp + M_er

   581  M_ameriflux_nee = sprintf("%.2f", M_nee)

   582  M_ameriflux_rad = sprintf("%.2f", M_rad)

   583  M_ameriflux_lh  = sprintf("%.2f", M_lh )

   584  M_ameriflux_sh  = sprintf("%.2f", M_sh )

   585  M_ameriflux_gpp = sprintf("%.2f", M_gpp)

   586  M_ameriflux_er  = sprintf("%.2f", M_er )

   587  M_ameriflux_all = sprintf("%.2f", M_all)

   589 ;*******************************************************************

   590 ; for station line plot

   591 ;*******************************************************************

   593 ; for x-axis in xyplot

   594   mon = ispan(1,12,1)

   595   mon@long_name = "month"

   597 ;-------------------------------------------------------------------

   599   plot_data   = new((/2,nmonth/),float)

   600   plot_data!0 = "case"

   601   plot_data!1 = "month"

   603   do n = 0,nstation-1

   604 ;----------------------------

   605 ; for observed

   607     plot_name = station(n)+"_ob"

   608     title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"

   609     res@tiMainString = title

   611     wks = gsn_open_wks (plot_type,plot_name)

   612     plot=new(nfield,graphic)                        ; create graphic array

   614     do i = 0,nfield-1

   615        plot_data(0,:) = (/data_ob_ann(i,:,n)/)

   616        plot_data@long_name = field(i)+" ("+field_unit(i)+")"

   617        plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res)    ; create plot

   618     end do

   620     gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot

   622     delete (wks)

   623     delete (plot)

   625     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   626            "rm "+plot_name+"."+plot_type)

   628 ;----------------------------

   629 ; for model_vs_ob

   631     plot_name = station(n)+"_model_vs_ob"

   632     title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"

   633     res@tiMainString = title

   635     wks = gsn_open_wks (plot_type,plot_name)

   636     plot=new(nfield,graphic)                         ; create graphic array

   638     do i = 0,nfield-1

   639        plot_data(0,:) = (/data_ob_ann(i,:,n)/)

   640        plot_data(1,:) = (/data_mod_ann(i,:,n)/)

   641        plot_data@long_name = field(i)+" ("+field_unit(i)+")"

   642        plot(i)=gsn_csm_xy(wks,mon,plot_data,res)     ; create plot

   643     end do

   645     gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot

   647     delete (wks)

   648     delete (plot)

   650     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   651            "rm "+plot_name+"."+plot_type)

   652  end do

   654 ;###################################################################

   655 ; for the following tables,

   656 ; sort by latitude in decending order (N->S)

   657 ; sort by lat in decending order (N->S)

   659   isort = dim_pqsort(lat_ob,-1)

   661   station_sort = station(isort)

   662   year_ob_sort = year_ob(isort)

   663   lat_ob_sort  = lat_ob(isort)

   664   lon_ob_sort  = lon_ob(isort)

   665   M_score_sort = M_score(isort,:)

   667 ; print(isort)

   668 ; print(lat_ob_sort)

   669 ;###################################################################

   670 ;*******************************************************************

   671 ; html table of site: observed

   672 ;*******************************************************************

   673   output_html = "line_ob.html"

   675   header = (/"<HTML>" \

   676             ,"<HEAD>" \

   677             ,"<TITLE>CLAMP metrics</TITLE>" \

   678             ,"</HEAD>" \

   679             ,"<H1>Energy at Site: Observation</H1>" \

   680             /)

   681   footer = "</HTML>"

   683   table_header = (/ \

   684         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \

   685        ,"<tr>" \

   686        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   687        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   688        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   689        ,"   <th bgcolor=DDDDDD >Observed</th>" \

   690        ,"</tr>" \

   691        /)

   692   table_footer = "</table>"

   693   row_header = "<tr>"

   694   row_footer = "</tr>"

   696   lines = new(50000,string)

   697   nline = 0

   699   set_line(lines,nline,header)

   700   set_line(lines,nline,table_header)

   701 ;-----------------------------------------------

   702 ; row of table

   704   do n = 0,nstation-1

   705      set_line(lines,nline,row_header)

   707      txt0 = station_sort(n)

   708      txt1 = sprintf("%5.2f", lat_ob_sort(n))

   709      txt2 = sprintf("%5.2f", lon_ob_sort(n))

   710      txt3 = year_ob_sort(n)

   712      set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")

   713      set_line(lines,nline,"<th>"+txt1+"</th>")

   714      set_line(lines,nline,"<th>"+txt2+"</th>")

   715      set_line(lines,nline,"<th>"+txt3+"</th>")

   717      set_line(lines,nline,row_footer)

   718   end do

   719 ;-----------------------------------------------

   720   set_line(lines,nline,table_footer)

   721   set_line(lines,nline,footer)

   723 ; Now write to an HTML file.

   724   idx = ind(.not.ismissing(lines))

   725   if(.not.any(ismissing(idx))) then

   726     asciiwrite(output_html,lines(idx))

   727   else

   728    print ("error?")

   729   end if

   730   delete (idx)

   732 ;*******************************************************************

   733 ; score and line table : model vs observed

   734 ;*******************************************************************

   735   output_html = "score+line_vs_ob.html"

   737   header = (/"<HTML>" \

   738             ,"<HEAD>" \

   739             ,"<TITLE>CLAMP metrics</TITLE>" \

   740             ,"</HEAD>" \

   741             ,"<H1>Energy at Site: Model "+model_name+"</H1>" \

   742             /)

   743   footer = "</HTML>"

   745   delete (table_header)

   746   table_header = (/ \

   747         "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \

   748        ,"<tr>" \

   749        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   750        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   751        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   752        ,"   <th bgcolor=DDDDDD >Observed</th>" \

   753        ,"   <th bgcolor=DDDDDD >NEE Flux</th>" \

   754        ,"   <th bgcolor=DDDDDD >Shortwave <br> Incoming</th>" \

   755        ,"   <th bgcolor=DDDDDD >Latent Heat</th>" \

   756        ,"   <th bgcolor=DDDDDD >Sensible Heat</th>" \

   757        ,"   <th bgcolor=DDDDDD >GPP Flux</th>" \

   758        ,"   <th bgcolor=DDDDDD >Respiration</th>" \

   759        ,"   <th bgcolor=DDDDDD >Average</th>" \

   760        ,"</tr>" \

   761        /)

   762   table_footer = "</table>"

   763   row_header = "<tr>"

   764   row_footer = "</tr>"

   766   lines = new(50000,string)

   767   nline = 0

   769   set_line(lines,nline,header)

   770   set_line(lines,nline,table_header)

   771 ;-----------------------------------------------

   772 ; row of table

   774   do n = 0,nstation-1

   775      set_line(lines,nline,row_header)

   777      txt0  = station_sort(n)

   778      txt1  = sprintf("%5.2f", lat_ob_sort(n))

   779      txt2  = sprintf("%5.2f", lon_ob_sort(n))

   780      txt3  = year_ob_sort(n)

   781      txt4  = sprintf("%5.2f", M_score_sort(n,0))

   782      txt5  = sprintf("%5.2f", M_score_sort(n,1))

   783      txt6  = sprintf("%5.2f", M_score_sort(n,2))

   784      txt7  = sprintf("%5.2f", M_score_sort(n,3))

   785      txt8  = sprintf("%5.2f", M_score_sort(n,4))

   786      txt9  = sprintf("%5.2f", M_score_sort(n,5))

   787      txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))

   789      set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")

   790      set_line(lines,nline,"<th>"+txt1+"</th>")

   791      set_line(lines,nline,"<th>"+txt2+"</th>")

   792      set_line(lines,nline,"<th>"+txt3+"</th>")

   793      set_line(lines,nline,"<th>"+txt4+"</th>")

   794      set_line(lines,nline,"<th>"+txt5+"</th>")

   795      set_line(lines,nline,"<th>"+txt6+"</th>")

   796      set_line(lines,nline,"<th>"+txt7+"</th>")

   797      set_line(lines,nline,"<th>"+txt8+"</th>")

   798      set_line(lines,nline,"<th>"+txt9+"</th>")

   799      set_line(lines,nline,"<th>"+txt10+"</th>")

   801      set_line(lines,nline,row_footer)

   802   end do

   804 ; last row, summary

   805   set_line(lines,nline,row_header)

   807   txt0  = "All_"+sprintf("%.0f", nstation)

   808   txt1  = "-"

   809   txt2  = "-"

   810   txt3  = "-"

   811   txt4  = M_ameriflux_nee

   812   txt5  = M_ameriflux_rad

   813   txt6  = M_ameriflux_lh

   814   txt7  = M_ameriflux_sh

   815   txt8  = M_ameriflux_gpp

   816   txt9  = M_ameriflux_er

   817   txt10 = M_ameriflux_all

   819   set_line(lines,nline,"<th>"+txt0+"</th>")

   820   set_line(lines,nline,"<th>"+txt1+"</th>")

   821   set_line(lines,nline,"<th>"+txt2+"</th>")

   822   set_line(lines,nline,"<th>"+txt3+"</th>")

   823   set_line(lines,nline,"<th>"+txt4+"</th>")

   824   set_line(lines,nline,"<th>"+txt5+"</th>")

   825   set_line(lines,nline,"<th>"+txt6+"</th>")

   826   set_line(lines,nline,"<th>"+txt7+"</th>")

   827   set_line(lines,nline,"<th>"+txt8+"</th>")

   828   set_line(lines,nline,"<th>"+txt9+"</th>")

   829   set_line(lines,nline,"<th>"+txt10+"</th>")

   831   set_line(lines,nline,row_footer)

   832 ;-----------------------------------------------

   833   set_line(lines,nline,table_footer)

   834   set_line(lines,nline,footer)

   836 ; Now write to an HTML file.

   837   idx = ind(.not.ismissing(lines))

   838   if(.not.any(ismissing(idx))) then

   839     asciiwrite(output_html,lines(idx))

   840   else

   841    print ("error?")

   842   end if

   843   delete (idx)

   845 ;**************************************************************************************

   846 ; update score

   847 ;**************************************************************************************

   849   if (isvar("compare")) then

   850      system("sed -e '1,/M_ameriflux_nee/s/M_ameriflux_nee/"+M_ameriflux_nee+"/' "+html_name2+" > "+html_new2+";"+ \

   851             "mv -f "+html_new2+" "+html_name2+";"+ \

   852             "sed -e '1,/M_ameriflux_rad/s/M_ameriflux_rad/"+M_ameriflux_rad+"/' "+html_name2+" > "+html_new2+";"+ \

   853             "mv -f "+html_new2+" "+html_name2+";"+ \

   854             "sed -e '1,/M_ameriflux_lh/s/M_ameriflux_lh/"+M_ameriflux_lh+"/' "+html_name2+" > "+html_new2+";"+ \

   855             "mv -f "+html_new2+" "+html_name2+";"+ \

   856             "sed -e '1,/M_ameriflux_sh/s/M_ameriflux_sh/"+M_ameriflux_sh+"/' "+html_name2+" > "+html_new2+";"+ \

   857             "mv -f "+html_new2+" "+html_name2+";"+ \

   858             "sed -e '1,/M_ameriflux_gpp/s/M_ameriflux_gpp/"+M_ameriflux_gpp+"/' "+html_name2+" > "+html_new2+";"+ \

   859             "mv -f "+html_new2+" "+html_name2+";"+ \

   860             "sed -e '1,/M_ameriflux_er/s/M_ameriflux_er/"+M_ameriflux_er+"/' "+html_name2+" > "+html_new2+";"+ \

   861             "mv -f "+html_new2+" "+html_name2)

   862   end if

   864   system("sed s#M_ameriflux_nee#"+M_ameriflux_nee+"# "+html_name+" > "+html_new+";"+ \

   865          "mv -f "+html_new+" "+html_name+";"+ \

   866          "sed s#M_ameriflux_rad#"+M_ameriflux_rad+"# "+html_name+" > "+html_new+";"+ \

   867          "mv -f "+html_new+" "+html_name+";"+ \

   868          "sed s#M_ameriflux_lh#"+M_ameriflux_lh+"# "+html_name+" > "+html_new+";"+ \

   869          "mv -f "+html_new+" "+html_name+";"+ \

   870          "sed s#M_ameriflux_sh#"+M_ameriflux_sh+"# "+html_name+" > "+html_new+";"+ \

   871          "mv -f "+html_new+" "+html_name+";"+ \

   872          "sed s#M_ameriflux_gpp#"+M_ameriflux_gpp+"# "+html_name+" > "+html_new+";"+ \

   873          "mv -f "+html_new+" "+html_name+";"+ \

   874          "sed s#M_ameriflux_er#"+M_ameriflux_er+"# "+html_name+" > "+html_new+";"+ \

   875          "mv -f "+html_new+" "+html_name)

   877 ;***************************************************************************

   878 ; add total score and write to file

   879 ;***************************************************************************

   880   M_total = M_ameriflux_all

   882   asciiwrite("M_save.ameriflux", M_total)

   884 ;***************************************************************************

   885 ; output plot and html

   886 ;***************************************************************************

   887   output_dir = model_name+"/ameriflux"

   889   system("mv *.png *.html " + output_dir)

   890 ;***************************************************************************

   892 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
child 1	4be95183fbcd
permissions	-rw-r--r--