;********************************************************

; required command line input parameters:

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

;

; histogram normalized by rain and compute correleration

;**************************************************************

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

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

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

;**************************************************************

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

begin

; add line to ascci/html file

  nnewlines = dimsizes(newlines)

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

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

    return

  end if 

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

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

  nline = nline + nnewlines

  return 

end

;**************************************************************

; Main code.

begin

  nclass = 20

  plot_type     = "ps"

  plot_type_new = "png"

;************************************************

; read data: model       

;************************************************

 co2_i = 283.1878

 co2_f = 364.1252

 model_grid = "T42"

;model_name_i = "i01.07cn"

;model_name_f = "i01.10cn"

 model_name_i = "i01.07casa"

 model_name_f = "i01.10casa"

 model_name = model_name_f

 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

 film_i = model_name_i + "_1990-2004_ANN_climo.nc"

 film_f = model_name_f + "_1990-2004_ANN_climo.nc"

 fm_i   = addfile (dirm+film_i,"r")

 fm_f   = addfile (dirm+film_f,"r")

 npp_i  = fm_i->NPP

 npp_f  = fm_f->NPP

;************************************************

; read data: observed

;************************************************

  ob_name = "MODIS MOD 15A2 2000-2005"

  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"

  filo  = "land_class_"+model_grid+".nc"

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

  classob = tofloat(fo->LAND_CLASS)                

;*******************************************************************

; Calculate "nice" bins for binning the data in equally spaced ranges

;********************************************************************

  nclassn     = nclass + 1

  range       = fspan(0,nclassn-1,nclassn)

; print (range)

; Use this range information to grab all the values in a

; particular range, and then take an average.

  nr           = dimsizes(range)

  nx           = nr-1

  xvalues      = new((/2,nx/),float)

  xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

  dx           = xvalues(0,1) - xvalues(0,0)       ; range width

  dx4          = dx/4                              ; 1/4 of the range

  xvalues(1,:) = xvalues(0,:) - dx/5.

; get data

  DATA11_1D = ndtooned(classob)

  DATA12_1D = ndtooned(npp_i)

  DATA22_1D = ndtooned(npp_f)

  yvalues      = new((/2,nx/),float)

  mn_yvalues   = new((/2,nx/),float)

  mx_yvalues   = new((/2,nx/),float)

  do nd=0,1

; See if we are doing model or observational data.

    if(nd.eq.0) then

      data_ob  = DATA11_1D

      data_mod = DATA12_1D

    else

      data_ob  = DATA11_1D

      data_mod = DATA22_1D

    end if

; Loop through each range and check for values.

    do i=0,nr-2

      if (i.ne.(nr-2)) then

;        print("")

;        print("In range ["+range(i)+","+range(i+1)+")")

         idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))

      else

;        print("")

;        print("In range ["+range(i)+",)")

         idx = ind(data_ob.ge.range(i))

      end if

; Calculate average, and get min and max.

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

        yvalues(nd,i)    = avg(data_mod(idx))

        mn_yvalues(nd,i) = min(data_mod(idx))

        mx_yvalues(nd,i) = max(data_mod(idx))

        count = dimsizes(idx)

      else

        count            = 0

        yvalues(nd,i)    = yvalues@_FillValue

        mn_yvalues(nd,i) = yvalues@_FillValue

        mx_yvalues(nd,i) = yvalues@_FillValue

      end if

;     print(nd + ": " + count + " points, avg = " + yvalues(nd,i))

;     print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

; Clean up for next time in loop.

      delete(idx)

    end do

    delete(data_ob)

    delete(data_mod)

  end do

;============================

;compute beta

;============================

 nsec_per_year = 60*60*24*365

 u = yvalues(0,:)

 v = yvalues(1,:)

 good = ind(.not.ismissing(u) .and. .not.ismissing(v))

 uu = u(good)* nsec_per_year 

 vv = v(good)* nsec_per_year 

 n_biome = dimsizes(uu)

 beta_biome = new((/n_biome/),float)

 beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)

 beta_biome_avg = avg(beta_biome)

 print (beta_biome_avg)

;*******************************************************************

; for html table

;*******************************************************************

; column (not including header column)

  col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta"/)

  ncol = dimsizes(col_head)

; row (not including header row)

  row_head  = (/"Water Bodies" \

               ,"Evergreen Needleleaf Forests" \

               ,"Evergreen Broadleaf Forests" \

               ,"Deciduous Needleleaf Forest" \

               ,"Deciduous Broadleaf Forests" \

               ,"Mixed Forests" \                      

               ,"Closed Bushlands" \                   

               ,"Open Bushlands" \                     

               ,"Woody Savannas (S. Hem.)" \           

               ,"Savannas (S. Hem.)" \                 

               ,"Grasslands" \                         

               ,"Permanent Wetlands" \                 

               ,"Croplands" \                                           

               ,"Cropland/Natural Vegetation Mosaic" \ 

               ,"Permanent Snow and Ice" \             

               ,"Barren or Sparsely Vegetated" \                             

               ,"Woody Savannas (N. Hem.)" \           

               ,"Savannas (N. Hem.)" \

               ,"All Biome" \                

               /)  

  nrow = dimsizes(row_head)                  

; arrays to be passed to table. 

  text4 = new ((/nrow, ncol/),string )

 do i=0,nrow-2

  text4(i,0) = sprintf("%.2f",co2_i)

  text4(i,1) = sprintf("%.2f",co2_f)

  text4(i,2) = sprintf("%.2f",uu(i))

  text4(i,3) = sprintf("%.2f",vv(i))

  text4(i,4) = sprintf("%.2f",beta_biome(i))

 end do

  text4(nrow-1,0) = "-"

  text4(nrow-1,1) = "-"

  text4(nrow-1,2) = "-"

  text4(nrow-1,3) = "-"

  text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)

;**************************************************

; html table

;**************************************************

  output_html = "table_biome.html"

  header_text = "<H1>Beta Factor: Model "+model_name+"</H1>" 

  header = (/"<HTML>" \

            ,"<HEAD>" \

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

            ,"</HEAD>" \

            ,header_text \

            /) 

  footer = "</HTML>"

  table_header = (/ \

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

       ,"<tr>" \

       ,"   <th bgcolor=DDDDDD >Biome Class</th>" \

       ,"   <th bgcolor=DDDDDD >CO2_i</th>" \

       ,"   <th bgcolor=DDDDDD >CO2_f</th>" \

       ,"   <th bgcolor=DDDDDD >NPP_i</th>" \

       ,"   <th bgcolor=DDDDDD >NPP_f</th>" \

       ,"   <th bgcolor=DDDDDD >Beta</th>" \

       ,"</tr>" \

       /)

  table_footer = "</table>"

  row_header = "<tr>"

  row_footer = "</tr>"

  lines = new(50000,string)

  nline = 0

  set_line(lines,nline,header)

  set_line(lines,nline,table_header)

;-----------------------------------------------

;row of table

  do n = 0,nrow-1

     set_line(lines,nline,row_header)

     txt1  = row_head(n)

     txt2  = text4(n,0)

     txt3  = text4(n,1)

     txt4  = text4(n,2)

     txt5  = text4(n,3)

     txt6  = text4(n,4)

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

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

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

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

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

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

     set_line(lines,nline,row_footer)

  end do

;-----------------------------------------------

  set_line(lines,nline,table_footer)

  set_line(lines,nline,footer) 

; Now write to an HTML file.

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

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

    asciiwrite(output_html,lines(idx))

  else

   print ("error?")

  end if

end