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

; add table

; add time series plot

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

; find year_ob >= 4

; add year_ob_i and year_ob_f

; change long_name

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

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

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

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

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

begin

  plot_type     = "ps"

  plot_type_new = "png"

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

; model data info

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

  model = "cn"

; model = "casa"

  model_name = "i01.10" + model

  ENERGY ="new"

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

; observed data info

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

 station = (/"ARM_Oklahoma" \

            ,"ARM_Oklahoma_burn" \

            ,"ARM_Oklahoma_control" \

            ,"Atqasuk" \

            ,"Audubon" \

            ,"AustinCary" \

            ,"Bartlett" \

            ,"Bondville" \

            ,"Brookings" \

            ,"Donaldson" \

            ,"Duke_Forest_Hardwoods" \

            ,"Duke_Forest_Open_Field" \

            ,"Duke_Forest_Pine" \

            ,"Fermi_Ag" \

            ,"Fermi_Prairie" \

            ,"Flagstaff_Managed" \

            ,"Flagstaff_Unmanaged" \

            ,"Flagstaff_Wildfire" \

            ,"FortPeck" \

            ,"FreemanRanch_mesquite" \

            ,"Goodwin_Creek" \

            ,"HarvardForest" \

            ,"HarvardForestHemlock" \

            ,"HowlandForestMain" \

            ,"HowlandForestWest" \

            ,"Ivotuk" \

            ,"KendallGrasslands" \

            ,"KennedySpaceCenterPine" \

            ,"KennedySpaceCenterScrub" \

            ,"LittleProspect" \

            ,"LostCreek" \

            ,"Mead-irrigated" \

            ,"Mead-irrigated-rotation" \

            ,"Mead-rainfed" \

            ,"Metolius_2nd_YoungPonderosaPine" \

            ,"MetoliusEyerly" \

            ,"MetoliusIntermediatePine" \

            ,"MetoliusOldPonderosaPine" \

            ,"MissouriOzark" \

            ,"Mize" \

            ,"MorganMonroe" \

            ,"NiwotRidge" \

            ,"NorthCarolina_cc" \

            ,"NorthCarolina_lp" \

            ,"ParkFalls" \

            ,"Rayonier" \

            ,"SantaRita" \

            ,"SkyOaks_Old" \

            ,"SkyOaks_PostFire" \

            ,"SkyOaks_Young" \

            ,"SylvaniaWilderness" \

            ,"Toledo" \

            ,"Tonzi" \

            ,"UCI_1850" \

            ,"UCI_1930" \

            ,"UCI_1964" \

            ,"UCI_1964wet" \

            ,"UCI_1981" \

            ,"UCI_1989" \

            ,"UCI_1998" \

            ,"UMBS" \

            ,"Vaira" \

            ,"WalkerBranch" \

            ,"WillowCreek" \

            ,"WindRiver" \

            ,"Wisconsin_ihw" \

            ,"Wisconsin_irp" \

            ,"Wisconsin_mrp" \

            ,"Wisconsin_myjp" \

            ,"Wisconsin_pb" \

            ,"Wisconsin_rpcc" \

            ,"Wisconsin_yhw" \

            ,"Wisconsin_yjp" \

            ,"Wisconsin_yrp" \

            /)

 year_ob = (/"2003-2006" \

            ,"2005-2006" \

            ,"2005-2006" \

            ,"1999-2006" \

            ,"2002-2006" \

            ,"2000-2005" \

            ,"2004-2005" \

            ,"1996-2006" \

            ,"2004-2006" \

            ,"1999-2004" \

            ,"2003-2005" \

            ,"2001-2005" \

            ,"2001-2005" \

            ,"2005-2006" \

            ,"2004-2006" \

            ,"2005-2006" \

            ,"2005-2006" \

            ,"2005-2006" \

            ,"2000-2006" \

            ,"2004-2006" \

            ,"2002-2006" \

            ,"1991-2004" \

            ,"2004-2004" \

            ,"1996-2004" \

            ,"1999-2004" \

            ,"2003-2006" \

            ,"2004-2006" \

            ,"2002-2002" \

            ,"2000-2006" \

            ,"2002-2005" \

            ,"2001-2005" \

            ,"2001-2005" \

            ,"2001-2005" \

            ,"2001-2005" \

            ,"2004-2005" \

            ,"2004-2005" \

            ,"2003-2005" \

            ,"1996-2000" \

            ,"2004-2006" \

            ,"1998-2004" \

            ,"1999-2005" \

            ,"1999-2003" \

            ,"2005-2006" \

            ,"2005-2006" \

            ,"1996-2003" \

            ,"1998-1998" \

            ,"2004-2006" \

            ,"1997-2006" \

            ,"2004-2006" \

            ,"1997-2006" \

            ,"2002-2006" \

            ,"2004-2005" \

            ,"2001-2006" \

            ,"2004-2005" \

            ,"2001-2005" \

            ,"2001-2005" \

            ,"2002-2004" \

            ,"2001-2005" \

            ,"2001-2005" \

            ,"2002-2005" \

            ,"1999-2003" \

            ,"2001-2006" \

            ,"1995-1999" \

            ,"1999-2005" \

            ,"1999-2004" \

            ,"2003-2003" \

            ,"2003-2003" \

            ,"2002-2005" \

            ,"2004-2004" \

            ,"2002-2002" \

            ,"2005-2005" \

            ,"2002-2002" \

            ,"2004-2005" \

            ,"2002-2002" \

            /)

 field   = (/"NEE Flux" \

            ,"Net Radiation" \

            ,"Latent Heat" \

            ,"Sensible Heat" \

            ,"GPP Flux" \

            ,"Respiration" \

            /)

 field_unit = (/"u mol/m2/s" \

               ,"W/m2" \

               ,"W/m2" \

               ,"W/m2" \

               ,"u mol/m2/s" \

               ,"u mol/m2/s" \

               /)

 nstation = dimsizes(station)

 nfield   = dimsizes(field)

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

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

 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"

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

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

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

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

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

 do n = 0, nstation-1

    diro = dir_root + station(n)+"/"

    filo = year_ob(n)+"_L4_m.nc"

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

    lat_ob(n) = fo->lat

    lon_ob(n) = fo->lon

    year      = fo->year

    year_station(n) = dimsizes(year)

    year_ob_i(n)    = year(0)

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

    delete (year) 

 end do

 delete (fo)

;print (year_ob_i)

;print (year_ob_f)

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

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

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

 station_long      = station(i_long_ob)

 lat_ob_long       = lat_ob(i_long_ob)

 lon_ob_long       = lat_ob(i_long_ob)

 year_ob_long      = year_ob(i_long_ob)

 year_ob_i_long    = year_ob_i(i_long_ob)

 year_ob_f_long    = year_ob_f(i_long_ob)

 year_station_long = year_station(i_long_ob)

 nstation_long     = dimsizes(station_long)

; print (i_long_ob)

; print (nstation_long)

; print (station_long)

; print (year_ob_i(i_long_ob))

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

;   get model data at observed lat-lon

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

    film = model_name+"_ameriflux_1990-2004_monthly.nc"

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

    xm   = fm->lon

    ym   = fm->lat 

    date = fm->date

    date_dim = dimsizes(date)

    nyear  = date_dim(0)

    nmonth =12

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

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

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

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

;   Model_units [=] gC/m2/s

;   12. = molecular weight of C

;   u mol = 1e-6 mol

    factor = 1.e6 /12.

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

;   interpolate model data into observed station

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

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

    ym(0) = -90.

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

if (ENERGY .eq. "old") then

  data = fm->NEE

  data_mod(0,:,:,:) = data(:,:,:) * factor

; data  = fm->LATENT

  data1 = fm->FCEV

  data2 = fm->FCTR

  data3 = fm->FGEV

  data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:) 

  delete (data1)

  delete (data2)

  delete (data3)

; data = fm->SENSIBLE

  data  = fm->FSH

  data_mod(3,:,:,:) = data(:,:,:) 

; data  = fm->NETRAD

  data1 = fm->FSA

  data2 = fm->FIRA

  data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:) 

  delete (data1)

  delete (data2)

else

  data = fm->NEE

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

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

  data_mod(0,:,:,:) = yy(:,:,:) * factor

; printVarSummary(yy)

  data = fm->NETRAD

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

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

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

  data = fm->LATENT

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

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

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

; data = fm->SENSIBLE

  data = fm->FSH

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

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

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

end if

  data = fm->GPP

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

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

  data_mod(4,:,:,:) = yy(:,:,:) * factor 

  if (model .eq. "cn") then

     data = fm->ER

  else

     data1 = fm->AR

     data2 = fm->HR

     data  = data1 + data2

     delete (data1)

     delete (data2)

  end if

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

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

  data_mod(5,:,:,:) = yy(:,:,:) * factor

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

  delete (data_mod)

  delete (fm) 

  delete (data)

  delete (yy)

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

; for station line plot

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

; for x-axis in xyplot

  mon = ispan(1,12,1)

  mon@long_name = "month"

  res                   = True               ; plot mods desired

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

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

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

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

; Add a boxed legend using the more simple method

  res@pmLegendDisplayMode    = "Always"

; res@pmLegendWidthF         = 0.1

  res@pmLegendWidthF         = 0.08

  res@pmLegendHeightF        = 0.06

; res@pmLegendOrthogonalPosF = -1.17

; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

  res@pmLegendOrthogonalPosF = -0.30  ;(downward)

; res@pmLegendParallelPosF   =  0.18

  res@pmLegendParallelPosF   =  0.23  ;(rightward)

; res@lgPerimOn             = False

  res@lgLabelFontHeightF     = 0.015

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

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

; for panel plot

  res@gsnFrame     = False                   ; Do not draw plot 

  res@gsnDraw      = False                   ; Do not advance frame

  pres                            = True     ; panel plot mods desired

  pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

                                             ; indiv. plots in panel

  pres@gsnMaximize                = True     ; fill the page

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

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

; get ob data at each site with long observation

  dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"

  do n = 0,nstation_long-1

; do n = 0,0

    year_setback = 0

    nyear = year_station_long(n)

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

       year_setback = 2006 -2004

    end if

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

       year_setback = 2005 -2004

    end if

    ntime = (nyear - year_setback) * nmonth

    print (nyear)

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

    diro = dir_root + station_long(n)+"/"

    filo = year_ob_long(n)+"_L4_m.nc"

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

;   print (diro)

;   print (filo)

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

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

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

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

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

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

    delete (fo)

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

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

    timeI = ispan(1,ntime,1)

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

    timeF@long_name = "year" 

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

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

; for model_vs_ob

    plot_name = station_long(n)+"_tseries_vs_ob"

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

    res@tiMainString = title

    wks = gsn_open_wks (plot_type,plot_name)

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

    i_year_mod_i = year_ob_i_long(n) - 1990 

    i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback

    i_year_ob_f =  nyear - year_setback - 1 

;   print (nyear)

;   print (i_year_ob_f)

;   print (i_year_mod_i)

;   print (i_year_mod_f)

    do i = 0,nfield-1                           

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

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

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

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

    end do

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

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

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

    clear (wks)  

    delete (plot)

    delete (data_ob)

    delete (timeI)    

    delete (timeF)

    delete (plot_data)

 end do

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

; for the following tables,

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

  isort = dim_pqsort(lat_ob_long,-1)

  station_sort = station_long(isort)

  year_ob_sort = year_ob_long(isort)

  lat_ob_sort  = lat_ob_long(isort)

  lon_ob_sort  = lon_ob_long(isort)

; print(isort)

; print(lat_ob_sort)

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

; html table of site: observed

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

  output_html = "tseries_vs_ob.html"

  header = (/"<HTML>" \

            ,"<HEAD>" \

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

            ,"</HEAD>" \

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

            /) 

  footer = "</HTML>"

  table_header = (/ \

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

       ,"<tr>" \

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

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

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

       ,"   <th bgcolor=DDDDDD >Observed</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,nstation_long-1

     set_line(lines,nline,row_header)

     txt0 = station_sort(n)

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

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

     txt3 = year_ob_sort(n)

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

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

  delete (idx)

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

; output plots

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

  output_dir = model_name+"/ameriflux"

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

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

end