Mercurial Mercurial > repos > public_software > clamp / file revision
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
ameriflux/14.plot.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Mon, 26 Jan 2009 22:08:20 -0500
changeset 0 0c6405ab2ff4
permissions -rw-r--r--
Initial commit of C-LAMP Diagnostics from Jeff Lee 
     1 ;************************************************************

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

     3 ; add year_ob_i and year_ob_f

     4 ; change long_name

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

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

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

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

     9 ;************************************************************

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

    11 begin

    12 ; add line to ascci/html file

    13     

    14   nnewlines = dimsizes(newlines)

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

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

    17     return

    18   end if 

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

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

    21   nline = nline + nnewlines

    22   return 

    23 end

    24 ;*************************************************************

    25 begin

    26 

    27   plot_type     = "ps"

    28   plot_type_new = "png"

    29 

    30 ; for 6 fields, 12-monthly

    31   nmon      = 12

    32   nfield    = 6

    33 

    34 ;************************************************

    35 ; read model data

    36 ;************************************************

    37   model = "cn"

    38 ; model = "casa"

    39   

    40   model_name = "i01.10" + model

    41 

    42   dirm  = "/fis/cgd/cseg/people/jeff/surface_data/"

    43   film  = "lnd_T42.nc"

    44   fm    = addfile(dirm+film,"r")

    45   

    46   xm    = fm->lon

    47   ym    = fm->lat

    48 ;------------------------------------------------

    49   nlat = dimsizes(ym)

    50   nlon = dimsizes(xm)

    51 

    52   data_mod0 = new ((/nfield,nmon,nlat,nlon/),float)

    53 

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

    55 ; Model_units [=] gC/m2/s

    56 ; 12. = molecular weight of C

    57 ; u mol = 1e-6 mol

    58   factor = 1e6 /12.

    59 

    60   ENERGY ="new"

    61 

    62 ;************************************************

    63 ; read data: observed

    64 ;************************************************

    65 

    66  station = (/"ARM_Oklahoma" \

    67             ,"ARM_Oklahoma_burn" \

    68             ,"ARM_Oklahoma_control" \

    69             ,"Atqasuk" \

    70             ,"Audubon" \

    71             ,"AustinCary" \

    72             ,"Bartlett" \

    73             ,"Bondville" \

    74             ,"Brookings" \

    75             ,"Donaldson" \

    76             ,"Duke_Forest_Hardwoods" \

    77             ,"Duke_Forest_Open_Field" \

    78             ,"Duke_Forest_Pine" \

    79             ,"Fermi_Ag" \

    80             ,"Fermi_Prairie" \

    81             ,"Flagstaff_Managed" \

    82             ,"Flagstaff_Unmanaged" \

    83             ,"Flagstaff_Wildfire" \

    84             ,"FortPeck" \

    85             ,"FreemanRanch_mesquite" \

    86             ,"Goodwin_Creek" \

    87             ,"HarvardForest" \

    88             ,"HarvardForestHemlock" \

    89             ,"HowlandForestMain" \

    90             ,"HowlandForestWest" \

    91             ,"Ivotuk" \

    92             ,"KendallGrasslands" \

    93             ,"KennedySpaceCenterPine" \

    94             ,"KennedySpaceCenterScrub" \

    95             ,"LittleProspect" \

    96             ,"LostCreek" \

    97             ,"Mead-irrigated" \

    98             ,"Mead-irrigated-rotation" \

    99             ,"Mead-rainfed" \

   100             ,"Metolius_2nd_YoungPonderosaPine" \

   101             ,"MetoliusEyerly" \

   102             ,"MetoliusIntermediatePine" \

   103             ,"MetoliusOldPonderosaPine" \

   104             ,"MissouriOzark" \

   105             ,"Mize" \

   106             ,"MorganMonroe" \

   107             ,"NiwotRidge" \

   108             ,"NorthCarolina_cc" \

   109             ,"NorthCarolina_lp" \

   110             ,"ParkFalls" \

   111             ,"Rayonier" \

   112             ,"SantaRita" \

   113             ,"SkyOaks_Old" \

   114             ,"SkyOaks_PostFire" \

   115             ,"SkyOaks_Young" \

   116             ,"SylvaniaWilderness" \

   117             ,"Toledo" \

   118             ,"Tonzi" \

   119             ,"UCI_1850" \

   120             ,"UCI_1930" \

   121             ,"UCI_1964" \

   122             ,"UCI_1964wet" \

   123             ,"UCI_1981" \

   124             ,"UCI_1989" \

   125             ,"UCI_1998" \

   126             ,"UMBS" \

   127             ,"Vaira" \

   128             ,"WalkerBranch" \

   129             ,"WillowCreek" \

   130             ,"WindRiver" \

   131             ,"Wisconsin_ihw" \

   132             ,"Wisconsin_irp" \

   133             ,"Wisconsin_mrp" \

   134             ,"Wisconsin_myjp" \

   135             ,"Wisconsin_pb" \

   136             ,"Wisconsin_rpcc" \

   137             ,"Wisconsin_yhw" \

   138             ,"Wisconsin_yjp" \

   139             ,"Wisconsin_yrp" \

   140             /)

   141 

   142  year_ob = (/"2003-2006" \

   143             ,"2005-2006" \

   144             ,"2005-2006" \

   145             ,"1999-2006" \

   146             ,"2002-2006" \

   147             ,"2000-2005" \

   148             ,"2004-2005" \

   149             ,"1996-2006" \

   150             ,"2004-2006" \

   151             ,"1999-2004" \

   152             ,"2003-2005" \

   153             ,"2001-2005" \

   154             ,"2001-2005" \

   155             ,"2005-2006" \

   156             ,"2004-2006" \

   157             ,"2005-2006" \

   158             ,"2005-2006" \

   159             ,"2005-2006" \

   160             ,"2000-2006" \

   161             ,"2004-2006" \

   162             ,"2002-2006" \

   163             ,"1991-2004" \

   164             ,"2004-2004" \

   165             ,"1996-2004" \

   166             ,"1999-2004" \

   167             ,"2003-2006" \

   168             ,"2004-2006" \

   169             ,"2002-2002" \

   170             ,"2000-2006" \

   171             ,"2002-2005" \

   172             ,"2001-2005" \

   173             ,"2001-2005" \

   174             ,"2001-2005" \

   175             ,"2001-2005" \

   176             ,"2004-2005" \

   177             ,"2004-2005" \

   178             ,"2003-2005" \

   179             ,"1996-2000" \

   180             ,"2004-2006" \

   181             ,"1998-2004" \

   182             ,"1999-2005" \

   183             ,"1999-2003" \

   184             ,"2005-2006" \

   185             ,"2005-2006" \

   186             ,"1996-2003" \

   187             ,"1998-1998" \

   188             ,"2004-2006" \

   189             ,"1997-2006" \

   190             ,"2004-2006" \

   191             ,"1997-2006" \

   192             ,"2002-2006" \

   193             ,"2004-2005" \

   194             ,"2001-2006" \

   195             ,"2004-2005" \

   196             ,"2001-2005" \

   197             ,"2001-2005" \

   198             ,"2002-2004" \

   199             ,"2001-2005" \

   200             ,"2001-2005" \

   201             ,"2002-2005" \

   202             ,"1999-2003" \

   203             ,"2001-2006" \

   204             ,"1995-1999" \

   205             ,"1999-2005" \

   206             ,"1999-2004" \

   207             ,"2003-2003" \

   208             ,"2003-2003" \

   209             ,"2002-2005" \

   210             ,"2004-2004" \

   211             ,"2002-2002" \

   212             ,"2005-2005" \

   213             ,"2002-2002" \

   214             ,"2004-2005" \

   215             ,"2002-2002" \

   216             /)

   217 

   218  year_ob_i = (/2003 \

   219             ,2005 \

   220             ,2005 \

   221             ,1999 \

   222             ,2002 \

   223             ,2000 \

   224             ,2004 \

   225             ,1996 \

   226             ,2004 \

   227             ,1999 \

   228             ,2003 \

   229             ,2001 \

   230             ,2001 \

   231             ,2005 \

   232             ,2004 \

   233             ,2005 \

   234             ,2005 \

   235             ,2005 \

   236             ,2000 \

   237             ,2004 \

   238             ,2002 \

   239             ,1991 \

   240             ,2004 \

   241             ,1996 \

   242             ,1999 \

   243             ,2003 \

   244             ,2004 \

   245             ,2002 \

   246             ,2000 \

   247             ,2002 \

   248             ,2001 \

   249             ,2001 \

   250             ,2001 \

   251             ,2001 \

   252             ,2004 \

   253             ,2004 \

   254             ,2003 \

   255             ,1996 \

   256             ,2004 \

   257             ,1998 \

   258             ,1999 \

   259             ,1999 \

   260             ,2005 \

   261             ,2005 \

   262             ,1996 \

   263             ,1998 \

   264             ,2004 \

   265             ,1997 \

   266             ,2004 \

   267             ,1997 \

   268             ,2002 \

   269             ,2004 \

   270             ,2001 \

   271             ,2004 \

   272             ,2001 \

   273             ,2001 \

   274             ,2002 \

   275             ,2001 \

   276             ,2001 \

   277             ,2002 \

   278             ,1999 \

   279             ,2001 \

   280             ,1995 \

   281             ,1999 \

   282             ,1999 \

   283             ,2003 \

   284             ,2003 \

   285             ,2002 \

   286             ,2004 \

   287             ,2002 \

   288             ,2005 \

   289             ,2002 \

   290             ,2004 \

   291             ,2002 \

   292             /)

   293 

   294  year_ob_f = (/2006 \

   295             ,2006 \

   296             ,2006 \

   297             ,2006 \

   298             ,2006 \

   299             ,2005 \

   300             ,2005 \

   301             ,2006 \

   302             ,2006 \

   303             ,2004 \

   304             ,2005 \

   305             ,2005 \

   306             ,2005 \

   307             ,2006 \

   308             ,2006 \

   309             ,2006 \

   310             ,2006 \

   311             ,2006 \

   312             ,2006 \

   313             ,2006 \

   314             ,2006 \

   315             ,2004 \

   316             ,2004 \

   317             ,2004 \

   318             ,2004 \

   319             ,2006 \

   320             ,2006 \

   321             ,2002 \

   322             ,2006 \

   323             ,2005 \

   324             ,2005 \

   325             ,2005 \

   326             ,2005 \

   327             ,2005 \

   328             ,2005 \

   329             ,2005 \

   330             ,2005 \

   331             ,2000 \

   332             ,2006 \

   333             ,2004 \

   334             ,2005 \

   335             ,2003 \

   336             ,2006 \

   337             ,2006 \

   338             ,2003 \

   339             ,1998 \

   340             ,2006 \

   341             ,2006 \

   342             ,2006 \

   343             ,2006 \

   344             ,2006 \

   345             ,2005 \

   346             ,2006 \

   347             ,2005 \

   348             ,2005 \

   349             ,2005 \

   350             ,2004 \

   351             ,2005 \

   352             ,2005 \

   353             ,2005 \

   354             ,2003 \

   355             ,2006 \

   356             ,1999 \

   357             ,2005 \

   358             ,2004 \

   359             ,2003 \

   360             ,2003 \

   361             ,2005 \

   362             ,2004 \

   363             ,2002 \

   364             ,2005 \

   365             ,2002 \

   366             ,2005 \

   367             ,2002 \

   368             /)

   369 

   370  field   = (/"NEE Flux" \

   371             ,"Net Radiation" \

   372             ,"Latent Heat" \

   373             ,"Sensible Heat" \

   374             ,"GPP Flux" \

   375             ,"Respiration" \

   376             /)

   377 

   378  field_unit = (/"u mol/m2/s" \

   379                ,"W/m2" \

   380                ,"W/m2" \

   381                ,"W/m2" \

   382                ,"u mol/m2/s" \

   383                ,"u mol/m2/s" \

   384                /)

   385 

   386  nstation  = dimsizes(station)

   387 

   388  data_mod  = new ((/nstation, nfield, nmon/),float)

   389  data_ob   = new ((/nstation, nfield, nmon/),float)

   390  lat_ob    = new ((/nstation/),float)

   391  lon_ob    = new ((/nstation/),float)

   392  unit_ob   = new ((/nfield/),string)

   393 

   394  diro_root  = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"

   395  dirm_root  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

   396 

   397  do n = 0,nstation-1

   398 

   399 ;-------------------------------------------------

   400 ;   get ob data

   401 

   402     diro = diro_root + station(n)+"/"

   403     filo = year_ob_i(n)+"-"+year_ob_f(n)+"_L4_m.nc"

   404     fo   = addfile (diro+filo,"r")

   405 

   406     print (filo)

   407  

   408     lon_ob(n) = fo->lon 

   409     lat_ob(n) = fo->lat

   410 

   411     data      = fo->NEE_or_fMDS

   412     data_ob(n,0,:) = dim_avg(data(month|:,year|:))

   413     delete (data)

   414 

   415     data      = fo->Rg_f

   416     data_ob(n,1,:) = dim_avg(data(month|:,year|:))

   417     delete (data)

   418 

   419     data      = fo->LE_f

   420     data_ob(n,2,:) = dim_avg(data(month|:,year|:))

   421     delete (data)

   422 

   423     data      = fo->H_f

   424     data_ob(n,3,:) = dim_avg(data(month|:,year|:))

   425     delete (data)

   426 

   427     data      = fo->GPP_or_MDS

   428     data_ob(n,4,:) = dim_avg(data(month|:,year|:))

   429     delete (data)

   430 

   431     data      = fo->Reco_or

   432     data_ob(n,5,:) = dim_avg(data(month|:,year|:))

   433     delete (data)

   434 

   435     delete (fo)

   436 ;---------------------------------------------------

   437 ;   get model data

   438 

   439 ;   film = model_name+"_"+ year_ob(n)+"_MONS_climo.nc"

   440     film = model_name+"_"+"1990-2004"+"_MONS_climo.nc"

   441     fm   = addfile (dirm_root+film,"r")

   442 

   443 ;   print (film)

   444 

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

   446 

   447   data = fm->NEE

   448   data_mod0(0,:,:,:) = data(:,:,:) * factor

   449   delete (data)

   450 

   451 ; data  = fm->LATENT

   452   data1 = fm->FCEV

   453   data2 = fm->FCTR

   454   data3 = fm->FGEV

   455   data_mod0(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:) 

   456   delete (data1)

   457   delete (data2)

   458   delete (data3)

   459 

   460 ; data = fm->SENSIBLE

   461   data  = fm->FSH

   462   data_mod0(3,:,:,:) = data(:,:,:) 

   463   delete (data)

   464 

   465 ; data  = fm->NETRAD

   466   data1 = fm->FSA

   467   data2 = fm->FIRA

   468   data_mod0(1,:,:,:) = data1(:,:,:)-data2(:,:,:)-data_mod0(2,:,:,:)-data_mod0(3,:,:,:) 

   469   delete (data1)

   470   delete (data2)

   471 

   472 else

   473 

   474   data = fm->NEE

   475   data_mod0(0,:,:,:) = data(:,:,:) * factor

   476   delete (data)

   477 

   478   data = fm->NETRAD

   479   data_mod0(1,:,:,:) = data(:,:,:) 

   480   delete (data)

   481 

   482   data = fm->LATENT

   483   data_mod0(2,:,:,:) = data(:,:,:) 

   484   delete (data)

   485 

   486 ; data = fm->SENSIBLE

   487   data = fm->FSH

   488   data_mod0(3,:,:,:) = data(:,:,:) 

   489   delete (data)

   490 end if

   491 

   492   data = fm->GPP

   493   data_mod0(4,:,:,:) = data(:,:,:) * factor 

   494   delete (data)

   495 

   496   if (model .eq. "cn") then

   497      data = fm->ER

   498   else

   499      data1 = fm->AR

   500      data2 = fm->HR

   501      data  = data1 + data2

   502     

   503      delete (data1)

   504      delete (data2)

   505   end if

   506 

   507   data_mod0(5,:,:,:) = data(:,:,:) * factor

   508   delete (data)

   509 

   510   delete (fm) 

   511 

   512 ;************************************************************

   513 ; interpolate model data into observed station

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

   515 ;************************************************************

   516 

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

   518   ym(0) = -90.  

   519 

   520   yy = linint2_points_Wrap(xm,ym,data_mod0,True,lon_ob(n),lat_ob(n),0)

   521 

   522 ; print (yy)

   523  

   524   data_mod(n,:,:) = yy(:,:,0)

   525 

   526   delete (yy)

   527 

   528  end do

   529 

   530 ;************************************************************

   531 ; compute correlation coef and M score

   532 ;************************************************************

   533 

   534  score_max = 5.

   535 

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

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

   538 

   539  do n=0,nstation-1

   540  do m=0,nfield-1   

   541     ccr(n,m) = esccr(data_ob(n,m,:),data_mod(n,m,:),0)

   542     bias = sum(abs(data_mod(n,m,:)-data_ob(n,m,:))/(abs(data_mod(n,m,:))+abs(data_ob(n,m,:))))

   543     M_score(n,m) = (1. -(bias/nmon)) * score_max

   544  end do

   545  end do

   546 

   547  M_co2 = avg(M_score(:,0))

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

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

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

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

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

   553  M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er

   554 

   555  M_energy_co2 = sprintf("%.2f", M_co2)

   556  M_energy_rad = sprintf("%.2f", M_rad)

   557  M_energy_lh  = sprintf("%.2f", M_lh )

   558  M_energy_sh  = sprintf("%.2f", M_sh )

   559  M_energy_gpp = sprintf("%.2f", M_gpp)

   560  M_energy_er  = sprintf("%.2f", M_er )

   561  M_energy_all = sprintf("%.2f", M_all)

   562 

   563 ;*******************************************************************

   564 ; for station line plot

   565 ;*******************************************************************

   566 

   567 ; for x-axis in xyplot

   568   mon = ispan(1,12,1)

   569   mon@long_name = "month"

   570 

   571   res                   = True               ; plot mods desired

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

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

   574 ;-------------------------------------------------------------------------

   575 ; Add a boxed legend using the more simple method

   576 

   577   res@pmLegendDisplayMode    = "Always"

   578 ; res@pmLegendWidthF         = 0.1

   579   res@pmLegendWidthF         = 0.08

   580   res@pmLegendHeightF        = 0.06

   581 ; res@pmLegendOrthogonalPosF = -1.17

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

   583   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   584 

   585 ; res@pmLegendParallelPosF   =  0.18

   586   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   587 

   588 ; res@lgPerimOn             = False

   589   res@lgLabelFontHeightF     = 0.015

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

   591 ;-------------------------------------------------------------------

   592 ; for panel plot

   593   res@gsnFrame     = False                   ; Do not draw plot 

   594   res@gsnDraw      = False                   ; Do not advance frame

   595 

   596   pres                            = True     ; panel plot mods desired

   597   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   598                                              ; indiv. plots in panel

   599   pres@gsnMaximize                = True     ; fill the page

   600 ;-------------------------------------------------------------------

   601 

   602   plot_data   = new((/2,12/),float)

   603   plot_data!0 = "case"

   604   plot_data!1 = "month"

   605 

   606   do n = 0,nstation-1

   607 ;----------------------------

   608 ; for observed

   609 

   610     plot_name = station(n)+"_ob"    

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

   612     res@tiMainString = title

   613 

   614     wks = gsn_open_wks (plot_type,plot_name)

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

   616 

   617     do i = 0,nfield-1                           

   618        plot_data(0,:) = (/data_ob (n,i,:)/)

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

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

   621     end do

   622 

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

   624 

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

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

   627 

   628     clear (wks)  

   629     delete (plot)

   630 ;----------------------------

   631 ; for model_vs_ob

   632 

   633     plot_name = station(n)+"_model_vs_ob"

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

   635     res@tiMainString = title

   636 

   637     wks = gsn_open_wks (plot_type,plot_name)

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

   639 

   640     do i = 0,nfield-1                           

   641        plot_data(0,:) = (/data_ob (n,i,:)/)

   642        plot_data(1,:) = (/data_mod(n,i,:)/)

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

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

   645     end do

   646    

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

   648 

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

   650            "rm "+plot_name+"."+plot_type)

   651 

   652     clear (wks)  

   653     delete (plot)

   654  end do

   655 

   656 ;###################################################################

   657 ; for the following tables,

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

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

   660 

   661   isort = dim_pqsort(lat_ob,-1)

   662 

   663   station_sort = station(isort)

   664   lat_ob_sort  = lat_ob(isort)

   665   lon_ob_sort  = lon_ob(isort)

   666   year_ob_sort = year_ob(isort)

   667   M_score_sort = M_score(isort,:)

   668  

   669   print(isort)

   670   print(lat_ob_sort)

   671 ;###################################################################

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

   673 ; html table of site: observed

   674 ;*******************************************************************

   675   output_html = "line_ob.html"

   676 

   677   header = (/"<HTML>" \

   678             ,"<HEAD>" \

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

   680             ,"</HEAD>" \

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

   682             /) 

   683   footer = "</HTML>"

   684 

   685   table_header = (/ \

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

   687        ,"<tr>" \

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

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

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

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

   692        ,"</tr>" \

   693        /)

   694   table_footer = "</table>"

   695   row_header = "<tr>"

   696   row_footer = "</tr>"

   697 

   698   lines = new(50000,string)

   699   nline = 0

   700 

   701   set_line(lines,nline,header)

   702   set_line(lines,nline,table_header)

   703 ;-----------------------------------------------

   704 ; row of table

   705   

   706   do n = 0,nstation-1

   707      set_line(lines,nline,row_header)

   708 

   709      txt0 = station_sort(n)

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

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

   712      txt3 = year_ob_sort(n)

   713 

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

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

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

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

   718 

   719      set_line(lines,nline,row_footer)

   720   end do

   721 ;-----------------------------------------------

   722   set_line(lines,nline,table_footer)

   723   set_line(lines,nline,footer) 

   724 

   725 ; Now write to an HTML file.

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

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

   728     asciiwrite(output_html,lines(idx))

   729   else

   730    print ("error?")

   731   end if

   732   delete (idx)

   733 

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

   735 ; score and line table : model vs observed

   736 ;*******************************************************************

   737   output_html = "score+line_vs_ob.html"

   738 

   739   header = (/"<HTML>" \

   740             ,"<HEAD>" \

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

   742             ,"</HEAD>" \

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

   744             /) 

   745   footer = "</HTML>"

   746 

   747   delete (table_header)

   748   table_header = (/ \

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

   750        ,"<tr>" \

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

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

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

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

   755        ,"   <th bgcolor=DDDDDD >CO2 Flux</th>" \

   756        ,"   <th bgcolor=DDDDDD >Net Radiation</th>" \

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

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

   759        ,"   <th bgcolor=DDDDDD >GPP Glux</th>" \

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

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

   762        ,"</tr>" \

   763        /)

   764   table_footer = "</table>"

   765   row_header = "<tr>"

   766   row_footer = "</tr>"

   767 

   768   lines = new(50000,string)

   769   nline = 0

   770 

   771   set_line(lines,nline,header)

   772   set_line(lines,nline,table_header)

   773 ;-----------------------------------------------

   774 ; row of table

   775   

   776   do n = 0,nstation-1

   777      set_line(lines,nline,row_header)

   778 

   779      txt0  = station_sort(n)

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

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

   782      txt3  = year_ob_sort(n)

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

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

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

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

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

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

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

   790 

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

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

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

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

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

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

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

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

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

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

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

   802 

   803      set_line(lines,nline,row_footer)

   804   end do

   805 

   806 ; last row, summary

   807   set_line(lines,nline,row_header)

   808 

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

   810   txt1  = "-"

   811   txt2  = "-"

   812   txt3  = "-"

   813   txt4  = M_energy_co2

   814   txt5  = M_energy_rad

   815   txt6  = M_energy_lh

   816   txt7  = M_energy_sh

   817   txt8  = M_energy_gpp

   818   txt9  = M_energy_er

   819   txt10 = M_energy_all

   820 

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

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

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

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

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

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

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

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

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

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

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

   832 

   833   set_line(lines,nline,row_footer)

   834 ;-----------------------------------------------

   835   set_line(lines,nline,table_footer)

   836   set_line(lines,nline,footer) 

   837 

   838 ; Now write to an HTML file.

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

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

   841     asciiwrite(output_html,lines(idx))

   842   else

   843    print ("error?")

   844   end if

   845   delete (idx)

   846 

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

   848 ; output plots

   849 ;***************************************************************************

   850   output_dir = model_name+"/ameriflux"

   851 

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

   853 ;***************************************************************************

   854 end
clamp