clamp: all/06.fluxnet.ncl@4be95183fbcd

     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

    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

    23   plot_type     = "ps"

    24   plot_type_new = "png"

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

    27 ; edit table.html of current model for movel1_vs_model2

    29  if (isvar("compare")) then

    30     html_name2 = compare+"/table.html"

    31     html_new2  = html_name2 +".new"

    32  end if

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

    35 ; edit table.html for current model

    37  html_name = model_name+"/table.html"

    38  html_new  = html_name +".new"

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

    41 ; read model data

    43   fm    = addfile(dirm+film2,"r")

    45   xm    = fm->lon

    46   ym    = fm->lat

    48   nlat = dimsizes(ym)

    49   nlon = dimsizes(xm)

    51 ; for 4 fields, 12-monthly

    52   nmon      = 12

    53   nfield    = 4

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

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

    58 ; Model_units [=] gC/m2/s

    59 ; 12. = molecular weight of C

    60 ; u mol = 1e-6 mol

    61   factor = 1e6 /12.

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

    65   data = fm->NEE

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

    67   delete (data)

    69 ; data  = fm->LATENT

    70   data1 = fm->FCEV

    71   data2 = fm->FCTR

    72   data3 = fm->FGEV

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

    74   delete (data1)

    75   delete (data2)

    76   delete (data3)

    78 ; data = fm->SENSIBLE

    79   data  = fm->FSH

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

    81   delete (data)

    83 ; data  = fm->NETRAD

    84   data1 = fm->FSA

    85   data2 = fm->FIRA

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

    87   delete (data1)

    88   delete (data2)

    90 else

    92   data = fm->NEE

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

    94   delete (data)

    96   data = fm->NETRAD

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

    98   delete (data)

   100   data = fm->LATENT

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

   102   delete (data)

   104   data = fm->FSH

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

   106   delete (data)

   107 end if

   109  delete (fm)

   111 ;************************************************

   112 ; read data: observed

   113 ;************************************************

   115  station = (/"BOREAS_NSA_OBS" \

   116             ,"CastelPorziano" \

   117             ,"Hyytiala" \

   118             ,"Kaamanen" \

   119             ,"LBA_Tapajos_KM67" \

   120             ,"Lethbridge" \

   121             ,"Tharandt" \

   122             ,"Vielsalm" \

   123             /)

   125  year_ob = (/"1994-2004" \

   126             ,"1997-2003" \

   127             ,"1996-2003" \

   128             ,"2000-2003" \

   129             ,"2002-2005" \

   130             ,"1999-2004" \

   131             ,"1996-2003" \

   132             ,"1998-2003" \

   133             /)

   135  field   = (/"NEE" \

   136             ,"Net Radiation" \

   137             ,"Latent Heat" \

   138             ,"Sensible Heat" \

   139             /)

   141  nstation  = dimsizes(station)

   142  nmon      = 12

   143  nfield    = dimsizes(field)

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

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

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

   149  diri_root  = diro + "fluxnet/"

   151  do n = 0,nstation-1

   152     diri = diri_root + station(n)+"/"

   153     fili = station(n)+"_"+year_ob(n)+"_monthly.nc"

   154     g     = addfile (diri+fili,"r")

   156     lon_ob(n) = g->lon

   157     lat_ob(n) = g->lat

   159     data      = g->CO2_FLUX

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

   161     delete (data)

   163     data      = g->RAD_FLUX

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

   165     delete (data)

   167     data      = g->LH_FLUX

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

   169     delete (data)

   171     data      = g->SH_FLUX

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

   173     delete (data)

   175     delete (g)

   176  end do

   178 ;************************************************************

   179 ; interpolate model data into observed station

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

   181 ;************************************************************

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

   184   ym(0) = -90.

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

   188   delete (data_mod0)

   189   yy!0 = "field"

   190   data_mod = yy(pts|:,field|:,time|:)

   192 ;************************************************************

   193 ; compute correlation coef and M score

   194 ;************************************************************

   196  score_max = 5.

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

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

   201  do n=0,nstation-1

   202  do m=0,nfield-1

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

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

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

   206  end do

   207  end do

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

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

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

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

   213  M_all = M_nee+ M_rad +M_lh + M_sh

   215  M_fluxnet_nee = sprintf("%.2f", M_nee)

   216  M_fluxnet_rad = sprintf("%.2f", M_rad)

   217  M_fluxnet_lh  = sprintf("%.2f", M_lh )

   218  M_fluxnet_sh  = sprintf("%.2f", M_sh )

   219  M_fluxnet_all = sprintf("%.2f", M_all)

   221 ;*******************************************************************

   222 ; for station line plot

   223 ;*******************************************************************

   225 ; for x-axis in xyplot

   226   mon = ispan(1,12,1)

   227   mon@long_name = "month"

   229   res                   = True               ; plot mods desired

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

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

   232 ;-------------------------------------------------------------------------

   233 ; Add a boxed legend using the more simple method

   235   res@pmLegendDisplayMode    = "Always"

   236 ; res@pmLegendWidthF         = 0.1

   237   res@pmLegendWidthF         = 0.08

   238   res@pmLegendHeightF        = 0.06

   239 ; res@pmLegendOrthogonalPosF = -1.17

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

   241   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   243 ; res@pmLegendParallelPosF   =  0.18

   244   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   246 ; res@lgPerimOn             = False

   247   res@lgLabelFontHeightF     = 0.015

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

   249 ;-------------------------------------------------------------------

   250 ; for panel plot

   251   res@gsnFrame     = False                   ; Do not draw plot

   252   res@gsnDraw      = False                   ; Do not advance frame

   254   pres                            = True     ; panel plot mods desired

   255   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   256                                              ; indiv. plots in panel

   257   pres@gsnMaximize                = True     ; fill the page

   258 ;-------------------------------------------------------------------

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

   261   plot_data!0 = "case"

   262   plot_data!1 = "month"

   264 ; change longitude from 0-360 to -180-180

   265   lon_ob = where(lon_ob .gt. 180.,lon_ob-360., lon_ob)

   267   do n = 0,nstation-1

   268 ;----------------------------

   269 ; for observed

   271     plot_name = station(n)+"_ob"

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

   273     res@tiMainString = title

   275     wks = gsn_open_wks (plot_type,plot_name)

   276     plot=new(4,graphic)                        ; create graphic array

   278     plot_data(0,:) = (/data_ob (n,0,:)/)

   279     plot_data@long_name = field(0)

   280     plot(0)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 1

   282     plot_data(0,:) = (/data_ob (n,1,:)/)

   283     plot_data@long_name = field(1)

   284     plot(1)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 2

   286     plot_data(0,:) = (/data_ob (n,2,:)/)

   287     plot_data@long_name = field(2)

   288     plot(2)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 3

   290     plot_data(0,:) = (/data_ob (n,3,:)/)

   291     plot_data@long_name = field(3)

   292     plot(3)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 4

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

   296     delete (wks)

   297     delete (plot)

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

   300            "rm "+plot_name+"."+plot_type)

   302 ;----------------------------

   303 ; for model_vs_ob

   305     plot_name = station(n)+"_model_vs_ob"

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

   307     res@tiMainString = title

   309     wks = gsn_open_wks (plot_type,plot_name)

   310     plot=new(4,graphic)                        ; create graphic array

   312     plot_data(0,:) = (/data_ob (n,0,:)/)

   313     plot_data(1,:) = (/data_mod(n,0,:)/)

   314     plot_data@long_name = field(0)

   315     plot(0)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 1

   317     plot_data(0,:) = (/data_ob (n,1,:)/)

   318     plot_data(1,:) = (/data_mod(n,1,:)/)

   319     plot_data@long_name = field(1)

   320     plot(1)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 2

   322     plot_data(0,:) = (/data_ob (n,2,:)/)

   323     plot_data(1,:) = (/data_mod(n,2,:)/)

   324     plot_data@long_name = field(2)

   325     plot(2)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 3

   327     plot_data(0,:) = (/data_ob (n,3,:)/)

   328     plot_data(1,:) = (/data_mod(n,3,:)/)

   329     plot_data@long_name = field(3)

   330     plot(3)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 4

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

   334     delete (wks)

   335     delete (plot)

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

   338            "rm "+plot_name+"."+plot_type)

   339  end do

   341 ;*******************************************************************

   342 ; html table of site: observed

   343 ;*******************************************************************

   344   output_html = "line_ob.html"

   346   header = (/"<HTML>" \

   347             ,"<HEAD>" \

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

   349             ,"</HEAD>" \

   350             ,"<H1>Fluxnet at Site: Observation</H1>" \

   351             /)

   352   footer = "</HTML>"

   354   table_header = (/ \

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

   356        ,"<tr>" \

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

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

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

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

   361        ,"</tr>" \

   362        /)

   363   table_footer = "</table>"

   364   row_header = "<tr>"

   365   row_footer = "</tr>"

   367   lines = new(50000,string)

   368   nline = 0

   370   set_line(lines,nline,header)

   371   set_line(lines,nline,table_header)

   372 ;-----------------------------------------------

   373 ; row of table

   375   do n = 0,nstation-1

   376      set_line(lines,nline,row_header)

   378      txt0 = station(n)

   379      txt1 = sprintf("%5.2f", lat_ob(n))

   380      txt2 = sprintf("%5.2f", lon_ob(n))

   381      txt3 = year_ob(n)

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

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

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

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

   388      set_line(lines,nline,row_footer)

   389   end do

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

   391   set_line(lines,nline,table_footer)

   392   set_line(lines,nline,footer)

   394 ; Now write to an HTML file.

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

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

   397     asciiwrite(output_html,lines(idx))

   398   else

   399    print ("error?")

   400   end if

   401   delete (idx)

   403 ;*******************************************************************

   404 ; score and line table : model vs observed

   405 ;*******************************************************************

   406   output_html = "score+line_vs_ob.html"

   408   header = (/"<HTML>" \

   409             ,"<HEAD>" \

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

   411             ,"</HEAD>" \

   412             ,"<H1>Fluxnet at Site: Model "+model_name+"</H1>" \

   413             /)

   414   footer = "</HTML>"

   416   delete (table_header)

   417   table_header = (/ \

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

   419        ,"<tr>" \

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

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

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

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

   424        ,"   <th bgcolor=DDDDDD >NEE</th>" \

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

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

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

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

   429        ,"</tr>" \

   430        /)

   431   table_footer = "</table>"

   432   row_header = "<tr>"

   433   row_footer = "</tr>"

   435   lines = new(50000,string)

   436   nline = 0

   438   set_line(lines,nline,header)

   439   set_line(lines,nline,table_header)

   440 ;-----------------------------------------------

   441 ; row of table

   443   do n = 0,nstation-1

   444      set_line(lines,nline,row_header)

   446      txt0 = station(n)

   447      txt1 = sprintf("%5.2f", lat_ob(n))

   448      txt2 = sprintf("%5.2f", lon_ob(n))

   449      txt3 = year_ob(n)

   450      txt4 = sprintf("%5.2f", M_score(n,0))

   451      txt5 = sprintf("%5.2f", M_score(n,1))

   452      txt6 = sprintf("%5.2f", M_score(n,2))

   453      txt7 = sprintf("%5.2f", M_score(n,3))

   454      txt8 = sprintf("%5.2f", avg(M_score(n,:)))

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

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

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

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

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

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

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

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

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

   466      set_line(lines,nline,row_footer)

   467   end do

   469 ; last row, summary

   470   set_line(lines,nline,row_header)

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

   473   txt1 = "-"

   474   txt2 = "-"

   475   txt3 = "-"

   476   txt4 = M_fluxnet_nee

   477   txt5 = M_fluxnet_rad

   478   txt6 = M_fluxnet_lh

   479   txt7 = M_fluxnet_sh

   480   txt8 = M_fluxnet_all

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

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

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

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

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

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

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

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

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

   492   set_line(lines,nline,row_footer)

   493 ;-----------------------------------------------

   494   set_line(lines,nline,table_footer)

   495   set_line(lines,nline,footer)

   497 ; Now write to an HTML file.

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

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

   500     asciiwrite(output_html,lines(idx))

   501   else

   502    print ("error?")

   503   end if

   504   delete (idx)

   506 ;**************************************************************************************

   507 ; update score

   508 ;**************************************************************************************

   510   if (isvar("compare")) then

   511      system("sed -e '1,/M_fluxnet_nee/s/M_fluxnet_nee/"+M_fluxnet_nee+"/' "+html_name2+" > "+html_new2+";"+ \

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

   513             "sed -e '1,/M_fluxnet_rad/s/M_fluxnet_rad/"+M_fluxnet_rad+"/' "+html_name2+" > "+html_new2+";"+ \

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

   515             "sed -e '1,/M_fluxnet_lh/s/M_fluxnet_lh/"+M_fluxnet_lh+"/' "+html_name2+" > "+html_new2+";"+ \

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

   517             "sed -e '1,/M_fluxnet_sh/s/M_fluxnet_sh/"+M_fluxnet_sh+"/' "+html_name2+" > "+html_new2+";"+ \

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

   519   end if

   521   system("sed s#M_fluxnet_nee#"+M_fluxnet_nee+"# "+html_name+" > "+html_new+";"+ \

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

   523          "sed s#M_fluxnet_rad#"+M_fluxnet_rad+"# "+html_name+" > "+html_new+";"+ \

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

   525          "sed s#M_fluxnet_lh#"+M_fluxnet_lh+"# "+html_name+" > "+html_new+";"+ \

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

   527          "sed s#M_fluxnet_sh#"+M_fluxnet_sh+"# "+html_name+" > "+html_new+";"+ \

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

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

   531 ; add total score and write to file

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

   533   ;M_total = M_fluxnet_all

   534   M_total = 0.00 ; we don't actually use these scores

   536   asciiwrite("M_save.fluxnet", M_total)

   538 ;***************************************************************************

   539 ; output plots

   540 ;***************************************************************************

   541   output_dir = model_name+"/fluxnet"

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

   544 ;***************************************************************************

   545 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
parent 0	0c6405ab2ff4
permissions	-rw-r--r--