ameriflux/21.plot.ncl
changeset 0 0c6405ab2ff4
equal deleted inserted replaced
-1:000000000000 0:749db502a8ea
       
     1 ;************************************************************
       
     2 ; add table
       
     3 ; add time series plot
       
     4 ; sort by latitude in decending order (N->S)
       
     5 ; find year_ob >= 4
       
     6 ; add year_ob_i and year_ob_f
       
     7 ; change long_name
       
     8 ;************************************************************
       
     9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
    10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
    11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
    12 ;************************************************************
       
    13 procedure set_line(lines:string,nline:integer,newlines:string) 
       
    14 begin
       
    15 ; add line to ascci/html file
       
    16     
       
    17   nnewlines = dimsizes(newlines)
       
    18   if(nline+nnewlines-1.ge.dimsizes(lines))
       
    19     print("set_line: bad index, not setting anything.") 
       
    20     return
       
    21   end if 
       
    22   lines(nline:nline+nnewlines-1) = newlines
       
    23 ;  print ("lines = " + lines(nline:nline+nnewlines-1))
       
    24   nline = nline + nnewlines
       
    25   return 
       
    26 end
       
    27 ;*************************************************************
       
    28 begin
       
    29 
       
    30   plot_type     = "ps"
       
    31   plot_type_new = "png"
       
    32 
       
    33 ;************************************************
       
    34 ; model data info
       
    35 ;************************************************
       
    36   model = "cn"
       
    37 ; model = "casa"
       
    38   
       
    39   model_name = "i01.10" + model
       
    40   
       
    41   ENERGY ="new"
       
    42 
       
    43 ;************************************************
       
    44 ; observed data info
       
    45 ;************************************************
       
    46 
       
    47  station = (/"ARM_Oklahoma" \
       
    48             ,"ARM_Oklahoma_burn" \
       
    49             ,"ARM_Oklahoma_control" \
       
    50             ,"Atqasuk" \
       
    51             ,"Audubon" \
       
    52             ,"AustinCary" \
       
    53             ,"Bartlett" \
       
    54             ,"Bondville" \
       
    55             ,"Brookings" \
       
    56             ,"Donaldson" \
       
    57             ,"Duke_Forest_Hardwoods" \
       
    58             ,"Duke_Forest_Open_Field" \
       
    59             ,"Duke_Forest_Pine" \
       
    60             ,"Fermi_Ag" \
       
    61             ,"Fermi_Prairie" \
       
    62             ,"Flagstaff_Managed" \
       
    63             ,"Flagstaff_Unmanaged" \
       
    64             ,"Flagstaff_Wildfire" \
       
    65             ,"FortPeck" \
       
    66             ,"FreemanRanch_mesquite" \
       
    67             ,"Goodwin_Creek" \
       
    68             ,"HarvardForest" \
       
    69             ,"HarvardForestHemlock" \
       
    70             ,"HowlandForestMain" \
       
    71             ,"HowlandForestWest" \
       
    72             ,"Ivotuk" \
       
    73             ,"KendallGrasslands" \
       
    74             ,"KennedySpaceCenterPine" \
       
    75             ,"KennedySpaceCenterScrub" \
       
    76             ,"LittleProspect" \
       
    77             ,"LostCreek" \
       
    78             ,"Mead-irrigated" \
       
    79             ,"Mead-irrigated-rotation" \
       
    80             ,"Mead-rainfed" \
       
    81             ,"Metolius_2nd_YoungPonderosaPine" \
       
    82             ,"MetoliusEyerly" \
       
    83             ,"MetoliusIntermediatePine" \
       
    84             ,"MetoliusOldPonderosaPine" \
       
    85             ,"MissouriOzark" \
       
    86             ,"Mize" \
       
    87             ,"MorganMonroe" \
       
    88             ,"NiwotRidge" \
       
    89             ,"NorthCarolina_cc" \
       
    90             ,"NorthCarolina_lp" \
       
    91             ,"ParkFalls" \
       
    92             ,"Rayonier" \
       
    93             ,"SantaRita" \
       
    94             ,"SkyOaks_Old" \
       
    95             ,"SkyOaks_PostFire" \
       
    96             ,"SkyOaks_Young" \
       
    97             ,"SylvaniaWilderness" \
       
    98             ,"Toledo" \
       
    99             ,"Tonzi" \
       
   100             ,"UCI_1850" \
       
   101             ,"UCI_1930" \
       
   102             ,"UCI_1964" \
       
   103             ,"UCI_1964wet" \
       
   104             ,"UCI_1981" \
       
   105             ,"UCI_1989" \
       
   106             ,"UCI_1998" \
       
   107             ,"UMBS" \
       
   108             ,"Vaira" \
       
   109             ,"WalkerBranch" \
       
   110             ,"WillowCreek" \
       
   111             ,"WindRiver" \
       
   112             ,"Wisconsin_ihw" \
       
   113             ,"Wisconsin_irp" \
       
   114             ,"Wisconsin_mrp" \
       
   115             ,"Wisconsin_myjp" \
       
   116             ,"Wisconsin_pb" \
       
   117             ,"Wisconsin_rpcc" \
       
   118             ,"Wisconsin_yhw" \
       
   119             ,"Wisconsin_yjp" \
       
   120             ,"Wisconsin_yrp" \
       
   121             /)
       
   122 
       
   123  year_ob = (/"2003-2006" \
       
   124             ,"2005-2006" \
       
   125             ,"2005-2006" \
       
   126             ,"1999-2006" \
       
   127             ,"2002-2006" \
       
   128             ,"2000-2005" \
       
   129             ,"2004-2005" \
       
   130             ,"1996-2006" \
       
   131             ,"2004-2006" \
       
   132             ,"1999-2004" \
       
   133             ,"2003-2005" \
       
   134             ,"2001-2005" \
       
   135             ,"2001-2005" \
       
   136             ,"2005-2006" \
       
   137             ,"2004-2006" \
       
   138             ,"2005-2006" \
       
   139             ,"2005-2006" \
       
   140             ,"2005-2006" \
       
   141             ,"2000-2006" \
       
   142             ,"2004-2006" \
       
   143             ,"2002-2006" \
       
   144             ,"1991-2004" \
       
   145             ,"2004-2004" \
       
   146             ,"1996-2004" \
       
   147             ,"1999-2004" \
       
   148             ,"2003-2006" \
       
   149             ,"2004-2006" \
       
   150             ,"2002-2002" \
       
   151             ,"2000-2006" \
       
   152             ,"2002-2005" \
       
   153             ,"2001-2005" \
       
   154             ,"2001-2005" \
       
   155             ,"2001-2005" \
       
   156             ,"2001-2005" \
       
   157             ,"2004-2005" \
       
   158             ,"2004-2005" \
       
   159             ,"2003-2005" \
       
   160             ,"1996-2000" \
       
   161             ,"2004-2006" \
       
   162             ,"1998-2004" \
       
   163             ,"1999-2005" \
       
   164             ,"1999-2003" \
       
   165             ,"2005-2006" \
       
   166             ,"2005-2006" \
       
   167             ,"1996-2003" \
       
   168             ,"1998-1998" \
       
   169             ,"2004-2006" \
       
   170             ,"1997-2006" \
       
   171             ,"2004-2006" \
       
   172             ,"1997-2006" \
       
   173             ,"2002-2006" \
       
   174             ,"2004-2005" \
       
   175             ,"2001-2006" \
       
   176             ,"2004-2005" \
       
   177             ,"2001-2005" \
       
   178             ,"2001-2005" \
       
   179             ,"2002-2004" \
       
   180             ,"2001-2005" \
       
   181             ,"2001-2005" \
       
   182             ,"2002-2005" \
       
   183             ,"1999-2003" \
       
   184             ,"2001-2006" \
       
   185             ,"1995-1999" \
       
   186             ,"1999-2005" \
       
   187             ,"1999-2004" \
       
   188             ,"2003-2003" \
       
   189             ,"2003-2003" \
       
   190             ,"2002-2005" \
       
   191             ,"2004-2004" \
       
   192             ,"2002-2002" \
       
   193             ,"2005-2005" \
       
   194             ,"2002-2002" \
       
   195             ,"2004-2005" \
       
   196             ,"2002-2002" \
       
   197             /)
       
   198 
       
   199  field   = (/"NEE Flux" \
       
   200             ,"Net Radiation" \
       
   201             ,"Latent Heat" \
       
   202             ,"Sensible Heat" \
       
   203             ,"GPP Flux" \
       
   204             ,"Respiration" \
       
   205             /)
       
   206 
       
   207  field_unit = (/"u mol/m2/s" \
       
   208                ,"W/m2" \
       
   209                ,"W/m2" \
       
   210                ,"W/m2" \
       
   211                ,"u mol/m2/s" \
       
   212                ,"u mol/m2/s" \
       
   213                /)
       
   214 
       
   215  nstation = dimsizes(station)
       
   216  nfield   = dimsizes(field)
       
   217 
       
   218 ;========================================================================
       
   219 ; get observed info: number of year, first/last year, lat, lon
       
   220 
       
   221  dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
       
   222 
       
   223  year_station = new ((/nstation/),integer) ; number of year
       
   224  year_ob_i    = new ((/nstation/),integer) ; first year
       
   225  year_ob_f    = new ((/nstation/),integer) ; last year
       
   226  lat_ob       = new ((/nstation/),float)   ; latitude
       
   227  lon_ob       = new ((/nstation/),float)   ; longitude 
       
   228 
       
   229  do n = 0, nstation-1
       
   230 
       
   231     diro = dir_root + station(n)+"/"
       
   232     filo = year_ob(n)+"_L4_m.nc"
       
   233     fo   = addfile (diro+filo,"r")
       
   234  
       
   235     lat_ob(n) = fo->lat
       
   236     lon_ob(n) = fo->lon
       
   237 
       
   238     year      = fo->year
       
   239 
       
   240     year_station(n) = dimsizes(year)
       
   241     year_ob_i(n)    = year(0)
       
   242     year_ob_f(n)    = year(year_station(n)-1)
       
   243  
       
   244     delete (year) 
       
   245     
       
   246  end do
       
   247 
       
   248  delete (fo)
       
   249 ;print (year_ob_i)
       
   250 ;print (year_ob_f)
       
   251 
       
   252 ;--------------------------------------------------------------
       
   253 ; find (# of year observed) >=4 and year_ob_i <= 2001
       
   254 
       
   255  i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
       
   256  
       
   257  station_long      = station(i_long_ob)
       
   258  lat_ob_long       = lat_ob(i_long_ob)
       
   259  lon_ob_long       = lat_ob(i_long_ob)
       
   260  year_ob_long      = year_ob(i_long_ob)
       
   261  year_ob_i_long    = year_ob_i(i_long_ob)
       
   262  year_ob_f_long    = year_ob_f(i_long_ob)
       
   263  year_station_long = year_station(i_long_ob)
       
   264 
       
   265  nstation_long     = dimsizes(station_long)
       
   266 
       
   267 ; print (i_long_ob)
       
   268 ; print (nstation_long)
       
   269 ; print (station_long)
       
   270 ; print (year_ob_i(i_long_ob))
       
   271 
       
   272 ;=========================================================
       
   273 ;   get model data at observed lat-lon
       
   274 
       
   275     dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"    
       
   276     film = model_name+"_ameriflux_1990-2004_monthly.nc"
       
   277     fm   = addfile (dirm+film,"r")
       
   278 
       
   279     xm   = fm->lon
       
   280     ym   = fm->lat 
       
   281     date = fm->date
       
   282 
       
   283     date_dim = dimsizes(date)
       
   284     nyear  = date_dim(0)
       
   285 
       
   286     nmonth =12
       
   287 
       
   288     data_mod      = new ((/nfield,nyear,nmonth,nstation/),float)
       
   289     data_mod_ann  = new ((/nfield,nmonth,nstation/),float)
       
   290     data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)
       
   291 
       
   292 ;   change to unit of observed (u mol/m2/s)
       
   293 ;   Model_units [=] gC/m2/s
       
   294 ;   12. = molecular weight of C
       
   295 ;   u mol = 1e-6 mol
       
   296 
       
   297     factor = 1.e6 /12.
       
   298 
       
   299 ;------------------------------------------------------------
       
   300 ;   interpolate model data into observed station
       
   301 ;   note: model is 0-360E, 90S-90N
       
   302 
       
   303 ;   to be able to handle observation at (-89.98,-24.80)
       
   304     ym(0) = -90.
       
   305 ;------------------------------------------------------------  
       
   306 
       
   307 if (ENERGY .eq. "old") then
       
   308 
       
   309   data = fm->NEE
       
   310   data_mod(0,:,:,:) = data(:,:,:) * factor
       
   311 
       
   312 ; data  = fm->LATENT
       
   313   data1 = fm->FCEV
       
   314   data2 = fm->FCTR
       
   315   data3 = fm->FGEV
       
   316   data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:) 
       
   317   delete (data1)
       
   318   delete (data2)
       
   319   delete (data3)
       
   320 
       
   321 ; data = fm->SENSIBLE
       
   322   data  = fm->FSH
       
   323   data_mod(3,:,:,:) = data(:,:,:) 
       
   324 
       
   325 ; data  = fm->NETRAD
       
   326   data1 = fm->FSA
       
   327   data2 = fm->FIRA
       
   328   data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:) 
       
   329   delete (data1)
       
   330   delete (data2)
       
   331 
       
   332 else
       
   333 
       
   334   data = fm->NEE
       
   335   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
       
   336   data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
       
   337   data_mod(0,:,:,:) = yy(:,:,:) * factor
       
   338 
       
   339 ; printVarSummary(yy)
       
   340 
       
   341   data = fm->NETRAD
       
   342   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
       
   343   data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
       
   344   data_mod(1,:,:,:) = yy(:,:,:) 
       
   345 
       
   346   data = fm->LATENT
       
   347   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
       
   348   data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
       
   349   data_mod(2,:,:,:) = yy(:,:,:) 
       
   350 
       
   351 ; data = fm->SENSIBLE
       
   352   data = fm->FSH
       
   353   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
       
   354   data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
       
   355   data_mod(3,:,:,:) = yy(:,:,:) 
       
   356 
       
   357 end if
       
   358 
       
   359   data = fm->GPP
       
   360   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
       
   361   data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
       
   362   data_mod(4,:,:,:) = yy(:,:,:) * factor 
       
   363 
       
   364   if (model .eq. "cn") then
       
   365      data = fm->ER
       
   366   else
       
   367      data1 = fm->AR
       
   368      data2 = fm->HR
       
   369      data  = data1 + data2
       
   370     
       
   371      delete (data1)
       
   372      delete (data2)
       
   373   end if
       
   374 
       
   375   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
       
   376   data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor
       
   377   data_mod(5,:,:,:) = yy(:,:,:) * factor
       
   378 
       
   379   data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)
       
   380 
       
   381   delete (data_mod)
       
   382   delete (fm) 
       
   383   delete (data)
       
   384   delete (yy)
       
   385 
       
   386 ;*******************************************************************
       
   387 ; for station line plot
       
   388 ;*******************************************************************
       
   389 
       
   390 ; for x-axis in xyplot
       
   391   mon = ispan(1,12,1)
       
   392   mon@long_name = "month"
       
   393 
       
   394   res                   = True               ; plot mods desired
       
   395   res@xyLineThicknesses = (/2.0,2.0/)        ; make 2nd lines thicker
       
   396   res@xyLineColors      = (/"blue","red"/)   ; line color (ob,model)
       
   397 
       
   398   res@tmXBFormat  = "f"                      ; not to add trailing zeros
       
   399 
       
   400 ;-------------------------------------------------------------------------
       
   401 ; Add a boxed legend using the more simple method
       
   402 
       
   403   res@pmLegendDisplayMode    = "Always"
       
   404 ; res@pmLegendWidthF         = 0.1
       
   405   res@pmLegendWidthF         = 0.08
       
   406   res@pmLegendHeightF        = 0.06
       
   407 ; res@pmLegendOrthogonalPosF = -1.17
       
   408 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
       
   409   res@pmLegendOrthogonalPosF = -0.30  ;(downward)
       
   410 
       
   411 ; res@pmLegendParallelPosF   =  0.18
       
   412   res@pmLegendParallelPosF   =  0.23  ;(rightward)
       
   413 
       
   414 ; res@lgPerimOn             = False
       
   415   res@lgLabelFontHeightF     = 0.015
       
   416   res@xyExplicitLegendLabels = (/"observed",model_name/)
       
   417 ;-------------------------------------------------------------------
       
   418 ; for panel plot
       
   419   res@gsnFrame     = False                   ; Do not draw plot 
       
   420   res@gsnDraw      = False                   ; Do not advance frame
       
   421 
       
   422   pres                            = True     ; panel plot mods desired
       
   423   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
       
   424                                              ; indiv. plots in panel
       
   425   pres@gsnMaximize                = True     ; fill the page
       
   426 ;-------------------------------------------------------------------
       
   427 
       
   428 ;==============================================================
       
   429 ; get ob data at each site with long observation
       
   430 
       
   431   dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
       
   432 
       
   433   do n = 0,nstation_long-1
       
   434 ; do n = 0,0
       
   435 
       
   436     year_setback = 0
       
   437 
       
   438     nyear = year_station_long(n)
       
   439 
       
   440     if (year_ob_f_long(n).eq. 2006) then
       
   441        year_setback = 2006 -2004
       
   442     end if
       
   443     if (year_ob_f_long(n).eq. 2005) then
       
   444        year_setback = 2005 -2004
       
   445     end if
       
   446     
       
   447     ntime = (nyear - year_setback) * nmonth
       
   448     print (nyear)
       
   449 
       
   450     data_ob   = new ((/nfield, nyear, nmonth/),float)
       
   451 
       
   452     diro = dir_root + station_long(n)+"/"
       
   453     filo = year_ob_long(n)+"_L4_m.nc"
       
   454     fo   = addfile (diro+filo,"r")
       
   455 
       
   456 ;   print (diro)
       
   457 ;   print (filo)
       
   458 
       
   459     data_ob(0,:,:) = fo->NEE_or_fMDS
       
   460     data_ob(1,:,:) = fo->Rg_f
       
   461     data_ob(2,:,:) = fo->LE_f
       
   462     data_ob(3,:,:) = fo->H_f
       
   463     data_ob(4,:,:) = fo->GPP_or_MDS
       
   464     data_ob(5,:,:) = fo->Reco_or
       
   465 
       
   466     delete (fo)
       
   467 
       
   468     timeI = new((/ntime/),integer)
       
   469     timeF = new((/ntime/),float)
       
   470     timeI = ispan(1,ntime,1)
       
   471     timeF = year_ob_i_long(n) + (timeI-1)/12.
       
   472     timeF@long_name = "year" 
       
   473 
       
   474     plot_data = new((/2,ntime/),float)
       
   475  
       
   476 ;----------------------------
       
   477 ; for model_vs_ob
       
   478 
       
   479     plot_name = station_long(n)+"_tseries_vs_ob"
       
   480     title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"    
       
   481     res@tiMainString = title
       
   482 
       
   483     wks = gsn_open_wks (plot_type,plot_name)
       
   484     plot=new(nfield,graphic)                         ; create graphic array   
       
   485 
       
   486     i_year_mod_i = year_ob_i_long(n) - 1990 
       
   487     i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
       
   488 
       
   489     i_year_ob_f =  nyear - year_setback - 1 
       
   490 
       
   491 ;   print (nyear)
       
   492 ;   print (i_year_ob_f)
       
   493 ;   print (i_year_mod_i)
       
   494 ;   print (i_year_mod_f)
       
   495 
       
   496     do i = 0,nfield-1                           
       
   497        plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
       
   498        plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))
       
   499        plot_data@long_name = field(i)+" ("+field_unit(i)+")"   
       
   500        plot(i)=gsn_csm_xy(wks,timeF,plot_data,res)            ; create plot 
       
   501     end do
       
   502    
       
   503     gsn_panel(wks,plot,(/3,2/),pres)                    ; create panel plot
       
   504 
       
   505     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
       
   506            "rm "+plot_name+"."+plot_type)
       
   507 
       
   508     clear (wks)  
       
   509     delete (plot)
       
   510 
       
   511     delete (data_ob)
       
   512     delete (timeI)    
       
   513     delete (timeF)
       
   514     delete (plot_data)
       
   515 
       
   516  end do
       
   517 
       
   518 ;###################################################################
       
   519 ; for the following tables,
       
   520 ; sort by latitude in decending order (N->S)
       
   521 
       
   522   isort = dim_pqsort(lat_ob_long,-1)
       
   523 
       
   524   station_sort = station_long(isort)
       
   525   year_ob_sort = year_ob_long(isort)
       
   526   lat_ob_sort  = lat_ob_long(isort)
       
   527   lon_ob_sort  = lon_ob_long(isort)
       
   528  
       
   529 ; print(isort)
       
   530 ; print(lat_ob_sort)
       
   531 
       
   532 ;*******************************************************************
       
   533 ; html table of site: observed
       
   534 ;*******************************************************************
       
   535   output_html = "tseries_vs_ob.html"
       
   536 
       
   537   header = (/"<HTML>" \
       
   538             ,"<HEAD>" \
       
   539             ,"<TITLE>CLAMP metrics</TITLE>" \
       
   540             ,"</HEAD>" \
       
   541             ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \
       
   542             /) 
       
   543   footer = "</HTML>"
       
   544 
       
   545   table_header = (/ \
       
   546         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
       
   547        ,"<tr>" \
       
   548        ,"   <th bgcolor=DDDDDD >Site Name</th>" \
       
   549        ,"   <th bgcolor=DDDDDD >Latitude</th>" \
       
   550        ,"   <th bgcolor=DDDDDD >Longitude</th>" \
       
   551        ,"   <th bgcolor=DDDDDD >Observed</th>" \ 
       
   552        ,"</tr>" \
       
   553        /)
       
   554   table_footer = "</table>"
       
   555   row_header = "<tr>"
       
   556   row_footer = "</tr>"
       
   557 
       
   558   lines = new(50000,string)
       
   559   nline = 0
       
   560 
       
   561   set_line(lines,nline,header)
       
   562   set_line(lines,nline,table_header)
       
   563 ;-----------------------------------------------
       
   564 ; row of table
       
   565   
       
   566   do n = 0,nstation_long-1
       
   567 
       
   568      set_line(lines,nline,row_header)
       
   569 
       
   570      txt0 = station_sort(n)
       
   571      txt1 = sprintf("%5.2f", lat_ob_sort(n))
       
   572      txt2 = sprintf("%5.2f", lon_ob_sort(n))
       
   573      txt3 = year_ob_sort(n)
       
   574 
       
   575      set_line(lines,nline,"<th><a href="+txt0+"_tseries_vs_ob.png>"+txt0+"</a></th>")
       
   576      set_line(lines,nline,"<th>"+txt1+"</th>")
       
   577      set_line(lines,nline,"<th>"+txt2+"</th>")
       
   578      set_line(lines,nline,"<th>"+txt3+"</th>")
       
   579 
       
   580      set_line(lines,nline,row_footer)
       
   581   end do
       
   582 ;-----------------------------------------------
       
   583   set_line(lines,nline,table_footer)
       
   584   set_line(lines,nline,footer) 
       
   585 
       
   586 ; Now write to an HTML file.
       
   587   idx = ind(.not.ismissing(lines))
       
   588   if(.not.any(ismissing(idx))) then
       
   589     asciiwrite(output_html,lines(idx))
       
   590   else
       
   591    print ("error?")
       
   592   end if
       
   593   delete (idx)
       
   594 
       
   595 ;***************************************************************************
       
   596 ; output plots
       
   597 ;***************************************************************************
       
   598   output_dir = model_name+"/ameriflux"
       
   599 
       
   600   system("mv *.png *.html " + output_dir) 
       
   601 ;***************************************************************************
       
   602 end