clamp: lai/99.all.ncl.1@4be95183fbcd

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

     2 ; histogram normalized by rain and compute correleration

     3 ;********************************************************

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

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

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

     8 procedure pminmax(data:numeric,name:string)

     9 begin

    10   print ("min/max " + name + " = " + min(data) + "/" + max(data))

    11   if(isatt(data,"units")) then

    12     print (name + " units = " + data@units)

    13   end if

    14 end

    16 ; Main code.

    17 begin

    19  nclass = 20

    21  plot_type     = "ps"

    22  plot_type_new = "png"

    24 ;************************************************

    25 ; read in data: model

    26 ;************************************************

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

    28 ;model_name = "b30.061n"

    29 ;model_grid = "T31"

    31 ;film  = "newcn05_ncep_1i_MONS_climo_lnd.nc"

    32 ;model_name = "newcn"

    33 ;model_grid = "1.9"

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

    36 ;model_name = "06cn"

    37 ;model_grid = "T42"

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

    40 ;model_name = "06casa"

    41 ;model_grid = "T42"

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

    44  model_name = "10cn"

    45  model_grid = "T42"

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

    48 ;model_name = "10casa"

    49 ;model_grid = "T42"

    51  system("sed s#model_name#"+model_name+"# table.html > table.html.new")

    52  system("mv -f table.html.new table.html")

    53 ;------------------------------------------------

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

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

    57  laimod  = fm->TLAI

    59 ;************************************************

    60 ; read in data: observed

    61 ;************************************************

    63  ob_name = "MODIS MOD 15A2 2000-2005"

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

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

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

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

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

    72  classob    = tofloat(fo1->LAND_CLASS)

    73  laiob      = fo2->LAI

    75 ;*******************************************************************

    76 ; for plotting table

    77 ;*******************************************************************

    78 ; table header name

    79   table_header_name = "LAI"

    81 ; column (not including header column)

    83   col_header1 = (/"Mean","Max","Phase","Growth"/)

    84   col_header2 = (/"ob","model","M" \

    85                  ,"ob","model","M" \

    86                  ,"ob","model","M" \

    87                  ,"ob","model","M" \

    88                  /)

    90   ncol1 = dimsizes(col_header1)

    91   ncol2 = dimsizes(col_header2)

    92   ncol  = ncol2

    94 ; row (not including header row)

    95   row_header = (/"Water Bodies" \

    96                 ,"Evergreen Needleleaf Forests" \

    97                 ,"Evergreen Broadleaf Forests" \

    98                 ,"Deciduous Needleleaf Forest" \

    99                 ,"Deciduous Broadleaf Forests" \

   100                 ,"Mixed Forests" \

   101                 ,"Closed Bushlands" \

   102                 ,"Open Bushlands" \

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

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

   105                 ,"Grasslands" \

   106                 ,"Permanent Wetlands" \

   107                 ,"Croplands" \

   108                 ,"Urban and Built-Up" \

   109                 ,"Cropland/Natural Vegetation Mosaic" \

   110                 ,"Permanent Snow and Ice" \

   111                 ,"Barren or Sparsely Vegetated" \

   112                 ,"Unclassified" \

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

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

   115                 ,"All biome average" \

   116                 /)

   117   nrow = dimsizes(row_header)

   119 ; arrays to be passed to table.

   120   value = new ((/nrow, ncol/),string )

   122   table_length = 0.995

   124 ; Table header

   125   ncr1  = (/1,1/)               ; 1 row, 1 column

   126   xx1    = (/0.005,0.25/)        ; Start and end X

   127   yy1    = (/0.900,0.995/)       ; Start and end Y

   128   text1 = table_header_name

   129   res1               = True

   130   res1@txFontHeightF = 0.03

   131   res1@gsFillColor   = "CornFlowerBlue"

   133 ; Column header (equally space in x2)

   134   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

   135   xx21    = (/xx1(1),0.995/)        ; start from end of x1

   136   yy21    = (/0.9475,0.995/)        ; half of y1

   137   text21 = col_header1

   138   res21               = True

   139   res21@txFontHeightF = 0.015

   140   res21@gsFillColor   = "Gray"

   142   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   143   xx22    = xx21                    ; start from end of x1

   144   yy22    = (/0.900,0.9475/)       ; half of y1

   145   text22 = col_header2

   146   res22               = True

   147   res22@txFontHeightF = 0.015

   148   res22@gsFillColor   = "Gray"

   150 ; Row header (equally space in y2)

   151   ncr3  = (/nrow,1/)              ; 20 rows, 1 columns

   152   xx3    = xx1                      ; same as x1

   153   yy3    = (/1.0-table_length,0.900/) ; end at start of y1

   154   text3 = row_header

   155   res3               = True

   156   res3@txFontHeightF = 0.01

   157   res3@gsFillColor   = "Gray"

   159 ; Main table body

   160   ncr4  = (/nrow,ncol/) ; 5 rows, 5 columns

   161   xx4    = xx21                      ; Start and end x

   162   yy4    = yy3                      ; Start and end Y

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

   165   color_fill4      = new((/nrow,ncol/),string)

   166   color_fill4      = "white"

   167   color_fill4(nrow-1,:) = "CornFlowerBlue"

   169   res4               = True       ; Set up resource list

   170 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   171   res4@txFontHeightF = 0.015

   172   res4@gsFillColor   = color_fill4

   174   delete (color_fill4)

   176 ;************************************************

   177 ; plot global land class: observed

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

   179 ;global res

   181   resg                     = True             ; Use plot options

   182   resg@cnFillOn            = True             ; Turn on color fill

   183   resg@gsnSpreadColors     = True             ; use full colormap

   184 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

   185 ; resg@lbLabelAutoStride   = True

   186   resg@cnLinesOn           = False            ; Turn off contourn lines

   187   resg@mpFillOn            = False            ; Turn off map fill

   189   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   190   resg@cnMinLevelValF       = 1.              ; Min level

   191   resg@cnMaxLevelValF       = 19.             ; Max level

   192   resg@cnLevelSpacingF      = 1.              ; interval

   194 ;global contour ob

   195   classob@_FillValue = 1.e+36

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

   198   plot_name = "global_class_ob"

   199   title     = ob_name

   200   resg@tiMainString  = title

   202   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   203   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   205   plot = gsn_csm_contour_map_ce(wks,classob,resg)

   206   frame(wks)

   208   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   209   system("rm "+plot_name+"."+plot_type)

   210   system("rm "+plot_name+"-1."+plot_type_new)

   211   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   213   clear (wks)

   214 ;*******************************************************************

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

   216 ;********************************************************************

   217   nclassn     = nclass + 1

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

   219 ; print (range)

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

   222 ; particular range, and then take an average.

   224   nr           = dimsizes(range)

   225   nx           = nr-1

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

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

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

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

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

   231 ;-----------------------------------------------------------------

   232 ;(A) mean

   233 ;--------------------------------------------------------------------

   234 ; get data

   236   laiob_mean  = dim_avg_Wrap(laiob(lat|:,lon|:,time|:))

   237   laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))

   239   DATA11_1D = ndtooned(classob)

   240   DATA12_1D = ndtooned(laiob_mean)

   241   DATA22_1D = ndtooned(laimod_mean)

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

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

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

   247   do nd=0,1

   249 ; See if we are doing model or observational data.

   251     if(nd.eq.0) then

   252       data_ob  = DATA11_1D

   253       data_mod = DATA12_1D

   254     else

   255       data_ob  = DATA11_1D

   256       data_mod = DATA22_1D

   257     end if

   259 ; Loop through each range and check for values.

   261     do i=0,nr-2

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

   263 ;        print("")

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

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

   266       else

   267 ;        print("")

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

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

   270       end if

   272 ; Calculate average, and get min and max.

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

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

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

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

   278         count = dimsizes(idx)

   279       else

   280         count            = 0

   281         yvalues(nd,i)    = yvalues@_FillValue

   282         mn_yvalues(nd,i) = yvalues@_FillValue

   283         mx_yvalues(nd,i) = yvalues@_FillValue

   284       end if

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

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

   289 ; Clean up for next time in loop.

   291       delete(idx)

   292     end do

   293     delete(data_ob)

   294     delete(data_mod)

   295   end do

   296 ;-----------------------------------------------------------------

   297 ; compute correlation coef and M score

   299   u = yvalues(0,:)

   300   v = yvalues(1,:)

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

   303   uu = u(good)

   304   vv = v(good)

   306   ccrMean = esccr(uu,vv,0)

   307 ; print (ccrMean)

   309 ; new eq

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

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

   313   M_lai_mean = sprintf("%.2f", Mmean)

   314   system("sed s#M_lai_mean#"+M_lai_mean+"# table.html > table.html.new")

   315   system("mv -f table.html.new table.html")

   316   print (M_lai_mean)

   318  do i=0,nrow-2

   319   text4(i,0) = sprintf("%5.2f",u(i))

   320   text4(i,1) = sprintf("%5.2f",v(i))

   321   text4(i,2) = "-"

   322  end do

   323   text4(nrow-1,0) = sprintf("%5.2f",avg(u))

   324   text4(nrow-1,1) = sprintf("%5.2f",avg(v))

   325   text4(nrow-1,2) = sprintf("%5.2f",Mmean)

   327  delete (u)

   328  delete (v)

   329  delete (uu)

   330  delete (vv)

   331 ;--------------------------------------------------------------------

   332 ; histogram res

   334   resm                = True

   335   resm@gsnMaximize    = True

   336   resm@gsnDraw        = False

   337   resm@gsnFrame       = False

   338   resm@xyMarkLineMode = "Markers"

   339   resm@xyMarkerSizeF  = 0.014

   340   resm@xyMarker       = 16

   341   resm@xyMarkerColors = (/"Brown","Blue"/)

   342 ; resm@trYMinF        = min(mn_yvalues) - 10.

   343 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   344   resm@trYMinF        = min(mn_yvalues) - 2

   345   resm@trYMaxF        = max(mx_yvalues) + 4

   347   resm@tiYAxisString  = "Mean LAI (Leaf Area Index)"

   348   resm@tiXAxisString  = "Land Cover Type"

   350   max_bar = new((/2,nx/),graphic)

   351   min_bar = new((/2,nx/),graphic)

   352   max_cap = new((/2,nx/),graphic)

   353   min_cap = new((/2,nx/),graphic)

   355   lnres = True

   356   line_colors = (/"brown","blue"/)

   357 ;------------------------------------------------------------------

   358 ; Start the graphics.

   360   plot_name = "histogram_mean"

   361   title     = model_name + " vs Observed"

   362   resm@tiMainString  = title

   364   wks   = gsn_open_wks (plot_type,plot_name)

   365 ;-----------------------------

   366 ; Add a boxed legend using the more simple method

   368   resm@pmLegendDisplayMode    = "Always"

   369 ; resm@pmLegendWidthF         = 0.1

   370   resm@pmLegendWidthF         = 0.08

   371   resm@pmLegendHeightF        = 0.05

   372   resm@pmLegendOrthogonalPosF = -1.17

   373 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

   374 ; resm@pmLegendParallelPosF   =  0.18

   375   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   377 ; resm@lgPerimOn              = False

   378   resm@lgLabelFontHeightF     = 0.015

   379   resm@xyExplicitLegendLabels = (/"observed",model_name/)

   380 ;-----------------------------

   381   tRes  = True

   382   tRes@txFontHeightF = 0.025

   384   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"

   386   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   388   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

   389 ;-------------------------------

   390 ;Attach the vertical bar and the horizontal cap line

   392   do nd=0,1

   393     lnres@gsLineColor = line_colors(nd)

   394     do i=0,nx-1

   396       if(.not.ismissing(mn_yvalues(nd,i)).and. \

   397          .not.ismissing(mx_yvalues(nd,i))) then

   399 ; Attach the vertical bar, both above and below the marker.

   401         x1 = xvalues(nd,i)

   402         y1 = yvalues(nd,i)

   403         y2 = mn_yvalues(nd,i)

   404         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   406         y2 = mx_yvalues(nd,i)

   407         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   409 ; Attach the horizontal cap line, both above and below the marker.

   411         x1 = xvalues(nd,i) - dx4

   412         x2 = xvalues(nd,i) + dx4

   413         y1 = mn_yvalues(nd,i)

   414         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   416         y1 = mx_yvalues(nd,i)

   417         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   418       end if

   419     end do

   420   end do

   422   draw(xy)

   423   frame(wks)

   425   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   426   system("rm "+plot_name+"."+plot_type)

   427 ; system("rm "+plot_name+"-1."+plot_type_new)

   428 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   430   clear (wks)

   432  delete (DATA11_1D)

   433  delete (DATA12_1D)

   434  delete (DATA22_1D)

   435 ;delete (range)

   436 ;delete (xvalues)

   437  delete (yvalues)

   438  delete (mn_yvalues)

   439  delete (mx_yvalues)

   440  delete (good)

   441  delete (max_bar)

   442  delete (min_bar)

   443  delete (max_cap)

   444  delete (min_cap)

   445 ;-----------------------------------------------------------------

   446 ;global res

   448   resg                     = True             ; Use plot options

   449   resg@cnFillOn            = True             ; Turn on color fill

   450   resg@gsnSpreadColors     = True             ; use full colormap

   451 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

   452 ; resg@lbLabelAutoStride   = True

   453   resg@cnLinesOn           = False            ; Turn off contourn lines

   454   resg@mpFillOn            = False            ; Turn off map fill

   456   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   457   resg@cnMinLevelValF       = 0.              ; Min level

   458   resg@cnMaxLevelValF       = 10.             ; Max level

   459   resg@cnLevelSpacingF      = 1.              ; interval

   461 ;global contour ob

   463   delta = 0.00001

   464   laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean)

   466   plot_name = "global_mean_ob"

   467   title     = ob_name

   468   resg@tiMainString  = title

   470   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   471   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   473   plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg)

   474   frame(wks)

   476   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   477   system("rm "+plot_name+"."+plot_type)

   478   system("rm "+plot_name+"-1."+plot_type_new)

   479   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   481   clear (wks)

   482 ;------------------------------------------------------------------------

   483 ;global contour model

   485   plot_name = "global_mean_model"

   486   title     = "Model " + model_name

   487   resg@tiMainString  = title

   489   wks = gsn_open_wks (plot_type,plot_name)

   490   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   492   plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg)

   493   frame(wks)

   495   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   496   system("rm "+plot_name+"."+plot_type)

   497   system("rm "+plot_name+"-1."+plot_type_new)

   498   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   500   clear (wks)

   501 ;-----------------------------------------------------------------

   502 ;(B) max

   503 ;--------------------------------------------------------------------

   504 ; get data

   506 ; observed

   507   laiob_max = laiob(0,:,:)

   508   s         = laiob(:,0,0)

   509   laiob_max@long_name = "Leaf Area Index Max"

   511   dsizes_z = dimsizes(laiob)

   512   nlat     = dsizes_z(1)

   513   nlon     = dsizes_z(2)

   515   do j = 0,nlat-1

   516   do i = 0,nlon-1

   517      s = laiob(:,j,i)

   518      laiob_max(j,i) = max(s)

   519   end do

   520   end do

   522 ; print (min(y)+"/"+max(y))

   523   delete (s)

   524   delete (dsizes_z)

   525 ;-------------------------

   526 ; model

   527   laimod_max = laimod(0,:,:)

   528   s          = laimod(:,0,0)

   529   laimod_max@long_name = "Leaf Area Index Max"

   531   dsizes_z = dimsizes(laimod)

   532   nlat     = dsizes_z(1)

   533   nlon     = dsizes_z(2)

   535   do j = 0,nlat-1

   536   do i = 0,nlon-1

   537      s = laimod(:,j,i)

   538      laimod_max(j,i) = max(s)

   539   end do

   540   end do

   542 ; print (min(laimod_max)+"/"+max(laimod_max))

   543   delete (s)

   544   delete (dsizes_z)

   545 ;------------------------

   546   DATA11_1D = ndtooned(classob)

   547   DATA12_1D = ndtooned(laiob_max)

   548   DATA22_1D = ndtooned(laimod_max)

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

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

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

   554   do nd=0,1

   556 ; See if we are doing model or observational data.

   558     if(nd.eq.0) then

   559       data_ob  = DATA11_1D

   560       data_mod = DATA12_1D

   561     else

   562       data_ob  = DATA11_1D

   563       data_mod = DATA22_1D

   564     end if

   566 ; Loop through each range and check for values.

   568     do i=0,nr-2

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

   570 ;        print("")

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

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

   573       else

   574 ;        print("")

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

   576         idx = ind(range(i).le.data_ob)

   577       end if

   579 ; Calculate average, and get min and max.

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

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

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

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

   585         count = dimsizes(idx)

   586       else

   587         count            = 0

   588         yvalues(nd,i)    = yvalues@_FillValue

   589         mn_yvalues(nd,i) = yvalues@_FillValue

   590         mx_yvalues(nd,i) = yvalues@_FillValue

   591       end if

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

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

   596 ; Clean up for next time in loop.

   598       delete(idx)

   599     end do

   600     delete(data_ob)

   601     delete(data_mod)

   602   end do

   603 ;-----------------------------------------------------------------

   604 ; compute correlation coef and M score

   606   u = yvalues(0,:)

   607   v = yvalues(1,:)

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

   610   uu = u(good)

   611   vv = v(good)

   613   ccrMax = esccr(uu,vv,0)

   614 ; print (ccrMax)

   616 ; new eq

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

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

   620   M_lai_max = sprintf("%.2f", Mmax)

   621   system("sed s#M_lai_max#"+M_lai_max+"# table.html > table.html.new")

   622   system("mv -f table.html.new table.html")

   623   print (M_lai_max)

   625  do i=0,nrow-2

   626   text4(i,3) = sprintf("%5.2f",u(i))

   627   text4(i,4) = sprintf("%5.2f",v(i))

   628   text4(i,5) = "-"

   629  end do

   630   text4(nrow-1,3) = sprintf("%5.2f",avg(u))

   631   text4(nrow-1,4) = sprintf("%5.2f",avg(v))

   632   text4(nrow-1,5) = sprintf("%5.2f",Mmax)

   634  delete (u)

   635  delete (v)

   636  delete (uu)

   637  delete (vv)

   638 ;--------------------------------------------------------------------

   639 ; histogram res

   641   resm                = True

   642   resm@gsnMaximize    = True

   643   resm@gsnDraw        = False

   644   resm@gsnFrame       = False

   645   resm@xyMarkLineMode = "Markers"

   646   resm@xyMarkerSizeF  = 0.014

   647   resm@xyMarker       = 16

   648   resm@xyMarkerColors = (/"Brown","Blue"/)

   649 ; resm@trYMinF        = min(mn_yvalues) - 10.

   650 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   651   resm@trYMinF        = min(mn_yvalues) - 2

   652   resm@trYMaxF        = max(mx_yvalues) + 4

   654   resm@tiYAxisString  = "Max LAI (Leaf Area Index)"

   655   resm@tiXAxisString  = "Land Cover Type"

   657   max_bar = new((/2,nx/),graphic)

   658   min_bar = new((/2,nx/),graphic)

   659   max_cap = new((/2,nx/),graphic)

   660   min_cap = new((/2,nx/),graphic)

   662   lnres = True

   663   line_colors = (/"brown","blue"/)

   664 ;------------------------------------------------------------------

   665 ; Start the graphics.

   667   plot_name = "histogram_max"

   668   title     = model_name + " vs Observed"

   669   resm@tiMainString  = title

   671   wks   = gsn_open_wks (plot_type,plot_name)

   672 ;-----------------------------

   673 ; Add a boxed legend using the more simple method

   675   resm@pmLegendDisplayMode    = "Always"

   676 ; resm@pmLegendWidthF         = 0.1

   677   resm@pmLegendWidthF         = 0.08

   678   resm@pmLegendHeightF        = 0.05

   679   resm@pmLegendOrthogonalPosF = -1.17

   680 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

   681 ; resm@pmLegendParallelPosF   =  0.18

   682   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   684 ; resm@lgPerimOn              = False

   685   resm@lgLabelFontHeightF     = 0.015

   686   resm@xyExplicitLegendLabels = (/"observed",model_name/)

   687 ;-----------------------------

   688   tRes  = True

   689   tRes@txFontHeightF = 0.025

   691   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"

   693   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   695   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

   696 ;-------------------------------

   697 ;Attach the vertical bar and the horizontal cap line

   699   do nd=0,1

   700     lnres@gsLineColor = line_colors(nd)

   701     do i=0,nx-1

   703       if(.not.ismissing(mn_yvalues(nd,i)).and. \

   704          .not.ismissing(mx_yvalues(nd,i))) then

   706 ; Attach the vertical bar, both above and below the marker.

   708         x1 = xvalues(nd,i)

   709         y1 = yvalues(nd,i)

   710         y2 = mn_yvalues(nd,i)

   711         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   713         y2 = mx_yvalues(nd,i)

   714         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   716 ; Attach the horizontal cap line, both above and below the marker.

   718         x1 = xvalues(nd,i) - dx4

   719         x2 = xvalues(nd,i) + dx4

   720         y1 = mn_yvalues(nd,i)

   721         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   723         y1 = mx_yvalues(nd,i)

   724         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   725       end if

   726     end do

   727   end do

   729   draw(xy)

   730   frame(wks)

   732   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   733   system("rm "+plot_name+"."+plot_type)

   734 ; system("rm "+plot_name+"-1."+plot_type_new)

   735 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   737   clear (wks)

   739  delete (DATA11_1D)

   740  delete (DATA12_1D)

   741  delete (DATA22_1D)

   742 ;delete (range)

   743 ;delete (xvalues)

   744  delete (yvalues)

   745  delete (mn_yvalues)

   746  delete (mx_yvalues)

   747  delete (good)

   748  delete (max_bar)

   749  delete (min_bar)

   750  delete (max_cap)

   751  delete (min_cap)

   752 ;-----------------------------------------------------------------

   753 ;global res

   755   resg                     = True             ; Use plot options

   756   resg@cnFillOn            = True             ; Turn on color fill

   757   resg@gsnSpreadColors     = True             ; use full colormap

   758 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

   759 ; resg@lbLabelAutoStride   = True

   760   resg@cnLinesOn           = False            ; Turn off contourn lines

   761   resg@mpFillOn            = False            ; Turn off map fill

   763   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   764   resg@cnMinLevelValF       = 0.              ; Min level

   765   resg@cnMaxLevelValF       = 10.             ; Max level

   766   resg@cnLevelSpacingF      = 1.              ; interval

   768 ;global contour ob

   770   delta = 0.00001

   771   laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)

   773   plot_name = "global_max_ob"

   774   title     = ob_name

   775   resg@tiMainString  = title

   777   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   778   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   780   plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg)

   781   frame(wks)

   783   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   784   system("rm "+plot_name+"."+plot_type)

   785   system("rm "+plot_name+"-1."+plot_type_new)

   786   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   788   clear (wks)

   789 ;------------------------------------------------------------------------

   790 ;global contour model

   792   plot_name = "global_max_model"

   793   title     = "Model " + model_name

   794   resg@tiMainString  = title

   796   wks = gsn_open_wks (plot_type,plot_name)

   797   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   799   plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg)

   800   frame(wks)

   802   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   803   system("rm "+plot_name+"."+plot_type)

   804   system("rm "+plot_name+"-1."+plot_type_new)

   805   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   807   clear (wks)

   808 ;------------------------------------------------------------------------

   809 ;(C) phase

   810 ;--------------------------------------------------------------------

   811 ; get data

   813 ; observed

   814   laiob_phase = laiob(0,:,:)

   815   s           = laiob(:,0,0)

   816   laiob_phase@long_name = "Leaf Area Index Max Month"

   818   dsizes_z = dimsizes(laiob)

   819   nlat     = dsizes_z(1)

   820   nlon     = dsizes_z(2)

   822   do j = 0,nlat-1

   823   do i = 0,nlon-1

   824      s = laiob(:,j,i)

   825      laiob_phase(j,i) = maxind(s) + 1

   826   end do

   827   end do

   829 ; print (min(laiob_phase)+"/"+max(laiob_phase))

   830   delete (s)

   831   delete (dsizes_z)

   832 ;-------------------------

   833 ; model

   834   laimod_phase = laimod(0,:,:)

   835   s            = laimod(:,0,0)

   836   laimod_phase@long_name = "Leaf Area Index Max Month"

   838   dsizes_z = dimsizes(laimod)

   839   nlat     = dsizes_z(1)

   840   nlon     = dsizes_z(2)

   842   do j = 0,nlat-1

   843   do i = 0,nlon-1

   844      s = laimod(:,j,i)

   845      laimod_phase(j,i) = maxind(s) + 1

   846   end do

   847   end do

   849 ; print (min(laimod_phase)+"/"+max(laimod_phase))

   850   delete (s)

   851   delete (dsizes_z)

   852 ;------------------------

   853   DATA11_1D = ndtooned(classob)

   854   DATA12_1D = ndtooned(laiob_phase)

   855   DATA22_1D = ndtooned(laimod_phase)

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

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

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

   861   do nd=0,1

   863 ; See if we are doing model or observational data.

   865     if(nd.eq.0) then

   866       data_ob  = DATA11_1D

   867       data_mod = DATA12_1D

   868     else

   869       data_ob  = DATA11_1D

   870       data_mod = DATA22_1D

   871     end if

   873 ; Loop through each range and check for values.

   875     do i=0,nr-2

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

   877 ;        print("")

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

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

   880       else

   881 ;        print("")

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

   883         idx = ind(range(i).le.data_ob)

   884       end if

   886 ; Calculate average, and get min and max.

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

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

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

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

   892         count = dimsizes(idx)

   893       else

   894         count            = 0

   895         yvalues(nd,i)    = yvalues@_FillValue

   896         mn_yvalues(nd,i) = yvalues@_FillValue

   897         mx_yvalues(nd,i) = yvalues@_FillValue

   898       end if

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

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

   903 ; Clean up for next time in loop.

   905       delete(idx)

   906     end do

   907     delete(data_ob)

   908     delete(data_mod)

   909   end do

   910 ;-----------------------------------------------------------------

   911 ; compute correlation coef and M score

   913   u = yvalues(0,:)

   914   v = yvalues(1,:)

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

   917   uu = u(good)

   918   vv = v(good)

   920   ccrPhase = esccr(uu,vv,0)

   921 ; print (ccrPhase)

   923 ; old eq

   924 ; bias   = abs(avg(vv)-avg(uu))

   925 ; new eq

   926   bias   = avg(abs(vv-uu))

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

   929   Mphase = ((6. - bias)/6.)*5.

   931   M_lai_phase = sprintf("%.2f", Mphase)

   932   system("sed s#M_lai_phase#"+M_lai_phase+"# table.html > table.html.new")

   933   system("mv -f table.html.new table.html")

   934   print (M_lai_phase)

   936  do i=0,nrow-2

   937   text4(i,6) = sprintf("%5.2f",u(i))

   938   text4(i,7) = sprintf("%5.2f",v(i))

   939   text4(i,8) = "-"

   940  end do

   941   text4(nrow-1,6) = sprintf("%5.2f",avg(u))

   942   text4(nrow-1,7) = sprintf("%5.2f",avg(v))

   943   text4(nrow-1,8) = sprintf("%5.2f",Mphase)

   945  delete (u)

   946  delete (v)

   947  delete (uu)

   948  delete (vv)

   949 ;--------------------------------------------------------------------

   950 ; histogram res

   952   resm                = True

   953   resm@gsnMaximize    = True

   954   resm@gsnDraw        = False

   955   resm@gsnFrame       = False

   956   resm@xyMarkLineMode = "Markers"

   957   resm@xyMarkerSizeF  = 0.014

   958   resm@xyMarker       = 16

   959   resm@xyMarkerColors = (/"Brown","Blue"/)

   960 ; resm@trYMinF        = min(mn_yvalues) - 10.

   961 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   962   resm@trYMinF        = min(mn_yvalues) - 2

   963   resm@trYMaxF        = max(mx_yvalues) + 4

   965   resm@tiYAxisString  = "Max LAI (Leaf Area Index) Month"

   966   resm@tiXAxisString  = "Land Cover Type"

   968   max_bar = new((/2,nx/),graphic)

   969   min_bar = new((/2,nx/),graphic)

   970   max_cap = new((/2,nx/),graphic)

   971   min_cap = new((/2,nx/),graphic)

   973   lnres = True

   974   line_colors = (/"brown","blue"/)

   975 ;------------------------------------------------------------------

   976 ; Start the graphics.

   978   plot_name = "histogram_phase"

   979   title     = model_name + " vs Observed"

   980   resm@tiMainString  = title

   982   wks   = gsn_open_wks (plot_type,plot_name)

   983 ;-----------------------------

   984 ; Add a boxed legend using the more simple method

   986   resm@pmLegendDisplayMode    = "Always"

   987 ; resm@pmLegendWidthF         = 0.1

   988   resm@pmLegendWidthF         = 0.08

   989   resm@pmLegendHeightF        = 0.05

   990   resm@pmLegendOrthogonalPosF = -1.17

   991 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

   992 ; resm@pmLegendParallelPosF   =  0.18

   993   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   995 ; resm@lgPerimOn              = False

   996   resm@lgLabelFontHeightF     = 0.015

   997   resm@xyExplicitLegendLabels = (/"observed",model_name/)

   998 ;-----------------------------

   999   tRes  = True

  1000   tRes@txFontHeightF = 0.025

  1002   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"

  1004   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

  1006   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

  1007 ;-------------------------------

  1008 ;Attach the vertical bar and the horizontal cap line

  1010   do nd=0,1

  1011     lnres@gsLineColor = line_colors(nd)

  1012     do i=0,nx-1

  1014       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1015          .not.ismissing(mx_yvalues(nd,i))) then

  1017 ; Attach the vertical bar, both above and below the marker.

  1019         x1 = xvalues(nd,i)

  1020         y1 = yvalues(nd,i)

  1021         y2 = mn_yvalues(nd,i)

  1022         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1024         y2 = mx_yvalues(nd,i)

  1025         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1027 ; Attach the horizontal cap line, both above and below the marker.

  1029         x1 = xvalues(nd,i) - dx4

  1030         x2 = xvalues(nd,i) + dx4

  1031         y1 = mn_yvalues(nd,i)

  1032         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1034         y1 = mx_yvalues(nd,i)

  1035         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1036       end if

  1037     end do

  1038   end do

  1040   draw(xy)

  1041   frame(wks)

  1043   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1044   system("rm "+plot_name+"."+plot_type)

  1045 ; system("rm "+plot_name+"-1."+plot_type_new)

  1046 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1048   clear (wks)

  1050  delete (DATA11_1D)

  1051  delete (DATA12_1D)

  1052  delete (DATA22_1D)

  1053 ;delete (range)

  1054 ;delete (xvalues)

  1055  delete (yvalues)

  1056  delete (mn_yvalues)

  1057  delete (mx_yvalues)

  1058  delete (good)

  1059  delete (max_bar)

  1060  delete (min_bar)

  1061  delete (max_cap)

  1062  delete (min_cap)

  1063 ;-----------------------------------------------------------------

  1064 ;global res

  1066   resg                     = True             ; Use plot options

  1067   resg@cnFillOn            = True             ; Turn on color fill

  1068   resg@gsnSpreadColors     = True             ; use full colormap

  1069 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

  1070 ; resg@lbLabelAutoStride   = True

  1071   resg@cnLinesOn           = False            ; Turn off contourn lines

  1072   resg@mpFillOn            = False            ; Turn off map fill

  1074   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1075   resg@cnMinLevelValF       = 1.              ; Min level

  1076   resg@cnMaxLevelValF       = 12.             ; Max level

  1077   resg@cnLevelSpacingF      = 1.              ; interval

  1079 ;global contour ob

  1081   delta = 0.00001

  1082   laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)

  1084   plot_name = "global_phase_ob"

  1085   title     = ob_name

  1086   resg@tiMainString  = title

  1088   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1089   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1091   plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)

  1092   frame(wks)

  1094   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1095   system("rm "+plot_name+"."+plot_type)

  1096   system("rm "+plot_name+"-1."+plot_type_new)

  1097   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1099   clear (wks)

  1100 ;------------------------------------------------------------------------

  1101 ;global contour model

  1103   plot_name = "global_phase_model"

  1104   title     = "Model " + model_name

  1105   resg@tiMainString  = title

  1107   wks = gsn_open_wks (plot_type,plot_name)

  1108   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1110   delete (plot)

  1111   plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)

  1112   frame(wks)

  1114   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1115   system("rm "+plot_name+"."+plot_type)

  1116   system("rm "+plot_name+"-1."+plot_type_new)

  1117   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1119   clear (wks)

  1120 ;-----------------------------------------------------------------

  1121 ;(D) grow day

  1122 ;--------------------------------------------------------------------

  1123 ; get data

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

  1127 ; observed

  1128   laiob_grow = laiob(0,:,:)

  1129   laiob_grow@long_name = "Days of Growing Season"

  1131   dsizes_z = dimsizes(laiob)

  1132   ntime    = dsizes_z(0)

  1133   nlat     = dsizes_z(1)

  1134   nlon     = dsizes_z(2)

  1136   do j = 0,nlat-1

  1137   do i = 0,nlon-1

  1138      nday = 0.

  1139      do k = 0,ntime-1

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

  1141            nday = nday + day_of_data(k)

  1142         end if

  1143      end do

  1145      laiob_grow(j,i) = nday

  1146   end do

  1147   end do

  1149 ; print (min(laiob_grow)+"/"+max(laiob_grow))

  1150 ;-------------------------

  1151 ; model

  1152   laimod_grow = laimod(0,:,:)

  1153   laimod_grow@long_name = "Days of Growing Season"

  1155   dsizes_z = dimsizes(laimod)

  1156   ntime    = dsizes_z(0)

  1157   nlat     = dsizes_z(1)

  1158   nlon     = dsizes_z(2)

  1160   do j = 0,nlat-1

  1161   do i = 0,nlon-1

  1162      nday = 0.

  1163      do k = 0,ntime-1

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

  1165            nday = nday + day_of_data(k)

  1166         end if

  1167      end do

  1169      laimod_grow(j,i) = nday

  1170   end do

  1171   end do

  1173 ; print (min(laimod_grow)+"/"+max(laimod_grow))

  1174 ;------------------------

  1175   DATA11_1D = ndtooned(classob)

  1176   DATA12_1D = ndtooned(laiob_grow)

  1177   DATA22_1D = ndtooned(laimod_grow)

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

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

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

  1183   do nd=0,1

  1185 ; See if we are doing model or observational data.

  1187     if(nd.eq.0) then

  1188       data_ob  = DATA11_1D

  1189       data_mod = DATA12_1D

  1190     else

  1191       data_ob  = DATA11_1D

  1192       data_mod = DATA22_1D

  1193     end if

  1195 ; Loop through each range and check for values.

  1197     do i=0,nr-2

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

  1199 ;        print("")

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

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

  1202       else

  1203 ;        print("")

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

  1205         idx = ind(range(i).le.data_ob)

  1206       end if

  1208 ; Calculate average, and get min and max.

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

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

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

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

  1214         count = dimsizes(idx)

  1215       else

  1216         count            = 0

  1217         yvalues(nd,i)    = yvalues@_FillValue

  1218         mn_yvalues(nd,i) = yvalues@_FillValue

  1219         mx_yvalues(nd,i) = yvalues@_FillValue

  1220       end if

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

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

  1225 ; Clean up for next time in loop.

  1227       delete(idx)

  1228     end do

  1229     delete(data_ob)

  1230     delete(data_mod)

  1231   end do

  1232 ;-----------------------------------------------------------------

  1233 ; compute correlation coef and M score

  1235   u = yvalues(0,:)

  1236   v = yvalues(1,:)

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

  1239   uu = u(good)

  1240   vv = v(good)

  1242   ccrGrow = esccr(uu,vv,0)

  1243 ; print (ccrGrow)

  1245 ; new eq

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

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

  1249   M_lai_grow = sprintf("%.2f", Mgrow)

  1250   system("sed s#M_lai_grow#"+M_lai_grow+"# table.html > table.html.new")

  1251   system("mv -f table.html.new table.html")

  1252   print (M_lai_grow)

  1254  do i=0,nrow-2

  1255   text4(i,9)  = sprintf("%5.2f",u(i))

  1256   text4(i,10) = sprintf("%5.2f",v(i))

  1257   text4(i,11) = "-"

  1258  end do

  1259   text4(nrow-1,9)  = sprintf("%5.2f",avg(u))

  1260   text4(nrow-1,10) = sprintf("%5.2f",avg(v))

  1261   text4(nrow-1,11) = sprintf("%5.2f",Mgrow)

  1263  delete (u)

  1264  delete (v)

  1265  delete (uu)

  1266  delete (vv)

  1267 ;--------------------------------------------------------------------

  1268 ; histogram res

  1270   resm                = True

  1271   resm@gsnMaximize    = True

  1272   resm@gsnDraw        = False

  1273   resm@gsnFrame       = False

  1274   resm@xyMarkLineMode = "Markers"

  1275   resm@xyMarkerSizeF  = 0.014

  1276   resm@xyMarker       = 16

  1277   resm@xyMarkerColors = (/"Brown","Blue"/)

  1278 ; resm@trYMinF        = min(mn_yvalues) - 10.

  1279 ; resm@trYMaxF        = max(mx_yvalues) + 10.

  1280   resm@trYMinF        = min(mn_yvalues) - 2

  1281   resm@trYMaxF        = max(mx_yvalues) + 4

  1283   resm@tiYAxisString = "Days of Growing season"

  1284   resm@tiXAxisString = "Land Cover Type"

  1286   max_bar = new((/2,nx/),graphic)

  1287   min_bar = new((/2,nx/),graphic)

  1288   max_cap = new((/2,nx/),graphic)

  1289   min_cap = new((/2,nx/),graphic)

  1291   lnres = True

  1292   line_colors = (/"brown","blue"/)

  1293 ;------------------------------------------------------------------

  1294 ; Start the graphics.

  1296   plot_name = "histogram_grow"

  1297   title     = model_name + " vs Observed"

  1298   resm@tiMainString  = title

  1300   wks   = gsn_open_wks (plot_type,plot_name)

  1301 ;-----------------------------

  1302 ; Add a boxed legend using the more simple method

  1304   resm@pmLegendDisplayMode    = "Always"

  1305 ; resm@pmLegendWidthF         = 0.1

  1306   resm@pmLegendWidthF         = 0.08

  1307   resm@pmLegendHeightF        = 0.05

  1308   resm@pmLegendOrthogonalPosF = -1.17

  1309 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

  1310 ; resm@pmLegendParallelPosF   =  0.18

  1311   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

  1313 ; resm@lgPerimOn              = False

  1314   resm@lgLabelFontHeightF     = 0.015

  1315   resm@xyExplicitLegendLabels = (/"observed",model_name/)

  1316 ;-----------------------------

  1317   tRes  = True

  1318   tRes@txFontHeightF = 0.025

  1320   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"

  1322   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

  1324   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

  1325 ;-------------------------------

  1326 ;Attach the vertical bar and the horizontal cap line

  1328   do nd=0,1

  1329     lnres@gsLineColor = line_colors(nd)

  1330     do i=0,nx-1

  1332       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1333          .not.ismissing(mx_yvalues(nd,i))) then

  1335 ; Attach the vertical bar, both above and below the marker.

  1337         x1 = xvalues(nd,i)

  1338         y1 = yvalues(nd,i)

  1339         y2 = mn_yvalues(nd,i)

  1340         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1342         y2 = mx_yvalues(nd,i)

  1343         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1345 ; Attach the horizontal cap line, both above and below the marker.

  1347         x1 = xvalues(nd,i) - dx4

  1348         x2 = xvalues(nd,i) + dx4

  1349         y1 = mn_yvalues(nd,i)

  1350         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1352         y1 = mx_yvalues(nd,i)

  1353         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1354       end if

  1355     end do

  1356   end do

  1358   draw(xy)

  1359   frame(wks)

  1361   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1362   system("rm "+plot_name+"."+plot_type)

  1363 ; system("rm "+plot_name+"-1."+plot_type_new)

  1364 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1366   clear (wks)

  1368  delete (DATA11_1D)

  1369  delete (DATA12_1D)

  1370  delete (DATA22_1D)

  1371 ;delete (range)

  1372 ;delete (xvalues)

  1373  delete (yvalues)

  1374  delete (mn_yvalues)

  1375  delete (mx_yvalues)

  1376  delete (good)

  1377  delete (max_bar)

  1378  delete (min_bar)

  1379  delete (max_cap)

  1380  delete (min_cap)

  1381 ;-----------------------------------------------------------------

  1382 ;global res

  1384   resg                     = True             ; Use plot options

  1385   resg@cnFillOn            = True             ; Turn on color fill

  1386   resg@gsnSpreadColors     = True             ; use full colormap

  1387 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

  1388 ; resg@lbLabelAutoStride   = True

  1389   resg@cnLinesOn           = False            ; Turn off contourn lines

  1390   resg@mpFillOn            = False            ; Turn off map fill

  1392   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1393   resg@cnMinLevelValF       = 60.             ; Min level

  1394   resg@cnMaxLevelValF       = 360.            ; Max level

  1395   resg@cnLevelSpacingF      = 20.             ; interval

  1397 ;global contour ob

  1399   laiob_grow@_FillValue = 1.e+36

  1400   laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow)

  1402   plot_name = "global_grow_ob"

  1403   title     = ob_name

  1404   resg@tiMainString  = title

  1406   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1407   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1409   plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)

  1410   frame(wks)

  1412   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1413   system("rm "+plot_name+"."+plot_type)

  1414   system("rm "+plot_name+"-1."+plot_type_new)

  1415   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1417   clear (wks)

  1418 ;------------------------------------------------------------------------

  1419 ;global contour model

  1421   laimod_grow@_FillValue = 1.e+36

  1422   laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow)

  1424   plot_name = "global_grow_model"

  1425   title     = "Model " + model_name

  1426   resg@tiMainString  = title

  1428   wks = gsn_open_wks (plot_type,plot_name)

  1429   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1431   delete (plot)

  1432   plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)

  1433   frame(wks)

  1435   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1436   system("rm "+plot_name+"."+plot_type)

  1437   system("rm "+plot_name+"-1."+plot_type_new)

  1438   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1440   clear (wks)

  1441 ;------------------------------------------------------------------------

  1442 ;global contour model vs ob

  1444   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1445   resg@cnMinLevelValF       = 0.              ; Min level

  1446   resg@cnMaxLevelValF       = 10.             ; Max level

  1447   resg@cnLevelSpacingF      = 1.              ; interval

  1449   plot_name = "global_mean_model_vs_ob"

  1451   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1452   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1454   delete (plot)

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

  1457   resg@gsnFrame             = False          ; Do not draw plot

  1458   resg@gsnDraw              = False          ; Do not advance frame

  1460 ; plot correlation coef

  1462   gRes               = True

  1463   gRes@txFontHeightF = 0.02

  1464   gRes@txAngleF      = 90

  1466   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"

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

  1469 ;--------------------------------------------------------------------

  1471 ;(a) ob

  1473   title     = ob_name

  1474   resg@tiMainString  = title

  1476   plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg)

  1478 ;(b) model

  1480   title     = "Model "+ model_name

  1481   resg@tiMainString  = title

  1483   plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg)

  1485 ;(c) model-ob

  1487   zz = laimod_mean

  1488   zz = laimod_mean - laiob_mean

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

  1490   resg@tiMainString    = title

  1492   resg@cnMinLevelValF  = -2.           ; Min level

  1493   resg@cnMaxLevelValF  =  2.           ; Max level

  1494   resg@cnLevelSpacingF =  0.4           ; interval

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

  1499   pres                            = True        ; panel plot mods desired

  1500   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1501                                                 ; indiv. plots in panel

  1502   pres@gsnMaximize                = True        ; fill the page

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

  1506   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1507   system("rm "+plot_name+"."+plot_type)

  1508 ; system("rm "+plot_name+"-1."+plot_type_new)

  1509 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1511   frame (wks)

  1512   clear (wks)

  1514   delete (plot)

  1515 ;-----------------------------------------------------------------

  1516 ;global contour model vs ob

  1518   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1519   resg@cnMinLevelValF       = 0.              ; Min level

  1520   resg@cnMaxLevelValF       = 10.             ; Max level

  1521   resg@cnLevelSpacingF      = 1.              ; interval

  1523   plot_name = "global_max_model_vs_ob"

  1525   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1526   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1528 ; delete(plot)

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

  1531   resg@gsnFrame             = False          ; Do not draw plot

  1532   resg@gsnDraw              = False          ; Do not advance frame

  1534 ; plot correlation coef

  1536   gRes               = True

  1537   gRes@txFontHeightF = 0.02

  1538   gRes@txAngleF      = 90

  1540   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"

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

  1543 ;--------------------------------------------------------------------

  1544 ;(a) ob

  1546   title     = ob_name

  1547   resg@tiMainString  = title

  1549   plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg)

  1551 ;(b) model

  1553   title     = "Model "+ model_name

  1554   resg@tiMainString  = title

  1556   plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg)

  1558 ;(c) model-ob

  1560   delete (zz)

  1561   zz = laimod_max

  1562   zz = laimod_max - laiob_max

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

  1564   resg@tiMainString    = title

  1566   resg@cnMinLevelValF  = -6.           ; Min level

  1567   resg@cnMaxLevelValF  =  6.           ; Max level

  1568   resg@cnLevelSpacingF =  1.           ; interval

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

  1573   pres                            = True        ; panel plot mods desired

  1574   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1575                                                 ; indiv. plots in panel

  1576   pres@gsnMaximize                = True        ; fill the page

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

  1580   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1581   system("rm "+plot_name+"."+plot_type)

  1582 ; system("rm "+plot_name+"-1."+plot_type_new)

  1583 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1585   frame (wks)

  1586   clear (wks)

  1588   delete (plot)

  1589 ;-----------------------------------------------------------------

  1590 ;global contour model vs ob

  1592   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1593   resg@cnMinLevelValF       = 1.              ; Min level

  1594   resg@cnMaxLevelValF       = 12.             ; Max level

  1595   resg@cnLevelSpacingF      = 1.              ; interval

  1597   plot_name = "global_phase_model_vs_ob"

  1599   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1600   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1602 ; delete (plot)

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

  1605   resg@gsnFrame             = False          ; Do not draw plot

  1606   resg@gsnDraw              = False          ; Do not advance frame

  1608 ; plot correlation coef

  1610   gRes               = True

  1611   gRes@txFontHeightF = 0.02

  1612   gRes@txAngleF      = 90

  1614   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"

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

  1617 ;--------------------------------------------------------------------

  1618 ;(a) ob

  1620   title     = ob_name

  1621   resg@tiMainString  = title

  1623   plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg)

  1625 ;(b) model

  1627   title     = "Model "+ model_name

  1628   resg@tiMainString  = title

  1630   plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg)

  1632 ;(c) model-ob

  1634   delete (zz)

  1635   zz = laimod_phase

  1636   zz = laimod_phase - laiob_phase

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

  1638   resg@tiMainString    = title

  1640   resg@cnMinLevelValF  = -6.           ; Min level

  1641   resg@cnMaxLevelValF  =  6.           ; Max level

  1642   resg@cnLevelSpacingF =  1.           ; interval

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

  1647 ; pres                            = True        ; panel plot mods desired

  1648 ; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1649                                                 ; indiv. plots in panel

  1650 ; pres@gsnMaximize                = True        ; fill the page

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

  1654   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1655   system("rm "+plot_name+"."+plot_type)

  1656 ; system("rm "+plot_name+"-1."+plot_type_new)

  1657 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1659   frame (wks)

  1660   clear (wks)

  1662   delete (plot)

  1663 ;------------------------------------------------------------------

  1664 ;global contour model vs ob

  1666   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1667   resg@cnMinLevelValF       = 60.             ; Min level

  1668   resg@cnMaxLevelValF       = 360.            ; Max level

  1669   resg@cnLevelSpacingF      = 20.             ; interval

  1671   plot_name = "global_grow_model_vs_ob"

  1673   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1674   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1676 ; delete (plot)

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

  1679   resg@gsnFrame             = False          ; Do not draw plot

  1680   resg@gsnDraw              = False          ; Do not advance frame

  1682 ; plot correlation coef

  1684   gRes               = True

  1685   gRes@txFontHeightF = 0.02

  1686   gRes@txAngleF      = 90

  1688   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"

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

  1691 ;--------------------------------------------------------------------

  1692 ;(a) ob

  1694   title     = ob_name

  1695   resg@tiMainString  = title

  1697   plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg)

  1699 ;(b) model

  1701   title     = "Model "+ model_name

  1702   resg@tiMainString  = title

  1704   plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg)

  1706 ;(c) model-ob

  1708   delete (zz)

  1709   zz = laimod_grow

  1710   zz = laimod_grow - laiob_grow

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

  1712   resg@tiMainString    = title

  1714   resg@cnMinLevelValF  = -100.           ; Min level

  1715   resg@cnMaxLevelValF  =  100.           ; Max level

  1716   resg@cnLevelSpacingF =  20.            ; interval

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

  1720   pres                            = True        ; panel plot mods desired

  1721   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1722                                                 ; indiv. plots in panel

  1723   pres@gsnMaximize                = True        ; fill the page

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

  1727   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1728   system("rm "+plot_name+"."+plot_type)

  1729 ; system("rm "+plot_name+"-1."+plot_type_new)

  1730 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1732   frame (wks)

  1733   clear (wks)

  1735   delete (plot)

  1736 ;**************************************************

  1737 ; plot lai table

  1738 ;**************************************************

  1740   plot_name = "table_lai"

  1741   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1743   gsn_table(wks,ncr1,xx1,yy1,text1,res1)

  1744   gsn_table(wks,ncr21,xx21,yy21,text21,res21)

  1745   gsn_table(wks,ncr22,xx22,yy22,text22,res22)

  1746   gsn_table(wks,ncr3,xx3,yy3,text3,res3)

  1747   gsn_table(wks,ncr4,xx4,yy4,text4,res4)

  1749   frame(wks)

  1751   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1752 ; system("rm "+plot_name+"."+plot_type)

  1753 ; system("rm "+plot_name+"-1."+plot_type_new)

  1754 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1755 ;-------------------------------------------------------------------

  1756   temp_name = "lai." + model_name

  1757   system("mkdir -p " + temp_name)

  1758   system("cp table.html " + temp_name)

  1759   system("mv *.png " + temp_name)

  1760   system("tar cf "+ temp_name +".tar " + temp_name)

  1761 ;-------------------------------------------------------------------

  1762 end

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