ameriflux/13.plot.ncl
changeset 0 0c6405ab2ff4
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/ameriflux/13.plot.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,683 @@
     1.4 +;************************************************************
     1.5 +; change long_name
     1.6 +;************************************************************
     1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
     1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
     1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    1.10 +;************************************************************
    1.11 +procedure set_line(lines:string,nline:integer,newlines:string) 
    1.12 +begin
    1.13 +; add line to ascci/html file
    1.14 +    
    1.15 +  nnewlines = dimsizes(newlines)
    1.16 +  if(nline+nnewlines-1.ge.dimsizes(lines))
    1.17 +    print("set_line: bad index, not setting anything.") 
    1.18 +    return
    1.19 +  end if 
    1.20 +  lines(nline:nline+nnewlines-1) = newlines
    1.21 +;  print ("lines = " + lines(nline:nline+nnewlines-1))
    1.22 +  nline = nline + nnewlines
    1.23 +  return 
    1.24 +end
    1.25 +;*************************************************************
    1.26 +begin
    1.27 +
    1.28 +  plot_type     = "ps"
    1.29 +  plot_type_new = "png"
    1.30 +
    1.31 +; for 6 fields, 12-monthly
    1.32 +  nmon      = 12
    1.33 +  nfield    = 6
    1.34 +
    1.35 +;************************************************
    1.36 +; read model data
    1.37 +;************************************************
    1.38 +  model = "cn"
    1.39 +; model = "casa"
    1.40 +  
    1.41 +  model_name = "i01.10" + model
    1.42 +
    1.43 +  dirm  = "/fis/cgd/cseg/people/jeff/surface_data/"
    1.44 +  film  = "lnd_T42.nc"
    1.45 +  fm    = addfile(dirm+film,"r")
    1.46 +  
    1.47 +  xm    = fm->lon
    1.48 +  ym    = fm->lat
    1.49 +;------------------------------------------------
    1.50 +  nlat = dimsizes(ym)
    1.51 +  nlon = dimsizes(xm)
    1.52 +
    1.53 +  data_mod0 = new ((/nfield,nmon,nlat,nlon/),float)
    1.54 +
    1.55 +; change to unit of observed (u mol/m2/s)
    1.56 +; Model_units [=] gC/m2/s
    1.57 +; 12. = molecular weight of C
    1.58 +; u mol = 1e-6 mol
    1.59 +  factor = 1e6 /12.
    1.60 +
    1.61 +  ENERGY ="new"
    1.62 +
    1.63 +;************************************************
    1.64 +; read data: observed
    1.65 +;************************************************
    1.66 +
    1.67 + station = (/"ARM_Oklahoma" \
    1.68 +            ,"ARM_Oklahoma_burn" \
    1.69 +            ,"ARM_Oklahoma_control" \
    1.70 +            ,"Atqasuk" \
    1.71 +            ,"Audubon" \
    1.72 +            ,"AustinCary" \
    1.73 +            ,"Bartlett" \
    1.74 +            ,"Bondville" \
    1.75 +            ,"Brookings" \
    1.76 +            ,"Donaldson" \
    1.77 +            ,"Duke_Forest_Hardwoods" \
    1.78 +            ,"Duke_Forest_Open_Field" \
    1.79 +            ,"Duke_Forest_Pine" \
    1.80 +            ,"Fermi_Ag" \
    1.81 +            ,"Fermi_Prairie" \
    1.82 +            ,"Flagstaff_Managed" \
    1.83 +            ,"Flagstaff_Unmanaged" \
    1.84 +            ,"Flagstaff_Wildfire" \
    1.85 +            ,"FortPeck" \
    1.86 +            ,"FreemanRanch_mesquite" \
    1.87 +            ,"Goodwin_Creek" \
    1.88 +            ,"HarvardForest" \
    1.89 +            ,"HarvardForestHemlock" \
    1.90 +            ,"HowlandForestMain" \
    1.91 +            ,"HowlandForestWest" \
    1.92 +            ,"Ivotuk" \
    1.93 +            ,"KendallGrasslands" \
    1.94 +            ,"KennedySpaceCenterPine" \
    1.95 +            ,"KennedySpaceCenterScrub" \
    1.96 +            ,"LittleProspect" \
    1.97 +            ,"LostCreek" \
    1.98 +            ,"Mead-irrigated" \
    1.99 +            ,"Mead-irrigated-rotation" \
   1.100 +            ,"Mead-rainfed" \
   1.101 +            ,"Metolius_2nd_YoungPonderosaPine" \
   1.102 +            ,"MetoliusEyerly" \
   1.103 +            ,"MetoliusIntermediatePine" \
   1.104 +            ,"MetoliusOldPonderosaPine" \
   1.105 +            ,"MissouriOzark" \
   1.106 +            ,"Mize" \
   1.107 +            ,"MorganMonroe" \
   1.108 +            ,"NiwotRidge" \
   1.109 +            ,"NorthCarolina_cc" \
   1.110 +            ,"NorthCarolina_lp" \
   1.111 +            ,"ParkFalls" \
   1.112 +            ,"Rayonier" \
   1.113 +            ,"SantaRita" \
   1.114 +            ,"SkyOaks_Old" \
   1.115 +            ,"SkyOaks_PostFire" \
   1.116 +            ,"SkyOaks_Young" \
   1.117 +            ,"SylvaniaWilderness" \
   1.118 +            ,"Toledo" \
   1.119 +            ,"Tonzi" \
   1.120 +            ,"UCI_1850" \
   1.121 +            ,"UCI_1930" \
   1.122 +            ,"UCI_1964" \
   1.123 +            ,"UCI_1964wet" \
   1.124 +            ,"UCI_1981" \
   1.125 +            ,"UCI_1989" \
   1.126 +            ,"UCI_1998" \
   1.127 +            ,"UMBS" \
   1.128 +            ,"Vaira" \
   1.129 +            ,"WalkerBranch" \
   1.130 +            ,"WillowCreek" \
   1.131 +            ,"WindRiver" \
   1.132 +            ,"Wisconsin_ihw" \
   1.133 +            ,"Wisconsin_irp" \
   1.134 +            ,"Wisconsin_mrp" \
   1.135 +            ,"Wisconsin_myjp" \
   1.136 +            ,"Wisconsin_pb" \
   1.137 +            ,"Wisconsin_rpcc" \
   1.138 +            ,"Wisconsin_yhw" \
   1.139 +            ,"Wisconsin_yjp" \
   1.140 +            ,"Wisconsin_yrp" \
   1.141 +            /)
   1.142 +
   1.143 + year_ob = (/"2003-2006" \
   1.144 +            ,"2005-2006" \
   1.145 +            ,"2005-2006" \
   1.146 +            ,"1999-2006" \
   1.147 +            ,"2002-2006" \
   1.148 +            ,"2000-2005" \
   1.149 +            ,"2004-2005" \
   1.150 +            ,"1996-2006" \
   1.151 +            ,"2004-2006" \
   1.152 +            ,"1999-2004" \
   1.153 +            ,"2003-2005" \
   1.154 +            ,"2001-2005" \
   1.155 +            ,"2001-2005" \
   1.156 +            ,"2005-2006" \
   1.157 +            ,"2004-2006" \
   1.158 +            ,"2005-2006" \
   1.159 +            ,"2005-2006" \
   1.160 +            ,"2005-2006" \
   1.161 +            ,"2000-2006" \
   1.162 +            ,"2004-2006" \
   1.163 +            ,"2002-2006" \
   1.164 +            ,"1991-2004" \
   1.165 +            ,"2004-2004" \
   1.166 +            ,"1996-2004" \
   1.167 +            ,"1999-2004" \
   1.168 +            ,"2003-2006" \
   1.169 +            ,"2004-2006" \
   1.170 +            ,"2002-2002" \
   1.171 +            ,"2000-2006" \
   1.172 +            ,"2002-2005" \
   1.173 +            ,"2001-2005" \
   1.174 +            ,"2001-2005" \
   1.175 +            ,"2001-2005" \
   1.176 +            ,"2001-2005" \
   1.177 +            ,"2004-2005" \
   1.178 +            ,"2004-2005" \
   1.179 +            ,"2003-2005" \
   1.180 +            ,"1996-2000" \
   1.181 +            ,"2004-2006" \
   1.182 +            ,"1998-2004" \
   1.183 +            ,"1999-2005" \
   1.184 +            ,"1999-2003" \
   1.185 +            ,"2005-2006" \
   1.186 +            ,"2005-2006" \
   1.187 +            ,"1996-2003" \
   1.188 +            ,"1998-1998" \
   1.189 +            ,"2004-2006" \
   1.190 +            ,"1997-2006" \
   1.191 +            ,"2004-2006" \
   1.192 +            ,"1997-2006" \
   1.193 +            ,"2002-2006" \
   1.194 +            ,"2004-2005" \
   1.195 +            ,"2001-2006" \
   1.196 +            ,"2004-2005" \
   1.197 +            ,"2001-2005" \
   1.198 +            ,"2001-2005" \
   1.199 +            ,"2002-2004" \
   1.200 +            ,"2001-2005" \
   1.201 +            ,"2001-2005" \
   1.202 +            ,"2002-2005" \
   1.203 +            ,"1999-2003" \
   1.204 +            ,"2001-2006" \
   1.205 +            ,"1995-1999" \
   1.206 +            ,"1999-2005" \
   1.207 +            ,"1999-2004" \
   1.208 +            ,"2003-2003" \
   1.209 +            ,"2003-2003" \
   1.210 +            ,"2002-2005" \
   1.211 +            ,"2004-2004" \
   1.212 +            ,"2002-2002" \
   1.213 +            ,"2005-2005" \
   1.214 +            ,"2002-2002" \
   1.215 +            ,"2004-2005" \
   1.216 +            ,"2002-2002" \
   1.217 +            /)
   1.218 +
   1.219 + field   = (/"NEE Flux" \
   1.220 +            ,"Net Radiation" \
   1.221 +            ,"Latent Heat" \
   1.222 +            ,"Sensible Heat" \
   1.223 +            ,"GPP Flux" \
   1.224 +            ,"Respiration" \
   1.225 +            /)
   1.226 +
   1.227 + field_unit = (/"u mol/m2/s" \
   1.228 +               ,"W/m2" \
   1.229 +               ,"W/m2" \
   1.230 +               ,"W/m2" \
   1.231 +               ,"u mol/m2/s" \
   1.232 +               ,"u mol/m2/s" \
   1.233 +               /)
   1.234 +
   1.235 + nstation  = dimsizes(station)
   1.236 +
   1.237 + data_mod  = new ((/nstation, nfield, nmon/),float)
   1.238 + data_ob   = new ((/nstation, nfield, nmon/),float)
   1.239 + lat_ob    = new ((/nstation/),float)
   1.240 + lon_ob    = new ((/nstation/),float)
   1.241 + unit_ob   = new ((/nfield/),string)
   1.242 +
   1.243 + diro_root  = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
   1.244 + dirm_root  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
   1.245 +
   1.246 + do n = 0,nstation-1
   1.247 +
   1.248 +;-------------------------------------------------
   1.249 +;   get ob data
   1.250 +
   1.251 +    diro = diro_root + station(n)+"/"
   1.252 +    filo = year_ob(n)+"_L4_m.nc"
   1.253 +    fo   = addfile (diro+filo,"r")
   1.254 +
   1.255 +    print (filo)
   1.256 + 
   1.257 +    lon_ob(n) = fo->lon 
   1.258 +    lat_ob(n) = fo->lat
   1.259 +
   1.260 +    data      = fo->NEE_or_fMDS
   1.261 +    data_ob(n,0,:) = dim_avg(data(month|:,year|:))
   1.262 +    delete (data)
   1.263 +
   1.264 +    data      = fo->Rg_f
   1.265 +    data_ob(n,1,:) = dim_avg(data(month|:,year|:))
   1.266 +    delete (data)
   1.267 +
   1.268 +    data      = fo->LE_f
   1.269 +    data_ob(n,2,:) = dim_avg(data(month|:,year|:))
   1.270 +    delete (data)
   1.271 +
   1.272 +    data      = fo->H_f
   1.273 +    data_ob(n,3,:) = dim_avg(data(month|:,year|:))
   1.274 +    delete (data)
   1.275 +
   1.276 +    data      = fo->GPP_or_MDS
   1.277 +    data_ob(n,4,:) = dim_avg(data(month|:,year|:))
   1.278 +    delete (data)
   1.279 +
   1.280 +    data      = fo->Reco_or
   1.281 +    data_ob(n,5,:) = dim_avg(data(month|:,year|:))
   1.282 +    delete (data)
   1.283 +
   1.284 +    delete (fo)
   1.285 +;---------------------------------------------------
   1.286 +;   get model data
   1.287 +
   1.288 +;   film = model_name+"_"+ year_ob(n)+"_MONS_climo.nc"
   1.289 +    film = model_name+"_"+"1990-2004"+"_MONS_climo.nc"
   1.290 +    fm   = addfile (dirm_root+film,"r")
   1.291 +
   1.292 +    print (film)
   1.293 +
   1.294 +if (ENERGY .eq. "old") then
   1.295 +
   1.296 +  data = fm->NEE
   1.297 +  data_mod0(0,:,:,:) = data(:,:,:) * factor
   1.298 +  delete (data)
   1.299 +
   1.300 +; data  = fm->LATENT
   1.301 +  data1 = fm->FCEV
   1.302 +  data2 = fm->FCTR
   1.303 +  data3 = fm->FGEV
   1.304 +  data_mod0(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:) 
   1.305 +  delete (data1)
   1.306 +  delete (data2)
   1.307 +  delete (data3)
   1.308 +
   1.309 +; data = fm->SENSIBLE
   1.310 +  data  = fm->FSH
   1.311 +  data_mod0(3,:,:,:) = data(:,:,:) 
   1.312 +  delete (data)
   1.313 +
   1.314 +; data  = fm->NETRAD
   1.315 +  data1 = fm->FSA
   1.316 +  data2 = fm->FIRA
   1.317 +  data_mod0(1,:,:,:) = data1(:,:,:)-data2(:,:,:)-data_mod0(2,:,:,:)-data_mod0(3,:,:,:) 
   1.318 +  delete (data1)
   1.319 +  delete (data2)
   1.320 +
   1.321 +else
   1.322 +
   1.323 +  data = fm->NEE
   1.324 +  data_mod0(0,:,:,:) = data(:,:,:) * factor
   1.325 +  delete (data)
   1.326 +
   1.327 +  data = fm->NETRAD
   1.328 +  data_mod0(1,:,:,:) = data(:,:,:) 
   1.329 +  delete (data)
   1.330 +
   1.331 +  data = fm->LATENT
   1.332 +  data_mod0(2,:,:,:) = data(:,:,:) 
   1.333 +  delete (data)
   1.334 +
   1.335 +; data = fm->SENSIBLE
   1.336 +  data = fm->FSH
   1.337 +  data_mod0(3,:,:,:) = data(:,:,:) 
   1.338 +  delete (data)
   1.339 +end if
   1.340 +
   1.341 +  data = fm->GPP
   1.342 +  data_mod0(4,:,:,:) = data(:,:,:) * factor 
   1.343 +  delete (data)
   1.344 +
   1.345 +  if (model .eq. "cn") then
   1.346 +     data = fm->ER
   1.347 +  else
   1.348 +     data1 = fm->AR
   1.349 +     data2 = fm->HR
   1.350 +     data  = data1 + data2
   1.351 +    
   1.352 +     delete (data1)
   1.353 +     delete (data2)
   1.354 +  end if
   1.355 +
   1.356 +  data_mod0(5,:,:,:) = data(:,:,:) * factor
   1.357 +  delete (data)
   1.358 +
   1.359 +  delete (fm) 
   1.360 +
   1.361 +;************************************************************
   1.362 +; interpolate model data into observed station
   1.363 +; note: model is 0-360E, 90S-90N
   1.364 +;************************************************************
   1.365 +
   1.366 +; to be able to handle observation at (-89.98,-24.80)
   1.367 +  ym(0) = -90.  
   1.368 +
   1.369 +  yy = linint2_points_Wrap(xm,ym,data_mod0,True,lon_ob(n),lat_ob(n),0)
   1.370 +
   1.371 +; print (yy)
   1.372 + 
   1.373 +  data_mod(n,:,:) = yy(:,:,0)
   1.374 +
   1.375 +  delete (yy)
   1.376 +
   1.377 + end do
   1.378 +;************************************************************
   1.379 +; compute correlation coef and M score
   1.380 +;************************************************************
   1.381 +
   1.382 + score_max = 5.
   1.383 +
   1.384 + ccr     = new ((/nstation, nfield/),float)
   1.385 + M_score = new ((/nstation, nfield/),float) 
   1.386 +
   1.387 + do n=0,nstation-1
   1.388 + do m=0,nfield-1   
   1.389 +    ccr(n,m) = esccr(data_ob(n,m,:),data_mod(n,m,:),0)
   1.390 +    bias = sum(abs(data_mod(n,m,:)-data_ob(n,m,:))/(abs(data_mod(n,m,:))+abs(data_ob(n,m,:))))
   1.391 +    M_score(n,m) = (1. -(bias/nmon)) * score_max
   1.392 + end do
   1.393 + end do
   1.394 +
   1.395 + M_co2 = avg(M_score(:,0))
   1.396 + M_rad = avg(M_score(:,1))
   1.397 + M_lh  = avg(M_score(:,2))
   1.398 + M_sh  = avg(M_score(:,3))
   1.399 + M_gpp = avg(M_score(:,4))
   1.400 + M_er  = avg(M_score(:,5))
   1.401 + M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er
   1.402 +
   1.403 + M_energy_co2 = sprintf("%.2f", M_co2)
   1.404 + M_energy_rad = sprintf("%.2f", M_rad)
   1.405 + M_energy_lh  = sprintf("%.2f", M_lh )
   1.406 + M_energy_sh  = sprintf("%.2f", M_sh )
   1.407 + M_energy_gpp = sprintf("%.2f", M_gpp)
   1.408 + M_energy_er  = sprintf("%.2f", M_er )
   1.409 + M_energy_all = sprintf("%.2f", M_all)
   1.410 +
   1.411 +;*******************************************************************
   1.412 +; for station line plot
   1.413 +;*******************************************************************
   1.414 +
   1.415 +; for x-axis in xyplot
   1.416 +  mon = ispan(1,12,1)
   1.417 +  mon@long_name = "month"
   1.418 +
   1.419 +  res                   = True               ; plot mods desired
   1.420 +  res@xyLineThicknesses = (/2.0,2.0/)        ; make 2nd lines thicker
   1.421 +  res@xyLineColors      = (/"blue","red"/)   ; line color (ob,model)
   1.422 +;-------------------------------------------------------------------------
   1.423 +; Add a boxed legend using the more simple method
   1.424 +
   1.425 +  res@pmLegendDisplayMode    = "Always"
   1.426 +; res@pmLegendWidthF         = 0.1
   1.427 +  res@pmLegendWidthF         = 0.08
   1.428 +  res@pmLegendHeightF        = 0.06
   1.429 +; res@pmLegendOrthogonalPosF = -1.17
   1.430 +; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
   1.431 +  res@pmLegendOrthogonalPosF = -0.30  ;(downward)
   1.432 +
   1.433 +; res@pmLegendParallelPosF   =  0.18
   1.434 +  res@pmLegendParallelPosF   =  0.23  ;(rightward)
   1.435 +
   1.436 +; res@lgPerimOn             = False
   1.437 +  res@lgLabelFontHeightF     = 0.015
   1.438 +  res@xyExplicitLegendLabels = (/"observed",model_name/)
   1.439 +;-------------------------------------------------------------------
   1.440 +; for panel plot
   1.441 +  res@gsnFrame     = False                   ; Do not draw plot 
   1.442 +  res@gsnDraw      = False                   ; Do not advance frame
   1.443 +
   1.444 +  pres                            = True     ; panel plot mods desired
   1.445 +  pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
   1.446 +                                             ; indiv. plots in panel
   1.447 +  pres@gsnMaximize                = True     ; fill the page
   1.448 +;-------------------------------------------------------------------
   1.449 +
   1.450 +  plot_data   = new((/2,12/),float)
   1.451 +  plot_data!0 = "case"
   1.452 +  plot_data!1 = "month"
   1.453 +
   1.454 +  do n = 0,nstation-1
   1.455 +;----------------------------
   1.456 +; for observed
   1.457 +
   1.458 +    plot_name = station(n)+"_ob"    
   1.459 +    title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
   1.460 +    res@tiMainString = title
   1.461 +
   1.462 +    wks = gsn_open_wks (plot_type,plot_name)
   1.463 +    plot=new(nfield,graphic)                        ; create graphic array   
   1.464 +
   1.465 +    do i = 0,nfield-1                           
   1.466 +       plot_data(0,:) = (/data_ob (n,i,:)/)
   1.467 +       plot_data@long_name = field(i)+" ("+field_unit(i)+")"   
   1.468 +       plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res)    ; create plot 
   1.469 +    end do
   1.470 +
   1.471 +    gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot
   1.472 +
   1.473 +    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   1.474 +           "rm "+plot_name+"."+plot_type)
   1.475 +
   1.476 +    clear (wks)  
   1.477 +    delete (plot)
   1.478 +;----------------------------
   1.479 +; for model_vs_ob
   1.480 +
   1.481 +    plot_name = station(n)+"_model_vs_ob"
   1.482 +    title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"    
   1.483 +    res@tiMainString = title
   1.484 +
   1.485 +    wks = gsn_open_wks (plot_type,plot_name)
   1.486 +    plot=new(nfield,graphic)                         ; create graphic array   
   1.487 +
   1.488 +    do i = 0,nfield-1                           
   1.489 +       plot_data(0,:) = (/data_ob (n,i,:)/)
   1.490 +       plot_data(1,:) = (/data_mod(n,i,:)/)
   1.491 +       plot_data@long_name = field(i)+" ("+field_unit(i)+")"   
   1.492 +       plot(i)=gsn_csm_xy(wks,mon,plot_data,res)     ; create plot 
   1.493 +    end do
   1.494 +   
   1.495 +    gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot
   1.496 +
   1.497 +    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   1.498 +           "rm "+plot_name+"."+plot_type)
   1.499 +
   1.500 +    clear (wks)  
   1.501 +    delete (plot)
   1.502 + end do
   1.503 +
   1.504 +;*******************************************************************
   1.505 +; html table of site: observed
   1.506 +;*******************************************************************
   1.507 +  output_html = "line_ob.html"
   1.508 +
   1.509 +  header = (/"<HTML>" \
   1.510 +            ,"<HEAD>" \
   1.511 +            ,"<TITLE>CLAMP metrics</TITLE>" \
   1.512 +            ,"</HEAD>" \
   1.513 +            ,"<H1>Energy at Site: Observation</H1>" \
   1.514 +            /) 
   1.515 +  footer = "</HTML>"
   1.516 +
   1.517 +  table_header = (/ \
   1.518 +        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
   1.519 +       ,"<tr>" \
   1.520 +       ,"   <th bgcolor=DDDDDD >Site Name</th>" \
   1.521 +       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
   1.522 +       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
   1.523 +       ,"   <th bgcolor=DDDDDD >Observed</th>" \ 
   1.524 +       ,"</tr>" \
   1.525 +       /)
   1.526 +  table_footer = "</table>"
   1.527 +  row_header = "<tr>"
   1.528 +  row_footer = "</tr>"
   1.529 +
   1.530 +  lines = new(50000,string)
   1.531 +  nline = 0
   1.532 +
   1.533 +  set_line(lines,nline,header)
   1.534 +  set_line(lines,nline,table_header)
   1.535 +;-----------------------------------------------
   1.536 +; row of table
   1.537 +  
   1.538 +  do n = 0,nstation-1
   1.539 +     set_line(lines,nline,row_header)
   1.540 +
   1.541 +     txt0 = station(n)
   1.542 +     txt1 = sprintf("%5.2f", lat_ob(n))
   1.543 +     txt2 = sprintf("%5.2f", lon_ob(n))
   1.544 +     txt3 = year_ob(n)
   1.545 +
   1.546 +     set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")
   1.547 +     set_line(lines,nline,"<th>"+txt1+"</th>")
   1.548 +     set_line(lines,nline,"<th>"+txt2+"</th>")
   1.549 +     set_line(lines,nline,"<th>"+txt3+"</th>")
   1.550 +
   1.551 +     set_line(lines,nline,row_footer)
   1.552 +  end do
   1.553 +;-----------------------------------------------
   1.554 +  set_line(lines,nline,table_footer)
   1.555 +  set_line(lines,nline,footer) 
   1.556 +
   1.557 +; Now write to an HTML file.
   1.558 +  idx = ind(.not.ismissing(lines))
   1.559 +  if(.not.any(ismissing(idx))) then
   1.560 +    asciiwrite(output_html,lines(idx))
   1.561 +  else
   1.562 +   print ("error?")
   1.563 +  end if
   1.564 +  delete (idx)
   1.565 +
   1.566 +;*******************************************************************
   1.567 +; score and line table : model vs observed
   1.568 +;*******************************************************************
   1.569 +  output_html = "score+line_vs_ob.html"
   1.570 +
   1.571 +  header = (/"<HTML>" \
   1.572 +            ,"<HEAD>" \
   1.573 +            ,"<TITLE>CLAMP metrics</TITLE>" \
   1.574 +            ,"</HEAD>" \
   1.575 +            ,"<H1>Energy at Site: Model "+model_name+"</H1>" \
   1.576 +            /) 
   1.577 +  footer = "</HTML>"
   1.578 +
   1.579 +  delete (table_header)
   1.580 +  table_header = (/ \
   1.581 +        "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
   1.582 +       ,"<tr>" \
   1.583 +       ,"   <th bgcolor=DDDDDD >Site Name</th>" \
   1.584 +       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
   1.585 +       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
   1.586 +       ,"   <th bgcolor=DDDDDD >Observed</th>" \
   1.587 +       ,"   <th bgcolor=DDDDDD >CO2 Flux</th>" \
   1.588 +       ,"   <th bgcolor=DDDDDD >Net Radiation</th>" \
   1.589 +       ,"   <th bgcolor=DDDDDD >Latent Heat</th>" \
   1.590 +       ,"   <th bgcolor=DDDDDD >Sensible Heat</th>" \
   1.591 +       ,"   <th bgcolor=DDDDDD >GPP Glux</th>" \
   1.592 +       ,"   <th bgcolor=DDDDDD >Respiration</th>" \
   1.593 +       ,"   <th bgcolor=DDDDDD >Average</th>" \
   1.594 +       ,"</tr>" \
   1.595 +       /)
   1.596 +  table_footer = "</table>"
   1.597 +  row_header = "<tr>"
   1.598 +  row_footer = "</tr>"
   1.599 +
   1.600 +  lines = new(50000,string)
   1.601 +  nline = 0
   1.602 +
   1.603 +  set_line(lines,nline,header)
   1.604 +  set_line(lines,nline,table_header)
   1.605 +;-----------------------------------------------
   1.606 +; row of table
   1.607 +  
   1.608 +  do n = 0,nstation-1
   1.609 +     set_line(lines,nline,row_header)
   1.610 +
   1.611 +     txt0  = station(n)
   1.612 +     txt1  = sprintf("%5.2f", lat_ob(n))
   1.613 +     txt2  = sprintf("%5.2f", lon_ob(n))
   1.614 +     txt3  = year_ob(n)
   1.615 +     txt4  = sprintf("%5.2f", M_score(n,0))
   1.616 +     txt5  = sprintf("%5.2f", M_score(n,1))
   1.617 +     txt6  = sprintf("%5.2f", M_score(n,2))
   1.618 +     txt7  = sprintf("%5.2f", M_score(n,3))
   1.619 +     txt8  = sprintf("%5.2f", M_score(n,4))
   1.620 +     txt9  = sprintf("%5.2f", M_score(n,5))
   1.621 +     txt10 = sprintf("%5.2f", avg(M_score(n,:)))
   1.622 +
   1.623 +     set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")
   1.624 +     set_line(lines,nline,"<th>"+txt1+"</th>")
   1.625 +     set_line(lines,nline,"<th>"+txt2+"</th>")
   1.626 +     set_line(lines,nline,"<th>"+txt3+"</th>")
   1.627 +     set_line(lines,nline,"<th>"+txt4+"</th>")
   1.628 +     set_line(lines,nline,"<th>"+txt5+"</th>")
   1.629 +     set_line(lines,nline,"<th>"+txt6+"</th>")
   1.630 +     set_line(lines,nline,"<th>"+txt7+"</th>")
   1.631 +     set_line(lines,nline,"<th>"+txt8+"</th>")
   1.632 +     set_line(lines,nline,"<th>"+txt9+"</th>")
   1.633 +     set_line(lines,nline,"<th>"+txt10+"</th>")
   1.634 +
   1.635 +     set_line(lines,nline,row_footer)
   1.636 +  end do
   1.637 +
   1.638 +; last row, summary
   1.639 +  set_line(lines,nline,row_header)
   1.640 +
   1.641 +  txt0  = "All_"+sprintf("%.0f", nstation)
   1.642 +  txt1  = "-"
   1.643 +  txt2  = "-"
   1.644 +  txt3  = "-"
   1.645 +  txt4  = M_energy_co2
   1.646 +  txt5  = M_energy_rad
   1.647 +  txt6  = M_energy_lh
   1.648 +  txt7  = M_energy_sh
   1.649 +  txt8  = M_energy_gpp
   1.650 +  txt9  = M_energy_er
   1.651 +  txt10 = M_energy_all
   1.652 +
   1.653 +  set_line(lines,nline,"<th>"+txt0+"</th>")
   1.654 +  set_line(lines,nline,"<th>"+txt1+"</th>")
   1.655 +  set_line(lines,nline,"<th>"+txt2+"</th>")
   1.656 +  set_line(lines,nline,"<th>"+txt3+"</th>")
   1.657 +  set_line(lines,nline,"<th>"+txt4+"</th>")
   1.658 +  set_line(lines,nline,"<th>"+txt5+"</th>")
   1.659 +  set_line(lines,nline,"<th>"+txt6+"</th>")
   1.660 +  set_line(lines,nline,"<th>"+txt7+"</th>")
   1.661 +  set_line(lines,nline,"<th>"+txt8+"</th>")
   1.662 +  set_line(lines,nline,"<th>"+txt9+"</th>")
   1.663 +  set_line(lines,nline,"<th>"+txt10+"</th>")
   1.664 +
   1.665 +  set_line(lines,nline,row_footer)
   1.666 +;-----------------------------------------------
   1.667 +  set_line(lines,nline,table_footer)
   1.668 +  set_line(lines,nline,footer) 
   1.669 +
   1.670 +; Now write to an HTML file.
   1.671 +  idx = ind(.not.ismissing(lines))
   1.672 +  if(.not.any(ismissing(idx))) then
   1.673 +    asciiwrite(output_html,lines(idx))
   1.674 +  else
   1.675 +   print ("error?")
   1.676 +  end if
   1.677 +  delete (idx)
   1.678 +
   1.679 +;***************************************************************************
   1.680 +; output plots
   1.681 +;***************************************************************************
   1.682 +  output_dir = model_name+"/ameriflux"
   1.683 +
   1.684 +  system("mv *.png *.html " + output_dir) 
   1.685 +;***************************************************************************
   1.686 +end