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/16.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 ; find year_ob >= 4

     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   nfield    = 6

    34 

    35 ;************************************************

    36 ; read model data

    37 ;************************************************

    38   model = "cn"

    39 ; model = "casa"

    40   

    41   model_name = "i01.10" + model

    42 

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

    44   film  = "lnd_T42.nc"

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

    46   

    47   xm    = fm->lon

    48   ym    = fm->lat

    49 ;------------------------------------------------

    50   nlat = dimsizes(ym)

    51   nlon = dimsizes(xm)

    52 

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

    54   

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

    56 ; Model_units [=] gC/m2/s

    57 ; 12. = molecular weight of C

    58 ; u mol = 1e-6 mol

    59   factor = 1e6 /12.

    60 

    61   ENERGY ="new"

    62 

    63 ;************************************************

    64 ; read data: observed

    65 ;************************************************

    66 

    67  station = (/"ARM_Oklahoma" \

    68             ,"ARM_Oklahoma_burn" \

    69             ,"ARM_Oklahoma_control" \

    70             ,"Atqasuk" \

    71             ,"Audubon" \

    72             ,"AustinCary" \

    73             ,"Bartlett" \

    74             ,"Bondville" \

    75             ,"Brookings" \

    76             ,"Donaldson" \

    77             ,"Duke_Forest_Hardwoods" \

    78             ,"Duke_Forest_Open_Field" \

    79             ,"Duke_Forest_Pine" \

    80             ,"Fermi_Ag" \

    81             ,"Fermi_Prairie" \

    82             ,"Flagstaff_Managed" \

    83             ,"Flagstaff_Unmanaged" \

    84             ,"Flagstaff_Wildfire" \

    85             ,"FortPeck" \

    86             ,"FreemanRanch_mesquite" \

    87             ,"Goodwin_Creek" \

    88             ,"HarvardForest" \

    89             ,"HarvardForestHemlock" \

    90             ,"HowlandForestMain" \

    91             ,"HowlandForestWest" \

    92             ,"Ivotuk" \

    93             ,"KendallGrasslands" \

    94             ,"KennedySpaceCenterPine" \

    95             ,"KennedySpaceCenterScrub" \

    96             ,"LittleProspect" \

    97             ,"LostCreek" \

    98             ,"Mead-irrigated" \

    99             ,"Mead-irrigated-rotation" \

   100             ,"Mead-rainfed" \

   101             ,"Metolius_2nd_YoungPonderosaPine" \

   102             ,"MetoliusEyerly" \

   103             ,"MetoliusIntermediatePine" \

   104             ,"MetoliusOldPonderosaPine" \

   105             ,"MissouriOzark" \

   106             ,"Mize" \

   107             ,"MorganMonroe" \

   108             ,"NiwotRidge" \

   109             ,"NorthCarolina_cc" \

   110             ,"NorthCarolina_lp" \

   111             ,"ParkFalls" \

   112             ,"Rayonier" \

   113             ,"SantaRita" \

   114             ,"SkyOaks_Old" \

   115             ,"SkyOaks_PostFire" \

   116             ,"SkyOaks_Young" \

   117             ,"SylvaniaWilderness" \

   118             ,"Toledo" \

   119             ,"Tonzi" \

   120             ,"UCI_1850" \

   121             ,"UCI_1930" \

   122             ,"UCI_1964" \

   123             ,"UCI_1964wet" \

   124             ,"UCI_1981" \

   125             ,"UCI_1989" \

   126             ,"UCI_1998" \

   127             ,"UMBS" \

   128             ,"Vaira" \

   129             ,"WalkerBranch" \

   130             ,"WillowCreek" \

   131             ,"WindRiver" \

   132             ,"Wisconsin_ihw" \

   133             ,"Wisconsin_irp" \

   134             ,"Wisconsin_mrp" \

   135             ,"Wisconsin_myjp" \

   136             ,"Wisconsin_pb" \

   137             ,"Wisconsin_rpcc" \

   138             ,"Wisconsin_yhw" \

   139             ,"Wisconsin_yjp" \

   140             ,"Wisconsin_yrp" \

   141             /)

   142 

   143  year_ob = (/"2003-2006" \

   144             ,"2005-2006" \

   145             ,"2005-2006" \

   146             ,"1999-2006" \

   147             ,"2002-2006" \

   148             ,"2000-2005" \

   149             ,"2004-2005" \

   150             ,"1996-2006" \

   151             ,"2004-2006" \

   152             ,"1999-2004" \

   153             ,"2003-2005" \

   154             ,"2001-2005" \

   155             ,"2001-2005" \

   156             ,"2005-2006" \

   157             ,"2004-2006" \

   158             ,"2005-2006" \

   159             ,"2005-2006" \

   160             ,"2005-2006" \

   161             ,"2000-2006" \

   162             ,"2004-2006" \

   163             ,"2002-2006" \

   164             ,"1991-2004" \

   165             ,"2004-2004" \

   166             ,"1996-2004" \

   167             ,"1999-2004" \

   168             ,"2003-2006" \

   169             ,"2004-2006" \

   170             ,"2002-2002" \

   171             ,"2000-2006" \

   172             ,"2002-2005" \

   173             ,"2001-2005" \

   174             ,"2001-2005" \

   175             ,"2001-2005" \

   176             ,"2001-2005" \

   177             ,"2004-2005" \

   178             ,"2004-2005" \

   179             ,"2003-2005" \

   180             ,"1996-2000" \

   181             ,"2004-2006" \

   182             ,"1998-2004" \

   183             ,"1999-2005" \

   184             ,"1999-2003" \

   185             ,"2005-2006" \

   186             ,"2005-2006" \

   187             ,"1996-2003" \

   188             ,"1998-1998" \

   189             ,"2004-2006" \

   190             ,"1997-2006" \

   191             ,"2004-2006" \

   192             ,"1997-2006" \

   193             ,"2002-2006" \

   194             ,"2004-2005" \

   195             ,"2001-2006" \

   196             ,"2004-2005" \

   197             ,"2001-2005" \

   198             ,"2001-2005" \

   199             ,"2002-2004" \

   200             ,"2001-2005" \

   201             ,"2001-2005" \

   202             ,"2002-2005" \

   203             ,"1999-2003" \

   204             ,"2001-2006" \

   205             ,"1995-1999" \

   206             ,"1999-2005" \

   207             ,"1999-2004" \

   208             ,"2003-2003" \

   209             ,"2003-2003" \

   210             ,"2002-2005" \

   211             ,"2004-2004" \

   212             ,"2002-2002" \

   213             ,"2005-2005" \

   214             ,"2002-2002" \

   215             ,"2004-2005" \

   216             ,"2002-2002" \

   217             /)

   218 

   219  year_ob_i = (/2003 \

   220             ,2005 \

   221             ,2005 \

   222             ,1999 \

   223             ,2002 \

   224             ,2000 \

   225             ,2004 \

   226             ,1996 \

   227             ,2004 \

   228             ,1999 \

   229             ,2003 \

   230             ,2001 \

   231             ,2001 \

   232             ,2005 \

   233             ,2004 \

   234             ,2005 \

   235             ,2005 \

   236             ,2005 \

   237             ,2000 \

   238             ,2004 \

   239             ,2002 \

   240             ,1991 \

   241             ,2004 \

   242             ,1996 \

   243             ,1999 \

   244             ,2003 \

   245             ,2004 \

   246             ,2002 \

   247             ,2000 \

   248             ,2002 \

   249             ,2001 \

   250             ,2001 \

   251             ,2001 \

   252             ,2001 \

   253             ,2004 \

   254             ,2004 \

   255             ,2003 \

   256             ,1996 \

   257             ,2004 \

   258             ,1998 \

   259             ,1999 \

   260             ,1999 \

   261             ,2005 \

   262             ,2005 \

   263             ,1996 \

   264             ,1998 \

   265             ,2004 \

   266             ,1997 \

   267             ,2004 \

   268             ,1997 \

   269             ,2002 \

   270             ,2004 \

   271             ,2001 \

   272             ,2004 \

   273             ,2001 \

   274             ,2001 \

   275             ,2002 \

   276             ,2001 \

   277             ,2001 \

   278             ,2002 \

   279             ,1999 \

   280             ,2001 \

   281             ,1995 \

   282             ,1999 \

   283             ,1999 \

   284             ,2003 \

   285             ,2003 \

   286             ,2002 \

   287             ,2004 \

   288             ,2002 \

   289             ,2005 \

   290             ,2002 \

   291             ,2004 \

   292             ,2002 \

   293             /)

   294 

   295  year_ob_f = (/2006 \

   296             ,2006 \

   297             ,2006 \

   298             ,2006 \

   299             ,2006 \

   300             ,2005 \

   301             ,2005 \

   302             ,2006 \

   303             ,2006 \

   304             ,2004 \

   305             ,2005 \

   306             ,2005 \

   307             ,2005 \

   308             ,2006 \

   309             ,2006 \

   310             ,2006 \

   311             ,2006 \

   312             ,2006 \

   313             ,2006 \

   314             ,2006 \

   315             ,2006 \

   316             ,2004 \

   317             ,2004 \

   318             ,2004 \

   319             ,2004 \

   320             ,2006 \

   321             ,2006 \

   322             ,2002 \

   323             ,2006 \

   324             ,2005 \

   325             ,2005 \

   326             ,2005 \

   327             ,2005 \

   328             ,2005 \

   329             ,2005 \

   330             ,2005 \

   331             ,2005 \

   332             ,2000 \

   333             ,2006 \

   334             ,2004 \

   335             ,2005 \

   336             ,2003 \

   337             ,2006 \

   338             ,2006 \

   339             ,2003 \

   340             ,1998 \

   341             ,2006 \

   342             ,2006 \

   343             ,2006 \

   344             ,2006 \

   345             ,2006 \

   346             ,2005 \

   347             ,2006 \

   348             ,2005 \

   349             ,2005 \

   350             ,2005 \

   351             ,2004 \

   352             ,2005 \

   353             ,2005 \

   354             ,2005 \

   355             ,2003 \

   356             ,2006 \

   357             ,1999 \

   358             ,2005 \

   359             ,2004 \

   360             ,2003 \

   361             ,2003 \

   362             ,2005 \

   363             ,2004 \

   364             ,2002 \

   365             ,2005 \

   366             ,2002 \

   367             ,2005 \

   368             ,2002 \

   369             /)

   370 

   371  field   = (/"NEE Flux" \

   372             ,"Net Radiation" \

   373             ,"Latent Heat" \

   374             ,"Sensible Heat" \

   375             ,"GPP Flux" \

   376             ,"Respiration" \

   377             /)

   378 

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

   380                ,"W/m2" \

   381                ,"W/m2" \

   382                ,"W/m2" \

   383                ,"u mol/m2/s" \

   384                ,"u mol/m2/s" \

   385                /)

   386 

   387 ;========================================================================

   388 ; find # year observed >=4 and year_ob_i <= 2001

   389 

   390  i_long_ob = ind((year_ob_f - year_ob_i) .ge. 4 .and. year_ob_i .le. 2001)

   391  

   392  station_long = station(i_long_ob)

   393  year_ob_long = year_ob(i_long_ob)

   394  year_ob_i_long = year_ob_i(i_long_ob)

   395  year_ob_f_long = year_ob_f(i_long_ob)

   396 

   397 ;print (i_long_ob)

   398 ;print (year_ob_i(i_long_ob))

   399 ;print (station_long)

   400 ;========================================================================

   401 

   402  nstation  = dimsizes(station_long)

   403 

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

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

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

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

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

   409 

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

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

   412 

   413  do n = 0,nstation-1

   414 

   415 ;-------------------------------------------------

   416 ;   get ob data

   417 

   418     diro = diro_root + station_long(n)+"/"

   419     filo = year_ob_i_long(n)+"-"+year_ob_f_long(n)+"_L4_m.nc"

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

   421 

   422     print (filo)

   423  

   424     lon_ob(n) = fo->lon 

   425     lat_ob(n) = fo->lat

   426 

   427     data      = fo->NEE_or_fMDS

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

   429     delete (data)

   430 

   431     data      = fo->Rg_f

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

   433     delete (data)

   434 

   435     data      = fo->LE_f

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

   437     delete (data)

   438 

   439     data      = fo->H_f

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

   441     delete (data)

   442 

   443     data      = fo->GPP_or_MDS

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

   445     delete (data)

   446 

   447     data      = fo->Reco_or

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

   449     delete (data)

   450 

   451     delete (fo)

   452 ;---------------------------------------------------

   453 ;   get model data

   454 

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

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

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

   458 

   459 ;   print (film)

   460 

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

   462 

   463   data = fm->NEE

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

   465   delete (data)

   466 

   467 ; data  = fm->LATENT

   468   data1 = fm->FCEV

   469   data2 = fm->FCTR

   470   data3 = fm->FGEV

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

   472   delete (data1)

   473   delete (data2)

   474   delete (data3)

   475 

   476 ; data = fm->SENSIBLE

   477   data  = fm->FSH

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

   479   delete (data)

   480 

   481 ; data  = fm->NETRAD

   482   data1 = fm->FSA

   483   data2 = fm->FIRA

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

   485   delete (data1)

   486   delete (data2)

   487 

   488 else

   489 

   490   data = fm->NEE

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

   492   delete (data)

   493 

   494   data = fm->NETRAD

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

   496   delete (data)

   497 

   498   data = fm->LATENT

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

   500   delete (data)

   501 

   502 ; data = fm->SENSIBLE

   503   data = fm->FSH

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

   505   delete (data)

   506 end if

   507 

   508   data = fm->GPP

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

   510   delete (data)

   511 

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

   513      data = fm->ER

   514   else

   515      data1 = fm->AR

   516      data2 = fm->HR

   517      data  = data1 + data2

   518     

   519      delete (data1)

   520      delete (data2)

   521   end if

   522 

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

   524   delete (data)

   525 

   526   delete (fm) 

   527 

   528 ;************************************************************

   529 ; interpolate model data into observed station

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

   531 ;************************************************************

   532 

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

   534   ym(0) = -90.  

   535 

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

   537 

   538 ; print (yy)

   539  

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

   541 

   542   delete (yy)

   543 

   544  end do

   545 

   546 ;************************************************************

   547 ; compute correlation coef and M score

   548 ;************************************************************

   549 

   550  score_max = 5.

   551 

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

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

   554 

   555  do n=0,nstation-1

   556  do m=0,nfield-1   

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

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

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

   560  end do

   561  end do

   562 

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

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

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

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

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

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

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

   570 

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

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

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

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

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

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

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

   578 

   579 ;*******************************************************************

   580 ; for station line plot

   581 ;*******************************************************************

   582 

   583 ; for x-axis in xyplot

   584   mon = ispan(1,12,1)

   585   mon@long_name = "month"

   586 

   587   res                   = True               ; plot mods desired

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

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

   590 ;-------------------------------------------------------------------------

   591 ; Add a boxed legend using the more simple method

   592 

   593   res@pmLegendDisplayMode    = "Always"

   594 ; res@pmLegendWidthF         = 0.1

   595   res@pmLegendWidthF         = 0.08

   596   res@pmLegendHeightF        = 0.06

   597 ; res@pmLegendOrthogonalPosF = -1.17

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

   599   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   600 

   601 ; res@pmLegendParallelPosF   =  0.18

   602   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   603 

   604 ; res@lgPerimOn             = False

   605   res@lgLabelFontHeightF     = 0.015

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

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

   608 ; for panel plot

   609   res@gsnFrame     = False                   ; Do not draw plot 

   610   res@gsnDraw      = False                   ; Do not advance frame

   611 

   612   pres                            = True     ; panel plot mods desired

   613   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   614                                              ; indiv. plots in panel

   615   pres@gsnMaximize                = True     ; fill the page

   616 ;-------------------------------------------------------------------

   617 

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

   619   plot_data!0 = "case"

   620   plot_data!1 = "month"

   621 

   622   do n = 0,nstation-1

   623 

   624 ;----------------------------

   625 ; for observed

   626 

   627     plot_name = station_long(n)+"_ob"    

   628     title = station_long(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"

   629     res@tiMainString = title

   630 

   631     wks = gsn_open_wks (plot_type,plot_name)

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

   633 

   634     do i = 0,nfield-1                           

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

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

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

   638     end do

   639 

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

   641 

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

   643            "rm "+plot_name+"."+plot_type)

   644 

   645     clear (wks)  

   646     delete (plot)

   647 ;----------------------------

   648 ; for model_vs_ob

   649 

   650     plot_name = station_long(n)+"_model_vs_ob"

   651     title = station_long(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"    

   652     res@tiMainString = title

   653 

   654     wks = gsn_open_wks (plot_type,plot_name)

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

   656 

   657     do i = 0,nfield-1                           

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

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

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

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

   662     end do

   663    

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

   665 

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

   667            "rm "+plot_name+"."+plot_type)

   668 

   669     clear (wks)  

   670     delete (plot)

   671  end do

   672 

   673 ;###################################################################

   674 ; for the following tables,

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

   676 

   677   isort = dim_pqsort(lat_ob,-1)

   678 

   679   station_sort = station_long(isort)

   680   year_ob_sort = year_ob_long(isort)

   681   lat_ob_sort  = lat_ob(isort)

   682   lon_ob_sort  = lon_ob(isort)

   683   M_score_sort = M_score(isort,:)

   684  

   685   print(isort)

   686   print(lat_ob_sort)

   687 

   688 ;*******************************************************************

   689 ; html table of site: observed

   690 ;*******************************************************************

   691   output_html = "line_ob.html"

   692 

   693   header = (/"<HTML>" \

   694             ,"<HEAD>" \

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

   696             ,"</HEAD>" \

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

   698             /) 

   699   footer = "</HTML>"

   700 

   701   table_header = (/ \

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

   703        ,"<tr>" \

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

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

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

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

   708        ,"</tr>" \

   709        /)

   710   table_footer = "</table>"

   711   row_header = "<tr>"

   712   row_footer = "</tr>"

   713 

   714   lines = new(50000,string)

   715   nline = 0

   716 

   717   set_line(lines,nline,header)

   718   set_line(lines,nline,table_header)

   719 ;-----------------------------------------------

   720 ; row of table

   721   

   722   do n = 0,nstation-1

   723 

   724      set_line(lines,nline,row_header)

   725 

   726      txt0 = station_sort(n)

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

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

   729      txt3 = year_ob_sort(n)

   730 

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

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

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

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

   735 

   736      set_line(lines,nline,row_footer)

   737   end do

   738 ;-----------------------------------------------

   739   set_line(lines,nline,table_footer)

   740   set_line(lines,nline,footer) 

   741 

   742 ; Now write to an HTML file.

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

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

   745     asciiwrite(output_html,lines(idx))

   746   else

   747    print ("error?")

   748   end if

   749   delete (idx)

   750 

   751 ;*******************************************************************

   752 ; score and line table : model vs observed

   753 ;*******************************************************************

   754   output_html = "score+line_vs_ob.html"

   755 

   756   header = (/"<HTML>" \

   757             ,"<HEAD>" \

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

   759             ,"</HEAD>" \

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

   761             /) 

   762   footer = "</HTML>"

   763 

   764   delete (table_header)

   765   table_header = (/ \

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

   767        ,"<tr>" \

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

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

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

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

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

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

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

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

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

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

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

   779        ,"</tr>" \

   780        /)

   781   table_footer = "</table>"

   782   row_header = "<tr>"

   783   row_footer = "</tr>"

   784 

   785   lines = new(50000,string)

   786   nline = 0

   787 

   788   set_line(lines,nline,header)

   789   set_line(lines,nline,table_header)

   790 ;-----------------------------------------------

   791 ; row of table

   792   

   793   do n = 0,nstation-1

   794 

   795      set_line(lines,nline,row_header)

   796 

   797      txt0  = station_sort(n)

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

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

   800      txt3  = year_ob_sort(n)

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

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

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

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

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

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

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

   808 

   809      set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></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   end do

   823 

   824 ; last row, summary

   825   set_line(lines,nline,row_header)

   826 

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

   828   txt1  = "-"

   829   txt2  = "-"

   830   txt3  = "-"

   831   txt4  = M_energy_co2

   832   txt5  = M_energy_rad

   833   txt6  = M_energy_lh

   834   txt7  = M_energy_sh

   835   txt8  = M_energy_gpp

   836   txt9  = M_energy_er

   837   txt10 = M_energy_all

   838 

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

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

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

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

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

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

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

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

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

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

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

   850 

   851   set_line(lines,nline,row_footer)

   852 ;-----------------------------------------------

   853   set_line(lines,nline,table_footer)

   854   set_line(lines,nline,footer) 

   855 

   856 ; Now write to an HTML file.

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

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

   859     asciiwrite(output_html,lines(idx))

   860   else

   861    print ("error?")

   862   end if

   863   delete (idx)

   864 

   865 ;***************************************************************************

   866 ; output plots

   867 ;***************************************************************************

   868   output_dir = model_name+"/ameriflux"

   869 

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

   871 ;***************************************************************************

   872 end
clamp