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