clamp: lai/99.all.ncl.0@0c6405ab2ff4

     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"

    50 ;------------------------------------------------

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

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

    54  laimod  = fm->TLAI

    56 ;************************************************

    57 ; read in data: observed

    58 ;************************************************

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

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

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

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

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

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

    69  classob    = tofloat(fo1->LAND_CLASS)

    70  laiob      = fo2->LAI

    72 ;*******************************************************************

    73 ; for plotting table

    74 ;*******************************************************************

    75 ; table header name

    76   table_header_name = "LAI"

    78 ; column (not including header column)

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

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

    82                  ,"ob","model","M" \

    83                  ,"ob","model","M" \

    84                  ,"ob","model","M" \

    85                  /)

    87   ncol1 = dimsizes(col_header1)

    88   ncol2 = dimsizes(col_header2)

    89   ncol  = ncol2

    91 ; row (not including header row)

    92   row_header = (/"Water Bodies" \

    93                 ,"Evergreen Needleleaf Forests" \

    94                 ,"Evergreen Broadleaf Forests" \

    95                 ,"Deciduous Needleleaf Forest" \

    96                 ,"Deciduous Broadleaf Forests" \

    97                 ,"Mixed Forests" \

    98                 ,"Closed Bushlands" \

    99                 ,"Open Bushlands" \

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

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

   102                 ,"Grasslands" \

   103                 ,"Permanent Wetlands" \

   104                 ,"Croplands" \

   105                 ,"Urban and Built-Up" \

   106                 ,"Cropland/Natural Vegetation Mosaic" \

   107                 ,"Permanent Snow and Ice" \

   108                 ,"Barren or Sparsely Vegetated" \

   109                 ,"Unclassified" \

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

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

   112                 ,"All biome average" \

   113                 /)

   114   nrow = dimsizes(row_header)

   116 ; arrays to be passed to table.

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

   119   table_length = 0.995

   121 ; Table header

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

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

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

   125   text1 = table_header_name

   126   res1               = True

   127   res1@txFontHeightF = 0.03

   128   res1@gsFillColor   = "CornFlowerBlue"

   130 ; Column header (equally space in x2)

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

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

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

   134   text21 = col_header1

   135   res21               = True

   136   res21@txFontHeightF = 0.015

   137   res21@gsFillColor   = "Gray"

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

   140   xx22    = xx21                    ; start from end of x1

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

   142   text22 = col_header2

   143   res22               = True

   144   res22@txFontHeightF = 0.015

   145   res22@gsFillColor   = "Gray"

   147 ; Row header (equally space in y2)

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

   149   xx3    = xx1                      ; same as x1

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

   151   text3 = row_header

   152   res3               = True

   153   res3@txFontHeightF = 0.01

   154   res3@gsFillColor   = "Gray"

   156 ; Main table body

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

   158   xx4    = xx21                      ; Start and end x

   159   yy4    = yy3                      ; Start and end Y

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

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

   163   color_fill4      = "white"

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

   166   res4               = True       ; Set up resource list

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

   168   res4@txFontHeightF = 0.015

   169   res4@gsFillColor   = color_fill4

   171   delete (color_fill4)

   173 ;************************************************

   174 ; plot global land class: observed

   175 ;************************************************

   176 ;global res

   178   resg                     = True             ; Use plot options

   179   resg@cnFillOn            = True             ; Turn on color fill

   180   resg@gsnSpreadColors     = True             ; use full colormap

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

   182 ; resg@lbLabelAutoStride   = True

   183   resg@cnLinesOn           = False            ; Turn off contourn lines

   184   resg@mpFillOn            = False            ; Turn off map fill

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

   187   resg@cnMinLevelValF       = 1.              ; Min level

   188   resg@cnMaxLevelValF       = 19.             ; Max level

   189   resg@cnLevelSpacingF      = 1.              ; interval

   191 ;global contour ob

   192   classob@_FillValue = 1.e+36

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

   195   plot_name = "global_class_ob"

   196   title     = ob_name

   197   resg@tiMainString  = title

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

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

   202   plot = gsn_csm_contour_map_ce(wks,classob,resg)

   203   frame(wks)

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

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

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

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

   210   clear (wks)

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

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

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

   214   nclassn     = nclass + 1

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

   216 ; print (range)

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

   219 ; particular range, and then take an average.

   221   nr           = dimsizes(range)

   222   nx           = nr-1

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

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

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

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

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

   228 ;-----------------------------------------------------------------

   229 ;(A) mean

   230 ;--------------------------------------------------------------------

   231 ; get data

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

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

   236   DATA11_1D = ndtooned(classob)

   237   DATA12_1D = ndtooned(laiob_mean)

   238   DATA22_1D = ndtooned(laimod_mean)

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

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

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

   244   do nd=0,1

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

   248     if(nd.eq.0) then

   249       data_ob  = DATA11_1D

   250       data_mod = DATA12_1D

   251     else

   252       data_ob  = DATA11_1D

   253       data_mod = DATA22_1D

   254     end if

   256 ; Loop through each range and check for values.

   258     do i=0,nr-2

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

   260 ;        print("")

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

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

   263       else

   264 ;        print("")

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

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

   267       end if

   269 ; Calculate average, and get min and max.

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

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

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

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

   275         count = dimsizes(idx)

   276       else

   277         count            = 0

   278         yvalues(nd,i)    = yvalues@_FillValue

   279         mn_yvalues(nd,i) = yvalues@_FillValue

   280         mx_yvalues(nd,i) = yvalues@_FillValue

   281       end if

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

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

   286 ; Clean up for next time in loop.

   288       delete(idx)

   289     end do

   290     delete(data_ob)

   291     delete(data_mod)

   292   end do

   293 ;-----------------------------------------------------------------

   294 ; compute correlation coef and M score

   296   u = yvalues(0,:)

   297   v = yvalues(1,:)

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

   300   uu = u(good)

   301   vv = v(good)

   303   ccrMean = esccr(uu,vv,0)

   304 ; print (ccrMean)

   306 ; new eq

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

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

   310   print (Mmean)

   312  do i=0,nrow-2

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

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

   315   text4(i,2) = "-"

   316  end do

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

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

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

   321  delete (u)

   322  delete (v)

   323  delete (uu)

   324  delete (vv)

   325 ;--------------------------------------------------------------------

   326 ; histogram res

   328   resm                = True

   329   resm@gsnMaximize    = True

   330   resm@gsnDraw        = False

   331   resm@gsnFrame       = False

   332   resm@xyMarkLineMode = "Markers"

   333   resm@xyMarkerSizeF  = 0.014

   334   resm@xyMarker       = 16

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

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

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

   338   resm@trYMinF        = min(mn_yvalues) - 2

   339   resm@trYMaxF        = max(mx_yvalues) + 4

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

   342   resm@tiXAxisString  = "Land Cover Type"

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

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

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

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

   349   lnres = True

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

   351 ;------------------------------------------------------------------

   352 ; Start the graphics.

   354   plot_name = "histogram_mean"

   355   title     = model_name + " vs Observed"

   356   resm@tiMainString  = title

   358   wks   = gsn_open_wks (plot_type,plot_name)

   359 ;-----------------------------

   360 ; Add a boxed legend using the more simple method

   362   resm@pmLegendDisplayMode    = "Always"

   363 ; resm@pmLegendWidthF         = 0.1

   364   resm@pmLegendWidthF         = 0.08

   365   resm@pmLegendHeightF        = 0.05

   366   resm@pmLegendOrthogonalPosF = -1.17

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

   368 ; resm@pmLegendParallelPosF   =  0.18

   369   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   371 ; resm@lgPerimOn              = False

   372   resm@lgLabelFontHeightF     = 0.015

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

   374 ;-----------------------------

   375   tRes  = True

   376   tRes@txFontHeightF = 0.025

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

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

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

   383 ;-------------------------------

   384 ;Attach the vertical bar and the horizontal cap line

   386   do nd=0,1

   387     lnres@gsLineColor = line_colors(nd)

   388     do i=0,nx-1

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

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

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

   395         x1 = xvalues(nd,i)

   396         y1 = yvalues(nd,i)

   397         y2 = mn_yvalues(nd,i)

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

   400         y2 = mx_yvalues(nd,i)

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

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

   405         x1 = xvalues(nd,i) - dx4

   406         x2 = xvalues(nd,i) + dx4

   407         y1 = mn_yvalues(nd,i)

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

   410         y1 = mx_yvalues(nd,i)

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

   412       end if

   413     end do

   414   end do

   416   draw(xy)

   417   frame(wks)

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

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

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

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

   424   clear (wks)

   426  delete (DATA11_1D)

   427  delete (DATA12_1D)

   428  delete (DATA22_1D)

   429 ;delete (range)

   430 ;delete (xvalues)

   431  delete (yvalues)

   432  delete (mn_yvalues)

   433  delete (mx_yvalues)

   434  delete (good)

   435  delete (max_bar)

   436  delete (min_bar)

   437  delete (max_cap)

   438  delete (min_cap)

   439 ;-----------------------------------------------------------------

   440 ;global res

   442   resg                     = True             ; Use plot options

   443   resg@cnFillOn            = True             ; Turn on color fill

   444   resg@gsnSpreadColors     = True             ; use full colormap

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

   446 ; resg@lbLabelAutoStride   = True

   447   resg@cnLinesOn           = False            ; Turn off contourn lines

   448   resg@mpFillOn            = False            ; Turn off map fill

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

   451   resg@cnMinLevelValF       = 0.              ; Min level

   452   resg@cnMaxLevelValF       = 10.             ; Max level

   453   resg@cnLevelSpacingF      = 1.              ; interval

   455 ;global contour ob

   457   delta = 0.00001

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

   460   plot_name = "global_mean_ob"

   461   title     = ob_name

   462   resg@tiMainString  = title

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

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

   467   plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg)

   468   frame(wks)

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

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

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

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

   475   clear (wks)

   476 ;------------------------------------------------------------------------

   477 ;global contour model

   479   plot_name = "global_mean_model"

   480   title     = "Model " + model_name

   481   resg@tiMainString  = title

   483   wks = gsn_open_wks (plot_type,plot_name)

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

   486   plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg)

   487   frame(wks)

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

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

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

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

   494   clear (wks)

   495 ;-----------------------------------------------------------------

   496 ;(B) max

   497 ;--------------------------------------------------------------------

   498 ; get data

   500 ; observed

   501   laiob_max = laiob(0,:,:)

   502   s         = laiob(:,0,0)

   503   laiob_max@long_name = "Leaf Area Index Max"

   505   dsizes_z = dimsizes(laiob)

   506   nlat     = dsizes_z(1)

   507   nlon     = dsizes_z(2)

   509   do j = 0,nlat-1

   510   do i = 0,nlon-1

   511      s = laiob(:,j,i)

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

   513   end do

   514   end do

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

   517   delete (s)

   518   delete (dsizes_z)

   519 ;-------------------------

   520 ; model

   521   laimod_max = laimod(0,:,:)

   522   s          = laimod(:,0,0)

   523   laimod_max@long_name = "Leaf Area Index Max"

   525   dsizes_z = dimsizes(laimod)

   526   nlat     = dsizes_z(1)

   527   nlon     = dsizes_z(2)

   529   do j = 0,nlat-1

   530   do i = 0,nlon-1

   531      s = laimod(:,j,i)

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

   533   end do

   534   end do

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

   537   delete (s)

   538   delete (dsizes_z)

   539 ;------------------------

   540   DATA11_1D = ndtooned(classob)

   541   DATA12_1D = ndtooned(laiob_max)

   542   DATA22_1D = ndtooned(laimod_max)

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

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

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

   548   do nd=0,1

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

   552     if(nd.eq.0) then

   553       data_ob  = DATA11_1D

   554       data_mod = DATA12_1D

   555     else

   556       data_ob  = DATA11_1D

   557       data_mod = DATA22_1D

   558     end if

   560 ; Loop through each range and check for values.

   562     do i=0,nr-2

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

   564 ;        print("")

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

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

   567       else

   568 ;        print("")

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

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

   571       end if

   573 ; Calculate average, and get min and max.

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

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

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

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

   579         count = dimsizes(idx)

   580       else

   581         count            = 0

   582         yvalues(nd,i)    = yvalues@_FillValue

   583         mn_yvalues(nd,i) = yvalues@_FillValue

   584         mx_yvalues(nd,i) = yvalues@_FillValue

   585       end if

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

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

   590 ; Clean up for next time in loop.

   592       delete(idx)

   593     end do

   594     delete(data_ob)

   595     delete(data_mod)

   596   end do

   597 ;-----------------------------------------------------------------

   598 ; compute correlation coef and M score

   600   u = yvalues(0,:)

   601   v = yvalues(1,:)

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

   604   uu = u(good)

   605   vv = v(good)

   607   ccrMax = esccr(uu,vv,0)

   608 ; print (ccrMax)

   610 ; new eq

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

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

   614   print (Mmax)

   616  do i=0,nrow-2

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

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

   619   text4(i,5) = "-"

   620  end do

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

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

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

   625  delete (u)

   626  delete (v)

   627  delete (uu)

   628  delete (vv)

   629 ;--------------------------------------------------------------------

   630 ; histogram res

   632   resm                = True

   633   resm@gsnMaximize    = True

   634   resm@gsnDraw        = False

   635   resm@gsnFrame       = False

   636   resm@xyMarkLineMode = "Markers"

   637   resm@xyMarkerSizeF  = 0.014

   638   resm@xyMarker       = 16

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

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

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

   642   resm@trYMinF        = min(mn_yvalues) - 2

   643   resm@trYMaxF        = max(mx_yvalues) + 4

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

   646   resm@tiXAxisString  = "Land Cover Type"

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

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

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

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

   653   lnres = True

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

   655 ;------------------------------------------------------------------

   656 ; Start the graphics.

   658   plot_name = "histogram_max"

   659   title     = model_name + " vs Observed"

   660   resm@tiMainString  = title

   662   wks   = gsn_open_wks (plot_type,plot_name)

   663 ;-----------------------------

   664 ; Add a boxed legend using the more simple method

   666   resm@pmLegendDisplayMode    = "Always"

   667 ; resm@pmLegendWidthF         = 0.1

   668   resm@pmLegendWidthF         = 0.08

   669   resm@pmLegendHeightF        = 0.05

   670   resm@pmLegendOrthogonalPosF = -1.17

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

   672 ; resm@pmLegendParallelPosF   =  0.18

   673   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   675 ; resm@lgPerimOn              = False

   676   resm@lgLabelFontHeightF     = 0.015

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

   678 ;-----------------------------

   679   tRes  = True

   680   tRes@txFontHeightF = 0.025

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

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

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

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

   688 ;Attach the vertical bar and the horizontal cap line

   690   do nd=0,1

   691     lnres@gsLineColor = line_colors(nd)

   692     do i=0,nx-1

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

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

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

   699         x1 = xvalues(nd,i)

   700         y1 = yvalues(nd,i)

   701         y2 = mn_yvalues(nd,i)

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

   704         y2 = mx_yvalues(nd,i)

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

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

   709         x1 = xvalues(nd,i) - dx4

   710         x2 = xvalues(nd,i) + dx4

   711         y1 = mn_yvalues(nd,i)

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

   714         y1 = mx_yvalues(nd,i)

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

   716       end if

   717     end do

   718   end do

   720   draw(xy)

   721   frame(wks)

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

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

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

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

   728   clear (wks)

   730  delete (DATA11_1D)

   731  delete (DATA12_1D)

   732  delete (DATA22_1D)

   733 ;delete (range)

   734 ;delete (xvalues)

   735  delete (yvalues)

   736  delete (mn_yvalues)

   737  delete (mx_yvalues)

   738  delete (good)

   739  delete (max_bar)

   740  delete (min_bar)

   741  delete (max_cap)

   742  delete (min_cap)

   743 ;-----------------------------------------------------------------

   744 ;global res

   746   resg                     = True             ; Use plot options

   747   resg@cnFillOn            = True             ; Turn on color fill

   748   resg@gsnSpreadColors     = True             ; use full colormap

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

   750 ; resg@lbLabelAutoStride   = True

   751   resg@cnLinesOn           = False            ; Turn off contourn lines

   752   resg@mpFillOn            = False            ; Turn off map fill

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

   755   resg@cnMinLevelValF       = 0.              ; Min level

   756   resg@cnMaxLevelValF       = 10.             ; Max level

   757   resg@cnLevelSpacingF      = 1.              ; interval

   759 ;global contour ob

   761   delta = 0.00001

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

   764   plot_name = "global_max_ob"

   765   title     = ob_name

   766   resg@tiMainString  = title

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

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

   771   plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg)

   772   frame(wks)

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

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

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

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

   779   clear (wks)

   780 ;------------------------------------------------------------------------

   781 ;global contour model

   783   plot_name = "global_max_model"

   784   title     = "Model " + model_name

   785   resg@tiMainString  = title

   787   wks = gsn_open_wks (plot_type,plot_name)

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

   790   plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg)

   791   frame(wks)

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

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

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

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

   798   clear (wks)

   799 ;------------------------------------------------------------------------

   800 ;(C) phase

   801 ;--------------------------------------------------------------------

   802 ; get data

   804 ; observed

   805   laiob_phase = laiob(0,:,:)

   806   s           = laiob(:,0,0)

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

   809   dsizes_z = dimsizes(laiob)

   810   nlat     = dsizes_z(1)

   811   nlon     = dsizes_z(2)

   813   do j = 0,nlat-1

   814   do i = 0,nlon-1

   815      s = laiob(:,j,i)

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

   817   end do

   818   end do

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

   821   delete (s)

   822   delete (dsizes_z)

   823 ;-------------------------

   824 ; model

   825   laimod_phase = laimod(0,:,:)

   826   s            = laimod(:,0,0)

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

   829   dsizes_z = dimsizes(laimod)

   830   nlat     = dsizes_z(1)

   831   nlon     = dsizes_z(2)

   833   do j = 0,nlat-1

   834   do i = 0,nlon-1

   835      s = laimod(:,j,i)

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

   837   end do

   838   end do

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

   841   delete (s)

   842   delete (dsizes_z)

   843 ;------------------------

   844   DATA11_1D = ndtooned(classob)

   845   DATA12_1D = ndtooned(laiob_phase)

   846   DATA22_1D = ndtooned(laimod_phase)

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

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

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

   852   do nd=0,1

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

   856     if(nd.eq.0) then

   857       data_ob  = DATA11_1D

   858       data_mod = DATA12_1D

   859     else

   860       data_ob  = DATA11_1D

   861       data_mod = DATA22_1D

   862     end if

   864 ; Loop through each range and check for values.

   866     do i=0,nr-2

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

   868 ;        print("")

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

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

   871       else

   872 ;        print("")

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

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

   875       end if

   877 ; Calculate average, and get min and max.

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

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

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

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

   883         count = dimsizes(idx)

   884       else

   885         count            = 0

   886         yvalues(nd,i)    = yvalues@_FillValue

   887         mn_yvalues(nd,i) = yvalues@_FillValue

   888         mx_yvalues(nd,i) = yvalues@_FillValue

   889       end if

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

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

   894 ; Clean up for next time in loop.

   896       delete(idx)

   897     end do

   898     delete(data_ob)

   899     delete(data_mod)

   900   end do

   901 ;-----------------------------------------------------------------

   902 ; compute correlation coef and M score

   904   u = yvalues(0,:)

   905   v = yvalues(1,:)

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

   908   uu = u(good)

   909   vv = v(good)

   911   ccrPhase = esccr(uu,vv,0)

   912 ; print (ccrPhase)

   914 ; old eq

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

   916 ; new eq

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

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

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

   922   print (Mphase)

   924  do i=0,nrow-2

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

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

   927   text4(i,8) = "-"

   928  end do

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

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

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

   933  delete (u)

   934  delete (v)

   935  delete (uu)

   936  delete (vv)

   937 ;--------------------------------------------------------------------

   938 ; histogram res

   940   resm                = True

   941   resm@gsnMaximize    = True

   942   resm@gsnDraw        = False

   943   resm@gsnFrame       = False

   944   resm@xyMarkLineMode = "Markers"

   945   resm@xyMarkerSizeF  = 0.014

   946   resm@xyMarker       = 16

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

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

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

   950   resm@trYMinF        = min(mn_yvalues) - 2

   951   resm@trYMaxF        = max(mx_yvalues) + 4

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

   954   resm@tiXAxisString  = "Land Cover Type"

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

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

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

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

   961   lnres = True

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

   963 ;------------------------------------------------------------------

   964 ; Start the graphics.

   966   plot_name = "histogram_phase"

   967   title     = model_name + " vs Observed"

   968   resm@tiMainString  = title

   970   wks   = gsn_open_wks (plot_type,plot_name)

   971 ;-----------------------------

   972 ; Add a boxed legend using the more simple method

   974   resm@pmLegendDisplayMode    = "Always"

   975 ; resm@pmLegendWidthF         = 0.1

   976   resm@pmLegendWidthF         = 0.08

   977   resm@pmLegendHeightF        = 0.05

   978   resm@pmLegendOrthogonalPosF = -1.17

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

   980 ; resm@pmLegendParallelPosF   =  0.18

   981   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   983 ; resm@lgPerimOn              = False

   984   resm@lgLabelFontHeightF     = 0.015

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

   986 ;-----------------------------

   987   tRes  = True

   988   tRes@txFontHeightF = 0.025

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

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

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

   995 ;-------------------------------

   996 ;Attach the vertical bar and the horizontal cap line

   998   do nd=0,1

   999     lnres@gsLineColor = line_colors(nd)

  1000     do i=0,nx-1

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

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

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

  1007         x1 = xvalues(nd,i)

  1008         y1 = yvalues(nd,i)

  1009         y2 = mn_yvalues(nd,i)

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

  1012         y2 = mx_yvalues(nd,i)

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

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

  1017         x1 = xvalues(nd,i) - dx4

  1018         x2 = xvalues(nd,i) + dx4

  1019         y1 = mn_yvalues(nd,i)

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

  1022         y1 = mx_yvalues(nd,i)

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

  1024       end if

  1025     end do

  1026   end do

  1028   draw(xy)

  1029   frame(wks)

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

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

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

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

  1036   clear (wks)

  1038  delete (DATA11_1D)

  1039  delete (DATA12_1D)

  1040  delete (DATA22_1D)

  1041 ;delete (range)

  1042 ;delete (xvalues)

  1043  delete (yvalues)

  1044  delete (mn_yvalues)

  1045  delete (mx_yvalues)

  1046  delete (good)

  1047  delete (max_bar)

  1048  delete (min_bar)

  1049  delete (max_cap)

  1050  delete (min_cap)

  1051 ;-----------------------------------------------------------------

  1052 ;global res

  1054   resg                     = True             ; Use plot options

  1055   resg@cnFillOn            = True             ; Turn on color fill

  1056   resg@gsnSpreadColors     = True             ; use full colormap

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

  1058 ; resg@lbLabelAutoStride   = True

  1059   resg@cnLinesOn           = False            ; Turn off contourn lines

  1060   resg@mpFillOn            = False            ; Turn off map fill

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

  1063   resg@cnMinLevelValF       = 1.              ; Min level

  1064   resg@cnMaxLevelValF       = 12.             ; Max level

  1065   resg@cnLevelSpacingF      = 1.              ; interval

  1067 ;global contour ob

  1069   delta = 0.00001

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

  1072   plot_name = "global_phase_ob"

  1073   title     = ob_name

  1074   resg@tiMainString  = title

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

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

  1079   plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)

  1080   frame(wks)

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

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

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

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

  1087   clear (wks)

  1088 ;------------------------------------------------------------------------

  1089 ;global contour model

  1091   plot_name = "global_phase_model"

  1092   title     = "Model " + model_name

  1093   resg@tiMainString  = title

  1095   wks = gsn_open_wks (plot_type,plot_name)

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

  1098   delete (plot)

  1099   plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)

  1100   frame(wks)

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

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

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

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

  1107   clear (wks)

  1108 ;-----------------------------------------------------------------

  1109 ;(D) grow day

  1110 ;--------------------------------------------------------------------

  1111 ; get data

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

  1115 ; observed

  1116   laiob_grow = laiob(0,:,:)

  1117   laiob_grow@long_name = "Days of Growing Season"

  1119   dsizes_z = dimsizes(laiob)

  1120   ntime    = dsizes_z(0)

  1121   nlat     = dsizes_z(1)

  1122   nlon     = dsizes_z(2)

  1124   do j = 0,nlat-1

  1125   do i = 0,nlon-1

  1126      nday = 0.

  1127      do k = 0,ntime-1

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

  1129            nday = nday + day_of_data(k)

  1130         end if

  1131      end do

  1133      laiob_grow(j,i) = nday

  1134   end do

  1135   end do

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

  1138 ;-------------------------

  1139 ; model

  1140   laimod_grow = laimod(0,:,:)

  1141   laimod_grow@long_name = "Days of Growing Season"

  1143   dsizes_z = dimsizes(laimod)

  1144   ntime    = dsizes_z(0)

  1145   nlat     = dsizes_z(1)

  1146   nlon     = dsizes_z(2)

  1148   do j = 0,nlat-1

  1149   do i = 0,nlon-1

  1150      nday = 0.

  1151      do k = 0,ntime-1

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

  1153            nday = nday + day_of_data(k)

  1154         end if

  1155      end do

  1157      laimod_grow(j,i) = nday

  1158   end do

  1159   end do

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

  1162 ;------------------------

  1163   DATA11_1D = ndtooned(classob)

  1164   DATA12_1D = ndtooned(laiob_grow)

  1165   DATA22_1D = ndtooned(laimod_grow)

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

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

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

  1171   do nd=0,1

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

  1175     if(nd.eq.0) then

  1176       data_ob  = DATA11_1D

  1177       data_mod = DATA12_1D

  1178     else

  1179       data_ob  = DATA11_1D

  1180       data_mod = DATA22_1D

  1181     end if

  1183 ; Loop through each range and check for values.

  1185     do i=0,nr-2

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

  1187 ;        print("")

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

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

  1190       else

  1191 ;        print("")

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

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

  1194       end if

  1196 ; Calculate average, and get min and max.

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

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

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

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

  1202         count = dimsizes(idx)

  1203       else

  1204         count            = 0

  1205         yvalues(nd,i)    = yvalues@_FillValue

  1206         mn_yvalues(nd,i) = yvalues@_FillValue

  1207         mx_yvalues(nd,i) = yvalues@_FillValue

  1208       end if

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

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

  1213 ; Clean up for next time in loop.

  1215       delete(idx)

  1216     end do

  1217     delete(data_ob)

  1218     delete(data_mod)

  1219   end do

  1220 ;-----------------------------------------------------------------

  1221 ; compute correlation coef and M score

  1223   u = yvalues(0,:)

  1224   v = yvalues(1,:)

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

  1227   uu = u(good)

  1228   vv = v(good)

  1230   ccrGrow = esccr(uu,vv,0)

  1231 ; print (ccrGrow)

  1233 ; new eq

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

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

  1237   print (Mgrow)

  1239  do i=0,nrow-2

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

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

  1242   text4(i,11) = "-"

  1243  end do

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

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

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

  1248  delete (u)

  1249  delete (v)

  1250  delete (uu)

  1251  delete (vv)

  1252 ;--------------------------------------------------------------------

  1253 ; histogram res

  1255   resm                = True

  1256   resm@gsnMaximize    = True

  1257   resm@gsnDraw        = False

  1258   resm@gsnFrame       = False

  1259   resm@xyMarkLineMode = "Markers"

  1260   resm@xyMarkerSizeF  = 0.014

  1261   resm@xyMarker       = 16

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

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

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

  1265   resm@trYMinF        = min(mn_yvalues) - 2

  1266   resm@trYMaxF        = max(mx_yvalues) + 4

  1268   resm@tiYAxisString = "Days of Growing season"

  1269   resm@tiXAxisString = "Land Cover Type"

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

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

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

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

  1276   lnres = True

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

  1278 ;------------------------------------------------------------------

  1279 ; Start the graphics.

  1281   plot_name = "histogram_grow"

  1282   title     = model_name + " vs Observed"

  1283   resm@tiMainString  = title

  1285   wks   = gsn_open_wks (plot_type,plot_name)

  1286 ;-----------------------------

  1287 ; Add a boxed legend using the more simple method

  1289   resm@pmLegendDisplayMode    = "Always"

  1290 ; resm@pmLegendWidthF         = 0.1

  1291   resm@pmLegendWidthF         = 0.08

  1292   resm@pmLegendHeightF        = 0.05

  1293   resm@pmLegendOrthogonalPosF = -1.17

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

  1295 ; resm@pmLegendParallelPosF   =  0.18

  1296   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

  1298 ; resm@lgPerimOn              = False

  1299   resm@lgLabelFontHeightF     = 0.015

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

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

  1302   tRes  = True

  1303   tRes@txFontHeightF = 0.025

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

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

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

  1310 ;-------------------------------

  1311 ;Attach the vertical bar and the horizontal cap line

  1313   do nd=0,1

  1314     lnres@gsLineColor = line_colors(nd)

  1315     do i=0,nx-1

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

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

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

  1322         x1 = xvalues(nd,i)

  1323         y1 = yvalues(nd,i)

  1324         y2 = mn_yvalues(nd,i)

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

  1327         y2 = mx_yvalues(nd,i)

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

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

  1332         x1 = xvalues(nd,i) - dx4

  1333         x2 = xvalues(nd,i) + dx4

  1334         y1 = mn_yvalues(nd,i)

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

  1337         y1 = mx_yvalues(nd,i)

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

  1339       end if

  1340     end do

  1341   end do

  1343   draw(xy)

  1344   frame(wks)

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

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

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

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

  1351   clear (wks)

  1353  delete (DATA11_1D)

  1354  delete (DATA12_1D)

  1355  delete (DATA22_1D)

  1356 ;delete (range)

  1357 ;delete (xvalues)

  1358  delete (yvalues)

  1359  delete (mn_yvalues)

  1360  delete (mx_yvalues)

  1361  delete (good)

  1362  delete (max_bar)

  1363  delete (min_bar)

  1364  delete (max_cap)

  1365  delete (min_cap)

  1366 ;-----------------------------------------------------------------

  1367 ;global res

  1369   resg                     = True             ; Use plot options

  1370   resg@cnFillOn            = True             ; Turn on color fill

  1371   resg@gsnSpreadColors     = True             ; use full colormap

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

  1373 ; resg@lbLabelAutoStride   = True

  1374   resg@cnLinesOn           = False            ; Turn off contourn lines

  1375   resg@mpFillOn            = False            ; Turn off map fill

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

  1378   resg@cnMinLevelValF       = 60.             ; Min level

  1379   resg@cnMaxLevelValF       = 360.            ; Max level

  1380   resg@cnLevelSpacingF      = 20.             ; interval

  1382 ;global contour ob

  1384   laiob_grow@_FillValue = 1.e+36

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

  1387   plot_name = "global_grow_ob"

  1388   title     = ob_name

  1389   resg@tiMainString  = title

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

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

  1394   plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)

  1395   frame(wks)

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

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

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

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

  1402   clear (wks)

  1403 ;------------------------------------------------------------------------

  1404 ;global contour model

  1406   laimod_grow@_FillValue = 1.e+36

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

  1409   plot_name = "global_grow_model"

  1410   title     = "Model " + model_name

  1411   resg@tiMainString  = title

  1413   wks = gsn_open_wks (plot_type,plot_name)

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

  1416   delete (plot)

  1417   plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)

  1418   frame(wks)

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

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

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

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

  1425   clear (wks)

  1426 ;------------------------------------------------------------------------

  1427 ;global contour model vs ob

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

  1430   resg@cnMinLevelValF       = 0.              ; Min level

  1431   resg@cnMaxLevelValF       = 10.             ; Max level

  1432   resg@cnLevelSpacingF      = 1.              ; interval

  1434   plot_name = "global_mean_model_vs_ob"

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

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

  1439   delete (plot)

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

  1442   resg@gsnFrame             = False          ; Do not draw plot

  1443   resg@gsnDraw              = False          ; Do not advance frame

  1445 ; plot correlation coef

  1447   gRes               = True

  1448   gRes@txFontHeightF = 0.02

  1449   gRes@txAngleF      = 90

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

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

  1454 ;--------------------------------------------------------------------

  1456 ;(a) ob

  1458   title     = ob_name

  1459   resg@tiMainString  = title

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

  1463 ;(b) model

  1465   title     = "Model "+ model_name

  1466   resg@tiMainString  = title

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

  1470 ;(c) model-ob

  1472   zz = laimod_mean

  1473   zz = laimod_mean - laiob_mean

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

  1475   resg@tiMainString    = title

  1477   resg@cnMinLevelValF  = -2.           ; Min level

  1478   resg@cnMaxLevelValF  =  2.           ; Max level

  1479   resg@cnLevelSpacingF =  0.4           ; interval

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

  1484   pres                            = True        ; panel plot mods desired

  1485   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1486                                                 ; indiv. plots in panel

  1487   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

  1496   frame (wks)

  1497   clear (wks)

  1499   delete (plot)

  1500 ;-----------------------------------------------------------------

  1501 ;global contour model vs ob

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

  1504   resg@cnMinLevelValF       = 0.              ; Min level

  1505   resg@cnMaxLevelValF       = 10.             ; Max level

  1506   resg@cnLevelSpacingF      = 1.              ; interval

  1508   plot_name = "global_max_model_vs_ob"

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

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

  1513 ; delete(plot)

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

  1516   resg@gsnFrame             = False          ; Do not draw plot

  1517   resg@gsnDraw              = False          ; Do not advance frame

  1519 ; plot correlation coef

  1521   gRes               = True

  1522   gRes@txFontHeightF = 0.02

  1523   gRes@txAngleF      = 90

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

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

  1528 ;--------------------------------------------------------------------

  1529 ;(a) ob

  1531   title     = ob_name

  1532   resg@tiMainString  = title

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

  1536 ;(b) model

  1538   title     = "Model "+ model_name

  1539   resg@tiMainString  = title

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

  1543 ;(c) model-ob

  1545   delete (zz)

  1546   zz = laimod_max

  1547   zz = laimod_max - laiob_max

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

  1549   resg@tiMainString    = title

  1551   resg@cnMinLevelValF  = -6.           ; Min level

  1552   resg@cnMaxLevelValF  =  6.           ; Max level

  1553   resg@cnLevelSpacingF =  1.           ; interval

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

  1558   pres                            = True        ; panel plot mods desired

  1559   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1560                                                 ; indiv. plots in panel

  1561   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

  1570   frame (wks)

  1571   clear (wks)

  1573   delete (plot)

  1574 ;-----------------------------------------------------------------

  1575 ;global contour model vs ob

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

  1578   resg@cnMinLevelValF       = 1.              ; Min level

  1579   resg@cnMaxLevelValF       = 12.             ; Max level

  1580   resg@cnLevelSpacingF      = 1.              ; interval

  1582   plot_name = "global_phase_model_vs_ob"

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

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

  1587 ; delete (plot)

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

  1590   resg@gsnFrame             = False          ; Do not draw plot

  1591   resg@gsnDraw              = False          ; Do not advance frame

  1593 ; plot correlation coef

  1595   gRes               = True

  1596   gRes@txFontHeightF = 0.02

  1597   gRes@txAngleF      = 90

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

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

  1602 ;--------------------------------------------------------------------

  1603 ;(a) ob

  1605   title     = ob_name

  1606   resg@tiMainString  = title

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

  1610 ;(b) model

  1612   title     = "Model "+ model_name

  1613   resg@tiMainString  = title

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

  1617 ;(c) model-ob

  1619   delete (zz)

  1620   zz = laimod_phase

  1621   zz = laimod_phase - laiob_phase

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

  1623   resg@tiMainString    = title

  1625   resg@cnMinLevelValF  = -6.           ; Min level

  1626   resg@cnMaxLevelValF  =  6.           ; Max level

  1627   resg@cnLevelSpacingF =  1.           ; interval

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

  1632 ; pres                            = True        ; panel plot mods desired

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

  1634                                                 ; indiv. plots in panel

  1635 ; pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

  1644   frame (wks)

  1645   clear (wks)

  1647   delete (plot)

  1648 ;------------------------------------------------------------------

  1649 ;global contour model vs ob

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

  1652   resg@cnMinLevelValF       = 60.             ; Min level

  1653   resg@cnMaxLevelValF       = 360.            ; Max level

  1654   resg@cnLevelSpacingF      = 20.             ; interval

  1656   plot_name = "global_grow_model_vs_ob"

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

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

  1661 ; delete (plot)

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

  1664   resg@gsnFrame             = False          ; Do not draw plot

  1665   resg@gsnDraw              = False          ; Do not advance frame

  1667 ; plot correlation coef

  1669   gRes               = True

  1670   gRes@txFontHeightF = 0.02

  1671   gRes@txAngleF      = 90

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

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

  1676 ;--------------------------------------------------------------------

  1677 ;(a) ob

  1679   title     = ob_name

  1680   resg@tiMainString  = title

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

  1684 ;(b) model

  1686   title     = "Model "+ model_name

  1687   resg@tiMainString  = title

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

  1691 ;(c) model-ob

  1693   delete (zz)

  1694   zz = laimod_grow

  1695   zz = laimod_grow - laiob_grow

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

  1697   resg@tiMainString    = title

  1699   resg@cnMinLevelValF  = -100.           ; Min level

  1700   resg@cnMaxLevelValF  =  100.           ; Max level

  1701   resg@cnLevelSpacingF =  20.            ; interval

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

  1705   pres                            = True        ; panel plot mods desired

  1706   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1707                                                 ; indiv. plots in panel

  1708   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

  1717   frame (wks)

  1718   clear (wks)

  1720   delete (plot)

  1721 ;**************************************************

  1722 ; plot lai table

  1723 ;**************************************************

  1725   plot_name = "table_lai"

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

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

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

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

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

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

  1734   frame(wks)

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

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

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

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

  1740 ;-------------------------------------------------------------------

  1741   temp_name = "temp." + model_name

  1742   system("mkdir -p " + temp_name)

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

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

  1745 ;-------------------------------------------------------------------

  1746 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
permissions	-rw-r--r--