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

; remove histogram plots.

;

; 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

  plot_type     = "ps"

  plot_type_new = "png"

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

; read data: model       

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

;film  = "b30.061n_1995-2004_MONS_climo_lnd.nc"

;model_name = "b30.061n"

;model_grid = "T31"

;film  = "newcn05_ncep_1i_MONS_climo_lnd.nc"

;model_name = "newcn"

;model_grid = "1.9"

;film  = "i01.06cn_1798-2004_MONS_climo.nc"

;model_name = "06cn"

;model_grid = "T42"

;film  = "i01.06casa_1798-2004_MONS_climo.nc"

;model_name = "06casa"

;model_grid = "T42"

 film  = "i01.10cn_1948-2004_MONS_climo.nc"

 model_name = "10cn"

 model_grid = "T42"

;film  = "i01.10casa_1948-2004_MONS_climo.nc"

;model_name = "10casa"

;model_grid = "T42"

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

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

 laimod  = fm->TLAI

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

; read data: observed

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

  ob_name = "MODIS MOD 15A2 2000-2005"

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

  filo1  = "land_class_"+model_grid+".nc"

  filo2  = "LAI_2000-2005_MONS_"+model_grid+".nc"

  fo1 = addfile(diro+filo1,"r")

  fo2 = addfile(diro+filo2,"r")

  classob    = tofloat(fo1->LAND_CLASS)               

  laiob      = fo2->LAI

; input observed data has 20 land-type classes

  nclass = 20

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

; global res

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

  resg                      = True             ; Use plot options

  resg@cnFillOn             = True             ; Turn on color fill

  resg@gsnSpreadColors      = True             ; use full colormap

  resg@cnLinesOn            = False            ; Turn off contourn lines

  resg@mpFillOn             = False            ; Turn off map fill

  resg@cnLevelSelectionMode = "ManualLevels"   ; Manual contour invtervals

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

; plot global land class: observed

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

  resg@cnMinLevelValF       = 1.              ; Min level

  resg@cnMaxLevelValF       = 19.             ; Max level

  resg@cnLevelSpacingF      = 1.              ; interval

; global contour ob

  classob@_FillValue = 1.e+36

  classob = where(classob.eq.0,classob@_FillValue,classob)

  plot_name = "global_class_ob"

  title     = ob_name

  resg@tiMainString  = title

  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  gsn_define_colormap(wks,"gui_default")     ; choose colormap

  plot = gsn_csm_contour_map_ce(wks,classob,resg)   

  frame(wks)

  clear (wks)

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

; for html table : all 4 components (Mean, Max, Phase, Growth)

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

; column (not including header column)

  component = (/"Mean","Max","Phase","Growth"/)

  col_head2 = (/"observed",model_name,"M_score" \

               ,"observed",model_name,"M_score" \

               ,"observed",model_name,"M_score" \

               ,"observed",model_name,"M_score" \

               /)

  n_comp = dimsizes(component) 

  ncol   = dimsizes(col_head2)

; row (not including header row)

  row_head  = (/"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" \             

               ,"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 )

; M_comp

  M_comp = (/"M_lai_mean","M_lai_max","M_lai_phase","M_lai_grow"/) 

; total M_score

  M_total = 0.

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

; use land-type class to bin 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.

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

; go through all components

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

  do n = 0,n_comp-1

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

; get data:

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

; (A) Mean

  if (n .eq. 0) then

     data_ob  = dim_avg_Wrap(laiob (lat|:,lon|:,time|:))

     data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))

  end if

; (B) Max

  if (n .eq. 1) then

;    observed  

     data_ob = laiob(0,:,:)

     s       = laiob(:,0,0)

     data_ob@long_name = "Leaf Area Index Max"

     dsizes_z = dimsizes(laiob)

     nlat     = dsizes_z(1)

     nlon     = dsizes_z(2)

     do j = 0,nlat-1

     do i = 0,nlon-1

        s = laiob(:,j,i) 

        data_ob(j,i) = max(s)

     end do

     end do

     delete (s)

     delete (dsizes_z)          

;    model  

     data_mod = laimod(0,:,:)

     s        = laimod(:,0,0)

     data_mod@long_name = "Leaf Area Index Max"

     dsizes_z = dimsizes(laimod)

     nlat     = dsizes_z(1)

     nlon     = dsizes_z(2)

     do j = 0,nlat-1

     do i = 0,nlon-1

        s = laimod(:,j,i) 

        data_mod(j,i) = max(s)

     end do

     end do

     delete (s)

     delete (dsizes_z)          

  end if

; (C) phase

  if (n .eq. 2) then  

;    observed

     data_ob = laiob(0,:,:)

     s       = laiob(:,0,0)

     data_ob@long_name = "Leaf Area Index Max Month"

     dsizes_z = dimsizes(laiob)

     nlat     = dsizes_z(1)

     nlon     = dsizes_z(2)

     do j = 0,nlat-1

     do i = 0,nlon-1

        s = laiob(:,j,i) 

        data_ob(j,i) = maxind(s) + 1

     end do

     end do

     delete (s)

     delete (dsizes_z)          

;    model

     data_mod = laimod(0,:,:)

     s        = laimod(:,0,0)

     data_mod@long_name = "Leaf Area Index Max Month"

     dsizes_z = dimsizes(laimod)

     nlat     = dsizes_z(1)

     nlon     = dsizes_z(2)

     do j = 0,nlat-1

     do i = 0,nlon-1

        s = laimod(:,j,i) 

        data_mod(j,i) = maxind(s) + 1

     end do

     end do

     delete (s)

     delete (dsizes_z)          

  end if

; (D) grow day

  if (n .eq. 3) then   

     day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)

;    observed

     data_ob = laiob(0,:,:)

     data_ob@long_name = "Days of Growing Season"

     dsizes_z = dimsizes(laiob)

     ntime    = dsizes_z(0)

     nlat     = dsizes_z(1)

     nlon     = dsizes_z(2)

     do j = 0,nlat-1

     do i = 0,nlon-1

        nday = 0.

        do k = 0,ntime-1

           if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then

              nday = nday + day_of_data(k)

           end if

        end do

        data_ob(j,i) = nday

     end do

     end do

     delete (dsizes_z)     

;    model

     data_mod = laimod(0,:,:)

     data_mod@long_name = "Days of Growing Season"

     dsizes_z = dimsizes(laimod)

     ntime    = dsizes_z(0)

     nlat     = dsizes_z(1)

     nlon     = dsizes_z(2)

     do j = 0,nlat-1

     do i = 0,nlon-1

        nday = 0.

        do k = 0,ntime-1

           if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then

              nday = nday + day_of_data(k)

           end if

        end do

        data_mod(j,i) = nday

     end do

     end do

     delete (dsizes_z)

  end if

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

; put data into bins

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

  base_1D  = ndtooned(classob)

  data1_1D = ndtooned(data_ob)

  data2_1D = ndtooned(data_mod)

; output for data in bins

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

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

; put data into bins

  do nd=0,1

;   See if we are doing data1 (nd=0) or data2 (nd=1).

    base = base_1D

    if(nd.eq.0) then

      data = data1_1D

    else

      data = data2_1D

    end if

;   Loop through each range, using base.

    do i=0,nr-2

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

;        print("")

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

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

      else

;        print("")

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

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

      end if

;     Calculate average 

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

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

        count(nd,i)   = dimsizes(idx)

      else

        yvalues(nd,i) = yvalues@_FillValue

        count(nd,i)   = 0

      end if

;#############################################################

;     set the following 4 classes to _FillValue:

;     Water Bodies(0), Urban and Build-Up(13),

;     Permenant Snow and Ice(15), Unclassified(17)

      if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then

         yvalues(nd,i) = yvalues@_FillValue

         count(nd,i)   = 0

      end if

;############################################################# 

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

;     Clean up for next time in loop.

      delete(idx)

    end do

    delete(data)

  end do

  delete (base)

  delete (base_1D)

  delete (data1_1D)

  delete (data2_1D)

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

; compute correlation coef and M score 

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

  u = yvalues(0,:)

  v = yvalues(1,:)

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

  uu = u(good)

  vv = v(good)

; compute correlation coef

  cc = esccr(uu,vv,0)

  if (n .eq. 2) then

     bias   = avg(abs(vv-uu))

     bias   = where((bias.gt. 6.),12.-bias,bias)

     Mscore = ((6. - bias)/6.)*5.

     M_score = sprintf("%.2f", Mscore)

  else

     bias  = sum(abs(vv-uu)/abs(vv+uu))

     Mscore = (1.- (bias/dimsizes(uu)))*5.

     M_score = sprintf("%.2f", Mscore)

  end if

; compute M_total

  M_total = M_total + Mscore

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

; output M_score

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

  print (Mscore)

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

; output to html table

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

  nn = n*3

  do i=0,nrow-2

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

     text4(i,nn+1) = sprintf("%.2f",v(i))

     text4(i,nn+2) = "-"

  end do

  text4(nrow-1,nn)   = sprintf("%.2f",avg(u))

  text4(nrow-1,nn+1) = sprintf("%.2f",avg(v))

  text4(nrow-1,nn+2) = M_score

  delete (u)

  delete (v)

  delete (uu)

  delete (vv)

  delete (yvalues)

  delete (good)

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

; global res changes for each component

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

  delta = 0.00001  

  if (n .eq. 0) then

     resg@cnMinLevelValF       = 0.             

     resg@cnMaxLevelValF       = 10.             

     resg@cnLevelSpacingF      = 1.

     data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)

  end if

  if (n .eq. 1) then

     resg@cnMinLevelValF       = 0.             

     resg@cnMaxLevelValF       = 10.             

     resg@cnLevelSpacingF      = 1.

     data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)

  end if

  if (n .eq. 2) then

     resg@cnMinLevelValF       = 1.             

     resg@cnMaxLevelValF       = 12.             

     resg@cnLevelSpacingF      = 1.

     data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)

  end if

  if (n .eq. 3) then

     resg@cnMinLevelValF       = 60.             

     resg@cnMaxLevelValF       = 360.             

     resg@cnLevelSpacingF      = 20.

     data_ob@_FillValue = 1.e+36

     data_ob = where(data_ob .lt. 10.,data_ob@_FillValue,data_ob)

     data_mod@_FillValue = 1.e+36

     data_mod = where(data_mod .lt. 10.,data_mod@_FillValue,data_mod)        

  end if                                                        

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

; global contour : ob

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

  plot_name = "global_"+component(n)+"_ob"

  title     = ob_name

  resg@tiMainString  = title

  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  gsn_define_colormap(wks,"gui_default")     ; choose colormap

  plot = gsn_csm_contour_map_ce(wks,data_ob,resg)   

  frame(wks)

  clear (wks)

  delete (plot)

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

; global contour : model

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

  plot_name = "global_"+component(n)+"_model"  

  title     = "Model " + model_name

  resg@tiMainString  = title

  wks = gsn_open_wks (plot_type,plot_name)

  gsn_define_colormap(wks,"gui_default")     

  plot = gsn_csm_contour_map_ce(wks,data_mod,resg)   

  frame(wks)

  clear (wks)

  delete (plot)

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

; global contour: model vs ob

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

  plot_name = "global_"+component(n)+"_model_vs_ob"

  wks = gsn_open_wks (plot_type,plot_name)   

  gsn_define_colormap(wks,"gui_default")     

  plot=new(3,graphic)                        ; create graphic array

  resg@gsnFrame             = False          ; Do not draw plot 

  resg@gsnDraw              = False          ; Do not advance frame

; plot correlation coef

  gRes               = True

  gRes@txFontHeightF = 0.02

  gRes@txAngleF      = 90

  correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"

  gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

; plot ob

  title     = ob_name

  resg@tiMainString  = title

  plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg)       

; plot model

  title     = "Model "+ model_name

  resg@tiMainString  = title

  plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg) 

; plot model-ob

  if (n .eq. 0) then

     resg@cnMinLevelValF  = -2.           

     resg@cnMaxLevelValF  =  2.            

     resg@cnLevelSpacingF =  0.4

  end if

  if (n .eq. 1) then

     resg@cnMinLevelValF  = -6.           

     resg@cnMaxLevelValF  =  6.            

     resg@cnLevelSpacingF =  1.

  end if

  if (n .eq. 2) then

     resg@cnMinLevelValF  = -6.           

     resg@cnMaxLevelValF  =  6.            

     resg@cnLevelSpacingF =  1.

  end if

  if (n .eq. 3) then

     resg@cnMinLevelValF  = -100.           

     resg@cnMaxLevelValF  =  100.            

     resg@cnLevelSpacingF =  20.

  end if            

  zz = data_mod

  zz = data_mod - data_ob

  title = "Model_"+model_name+" - Observed"

  resg@tiMainString    = title

  plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 

; plot panel

  pres                            = True        ; panel plot mods desired

  pres@gsnMaximize                = True        ; fill the page

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

  clear (wks)

  delete (plot)

  end do

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

; html table

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

  output_html = "table_model_vs_ob.html"

  header_text = "<H1>LAI: Model "+model_name+" vs Observed</H1>" 

  header = (/"<HTML>" \

            ,"<HEAD>" \

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

            ,"</HEAD>" \

            ,header_text \

            /) 

  footer = "</HTML>"

  table_header = (/ \

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

       ,"<tr>" \

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

       ,"   <th bgcolor=DDDDDD colspan=3>"+component(0)+"</th>" \

       ,"   <th bgcolor=DDDDDD colspan=3>"+component(1)+"</th>" \

       ,"   <th bgcolor=DDDDDD colspan=3>"+component(2)+"</th>" \

       ,"   <th bgcolor=DDDDDD colspan=3>"+component(3)+"</th>" \

       ,"</tr>" \

       ,"<tr>" \

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

       ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

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

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

       ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

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

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

       ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

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

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

       ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

       ,"   <th bgcolor=DDDDDD >M_score</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)

     txt7  = text4(n,5)

     txt8  = text4(n,6)

     txt9  = text4(n,7)

     txt10 = text4(n,8)

     txt11 = text4(n,9)

     txt12 = text4(n,10)

     txt13 = text4(n,11) 

     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,"<th>"+txt7+"</th>")

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

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

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

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

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

     set_line(lines,nline,"<th>"+txt13+"</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

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

; write total score to file

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

  asciiwrite("M_save.lai", M_total)

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

end