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/11.ameriflux.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 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

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

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

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

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

     7 begin

     8 ; add line to ascci/html file

     9     

    10   nnewlines = dimsizes(newlines)

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

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

    13     return

    14   end if 

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

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

    17   nline = nline + nnewlines

    18   return 

    19 end

    20 ;*************************************************************

    21 begin

    22 

    23   plot_type     = "ps"

    24   plot_type_new = "png"

    25 

    26 ;------------------------------------------------------

    27 ; edit table.html of current model for movel1_vs_model2

    28 

    29  if (isvar("compare")) then

    30     html_name2 = compare+"/table.html"  

    31     html_new2  = html_name2 +".new"

    32  end if

    33 

    34 ;------------------------------------------------------

    35 ; edit table.html for current model

    36 

    37  html_name = model_name+"/table.html"  

    38  html_new  = html_name +".new"

    39 

    40 ;------------------------------------------------------

    41 

    42   nmonth = 12

    43 

    44 ; for nee, gpp, and ar

    45 ; observed unit is gC/m2/day

    46 ; model    unit is gC/m2/s

    47 ; to change to observed unit,

    48 

    49   factor_flux = 86400.

    50   

    51 ; for incident solar radiation,

    52 ; observed Rg_f unit is MJ/m2/day

    53 ; model (FSDS)  unit is  W/m2

    54 ; to change to model unit,

    55 

    56   factor_rad = 1.e6/86400.

    57 

    58 ;************************************************

    59 ; observed data info

    60 ;************************************************

    61 

    62  station = (/"ARM_Oklahoma" \

    63             ,"ARM_Oklahoma_burn" \

    64             ,"ARM_Oklahoma_control" \

    65             ,"Atqasuk" \

    66             ,"Audubon" \

    67             ,"AustinCary" \

    68             ,"Bartlett" \

    69             ,"Bondville" \

    70             ,"Brookings" \

    71             ,"Donaldson" \

    72             ,"Duke_Forest_Hardwoods" \

    73             ,"Duke_Forest_Open_Field" \

    74             ,"Duke_Forest_Pine" \

    75             ,"Fermi_Ag" \

    76             ,"Fermi_Prairie" \

    77             ,"Flagstaff_Managed" \

    78             ,"Flagstaff_Unmanaged" \

    79             ,"Flagstaff_Wildfire" \

    80             ,"FortPeck" \

    81             ,"FreemanRanch_mesquite" \

    82             ,"Goodwin_Creek" \

    83             ,"HarvardForest" \

    84             ,"HarvardForestHemlock" \

    85             ,"HowlandForestMain" \

    86             ,"HowlandForestWest" \

    87             ,"Ivotuk" \

    88             ,"KendallGrasslands" \

    89             ,"KennedySpaceCenterPine" \

    90             ,"KennedySpaceCenterScrub" \

    91             ,"LittleProspect" \

    92             ,"LostCreek" \

    93             ,"Mead-irrigated" \

    94             ,"Mead-irrigated-rotation" \

    95             ,"Mead-rainfed" \

    96             ,"Metolius_2nd_YoungPonderosaPine" \

    97             ,"MetoliusEyerly" \

    98             ,"MetoliusIntermediatePine" \

    99             ,"MetoliusOldPonderosaPine" \

   100             ,"MissouriOzark" \

   101             ,"Mize" \

   102             ,"MorganMonroe" \

   103             ,"NiwotRidge" \

   104             ,"NorthCarolina_cc" \

   105             ,"NorthCarolina_lp" \

   106             ,"ParkFalls" \

   107             ,"Rayonier" \

   108             ,"SantaRita" \

   109             ,"SkyOaks_Old" \

   110             ,"SkyOaks_PostFire" \

   111             ,"SkyOaks_Young" \

   112             ,"SylvaniaWilderness" \

   113             ,"Toledo" \

   114             ,"Tonzi" \

   115             ,"UCI_1850" \

   116             ,"UCI_1930" \

   117             ,"UCI_1964" \

   118             ,"UCI_1964wet" \

   119             ,"UCI_1981" \

   120             ,"UCI_1989" \

   121             ,"UCI_1998" \

   122             ,"UMBS" \

   123             ,"Vaira" \

   124             ,"WalkerBranch" \

   125             ,"WillowCreek" \

   126             ,"WindRiver" \

   127             ,"Wisconsin_ihw" \

   128             ,"Wisconsin_irp" \

   129             ,"Wisconsin_mrp" \

   130             ,"Wisconsin_myjp" \

   131             ,"Wisconsin_pb" \

   132             ,"Wisconsin_rpcc" \

   133             ,"Wisconsin_yhw" \

   134             ,"Wisconsin_yjp" \

   135             ,"Wisconsin_yrp" \

   136             /)

   137 

   138  field   = (/"NEE Flux" \

   139             ,"Shortwave Incoming" \

   140             ,"Latent Heat" \

   141             ,"Sensible Heat" \

   142             ,"GPP Flux" \

   143             ,"Respiration" \

   144             /)

   145 

   146  field_unit = (/"gC/m2/day" \

   147                ,"W/m2" \

   148                ,"W/m2" \

   149                ,"W/m2" \

   150                ,"gC/m2/day" \

   151                ,"gC/m2/day" \

   152                /)

   153 

   154  nstation = dimsizes(station)

   155  nfield   = dimsizes(field)

   156 

   157 ;========================================================================

   158 ; get observed info: number of year, first/last year, lat, lon

   159 ; and annual data

   160 

   161  dir_root = diro + "ameriflux/"

   162 

   163  year_station = new ((/nstation/),integer) ; number of year

   164  year_ob      = new ((/nstation/),string)  ; observed year

   165  year_ob_i    = new ((/nstation/),integer) ; first year

   166  year_ob_f    = new ((/nstation/),integer) ; last year

   167  lat_ob       = new ((/nstation/),float)   ; latitude

   168  lon_ob       = new ((/nstation/),float)   ; longitude 

   169 

   170  data_ob_ann  = new ((/nfield, nmonth, nstation/),float)

   171 

   172  do n = 0, nstation-1

   173 

   174     dir_f = dir_root + station(n)+"/"

   175     fil_f = "timeseries_L4_m.nc"

   176     fo   = addfile (dir_f+fil_f,"r")

   177  

   178     lat_ob(n) = fo->lat

   179     lon_ob(n) = fo->lon

   180 

   181     year      = fo->year

   182 

   183     year_station(n) = dimsizes(year)

   184     year_ob_i(n)    = year(0)

   185     year_ob_f(n)    = year(year_station(n)-1)

   186     year_ob(n)      = year_ob_i(n) + "-" + year_ob_f(n)

   187  

   188     delete (year)

   189  

   190     data = fo->NEE_or_fMDS

   191     data_ob_ann(0,:,n) = dim_avg(data(month|:,year|:))

   192 

   193     data = fo->Rg_f

   194     data_ob_ann(1,:,n) = dim_avg(data(month|:,year|:)) * factor_rad

   195 

   196     data = fo->LE_f

   197     data_ob_ann(2,:,n) = dim_avg(data(month|:,year|:))

   198 

   199     data = fo->H_f

   200     data_ob_ann(3,:,n) = dim_avg(data(month|:,year|:))

   201 

   202     data = fo->GPP_or_MDS

   203     data_ob_ann(4,:,n) = dim_avg(data(month|:,year|:))

   204 

   205     data = fo->Reco_or

   206     data_ob_ann(5,:,n) = dim_avg(data(month|:,year|:))

   207 

   208     delete (data) 

   209     delete (fo)   

   210  end do

   211 

   212 ;--------------------------------------------------------------

   213 ; find (# of year observed) >=4 and year_ob_i <= 2001

   214 

   215  i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)

   216  

   217  station_long      = station(i_long_ob)

   218  lat_ob_long       = lat_ob(i_long_ob)

   219  lon_ob_long       = lat_ob(i_long_ob)

   220  year_ob_long      = year_ob(i_long_ob)

   221  year_ob_i_long    = year_ob_i(i_long_ob)

   222  year_ob_f_long    = year_ob_f(i_long_ob)

   223  year_station_long = year_station(i_long_ob)

   224 

   225  nstation_long     = dimsizes(station_long)

   226 

   227 ;=========================================================

   228 ;   get model data at observed lat-lon

   229 

   230     fm   = addfile (dirm+film8,"r")

   231 

   232     xm   = fm->lon

   233     ym   = fm->lat 

   234     date = fm->date

   235 

   236     date_dim = dimsizes(date)

   237     nyear  = date_dim(0)

   238 

   239     data_mod      = new ((/nfield,nyear,nmonth,nstation/),float)

   240     data_mod_ann  = new ((/nfield,nmonth,nstation/),float)

   241     data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)

   242 

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

   244 ;   Model_units [=] gC/m2/s

   245 ;   12. = molecular weight of C

   246 ;   u mol = 1e-6 mol

   247 

   248     factor = 1.e6 /12.

   249 

   250 ;------------------------------------------------------------

   251 ;   interpolate model data into observed station

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

   253 

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

   255     ym(0) = -90.

   256 ;------------------------------------------------------------  

   257 

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

   259 

   260   data = fm->NEE

   261   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   262   data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   263   data_mod(0,:,:,:) = yy(:,:,:) * factor_flux

   264 

   265 ;;data  = fm->NETRAD

   266 ; data  = fm->FSA

   267 ; data1 = fm->FIRA

   268 ; data  = data - data1

   269 ; delete (data1)

   270 

   271   data  = fm->FSDS

   272   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   273   data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   274   data_mod(1,:,:,:) = yy(:,:,:)  

   275 

   276 

   277 ; data  = fm->LATENT

   278   data  = fm->FCEV

   279   data1 = fm->FCTR

   280   data2 = fm->FGEV

   281   data  = data + data1 + data2

   282   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   283   data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))

   284   data_mod(2,:,:,:) = yy(:,:,:) 

   285   delete (data1)

   286   delete (data2)

   287  

   288 ; data = fm->SENSIBLE

   289   data  = fm->FSH

   290   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   291   data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   292   data_mod(3,:,:,:) = yy(:,:,:)  

   293 

   294 else

   295 

   296   data = fm->NEE

   297   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   298   data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   299   data_mod(0,:,:,:) = yy(:,:,:) * factor_flux

   300 

   301 ; data = fm->NETRAD

   302   data = fm->FSDS

   303   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   304   data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   305   data_mod(1,:,:,:) = yy(:,:,:) 

   306 

   307   data = fm->LATENT

   308   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   309   data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   310   data_mod(2,:,:,:) = yy(:,:,:) 

   311 

   312 ; data = fm->SENSIBLE

   313   data = fm->FSH

   314   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   315   data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   316   data_mod(3,:,:,:) = yy(:,:,:) 

   317 

   318 end if

   319 

   320   data = fm->GPP

   321   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   322   data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   323   data_mod(4,:,:,:) = yy(:,:,:) * factor_flux 

   324 

   325   data = fm->ER

   326   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   327   data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux

   328   data_mod(5,:,:,:) = yy(:,:,:) * factor_flux

   329 

   330   data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)

   331 

   332   delete (data_mod)

   333   delete (fm) 

   334   delete (data)

   335   delete (yy)

   336 

   337   asciiwrite(station(0)+"_ob.txt", data_ob_ann(:,:,0))

   338   asciiwrite(station(0)+"_"+model_name+".txt", data_mod_ann(:,:,0))

   339 

   340   asciiwrite(station(7)+"_ob.txt", data_ob_ann(:,:,7))

   341   asciiwrite(station(7)+"_"+model_name+".txt", data_mod_ann(:,:,7))

   342 

   343   asciiwrite(station(21)+"_ob.txt", data_ob_ann(:,:,21))

   344   asciiwrite(station(21)+"_"+model_name+".txt", data_mod_ann(:,:,21))

   345 

   346   asciiwrite(station(44)+"_ob.txt", data_ob_ann(:,:,44))

   347   asciiwrite(station(44)+"_"+model_name+".txt", data_mod_ann(:,:,44))

   348 

   349   asciiwrite(station(50)+"_ob.txt", data_ob_ann(:,:,50))

   350   asciiwrite(station(50)+"_"+model_name+".txt", data_mod_ann(:,:,50))

   351 

   352   asciiwrite(station(54)+"_ob.txt", data_ob_ann(:,:,54))

   353   asciiwrite(station(54)+"_"+model_name+".txt", data_mod_ann(:,:,54))

   354 

   355   asciiwrite(station(62)+"_ob.txt", data_ob_ann(:,:,62))

   356   asciiwrite(station(62)+"_"+model_name+".txt", data_mod_ann(:,:,62))

   357 exit

   358 

   359 ;*******************************************************************

   360 ; for station line plot

   361 ;*******************************************************************

   362 

   363 ; for x-axis in xyplot

   364   mon = ispan(1,12,1)

   365   mon@long_name = "month"

   366 

   367   res                   = True               ; plot mods desired

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

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

   370 

   371   res@tmXBFormat  = "f"                      ; not to add trailing zeros

   372 

   373 ;-------------------------------------------------------------------------

   374 ; Add a boxed legend using the more simple method

   375 

   376   res@pmLegendDisplayMode    = "Always"

   377 ; res@pmLegendWidthF         = 0.1

   378   res@pmLegendWidthF         = 0.08

   379   res@pmLegendHeightF        = 0.06

   380 ; res@pmLegendOrthogonalPosF = -1.17

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

   382   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   383 

   384 ; res@pmLegendParallelPosF   =  0.18

   385   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   386 

   387 ; res@lgPerimOn             = False

   388   res@lgLabelFontHeightF     = 0.015

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

   390 ;-------------------------------------------------------------------

   391 ; for panel plot

   392   res@gsnFrame     = False                   ; Do not draw plot 

   393   res@gsnDraw      = False                   ; Do not advance frame

   394 

   395   pres                            = True     ; panel plot mods desired

   396   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   397                                              ; indiv. plots in panel

   398   pres@gsnMaximize                = True     ; fill the page

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

   400 

   401 ;==============================================================

   402 ; get ob data at each site with long observation

   403 

   404   do n = 0,nstation_long-1

   405 

   406 ;##################################################################

   407 ; hardwired: model    up to year 2004

   408 ;            observed up to year 2006

   409 

   410     year_setback = 0

   411 

   412     nyear = year_station_long(n)

   413 

   414     if (year_ob_f_long(n).eq. 2006) then

   415        year_setback = 2006 -2004

   416     end if

   417     if (year_ob_f_long(n).eq. 2005) then

   418        year_setback = 2005 -2004

   419     end if

   420 ;##################################################################

   421     

   422     ntime = (nyear - year_setback) * nmonth

   423 

   424     data_ob   = new ((/nfield, nyear, nmonth/),float)

   425 

   426     dir_f = dir_root + station_long(n)+"/"

   427     fil_f = "timeseries_L4_m.nc"

   428     fo    = addfile (dir_f+fil_f,"r")

   429 

   430     data_ob(0,:,:) = fo->NEE_or_fMDS

   431     data_ob(1,:,:) = fo->Rg_f

   432     data_ob(2,:,:) = fo->LE_f

   433     data_ob(3,:,:) = fo->H_f

   434     data_ob(4,:,:) = fo->GPP_or_MDS

   435     data_ob(5,:,:) = fo->Reco_or

   436 

   437     data_ob(1,:,:) = data_ob(1,:,:) * factor_rad

   438 

   439     delete (fo)

   440 

   441     timeI = new((/ntime/),integer)

   442     timeF = new((/ntime/),float)

   443     timeI = ispan(1,ntime,1)

   444     timeF = year_ob_i_long(n) + (timeI-1)/12.

   445     timeF@long_name = "year" 

   446 

   447     plot_data = new((/2,ntime/),float)

   448  

   449 ;----------------------------

   450 ; for model_vs_ob

   451 

   452     plot_name = station_long(n)+"_tseries_vs_ob"

   453     title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"    

   454     res@tiMainString = title

   455 

   456     wks = gsn_open_wks (plot_type,plot_name)

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

   458 

   459     i_year_mod_i = year_ob_i_long(n) - 1990 

   460     i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback

   461 

   462     i_year_ob_f =  nyear - year_setback - 1 

   463 

   464 ;   print (nyear)

   465 ;   print (i_year_ob_f)

   466 ;   print (i_year_mod_i)

   467 ;   print (i_year_mod_f)

   468 

   469     do i = 0,nfield-1                           

   470        plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))

   471        plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))

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

   473        plot(i)=gsn_csm_xy(wks,timeF,plot_data,res)            ; create plot 

   474     end do

   475    

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

   477 

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

   479            "rm "+plot_name+"."+plot_type)

   480 

   481     clear (wks)  

   482     delete (plot)

   483 

   484     delete (data_ob)

   485     delete (timeI)    

   486     delete (timeF)

   487     delete (plot_data)

   488 

   489  end do

   490 

   491 ;###################################################################

   492 ; for the following tables,

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

   494 

   495   isort = dim_pqsort(lat_ob_long,-1)

   496 

   497   station_sort = station_long(isort)

   498   year_ob_sort = year_ob_long(isort)

   499   lat_ob_sort  = lat_ob_long(isort)

   500   lon_ob_sort  = lon_ob_long(isort)

   501  

   502 ; print(isort)

   503 ; print(lat_ob_sort)

   504 

   505 ;*******************************************************************

   506 ; html table of site: observed

   507 ;*******************************************************************

   508   output_html = "tseries_vs_ob.html"

   509 

   510   header = (/"<HTML>" \

   511             ,"<HEAD>" \

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

   513             ,"</HEAD>" \

   514             ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \

   515             /) 

   516   footer = "</HTML>"

   517 

   518   table_header = (/ \

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

   520        ,"<tr>" \

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

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

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

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

   525        ,"</tr>" \

   526        /)

   527   table_footer = "</table>"

   528   row_header = "<tr>"

   529   row_footer = "</tr>"

   530 

   531   lines = new(50000,string)

   532   nline = 0

   533 

   534   set_line(lines,nline,header)

   535   set_line(lines,nline,table_header)

   536 ;-----------------------------------------------

   537 ; row of table

   538   

   539   do n = 0,nstation_long-1

   540 

   541      set_line(lines,nline,row_header)

   542 

   543      txt0 = station_sort(n)

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

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

   546      txt3 = year_ob_sort(n)

   547 

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

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

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

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

   552 

   553      set_line(lines,nline,row_footer)

   554   end do

   555 ;-----------------------------------------------

   556   set_line(lines,nline,table_footer)

   557   set_line(lines,nline,footer) 

   558 

   559 ; Now write to an HTML file.

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

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

   562     asciiwrite(output_html,lines(idx))

   563   else

   564    print ("error?")

   565   end if

   566   delete (idx)

   567 

   568   delete (isort)

   569   delete (station_sort)

   570   delete (year_ob_sort)

   571   delete (lat_ob_sort)

   572   delete (lon_ob_sort)

   573 

   574 ;************************************************************

   575 ; compute annual cycle correlation coef and M score

   576 ;************************************************************

   577 

   578  score_max = 1.

   579 

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

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

   582 

   583  do n=0,nstation-1

   584  do m=0,nfield-1   

   585     ccr(n,m) = esccr(data_ob_ann(m,:,n),data_mod_ann(m,:,n),0)

   586     bias = sum(abs(data_mod_ann(m,:,n)-data_ob_ann(m,:,n))/(abs(data_mod_ann(m,:,n))+abs(data_ob_ann(m,:,n))))

   587     M_score(n,m) = (1. -(bias/nmonth)) * score_max

   588  end do

   589  end do

   590 

   591  M_nee = avg(M_score(:,0))

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

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

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

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

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

   597  M_all = M_nee+ M_rad +M_lh + M_sh + M_gpp + M_er

   598 

   599  M_ameriflux_nee = sprintf("%.2f", M_nee)

   600  M_ameriflux_rad = sprintf("%.2f", M_rad)

   601  M_ameriflux_lh  = sprintf("%.2f", M_lh )

   602  M_ameriflux_sh  = sprintf("%.2f", M_sh )

   603  M_ameriflux_gpp = sprintf("%.2f", M_gpp)

   604  M_ameriflux_er  = sprintf("%.2f", M_er )

   605  M_ameriflux_all = sprintf("%.2f", M_all)

   606 

   607 ;*******************************************************************

   608 ; for station line plot

   609 ;*******************************************************************

   610 

   611 ; for x-axis in xyplot

   612   mon = ispan(1,12,1)

   613   mon@long_name = "month"

   614 

   615 ;-------------------------------------------------------------------

   616 

   617   plot_data   = new((/2,nmonth/),float)

   618   plot_data!0 = "case"

   619   plot_data!1 = "month"

   620 

   621   do n = 0,nstation-1

   622 ;----------------------------

   623 ; for observed

   624 

   625     plot_name = station(n)+"_ob"    

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

   627     res@tiMainString = title

   628 

   629     wks = gsn_open_wks (plot_type,plot_name)

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

   631 

   632     do i = 0,nfield-1                           

   633        plot_data(0,:) = (/data_ob_ann(i,:,n)/)

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

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

   636     end do

   637 

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

   639 

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

   641            "rm "+plot_name+"."+plot_type)

   642 

   643     clear (wks)  

   644     delete (plot)

   645 ;----------------------------

   646 ; for model_vs_ob

   647 

   648     plot_name = station(n)+"_model_vs_ob"

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

   650     res@tiMainString = title

   651 

   652     wks = gsn_open_wks (plot_type,plot_name)

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

   654 

   655     do i = 0,nfield-1                           

   656        plot_data(0,:) = (/data_ob_ann(i,:,n)/)

   657        plot_data(1,:) = (/data_mod_ann(i,:,n)/)

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

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

   660     end do

   661    

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

   663 

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

   665            "rm "+plot_name+"."+plot_type)

   666 

   667     clear (wks)  

   668     delete (plot)

   669  end do

   670 

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

   672 ; for the following tables,

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

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

   675 

   676   isort = dim_pqsort(lat_ob,-1)

   677 

   678   station_sort = station(isort)

   679   year_ob_sort = year_ob(isort)

   680   lat_ob_sort  = lat_ob(isort)

   681   lon_ob_sort  = lon_ob(isort)

   682   M_score_sort = M_score(isort,:)

   683  

   684 ; print(isort)

   685 ; print(lat_ob_sort)

   686 ;###################################################################

   687 ;*******************************************************************

   688 ; html table of site: observed

   689 ;*******************************************************************

   690   output_html = "line_ob.html"

   691 

   692   header = (/"<HTML>" \

   693             ,"<HEAD>" \

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

   695             ,"</HEAD>" \

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

   697             /) 

   698   footer = "</HTML>"

   699 

   700   table_header = (/ \

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

   702        ,"<tr>" \

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

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

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

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

   707        ,"</tr>" \

   708        /)

   709   table_footer = "</table>"

   710   row_header = "<tr>"

   711   row_footer = "</tr>"

   712 

   713   lines = new(50000,string)

   714   nline = 0

   715 

   716   set_line(lines,nline,header)

   717   set_line(lines,nline,table_header)

   718 ;-----------------------------------------------

   719 ; row of table

   720   

   721   do n = 0,nstation-1

   722      set_line(lines,nline,row_header)

   723 

   724      txt0 = station_sort(n)

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

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

   727      txt3 = year_ob_sort(n)

   728 

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

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

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

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

   733 

   734      set_line(lines,nline,row_footer)

   735   end do

   736 ;-----------------------------------------------

   737   set_line(lines,nline,table_footer)

   738   set_line(lines,nline,footer) 

   739 

   740 ; Now write to an HTML file.

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

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

   743     asciiwrite(output_html,lines(idx))

   744   else

   745    print ("error?")

   746   end if

   747   delete (idx)

   748 

   749 ;*******************************************************************

   750 ; score and line table : model vs observed

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

   752   output_html = "score+line_vs_ob.html"

   753 

   754   header = (/"<HTML>" \

   755             ,"<HEAD>" \

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

   757             ,"</HEAD>" \

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

   759             /) 

   760   footer = "</HTML>"

   761 

   762   delete (table_header)

   763   table_header = (/ \

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

   765        ,"<tr>" \

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

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

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

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

   770        ,"   <th bgcolor=DDDDDD >NEE Flux</th>" \

   771        ,"   <th bgcolor=DDDDDD >Shortwave <br> Incoming</th>" \

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

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

   774        ,"   <th bgcolor=DDDDDD >GPP Flux</th>" \

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

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

   777        ,"</tr>" \

   778        /)

   779   table_footer = "</table>"

   780   row_header = "<tr>"

   781   row_footer = "</tr>"

   782 

   783   lines = new(50000,string)

   784   nline = 0

   785 

   786   set_line(lines,nline,header)

   787   set_line(lines,nline,table_header)

   788 ;-----------------------------------------------

   789 ; row of table

   790   

   791   do n = 0,nstation-1

   792      set_line(lines,nline,row_header)

   793 

   794      txt0  = station_sort(n)

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

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

   797      txt3  = year_ob_sort(n)

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

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

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

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

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

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

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

   805 

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

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

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

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

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

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

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

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

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

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

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

   817 

   818      set_line(lines,nline,row_footer)

   819   end do

   820 

   821 ; last row, summary

   822   set_line(lines,nline,row_header)

   823 

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

   825   txt1  = "-"

   826   txt2  = "-"

   827   txt3  = "-"

   828   txt4  = M_ameriflux_nee

   829   txt5  = M_ameriflux_rad

   830   txt6  = M_ameriflux_lh

   831   txt7  = M_ameriflux_sh

   832   txt8  = M_ameriflux_gpp

   833   txt9  = M_ameriflux_er

   834   txt10 = M_ameriflux_all

   835 

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

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

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

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

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

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

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

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

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

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

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

   847 

   848   set_line(lines,nline,row_footer)

   849 ;-----------------------------------------------

   850   set_line(lines,nline,table_footer)

   851   set_line(lines,nline,footer) 

   852 

   853 ; Now write to an HTML file.

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

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

   856     asciiwrite(output_html,lines(idx))

   857   else

   858    print ("error?")

   859   end if

   860   delete (idx)

   861 

   862 ;**************************************************************************************

   863 ; update score

   864 ;**************************************************************************************

   865  

   866   if (isvar("compare")) then

   867      system("sed 1,/M_ameriflux_nee/s//"+M_ameriflux_nee+"/ "+html_name2+" > "+html_new2+";"+ \

   868             "mv -f "+html_new2+" "+html_name2+";"+ \

   869             "sed 1,/M_ameriflux_rad/s//"+M_ameriflux_rad+"/ "+html_name2+" > "+html_new2+";"+ \

   870             "mv -f "+html_new2+" "+html_name2+";"+ \

   871             "sed 1,/M_ameriflux_lh/s//"+M_ameriflux_lh+"/ "+html_name2+" > "+html_new2+";"+ \

   872             "mv -f "+html_new2+" "+html_name2+";"+ \

   873             "sed 1,/M_ameriflux_sh/s//"+M_ameriflux_sh+"/ "+html_name2+" > "+html_new2+";"+ \

   874             "mv -f "+html_new2+" "+html_name2+";"+ \

   875             "sed 1,/M_ameriflux_gpp/s//"+M_ameriflux_gpp+"/ "+html_name2+" > "+html_new2+";"+ \

   876             "mv -f "+html_new2+" "+html_name2+";"+ \

   877             "sed 1,/M_ameriflux_er/s//"+M_ameriflux_er+"/ "+html_name2+" > "+html_new2+";"+ \

   878             "mv -f "+html_new2+" "+html_name2)

   879   end if

   880 

   881   system("sed s#M_ameriflux_nee#"+M_ameriflux_nee+"# "+html_name+" > "+html_new+";"+ \

   882          "mv -f "+html_new+" "+html_name+";"+ \

   883          "sed s#M_ameriflux_rad#"+M_ameriflux_rad+"# "+html_name+" > "+html_new+";"+ \

   884          "mv -f "+html_new+" "+html_name+";"+ \

   885          "sed s#M_ameriflux_lh#"+M_ameriflux_lh+"# "+html_name+" > "+html_new+";"+ \

   886          "mv -f "+html_new+" "+html_name+";"+ \

   887          "sed s#M_ameriflux_sh#"+M_ameriflux_sh+"# "+html_name+" > "+html_new+";"+ \

   888          "mv -f "+html_new+" "+html_name+";"+ \

   889          "sed s#M_ameriflux_gpp#"+M_ameriflux_gpp+"# "+html_name+" > "+html_new+";"+ \

   890          "mv -f "+html_new+" "+html_name+";"+ \

   891          "sed s#M_ameriflux_er#"+M_ameriflux_er+"# "+html_name+" > "+html_new+";"+ \

   892          "mv -f "+html_new+" "+html_name) 

   893 

   894 ;***************************************************************************

   895 ; add total score and write to file

   896 ;***************************************************************************

   897   M_total = M_ameriflux_all

   898 

   899   asciiwrite("M_save.ameriflux", M_total)

   900 

   901 ;***************************************************************************

   902 ; output plot and html

   903 ;***************************************************************************

   904   output_dir = model_name+"/ameriflux"

   905 

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

   907 ;***************************************************************************

   908 

   909 end

   910
clamp