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/20.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 ; merge 14 and 18

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

     4 ; add year_ob_i and year_ob_f

     5 ; change long_name

     6 ;************************************************************

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

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

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

    10 ;************************************************************

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

    12 begin

    13 ; add line to ascci/html file

    14     

    15   nnewlines = dimsizes(newlines)

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

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

    18     return

    19   end if 

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

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

    22   nline = nline + nnewlines

    23   return 

    24 end

    25 ;*************************************************************

    26 begin

    27 

    28   plot_type     = "ps"

    29   plot_type_new = "png"

    30 

    31 ; for 6 fields, 12-monthly

    32   nmon      = 12

    33 

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

    35 ; read model 

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

    37   model = "cn"

    38 ; model = "casa"

    39   

    40   model_name = "i01.10" + model

    41 

    42 ;************************************************

    43 ; read data: observed

    44 ;************************************************

    45 

    46  station = (/"ARM_Oklahoma" \

    47             ,"ARM_Oklahoma_burn" \

    48             ,"ARM_Oklahoma_control" \

    49             ,"Atqasuk" \

    50             ,"Audubon" \

    51             ,"AustinCary" \

    52             ,"Bartlett" \

    53             ,"Bondville" \

    54             ,"Brookings" \

    55             ,"Donaldson" \

    56             ,"Duke_Forest_Hardwoods" \

    57             ,"Duke_Forest_Open_Field" \

    58             ,"Duke_Forest_Pine" \

    59             ,"Fermi_Ag" \

    60             ,"Fermi_Prairie" \

    61             ,"Flagstaff_Managed" \

    62             ,"Flagstaff_Unmanaged" \

    63             ,"Flagstaff_Wildfire" \

    64             ,"FortPeck" \

    65             ,"FreemanRanch_mesquite" \

    66             ,"Goodwin_Creek" \

    67             ,"HarvardForest" \

    68             ,"HarvardForestHemlock" \

    69             ,"HowlandForestMain" \

    70             ,"HowlandForestWest" \

    71             ,"Ivotuk" \

    72             ,"KendallGrasslands" \

    73             ,"KennedySpaceCenterPine" \

    74             ,"KennedySpaceCenterScrub" \

    75             ,"LittleProspect" \

    76             ,"LostCreek" \

    77             ,"Mead-irrigated" \

    78             ,"Mead-irrigated-rotation" \

    79             ,"Mead-rainfed" \

    80             ,"Metolius_2nd_YoungPonderosaPine" \

    81             ,"MetoliusEyerly" \

    82             ,"MetoliusIntermediatePine" \

    83             ,"MetoliusOldPonderosaPine" \

    84             ,"MissouriOzark" \

    85             ,"Mize" \

    86             ,"MorganMonroe" \

    87             ,"NiwotRidge" \

    88             ,"NorthCarolina_cc" \

    89             ,"NorthCarolina_lp" \

    90             ,"ParkFalls" \

    91             ,"Rayonier" \

    92             ,"SantaRita" \

    93             ,"SkyOaks_Old" \

    94             ,"SkyOaks_PostFire" \

    95             ,"SkyOaks_Young" \

    96             ,"SylvaniaWilderness" \

    97             ,"Toledo" \

    98             ,"Tonzi" \

    99             ,"UCI_1850" \

   100             ,"UCI_1930" \

   101             ,"UCI_1964" \

   102             ,"UCI_1964wet" \

   103             ,"UCI_1981" \

   104             ,"UCI_1989" \

   105             ,"UCI_1998" \

   106             ,"UMBS" \

   107             ,"Vaira" \

   108             ,"WalkerBranch" \

   109             ,"WillowCreek" \

   110             ,"WindRiver" \

   111             ,"Wisconsin_ihw" \

   112             ,"Wisconsin_irp" \

   113             ,"Wisconsin_mrp" \

   114             ,"Wisconsin_myjp" \

   115             ,"Wisconsin_pb" \

   116             ,"Wisconsin_rpcc" \

   117             ,"Wisconsin_yhw" \

   118             ,"Wisconsin_yjp" \

   119             ,"Wisconsin_yrp" \

   120             /)

   121 

   122  year_ob = (/"2003-2006" \

   123             ,"2005-2006" \

   124             ,"2005-2006" \

   125             ,"1999-2006" \

   126             ,"2002-2006" \

   127             ,"2000-2005" \

   128             ,"2004-2005" \

   129             ,"1996-2006" \

   130             ,"2004-2006" \

   131             ,"1999-2004" \

   132             ,"2003-2005" \

   133             ,"2001-2005" \

   134             ,"2001-2005" \

   135             ,"2005-2006" \

   136             ,"2004-2006" \

   137             ,"2005-2006" \

   138             ,"2005-2006" \

   139             ,"2005-2006" \

   140             ,"2000-2006" \

   141             ,"2004-2006" \

   142             ,"2002-2006" \

   143             ,"1991-2004" \

   144             ,"2004-2004" \

   145             ,"1996-2004" \

   146             ,"1999-2004" \

   147             ,"2003-2006" \

   148             ,"2004-2006" \

   149             ,"2002-2002" \

   150             ,"2000-2006" \

   151             ,"2002-2005" \

   152             ,"2001-2005" \

   153             ,"2001-2005" \

   154             ,"2001-2005" \

   155             ,"2001-2005" \

   156             ,"2004-2005" \

   157             ,"2004-2005" \

   158             ,"2003-2005" \

   159             ,"1996-2000" \

   160             ,"2004-2006" \

   161             ,"1998-2004" \

   162             ,"1999-2005" \

   163             ,"1999-2003" \

   164             ,"2005-2006" \

   165             ,"2005-2006" \

   166             ,"1996-2003" \

   167             ,"1998-1998" \

   168             ,"2004-2006" \

   169             ,"1997-2006" \

   170             ,"2004-2006" \

   171             ,"1997-2006" \

   172             ,"2002-2006" \

   173             ,"2004-2005" \

   174             ,"2001-2006" \

   175             ,"2004-2005" \

   176             ,"2001-2005" \

   177             ,"2001-2005" \

   178             ,"2002-2004" \

   179             ,"2001-2005" \

   180             ,"2001-2005" \

   181             ,"2002-2005" \

   182             ,"1999-2003" \

   183             ,"2001-2006" \

   184             ,"1995-1999" \

   185             ,"1999-2005" \

   186             ,"1999-2004" \

   187             ,"2003-2003" \

   188             ,"2003-2003" \

   189             ,"2002-2005" \

   190             ,"2004-2004" \

   191             ,"2002-2002" \

   192             ,"2005-2005" \

   193             ,"2002-2002" \

   194             ,"2004-2005" \

   195             ,"2002-2002" \

   196             /)

   197 

   198  year_ob_i = (/2003 \

   199             ,2005 \

   200             ,2005 \

   201             ,1999 \

   202             ,2002 \

   203             ,2000 \

   204             ,2004 \

   205             ,1996 \

   206             ,2004 \

   207             ,1999 \

   208             ,2003 \

   209             ,2001 \

   210             ,2001 \

   211             ,2005 \

   212             ,2004 \

   213             ,2005 \

   214             ,2005 \

   215             ,2005 \

   216             ,2000 \

   217             ,2004 \

   218             ,2002 \

   219             ,1991 \

   220             ,2004 \

   221             ,1996 \

   222             ,1999 \

   223             ,2003 \

   224             ,2004 \

   225             ,2002 \

   226             ,2000 \

   227             ,2002 \

   228             ,2001 \

   229             ,2001 \

   230             ,2001 \

   231             ,2001 \

   232             ,2004 \

   233             ,2004 \

   234             ,2003 \

   235             ,1996 \

   236             ,2004 \

   237             ,1998 \

   238             ,1999 \

   239             ,1999 \

   240             ,2005 \

   241             ,2005 \

   242             ,1996 \

   243             ,1998 \

   244             ,2004 \

   245             ,1997 \

   246             ,2004 \

   247             ,1997 \

   248             ,2002 \

   249             ,2004 \

   250             ,2001 \

   251             ,2004 \

   252             ,2001 \

   253             ,2001 \

   254             ,2002 \

   255             ,2001 \

   256             ,2001 \

   257             ,2002 \

   258             ,1999 \

   259             ,2001 \

   260             ,1995 \

   261             ,1999 \

   262             ,1999 \

   263             ,2003 \

   264             ,2003 \

   265             ,2002 \

   266             ,2004 \

   267             ,2002 \

   268             ,2005 \

   269             ,2002 \

   270             ,2004 \

   271             ,2002 \

   272             /)

   273 

   274  year_ob_f = (/2006 \

   275             ,2006 \

   276             ,2006 \

   277             ,2006 \

   278             ,2006 \

   279             ,2005 \

   280             ,2005 \

   281             ,2006 \

   282             ,2006 \

   283             ,2004 \

   284             ,2005 \

   285             ,2005 \

   286             ,2005 \

   287             ,2006 \

   288             ,2006 \

   289             ,2006 \

   290             ,2006 \

   291             ,2006 \

   292             ,2006 \

   293             ,2006 \

   294             ,2006 \

   295             ,2004 \

   296             ,2004 \

   297             ,2004 \

   298             ,2004 \

   299             ,2006 \

   300             ,2006 \

   301             ,2002 \

   302             ,2006 \

   303             ,2005 \

   304             ,2005 \

   305             ,2005 \

   306             ,2005 \

   307             ,2005 \

   308             ,2005 \

   309             ,2005 \

   310             ,2005 \

   311             ,2000 \

   312             ,2006 \

   313             ,2004 \

   314             ,2005 \

   315             ,2003 \

   316             ,2006 \

   317             ,2006 \

   318             ,2003 \

   319             ,1998 \

   320             ,2006 \

   321             ,2006 \

   322             ,2006 \

   323             ,2006 \

   324             ,2006 \

   325             ,2005 \

   326             ,2006 \

   327             ,2005 \

   328             ,2005 \

   329             ,2005 \

   330             ,2004 \

   331             ,2005 \

   332             ,2005 \

   333             ,2005 \

   334             ,2003 \

   335             ,2006 \

   336             ,1999 \

   337             ,2005 \

   338             ,2004 \

   339             ,2003 \

   340             ,2003 \

   341             ,2005 \

   342             ,2004 \

   343             ,2002 \

   344             ,2005 \

   345             ,2002 \

   346             ,2005 \

   347             ,2002 \

   348             /)

   349 

   350  field   = (/"NEE Flux" \

   351             ,"Net Radiation" \

   352             ,"Latent Heat" \

   353             ,"Sensible Heat" \

   354             ,"GPP Flux" \

   355             ,"Respiration" \

   356             /)

   357 

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

   359                ,"W/m2" \

   360                ,"W/m2" \

   361                ,"W/m2" \

   362                ,"u mol/m2/s" \

   363                ,"u mol/m2/s" \

   364                /)

   365 

   366  nstation = dimsizes(station)

   367  nfield   = dimsizes(field)

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

   369 ;   get ob data

   370 

   371     nmon = 12

   372 

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

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

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

   376 

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

   378 

   379  do n = 0,nstation-1

   380 

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

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

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

   384 

   385     print (filo)

   386  

   387     lon_ob(n) = fo->lon 

   388     lat_ob(n) = fo->lat

   389 

   390     data      = fo->NEE_or_fMDS

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

   392 

   393     data      = fo->Rg_f

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

   395 

   396     data      = fo->LE_f

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

   398 

   399     data      = fo->H_f

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

   401 

   402     data      = fo->GPP_or_MDS

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

   404 

   405     data      = fo->Reco_or

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

   407 

   408     delete (data)

   409     delete (fo)

   410 

   411   end do

   412 

   413 ;---------------------------------------------------

   414 ;   get model data

   415 

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

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

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

   419 

   420 ;   print (film)

   421 

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

   423 ;   Model_units [=] gC/m2/s

   424 ;   12. = molecular weight of C

   425 ;   u mol = 1e-6 mol

   426     factor = 1e6 /12.

   427 

   428     ENERGY ="new"

   429 

   430     xm   = fm->lon

   431     ym   = fm->lat 

   432 

   433     nlat = dimsizes(ym)

   434     nlon = dimsizes(xm)

   435 

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

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

   438 

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

   440 

   441   data = fm->NEE

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

   443 

   444 ; data  = fm->LATENT

   445   data1 = fm->FCEV

   446   data2 = fm->FCTR

   447   data3 = fm->FGEV

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

   449   delete (data1)

   450   delete (data2)

   451   delete (data3)

   452 

   453 ; data = fm->SENSIBLE

   454   data  = fm->FSH

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

   456 

   457 ; data  = fm->NETRAD

   458   data1 = fm->FSA

   459   data2 = fm->FIRA

   460   data_mod0(1,:,:,:) = data1(:,:,:)-data2(:,:,:) 

   461   delete (data1)

   462   delete (data2)

   463 

   464 else

   465 

   466   data = fm->NEE

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

   468 

   469   data = fm->NETRAD

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

   471 

   472   data = fm->LATENT

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

   474 

   475 ; data = fm->SENSIBLE

   476   data = fm->FSH

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

   478 end if

   479 

   480   data = fm->GPP

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

   482 

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

   484      data = fm->ER

   485   else

   486      data1 = fm->AR

   487      data2 = fm->HR

   488      data  = data1 + data2

   489     

   490      delete (data1)

   491      delete (data2)

   492   end if

   493 

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

   495 

   496   delete (data)

   497   delete (fm) 

   498 

   499 ;************************************************************

   500 ; interpolate model data into observed station

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

   502 ;************************************************************

   503 

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

   505   ym(0) = -90.  

   506 

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

   508 

   509 ; print (yy)

   510  

   511   do n = 0, nstation-1

   512      data_mod(n,:,:) = yy(:,:,n)

   513   end do

   514 

   515   delete (data_mod0)

   516   delete (yy)

   517 

   518 ;************************************************************

   519 ; compute correlation coef and M score

   520 ;************************************************************

   521 

   522  score_max = 5.

   523 

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

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

   526 

   527  do n=0,nstation-1

   528  do m=0,nfield-1   

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

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

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

   532  end do

   533  end do

   534 

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

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

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

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

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

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

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

   542 

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

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

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

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

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

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

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

   550 

   551 ;*******************************************************************

   552 ; for station line plot

   553 ;*******************************************************************

   554 

   555 ; for x-axis in xyplot

   556   mon = ispan(1,12,1)

   557   mon@long_name = "month"

   558 

   559   res                   = True               ; plot mods desired

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

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

   562 ;-------------------------------------------------------------------------

   563 ; Add a boxed legend using the more simple method

   564 

   565   res@pmLegendDisplayMode    = "Always"

   566 ; res@pmLegendWidthF         = 0.1

   567   res@pmLegendWidthF         = 0.08

   568   res@pmLegendHeightF        = 0.06

   569 ; res@pmLegendOrthogonalPosF = -1.17

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

   571   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   572 

   573 ; res@pmLegendParallelPosF   =  0.18

   574   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   575 

   576 ; res@lgPerimOn             = False

   577   res@lgLabelFontHeightF     = 0.015

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

   579 ;-------------------------------------------------------------------

   580 ; for panel plot

   581   res@gsnFrame     = False                   ; Do not draw plot 

   582   res@gsnDraw      = False                   ; Do not advance frame

   583 

   584   pres                            = True     ; panel plot mods desired

   585   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   586                                              ; indiv. plots in panel

   587   pres@gsnMaximize                = True     ; fill the page

   588 ;-------------------------------------------------------------------

   589 

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

   591   plot_data!0 = "case"

   592   plot_data!1 = "month"

   593 

   594   do n = 0,nstation-1

   595 ;----------------------------

   596 ; for observed

   597 

   598     plot_name = station(n)+"_ob"    

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

   600     res@tiMainString = title

   601 

   602     wks = gsn_open_wks (plot_type,plot_name)

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

   604 

   605     do i = 0,nfield-1                           

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

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

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

   609     end do

   610 

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

   612 

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

   614            "rm "+plot_name+"."+plot_type)

   615 

   616     clear (wks)  

   617     delete (plot)

   618 ;----------------------------

   619 ; for model_vs_ob

   620 

   621     plot_name = station(n)+"_model_vs_ob"

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

   623     res@tiMainString = title

   624 

   625     wks = gsn_open_wks (plot_type,plot_name)

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

   627 

   628     do i = 0,nfield-1                           

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

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

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

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

   633     end do

   634    

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

   636 

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

   638            "rm "+plot_name+"."+plot_type)

   639 

   640     clear (wks)  

   641     delete (plot)

   642  end do

   643 

   644 ;###################################################################

   645 ; for the following tables,

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

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

   648 

   649   isort = dim_pqsort(lat_ob,-1)

   650 

   651   station_sort = station(isort)

   652   lat_ob_sort  = lat_ob(isort)

   653   lon_ob_sort  = lon_ob(isort)

   654   year_ob_sort = year_ob(isort)

   655   M_score_sort = M_score(isort,:)

   656  

   657   print(isort)

   658   print(lat_ob_sort)

   659 ;###################################################################

   660 ;*******************************************************************

   661 ; html table of site: observed

   662 ;*******************************************************************

   663   output_html = "line_ob.html"

   664 

   665   header = (/"<HTML>" \

   666             ,"<HEAD>" \

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

   668             ,"</HEAD>" \

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

   670             /) 

   671   footer = "</HTML>"

   672 

   673   table_header = (/ \

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

   675        ,"<tr>" \

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

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

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

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

   680        ,"</tr>" \

   681        /)

   682   table_footer = "</table>"

   683   row_header = "<tr>"

   684   row_footer = "</tr>"

   685 

   686   lines = new(50000,string)

   687   nline = 0

   688 

   689   set_line(lines,nline,header)

   690   set_line(lines,nline,table_header)

   691 ;-----------------------------------------------

   692 ; row of table

   693   

   694   do n = 0,nstation-1

   695      set_line(lines,nline,row_header)

   696 

   697      txt0 = station_sort(n)

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

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

   700      txt3 = year_ob_sort(n)

   701 

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

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

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

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

   706 

   707      set_line(lines,nline,row_footer)

   708   end do

   709 ;-----------------------------------------------

   710   set_line(lines,nline,table_footer)

   711   set_line(lines,nline,footer) 

   712 

   713 ; Now write to an HTML file.

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

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

   716     asciiwrite(output_html,lines(idx))

   717   else

   718    print ("error?")

   719   end if

   720   delete (idx)

   721 

   722 ;*******************************************************************

   723 ; score and line table : model vs observed

   724 ;*******************************************************************

   725   output_html = "score+line_vs_ob.html"

   726 

   727   header = (/"<HTML>" \

   728             ,"<HEAD>" \

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

   730             ,"</HEAD>" \

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

   732             /) 

   733   footer = "</HTML>"

   734 

   735   delete (table_header)

   736   table_header = (/ \

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

   738        ,"<tr>" \

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

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

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

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

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

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

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

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

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

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

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

   750        ,"</tr>" \

   751        /)

   752   table_footer = "</table>"

   753   row_header = "<tr>"

   754   row_footer = "</tr>"

   755 

   756   lines = new(50000,string)

   757   nline = 0

   758 

   759   set_line(lines,nline,header)

   760   set_line(lines,nline,table_header)

   761 ;-----------------------------------------------

   762 ; row of table

   763   

   764   do n = 0,nstation-1

   765      set_line(lines,nline,row_header)

   766 

   767      txt0  = station_sort(n)

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

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

   770      txt3  = year_ob_sort(n)

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

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

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

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

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

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

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

   778 

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

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

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

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

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

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

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

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

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

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

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

   790 

   791      set_line(lines,nline,row_footer)

   792   end do

   793 

   794 ; last row, summary

   795   set_line(lines,nline,row_header)

   796 

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

   798   txt1  = "-"

   799   txt2  = "-"

   800   txt3  = "-"

   801   txt4  = M_energy_co2

   802   txt5  = M_energy_rad

   803   txt6  = M_energy_lh

   804   txt7  = M_energy_sh

   805   txt8  = M_energy_gpp

   806   txt9  = M_energy_er

   807   txt10 = M_energy_all

   808 

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

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

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

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

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

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

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

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

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

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

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

   820 

   821   set_line(lines,nline,row_footer)

   822 ;-----------------------------------------------

   823   set_line(lines,nline,table_footer)

   824   set_line(lines,nline,footer) 

   825 

   826 ; Now write to an HTML file.

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

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

   829     asciiwrite(output_html,lines(idx))

   830   else

   831    print ("error?")

   832   end if

   833   delete (idx)

   834 

   835 ;***************************************************************************

   836 ; output plots

   837 ;***************************************************************************

   838   output_dir = model_name+"/ameriflux"

   839 

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

   841 ;***************************************************************************

   842 end
clamp