clamp: energy/99.new.ncl@4be95183fbcd

     1 ;************************************************************

     2 ; required command line input parameters:

     3 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

4 ;

     5 ; using gsn_table for all

     6 ; output: line plot for each site (4 fields)

     7 ;         table for M_score

     8 ;************************************************************

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

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

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

    12 ;************************************************************

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

    14 begin

    15 ; add line to ascci/html file

    17   nnewlines = dimsizes(newlines)

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

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

    20     return

    21   end if

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

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

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;*************************************************************

    28 begin

    30   plot_type     = "ps"

    31   plot_type_new = "png"

    33 ;************************************************

    34 ; read model data

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

    37 ; model_name = "i01.10cn"

    38   model_name = "i01.10casa"

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

    41   film  = "lnd_T42.nc"

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

    44   xm    = fm->lon

    45   ym    = fm->lat

    46 ;------------------------------------------------

    47   nlat = dimsizes(ym)

    48   nlon = dimsizes(xm)

    50 ; for 4 fields, 12-monthly

    51   nmon      = 12

    52   nfield    = 4

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

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

    57 ; Model_units [=] gC/m2/s

    58 ; 12. = molecular weight of C

    59 ; u mol = 1e-6 mol

    60   factor = 1e6 /12.

    62   ENERGY ="new"

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

    65 ; read data: observed

    66 ;************************************************

    68  station = (/"BOREAS_NSA_OBS" \

    69             ,"CastelPorziano" \

    70             ,"Hyytiala" \

    71             ,"Kaamanen" \

    72             ,"LBA_Tapajos_KM67" \

    73             ,"Lethbridge" \

    74             ,"Tharandt" \

    75             ,"Vielsalm" \

    76             /)

    78  year_ob = (/"1994-2004" \

    79             ,"1997-2003" \

    80             ,"1996-2003" \

    81             ,"2000-2003" \

    82             ,"2002-2005" \

    83             ,"1999-2004" \

    84             ,"1996-2003" \

    85             ,"1998-2003" \

    86             /)

    88  field   = (/"CO2 Flux" \

    89             ,"Net Radiation" \

    90             ,"Latent Heat" \

    91             ,"Sensible Heat" \

    92             /)

    94  nstation  = dimsizes(station)

    95  nmon      = 12

    96  nfield    = dimsizes(field)

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

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

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

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

   103  diro_root  = "/fis/cgd/cseg/people/jeff/clamp_data/fluxnet/"

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

   106  do n = 0,nstation-1

   108 ;-------------------------------------------------

   109 ;   get ob data

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

   112     filo = station(n)+"_"+year_ob(n)+"_monthly.nc"

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

   115     print (filo)

   117     lon_ob(n) = fo->lon

   118     lat_ob(n) = fo->lat

   120     data      = fo->CO2_FLUX

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

   122     delete (data)

   124     data      = fo->RAD_FLUX

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

   126     delete (data)

   128     data      = fo->LH_FLUX

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

   130     delete (data)

   132     data      = fo->SH_FLUX

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

   134     delete (data)

   136     delete (fo)

   137 ;---------------------------------------------------

   138 ;   get model data

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

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

   143     print (film)

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

   147   data = fm->NEE

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

   149   delete (data)

   151 ; data  = fm->LATENT

   152   data1 = fm->FCEV

   153   data2 = fm->FCTR

   154   data3 = fm->FGEV

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

   156   delete (data1)

   157   delete (data2)

   158   delete (data3)

   160 ; data = fm->SENSIBLE

   161   data  = fm->FSH

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

   163   delete (data)

   165 ; data  = fm->NETRAD

   166   data1 = fm->FSA

   167   data2 = fm->FIRA

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

   169   delete (data1)

   170   delete (data2)

   172 else

   174   data = fm->NEE

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

   176   delete (data)

   178   data = fm->NETRAD

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

   180   delete (data)

   182   data = fm->LATENT

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

   184   delete (data)

   186 ; data = fm->SENSIBLE

   187   data = fm->FSH

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

   189   delete (data)

   190 end if

   192     delete (fm)

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

   195 ; interpolate model data into observed station

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

   197 ;************************************************************

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

   200   ym(0) = -90.

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

   204 ; print (yy)

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

   208   delete (yy)

   210  end do

   211 ;************************************************************

   212 ; compute correlation coef and M score

   213 ;************************************************************

   215  score_max = 5.

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

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

   220  do n=0,nstation-1

   221  do m=0,nfield-1

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

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

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

   225  end do

   226  end do

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

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

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

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

   232  M_all = M_co2+ M_rad +M_lh + M_sh

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

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

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

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

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

   240 ;*******************************************************************

   241 ; for station line plot

   242 ;*******************************************************************

   244 ; for x-axis in xyplot

   245   mon = ispan(1,12,1)

   246   mon@long_name = "month"

   248   res                   = True               ; plot mods desired

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

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

   251 ;-------------------------------------------------------------------------

   252 ; Add a boxed legend using the more simple method

   254   res@pmLegendDisplayMode    = "Always"

   255 ; res@pmLegendWidthF         = 0.1

   256   res@pmLegendWidthF         = 0.08

   257   res@pmLegendHeightF        = 0.06

   258 ; res@pmLegendOrthogonalPosF = -1.17

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

   260   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   262 ; res@pmLegendParallelPosF   =  0.18

   263   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   265 ; res@lgPerimOn             = False

   266   res@lgLabelFontHeightF     = 0.015

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

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

   269 ; for panel plot

   270   res@gsnFrame     = False                   ; Do not draw plot

   271   res@gsnDraw      = False                   ; Do not advance frame

   273   pres                            = True     ; panel plot mods desired

   274   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   275                                              ; indiv. plots in panel

   276   pres@gsnMaximize                = True     ; fill the page

   277 ;-------------------------------------------------------------------

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

   280   plot_data!0 = "case"

   281   plot_data!1 = "month"

   283   do n = 0,nstation-1

   284 ;----------------------------

   285 ; for observed

   287     plot_name = station(n)+"_ob"

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

   289     res@tiMainString = title

   291     wks = gsn_open_wks (plot_type,plot_name)

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

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

   295     plot_data@long_name = field(0)

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

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

   299     plot_data@long_name = field(1)

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

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

   303     plot_data@long_name = field(2)

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

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

   307     plot_data@long_name = field(3)

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

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

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

   313            "rm "+plot_name+"."+plot_type)

   315     clear (wks)

   316     delete (plot)

   317 ;----------------------------

   318 ; for model_vs_ob

   320     plot_name = station(n)+"_model_vs_ob"

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

   322     res@tiMainString = title

   324     wks = gsn_open_wks (plot_type,plot_name)

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

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

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

   329     plot_data@long_name = field(0)

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

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

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

   334     plot_data@long_name = field(1)

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

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

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

   339     plot_data@long_name = field(2)

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

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

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

   344     plot_data@long_name = field(3)

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

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

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

   350            "rm "+plot_name+"."+plot_type)

   352     clear (wks)

   353     delete (plot)

   354  end do

   356 ;*******************************************************************

   357 ; html table of site: observed

   358 ;*******************************************************************

   359   output_html = "line_ob.html"

   361   header = (/"<HTML>" \

   362             ,"<HEAD>" \

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

   364             ,"</HEAD>" \

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

   366             /)

   367   footer = "</HTML>"

   369   table_header = (/ \

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

   371        ,"<tr>" \

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

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

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

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

   376        ,"</tr>" \

   377        /)

   378   table_footer = "</table>"

   379   row_header = "<tr>"

   380   row_footer = "</tr>"

   382   lines = new(50000,string)

   383   nline = 0

   385   set_line(lines,nline,header)

   386   set_line(lines,nline,table_header)

   387 ;-----------------------------------------------

   388 ; row of table

   390   do n = 0,nstation-1

   391      set_line(lines,nline,row_header)

   393      txt0 = station(n)

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

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

   396      txt3 = year_ob(n)

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

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

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

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

   403      set_line(lines,nline,row_footer)

   404   end do

   405 ;-----------------------------------------------

   406   set_line(lines,nline,table_footer)

   407   set_line(lines,nline,footer)

   409 ; Now write to an HTML file.

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

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

   412     asciiwrite(output_html,lines(idx))

   413   else

   414    print ("error?")

   415   end if

   416   delete (idx)

   418 ;*******************************************************************

   419 ; score and line table : model vs observed

   420 ;*******************************************************************

   421   output_html = "score+line_vs_ob.html"

   423   header = (/"<HTML>" \

   424             ,"<HEAD>" \

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

   426             ,"</HEAD>" \

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

   428             /)

   429   footer = "</HTML>"

   431   delete (table_header)

   432   table_header = (/ \

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

   434        ,"<tr>" \

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

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

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

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

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

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

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

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

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

   444        ,"</tr>" \

   445        /)

   446   table_footer = "</table>"

   447   row_header = "<tr>"

   448   row_footer = "</tr>"

   450   lines = new(50000,string)

   451   nline = 0

   453   set_line(lines,nline,header)

   454   set_line(lines,nline,table_header)

   455 ;-----------------------------------------------

   456 ; row of table

   458   do n = 0,nstation-1

   459      set_line(lines,nline,row_header)

   461      txt0 = station(n)

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

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

   464      txt3 = year_ob(n)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   481      set_line(lines,nline,row_footer)

   482   end do

   484 ; last row, summary

   485   set_line(lines,nline,row_header)

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

   488   txt1 = "-"

   489   txt2 = "-"

   490   txt3 = "-"

   491   txt4 = M_energy_co2

   492   txt5 = M_energy_rad

   493   txt6 = M_energy_lh

   494   txt7 = M_energy_sh

   495   txt8 = M_energy_all

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

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

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

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

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

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

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

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

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

   507   set_line(lines,nline,row_footer)

   508 ;-----------------------------------------------

   509   set_line(lines,nline,table_footer)

   510   set_line(lines,nline,footer)

   512 ; Now write to an HTML file.

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

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

   515     asciiwrite(output_html,lines(idx))

   516   else

   517    print ("error?")

   518   end if

   519   delete (idx)

   521 ;***************************************************************************

   522 ; output plots

   523 ;***************************************************************************

   524   output_dir = model_name+"/energy"

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

   527 ;***************************************************************************

   528 end

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