clamp: npp/99.all.ncl.0@4be95183fbcd

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

     2 ; combine scatter, histogram, global and zonal plots

     3 ; compute all correlation coef and M score

     4 ; add 1-to-1 line in scatter plots

     5 ; ****************************************************

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

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

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

     9 ;************************************************

    10 procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \

    11                  ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \

    12                  ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \

    13                  ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \

    14                  ,dx4[1]:numeric \

    15                   )

    16 begin

    17 ; Calculaee "nice" bins for binning the data in equally spaced ranges.

    18 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D

    19 ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4

    21   nbins       = 15     ; Number of bins to use.

    23   nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)

    24   nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)

    25   range       = fspan(nicevals(0),nicevals(1),nvals)

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

    28 ; particular range, and then take an average.

    30   nr           = dimsizes(range)

    31   nx           = nr-1

    32 ; print (nx)

    34 ; xvalues      = new((/2,nx/),typeof(RAIN1_1D))

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

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

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

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

    39 ; yvalues      = new((/2,nx/),typeof(RAIN1_1D))

    40 ; mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

    41 ; mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

    43   do nd=0,1

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

    47     if(nd.eq.0) then

    48       data     = RAIN1_1D

    49       npp_data = NPP1_1D

    50     else

    51       data     = RAIN2_1D

    52       npp_data = NPP2_1D

    53     end if

    55 ; Loop through each range and check for values.

    57     do i=0,nr-2

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

    59 ;        print("")

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

    61         idx = ind((range(i).le.data).and.(data.lt.range(i+1)))

    62       else

    63 ;        print("")

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

    65         idx = ind(range(i).le.data)

    66       end if

    68 ; Calculate average, and get min and max.

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

    71         yvalues(nd,i)    = avg(npp_data(idx))

    72         mn_yvalues(nd,i) = min(npp_data(idx))

    73         mx_yvalues(nd,i) = max(npp_data(idx))

    74         count = dimsizes(idx)

    75       else

    76         count            = 0

    77         yvalues(nd,i)    = yvalues@_FillValue

    78         mn_yvalues(nd,i) = yvalues@_FillValue

    79         mx_yvalues(nd,i) = yvalues@_FillValue

    80       end if

    82 ; Print out information.

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

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

    88 ; Clean up for next time in loop.

    90       delete(idx)

    91     end do

    92     delete(data)

    93     delete(npp_data)

    94   end do

    95 end

    97 ; Main code.

    98 begin

   100  plot_type     = "ps"

   101  plot_type_new = "png"

   102 ;************************************************

   103 ; read in model data

   104 ;************************************************

   105 ;film = "i01.06cn_1798-2004_ANN_climo.nc"

   106 ;model_name = "06cn"

   107 ;model_grid = "T42"

   109 ;film = "i01.06casa_1798-2004_ANN_climo.nc"

   110 ;model_name = "06casa"

   111 ;model_grid = "T42"

   113  film = "i01.10casa_1948-2004_ANN_climo.nc"

   114  model_name = "10casa"

   115  model_grid = "T42"

   117 ;film = "i01.10cn_1948-2004_ANN_climo.nc"

   118 ;model_name = "10cn"

   119 ;model_grid = "T42"

   120 ;--------------------------------------------------

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

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

   124  nppmod0  = fm->NPP

   125  rainmod0 = fm->RAIN

   126  xm       = fm->lon

   127  ym       = fm->lat

   129 ;************************************************

   130 ; read in ob data

   131 ;************************************************

   132 ;(1)

   133  dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

   134  fil81 = "data.81.nc"

   135  f81   = addfile (dir81+fil81,"r")

   137  id81   = f81->SITE_ID

   138  npp81  = f81->TNPP_C

   139  rain81 = tofloat(f81->PREC_ANN)

   140  x81    = f81->LONG_DD

   141  y81    = f81->LAT_DD

   143  id81@long_name  = "SITE_ID"

   145 ;change longitude from (-180,180) to (0,360)

   146 ;for model data interpolation

   148  do i= 0,dimsizes(x81)-1

   149     if (x81(i) .lt. 0.) then

   150         x81(i) = x81(i)+ 360.

   151     end if

   152  end do

   153 ;print (x81)

   154 ;-------------------------------------

   155 ;(2)

   156  dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

   157  fil933 = "data.933.nc"

   158  f933   = addfile (dir933+fil933,"r")

   160  id933   = f933->SITE_ID

   161  npp933  = f933->TNPP_C

   162  rain933 = f933->PREC

   163  x933    = f933->LONG_DD

   164  y933    = f933->LAT_DD

   166  id933@long_name  = "SITE_ID"

   168 ;change longitude from (-180,180) to (0,360)

   169 ;for model data interpolation

   171  do i= 0,dimsizes(x933)-1

   172     if (x933(i) .lt. 0.) then

   173         x933(i) = x933(i)+ 360.

   174     end if

   175  end do

   176 ;print (x933)

   177 ;----------------------------------------

   178 ;(3)

   179  dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

   180  filglobe = "Npp_"+model_grid+"_mean.nc"

   181  fglobe   = addfile (dirglobe+filglobe,"r")

   183  nppglobe0 = fglobe->NPP

   184  nppglobe  = nppglobe0

   185 ;***********************************************************************

   186 ; interpolate model data

   187 ;***********************************************************************

   188  nppmod   = nppmod0(0,:,:)

   189  rainmod  = rainmod0(0,:,:)

   190  delete (nppmod0)

   191  delete (rainmod0)

   193  nppmod81  =linint2_points(xm,ym,nppmod,True,x81,y81,0)

   195  nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)

   197  rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)

   199  rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)

   201 ;************************************************

   202 ; Units for these variables are:

   203 ;

   204 ; rain81  : mm/year

   205 ; rainmod : mm/s

   206 ; npp81   : g C/m^2/year

   207 ; nppmod81: g C/m^2/s

   208 ; nppglobe: g C/m^2/year

   209 ;

   210 ; We want to convert these to "m/year" and "g C/m^2/year".

   211 ;

   212   nsec_per_year = 60*60*24*365

   214   rain81    = rain81 / 1000.

   215   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   216   nppmod81  = nppmod81 * nsec_per_year

   218   rain933    = rain933 / 1000.

   219   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   220   nppmod933  = nppmod933 * nsec_per_year

   222   nppmod  = nppmod * nsec_per_year

   224   npp81@units      = "gC/m^2/yr"

   225   nppmod81@units   = "gC/m^2/yr"

   226   npp933@units     = "gC/m^2/yr"

   227   nppmod933@units  = "gC/m^2/yr"

   228   nppmod@units     = "gC/m^2/yr"

   229   nppglobe@units   = "gC/m^2/yr"

   230   rain81@units     = "m/yr"

   231   rainmod81@units  = "m/yr"

   232   rain933@units    = "m/yr"

   233   rainmod933@units = "m/yr"

   235   npp81@long_name      = "NPP (gC/m2/year)"

   236   npp933@long_name     = "NPP (gC/m2/year)"

   237   nppmod81@long_name   = "NPP (gC/m2/year)"

   238   nppmod933@long_name  = "NPP (gC/m2/year)"

   239   nppmod@long_name     = "NPP (gC/m2/year)"

   240   nppglobe@long_name   = "NPP (gC/m2/year)"

   241   rain81@long_name     = "PREC (m/year)"

   242   rain933@long_name    = "PREC (m/year)"

   243   rainmod81@long_name  = "PREC (m/year)"

   244   rainmod933@long_name = "PREC (m/year)"

   246 ;(A)-1 plot scatter ob 81

   248 ;plot_name = "scatter_ob_81"

   249 ;title     = "Observed 81 sites"

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

   252 ;res                   = True                  ; plot mods desired

   253 ;res@tiMainString      = title                 ; add title

   254 ;res@xyMarkLineModes   = "Markers"             ; choose which have markers

   255 ;res@xyMarkers         =  16                   ; choose type of marker

   256 ;res@xyMarkerColor     = "red"                 ; Marker color

   257 ;res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   258 ;res@tmLabelAutoStride = True                  ; nice tick mark labels

   260 ;plot  = gsn_csm_xy (wks,id81,npp81,res)       ; create plot

   261 ;frame(wks)

   263 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

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

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

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

   268 ;clear (wks)

   269 ;*******************************************************************

   270 ;(A)-2 for table of site(1) (the first 41 sites)

   271 ;*******************************************************************

   273   table_length = 0.95

   275 ; table header name

   276   table_header_name = "Site ID"

   278 ; column (not including header column)

   279   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)

   280   ncol = dimsizes(col_header)

   282 ; row (not including header row)

   283   nrow       = 41

   284   row_header = new ((/nrow/),string )

   285   row_header(0:nrow-1) = id81(0:nrow-1)

   287 ; Table header

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

   289   x1    = (/0.005,0.15/)        ; Start and end X

   290   y1    = (/0.900,0.95/)       ; Start and end Y

   291   text1 = table_header_name

   292   res1               = True

   293   res1@txFontHeightF = 0.015

   294   res1@gsFillColor   = "CornFlowerBlue"

   296 ; Column header (equally space in x2)

   297   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   298   x2    = (/x1(1),0.995/)       ; start from end of x1

   299   y2    = y1                    ; same as y1

   300   text2 = col_header

   301   res2               = True

   302   res2@txFontHeightF = 0.015

   303   res2@gsFillColor   = "Gray"

   305 ; Row header (equally space in y2)

   306   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   307   x3    = x1                         ; same as x1

   308   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   309   text3 = row_header

   310   res3               = True

   311   res3@txFontHeightF = 0.010

   312   res3@gsFillColor   = "Gray"

   314 ; Main table body

   315   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   316   x4    = x2                      ; Start and end x

   317   y4    = y3                      ; Start and end Y

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

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

   321   color_fill4           = "white"

   322 ; color_fill4(:,ncol-1) = "grey"

   323 ; color_fill4(nrow-1,:) = "green"

   325   res4               = True       ; Set up resource list

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

   327   res4@txFontHeightF = 0.015

   328   res4@gsFillColor   = color_fill4

   330   delete (color_fill4)

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

   332 ; for table value

   334   do n = 0,nrow-1

   335      text4(n,0) = sprintf("%5.2f", y81(n))

   336      text4(n,1) = sprintf("%5.2f", x81(n))

   337      text4(n,2) = sprintf("%5.2f", npp81(n))

   338      text4(n,3) = sprintf("%5.2f", rain81(n))

   339   end do

   340 ;---------------------------------------------------------------------------

   342   plot_name = "table_site81_ob1"

   344   wks = gsn_open_wks (plot_type,plot_name)

   345 ;------------------------------------------

   346 ; for table title

   348   gRes               = True

   349   gRes@txFontHeightF = 0.02

   350 ; gRes@txAngleF      = 90

   352   title_text = "Observation at 81 Sites (1)"

   354   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   355 ;------------------------------------------

   357   gsn_table(wks,ncr1,x1,y1,text1,res1)

   358   gsn_table(wks,ncr2,x2,y2,text2,res2)

   359   gsn_table(wks,ncr3,x3,y3,text3,res3)

   360   gsn_table(wks,ncr4,x4,y4,text4,res4)

   362   frame(wks)

   363   clear (wks)

   365   delete (row_header)

   366   delete (res3)

   367   delete (ncr3)

   368   delete (text3)

   369   delete (res4)

   370   delete (ncr4)

   371   delete (text4)

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

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

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

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

   378 ;*******************************************************************

   379 ;(A)-3 for table of site(2) (the second 40 sites)

   380 ;*******************************************************************

   382   table_length = 0.95

   384 ; table header name

   385   table_header_name = "Site ID"

   387 ; column (not including header column)

   388   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)

   389   ncol = dimsizes(col_header)

   391 ; row (not including header row)

   392   nrow       = 40

   393   row_header = new ((/nrow/),string )

   394   row_header(0:nrow-1) = id81(41:41+nrow-1)

   396 ; Table header

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

   398   x1    = (/0.005,0.15/)        ; Start and end X

   399   y1    = (/0.900,0.95/)       ; Start and end Y

   400   text1 = table_header_name

   401   res1               = True

   402   res1@txFontHeightF = 0.015

   403   res1@gsFillColor   = "CornFlowerBlue"

   405 ; Column header (equally space in x2)

   406   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   407   x2    = (/x1(1),0.995/)       ; start from end of x1

   408   y2    = y1                    ; same as y1

   409   text2 = col_header

   410   res2               = True

   411   res2@txFontHeightF = 0.015

   412   res2@gsFillColor   = "Gray"

   414 ; Row header (equally space in y2)

   415   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   416   x3    = x1                         ; same as x1

   417   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   418   text3 = row_header

   419   res3               = True

   420   res3@txFontHeightF = 0.010

   421   res3@gsFillColor   = "Gray"

   423 ; Main table body

   424   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   425   x4    = x2                      ; Start and end x

   426   y4    = y3                      ; Start and end Y

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

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

   430   color_fill4           = "white"

   431 ; color_fill4(:,ncol-1) = "grey"

   432 ; color_fill4(nrow-1,:) = "green"

   434   res4               = True       ; Set up resource list

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

   436   res4@txFontHeightF = 0.015

   437   res4@gsFillColor   = color_fill4

   439   delete (color_fill4)

   440 ;-------------------------------------------------------------------

   441 ; for table value

   443   do n = 0,nrow-1

   444      text4(n,0) = sprintf("%5.2f", y81(n+41))

   445      text4(n,1) = sprintf("%5.2f", x81(n+41))

   446      text4(n,2) = sprintf("%5.2f", npp81(n+41))

   447      text4(n,3) = sprintf("%5.2f", rain81(n+41))

   448   end do

   449 ;---------------------------------------------------------------------------

   451   plot_name = "table_site81_ob2"

   453   wks = gsn_open_wks (plot_type,plot_name)

   454 ;------------------------------------------

   455 ; for table title

   457   gRes               = True

   458   gRes@txFontHeightF = 0.02

   459 ; gRes@txAngleF      = 90

   461   title_text = "Observation at 81 Sites (2)"

   463   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   464 ;------------------------------------------

   466   gsn_table(wks,ncr1,x1,y1,text1,res1)

   467   gsn_table(wks,ncr2,x2,y2,text2,res2)

   468   gsn_table(wks,ncr3,x3,y3,text3,res3)

   469   gsn_table(wks,ncr4,x4,y4,text4,res4)

   471   frame(wks)

   472   clear (wks)

   474   delete (row_header)

   475   delete (res3)

   476   delete (ncr3)

   477   delete (text3)

   478   delete (res4)

   479   delete (ncr4)

   480   delete (text4)

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

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

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

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

   486 ;**************************************************************************

   487 ;(A)-4 plot scatter ob 933

   489  plot_name = "scatter_ob_933"

   490  title     = "Observed 933 sites"

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

   493  res                   = True                  ; plot mods desired

   494  res@tiMainString      = title                 ; add title

   495  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   496  res@xyMarkers         =  16                   ; choose type of marker

   497  res@xyMarkerColor     = "red"                 ; Marker color

   498  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   499  res@tmLabelAutoStride = True                  ; nice tick mark labels

   501  plot  = gsn_csm_xy (wks,id933,npp933,res)     ; create plot

   502  frame(wks)

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

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

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

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

   509  clear (wks)

   511 ;(A)-5 plot scatter model 81

   513 ;plot_name = "scatter_model_81"

   514 ;title     = "Model "+ model_name +" 81 sites"

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

   517 ;res                   = True                  ; plot mods desired

   518 ;res@tiMainString      = title                 ; add title

   519 ;res@xyMarkLineModes   = "Markers"             ; choose which have markers

   520 ;res@xyMarkers         =  16                   ; choose type of marker

   521 ;res@xyMarkerColor     = "red"                 ; Marker color

   522 ;res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   523 ;res@tmLabelAutoStride = True                  ; nice tick mark labels

   525 ;plot  = gsn_csm_xy (wks,id81,nppmod81,res)    ; create plot

   526 ;frame(wks)

   528 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

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

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

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

   533 ;clear (wks)

   535 ;(A)-6 plot scatter model 933

   537  plot_name = "scatter_model_933"

   538  title     = "Model "+ model_name +" 933 sites"

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

   541  res                   = True                  ; plot mods desired

   542  res@tiMainString      = title                 ; add title

   543  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   544  res@xyMarkers         =  16                   ; choose type of marker

   545  res@xyMarkerColor     = "red"                 ; Marker color

   546  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   547  res@tmLabelAutoStride = True                  ; nice tick mark labels

   549  plot  = gsn_csm_xy (wks,id933,nppmod933,res)    ; create plot

   550  frame(wks)

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

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

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

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

   557  clear (wks)

   559 ;(A)-7 scatter model vs ob 81

   561  plot_name = "scatter_model_vs_ob_81"

   562  title     = model_name +" vs ob 81 sites"

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

   565  res                   = True                  ; plot mods desired

   566  res@tiMainString      = title                 ; add title

   567  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   568  res@xyMarkers         =  16                   ; choose type of marker

   569  res@xyMarkerColor     = "red"                 ; Marker color

   570  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   571  res@tmLabelAutoStride = True                  ; nice tick mark labels

   573  res@gsnDraw           = False

   574  res@gsnFrame          = False                 ; don't advance frame yet

   575 ;-------------------------------

   576 ;compute correlation coef. and M

   577  ccr81 = esccr(nppmod81,npp81,0)

   578 ;print (ccr81)

   580 ;new eq

   581  bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))

   582  M81  = (1. - (bias/dimsizes(y81)))*5.

   583  print (M81)

   585  tRes  = True

   586  tRes@txFontHeightF = 0.025

   588  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"

   590  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

   591 ;--------------------------------

   592  plot  = gsn_csm_xy (wks,npp81,nppmod81,res)       ; create plot

   593 ;-------------------------------

   594 ; add polyline

   596  dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)

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

   598  draw (plot)

   599  frame(wks)

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

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

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

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

   606  clear (wks)

   607  delete (dum)

   609 ;(A)-8 scatter model vs ob 933

   611  plot_name = "scatter_model_vs_ob_933"

   612  title     = model_name +" vs ob 933 sites"

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

   615  res                   = True                  ; plot mods desired

   616  res@tiMainString      = title                 ; add title

   617  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   618  res@xyMarkers         =  16                   ; choose type of marker

   619  res@xyMarkerColor     = "red"                 ; Marker color

   620  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   621  res@tmLabelAutoStride = True                  ; nice tick mark labels

   623  res@gsnDraw           = False

   624  res@gsnFrame          = False                 ; don't advance frame yet

   625 ;-------------------------------

   626 ;compute correlation coef. and M

   627  ccr933 = esccr(nppmod933,npp933,0)

   628 ;print (ccr933)

   630 ;new eq

   631  bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))

   632  M933 = (1. - (bias/dimsizes(y933)))*5.

   633  print (M933)

   635  tRes  = True

   636  tRes@txFontHeightF = 0.025

   638  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"

   640  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

   641 ;--------------------------------

   642  plot  = gsn_csm_xy (wks,npp933,nppmod933,res)    ; create plot

   643 ;-------------------------------

   644 ; add polyline

   646  dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)

   647 ;-------------------------------

   648  draw (plot)

   649  frame(wks)

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

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

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

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

   656  clear (wks)

   657  delete (dum)

   659 ;*******************************************************************

   660 ;(A)-9 for table of site(1) (the first 41 sites), model vs ob

   661 ;*******************************************************************

   663   table_length = 0.95

   665 ; table header name

   666   table_header_name = "Site ID"

   668 ; row (not including header row)

   669   nrow       = 41

   670   row_header = new ((/nrow/),string )

   671   row_header(0:nrow-1) = id81(0:nrow-1)

   673 ; Table header

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

   675   x1    = (/0.005,0.15/)        ; Start and end X

   676   y1    = (/0.85 ,0.95/)        ; Start and end Y

   677   text1 = table_header_name

   678   res1               = True

   679   res1@txFontHeightF = 0.015

   680   res1@gsFillColor   = "CornFlowerBlue"

   682 ; Column header1 (equally space in x2)

   684   delete (col_header)

   685   delete (ncr2)

   686   delete (text2)

   688   col_header = (/"Latitude","Longitude"/)

   689   ncol = dimsizes(col_header)

   691   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   692   x2    = (/x1(1),0.35/)        ; start from end of x1

   693   y2    = y1                    ; same as y1

   694   text2 = col_header

   695   res2               = True

   696   res2@txFontHeightF = 0.015

   697   res2@gsFillColor   = "Gray"

   699 ; Column header2 (equally space in x2)

   701   col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   702   col_header2 = (/"ob",model_name \

   703                  ,"ob",model_name \

   704                  /)

   706   ncol1 = dimsizes(col_header1)

   707   ncol2 = dimsizes(col_header2)

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

   710   x21    = (/x2(1),0.995/)        ; start from end of x2

   711   yhalf  = (y1(0)+y1(1))*0.5

   712   y21    = (/yhalf,y1(1)/)        ; top half of y1

   713   text21 = col_header1

   714   res21               = True

   715   res21@txFontHeightF = 0.015

   716   res21@gsFillColor   = "Gray"

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

   719   x22    = x21                    ; start from end of x1

   720   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   721   text22 = col_header2

   722   res22               = True

   723   res22@txFontHeightF = 0.015

   724   res22@gsFillColor   = "Gray"

   726 ; Row header (equally space in y2)

   727   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   728   x3    = x1                         ; same as x1

   729   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   730   text3 = row_header

   731   res3               = True

   732   res3@txFontHeightF = 0.010

   733   res3@gsFillColor   = "Gray"

   735 ; Main table body-1

   736   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   737   x4    = x2                      ; Start and end x

   738   y4    = y3                      ; Start and end Y

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

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

   742   color_fill4           = "white"

   743 ; color_fill4(:,ncol-1) = "grey"

   744 ; color_fill4(nrow-1,:) = "green"

   746   res4               = True       ; Set up resource list

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

   748   res4@txFontHeightF = 0.015

   749   res4@gsFillColor   = color_fill4

   751   delete (color_fill4)

   753 ; Main table body-2

   754   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   755   x42    = x21                      ; Start and end x

   756   y42    = y3                       ; Start and end Y

   757   text42 = new((/nrow,ncol2/),string)

   759   color_fill42           = new((/nrow,ncol2/),string)

   760   color_fill42           = "white"

   761 ; color_fill42(:,ncol-1) = "grey"

   762 ; color_fill42(nrow-1,:) = "green"

   764   res42               = True       ; Set up resource list

   765 ; res42@gsnDebug      = True       ; Useful to print NDC row,col values used.

   766   res42@txFontHeightF = 0.015

   767   res42@gsFillColor   = color_fill42

   769   delete (color_fill42)

   770 ;-------------------------------------------------------------------

   771 ; for table value

   773   do n = 0,nrow-1

   774      text4(n,0) = sprintf("%5.2f", y81(n))

   775      text4(n,1) = sprintf("%5.2f", x81(n))

   776   end do

   778   do n = 0,nrow-1

   779      text42(n,0) = sprintf("%5.2f", npp81(n))

   780      text42(n,1) = sprintf("%5.2f", nppmod81(n))

   781      text42(n,2) = sprintf("%5.2f", rain81(n))

   782      text42(n,3) = sprintf("%5.2f", rainmod81(n))

   783   end do

   784 ;---------------------------------------------------------------------------

   786   plot_name = "table_site81_model_vs_ob1"

   788   wks = gsn_open_wks (plot_type,plot_name)

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

   790 ; for table title

   792   gRes               = True

   793   gRes@txFontHeightF = 0.02

   794 ; gRes@txAngleF      = 90

   796   title_text = "Model vs Observation at 81 Sites (1)"

   798   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

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

   801   gsn_table(wks,ncr1,x1,y1,text1,res1)

   802   gsn_table(wks,ncr2,x2,y2,text2,res2)

   803   gsn_table(wks,ncr21,x21,y21,text21,res21)

   804   gsn_table(wks,ncr22,x22,y22,text22,res22)

   805   gsn_table(wks,ncr3,x3,y3,text3,res3)

   806   gsn_table(wks,ncr4,x4,y4,text4,res4)

   807   gsn_table(wks,ncr42,x42,y42,text42,res42)

   809   frame(wks)

   810   clear (wks)

   812   delete (col_header)

   813   delete (row_header)

   814   delete (res1)

   815   delete (ncr1)

   816   delete (text1)

   817   delete (res2)

   818   delete (ncr2)

   819   delete (text2)

   820   delete (res21)

   821   delete (ncr21)

   822   delete (text21)

   823   delete (res22)

   824   delete (ncr22)

   825   delete (text22)

   826   delete (res3)

   827   delete (ncr3)

   828   delete (text3)

   829   delete (res4)

   830   delete (ncr4)

   831   delete (text4)

   832   delete (res42)

   833   delete (ncr42)

   834   delete (text42)

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

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

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

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

   841 ;*******************************************************************

   842 ;(A)-10 for table of site(2) (the last 40 sites), model vs ob

   843 ;*******************************************************************

   845   table_length = 0.95

   847 ; table header name

   848   table_header_name = "Site ID"

   850 ; row (not including header row)

   851   nrow       = 40

   852   row_header = new ((/nrow/),string )

   853   row_header(0:nrow-1) = id81(41:41+nrow-1)

   855 ; Table header

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

   857   x1    = (/0.005,0.15/)        ; Start and end X

   858   y1    = (/0.85 ,0.95/)        ; Start and end Y

   859   text1 = table_header_name

   860   res1               = True

   861   res1@txFontHeightF = 0.015

   862   res1@gsFillColor   = "CornFlowerBlue"

   864 ; Column header1 (equally space in x2)

   866   col_header = (/"Latitude","Longitude"/)

   867   ncol = dimsizes(col_header)

   869   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   870   x2    = (/x1(1),0.35/)        ; start from end of x1

   871   y2    = y1                    ; same as y1

   872   text2 = col_header

   873   res2               = True

   874   res2@txFontHeightF = 0.015

   875   res2@gsFillColor   = "Gray"

   877 ; Column header2 (equally space in x2)

   879 ; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   880   col_header2 = (/"ob",model_name \

   881                  ,"ob",model_name \

   882                  /)

   884   ncol1 = dimsizes(col_header1)

   885   ncol2 = dimsizes(col_header2)

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

   888   x21    = (/x2(1),0.995/)        ; start from end of x2

   889   yhalf  = (y1(0)+y1(1))*0.5

   890   y21    = (/yhalf,y1(1)/)        ; top half of y1

   891   text21 = col_header1

   892   res21               = True

   893   res21@txFontHeightF = 0.015

   894   res21@gsFillColor   = "Gray"

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

   897   x22    = x21                    ; start from end of x1

   898   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   899   text22 = col_header2

   900   res22               = True

   901   res22@txFontHeightF = 0.015

   902   res22@gsFillColor   = "Gray"

   904 ; Row header (equally space in y2)

   905   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   906   x3    = x1                         ; same as x1

   907   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   908   text3 = row_header

   909   res3               = True

   910   res3@txFontHeightF = 0.010

   911   res3@gsFillColor   = "Gray"

   913 ; Main table body-1

   914   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   915   x4    = x2                      ; Start and end x

   916   y4    = y3                      ; Start and end Y

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

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

   920   color_fill4           = "white"

   921 ; color_fill4(:,ncol-1) = "grey"

   922 ; color_fill4(nrow-1,:) = "green"

   924   res4               = True       ; Set up resource list

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

   926   res4@txFontHeightF = 0.015

   927   res4@gsFillColor   = color_fill4

   929   delete (color_fill4)

   931 ; Main table body-2

   932   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   933   x42    = x21                      ; Start and end x

   934   y42    = y3                       ; Start and end Y

   935   text42 = new((/nrow,ncol2/),string)

   937   color_fill42           = new((/nrow,ncol2/),string)

   938   color_fill42           = "white"

   939 ; color_fill42(:,ncol-1) = "grey"

   940 ; color_fill42(nrow-1,:) = "green"

   942   res42               = True       ; Set up resource list

   943 ; res42@gsnDebug      = True       ; Useful to print NDC row,col values used.

   944   res42@txFontHeightF = 0.015

   945   res42@gsFillColor   = color_fill42

   947   delete (color_fill42)

   948 ;-------------------------------------------------------------------

   949 ; for table value

   951   do n = 0,nrow-1

   952      text4(n,0) = sprintf("%5.2f", y81(n+41))

   953      text4(n,1) = sprintf("%5.2f", x81(n+41))

   954   end do

   956   do n = 0,nrow-1

   957      text42(n,0) = sprintf("%5.2f", npp81(n+41))

   958      text42(n,1) = sprintf("%5.2f", nppmod81(n+41))

   959      text42(n,2) = sprintf("%5.2f", rain81(n+41))

   960      text42(n,3) = sprintf("%5.2f", rainmod81(n+41))

   961   end do

   962 ;---------------------------------------------------------------------------

   964   plot_name = "table_site81_model_vs_ob2"

   966   wks = gsn_open_wks (plot_type,plot_name)

   967 ;------------------------------------------

   968 ; for table title

   970   gRes               = True

   971   gRes@txFontHeightF = 0.02

   972 ; gRes@txAngleF      = 90

   974   title_text = "Model vs Observation at 81 Sites (2)"

   976   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   977 ;------------------------------------------

   979   gsn_table(wks,ncr1,x1,y1,text1,res1)

   980   gsn_table(wks,ncr2,x2,y2,text2,res2)

   981   gsn_table(wks,ncr21,x21,y21,text21,res21)

   982   gsn_table(wks,ncr22,x22,y22,text22,res22)

   983   gsn_table(wks,ncr3,x3,y3,text3,res3)

   984   gsn_table(wks,ncr4,x4,y4,text4,res4)

   985   gsn_table(wks,ncr42,x42,y42,text42,res42)

   987   frame(wks)

   988   clear (wks)

   990   delete (col_header)

   991   delete (row_header)

   992   delete (res1)

   993   delete (ncr1)

   994   delete (text1)

   995   delete (res2)

   996   delete (ncr2)

   997   delete (text2)

   998   delete (res21)

   999   delete (ncr21)

  1000   delete (text21)

  1001   delete (res22)

  1002   delete (ncr22)

  1003   delete (text22)

  1004   delete (res3)

  1005   delete (ncr3)

  1006   delete (text3)

  1007   delete (res4)

  1008   delete (ncr4)

  1009   delete (text4)

  1010   delete (res42)

  1011   delete (ncr42)

  1012   delete (text42)

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

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

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

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

  1019 ;**************************************************************************

  1020 ;(B) histogram 81

  1021 ;***********************************************************************-

  1022 ;get data

  1024   RAIN1_1D = ndtooned(rain81)

  1025   RAIN2_1D = ndtooned(rainmod81)

  1026   NPP1_1D  = ndtooned(npp81)

  1027   NPP2_1D  = ndtooned(nppmod81)

  1029   nx = 8

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

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

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

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

  1035   dx4          = new((/1/),float)

  1037   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

  1038          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

  1040 ;----------------------------------------

  1041 ;compute correlation coeff and M score

  1043  u = yvalues(0,:)

  1044  v = yvalues(1,:)

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

  1047  uu = u(good)

  1048  vv = v(good)

  1050  ccr81h = esccr(uu,vv,0)

  1051 ;print (ccr81h)

  1053 ;new eq

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

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

  1056  print (M81h)

  1058  delete (u)

  1059  delete (v)

  1060  delete (uu)

  1061  delete (vv)

  1062 ;----------------------------------------------------------------------

  1063 ; histogram res

  1065   resh                = True

  1066   resh@gsnMaximize    = True

  1067   resh@gsnDraw        = False

  1068   resh@gsnFrame       = False

  1069   resh@xyMarkLineMode = "Markers"

  1070   resh@xyMarkerSizeF  = 0.014

  1071   resh@xyMarker       = 16

  1072   resh@xyMarkerColors = (/"Brown","Blue"/)

  1073   resh@trYMinF        = min(mn_yvalues) - 10.

  1074   resh@trYMaxF        = max(mx_yvalues) + 10.

  1076   resh@tiYAxisString  = "NPP (g C/m2/year)"

  1077   resh@tiXAxisString  = "Precipitation (m/year)"

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

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

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

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

  1084   lnres = True

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

  1086 ;=================================================================

  1087 ;(B)-1 histogram ob 81 site only

  1088 ;

  1089   plot_name = "histogram_ob_81"

  1090   title     = "Observed 81 site"

  1091   resh@tiMainString  = title

  1093   wks   = gsn_open_wks (plot_type,plot_name)

  1095   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

  1097 ;-------------------------------

  1098 ;Attach the vertical bar and the horizontal cap line

  1100   delete (x1)

  1101   delete (x2)

  1102   delete (y1)

  1103   delete (y2)

  1105   do nd=0,0

  1106     lnres@gsLineColor = line_colors(nd)

  1107     do i=0,nx-1

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

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

  1111 ;

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

  1113 ;

  1114         x1 = xvalues(nd,i)

  1115         y1 = yvalues(nd,i)

  1116         y2 = mn_yvalues(nd,i)

  1117         plx = (/x1,x1/)

  1118         ply = (/y1,y2/)

  1119         min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1121         y2 = mx_yvalues(nd,i)

  1122         plx = (/x1,x1/)

  1123         ply = (/y1,y2/)

  1124         max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1125 ;

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

  1127 ;

  1128         x1 = xvalues(nd,i) - dx4

  1129         x2 = xvalues(nd,i) + dx4

  1130         y1 = mn_yvalues(nd,i)

  1131         plx = (/x1,x2/)

  1132         ply = (/y1,y1/)

  1133         min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1135         y1 = mx_yvalues(nd,i)

  1136         plx = (/x1,x2/)

  1137         ply = (/y1,y1/)

  1138         max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1139       end if

  1140     end do

  1141   end do

  1143   draw(xy)

  1144   frame(wks)

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

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

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

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

  1151   clear (wks)

  1152 ;===========================================================================

  1153 ;(B)-2 histogram model vs ob 81 site

  1155   plot_name = "histogram_mod-ob_81"

  1156   title     = model_name+ " vs Observed 81 site"

  1157   resh@tiMainString  = title

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

  1161 ;-----------------------------

  1162 ; Add a boxed legend using the more simple method, which won't have

  1163 ; vertical lines going through the markers.

  1165   resh@pmLegendDisplayMode    = "Always"

  1166 ; resh@pmLegendWidthF         = 0.1

  1167   resh@pmLegendWidthF         = 0.08

  1168   resh@pmLegendHeightF        = 0.05

  1169   resh@pmLegendOrthogonalPosF = -1.17

  1170 ; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

  1171 ; resh@pmLegendParallelPosF   =  0.18

  1172   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1174 ; resh@lgPerimOn              = False

  1175   resh@lgLabelFontHeightF     = 0.015

  1176   resh@xyExplicitLegendLabels = (/"observed",model_name/)

  1177 ;-----------------------------

  1178   tRes  = True

  1179   tRes@txFontHeightF = 0.025

  1181   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"

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

  1185   xy = gsn_csm_xy(wks,xvalues,yvalues,resh)

  1186 ;-------------------------------

  1187 ;Attach the vertical bar and the horizontal cap line

  1189   do nd=0,1

  1190     lnres@gsLineColor = line_colors(nd)

  1191     do i=0,nx-1

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

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

  1195 ;

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

  1197 ;

  1198         x1 = xvalues(nd,i)

  1199         y1 = yvalues(nd,i)

  1200         y2 = mn_yvalues(nd,i)

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

  1203         y2 = mx_yvalues(nd,i)

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

  1205 ;

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

  1207 ;

  1208         x1 = xvalues(nd,i) - dx4

  1209         x2 = xvalues(nd,i) + dx4

  1210         y1 = mn_yvalues(nd,i)

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

  1213         y1 = mx_yvalues(nd,i)

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

  1215       end if

  1216     end do

  1217   end do

  1218   draw(xy)

  1219   frame(wks)

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

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

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

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

  1226  clear (wks)

  1228  delete (RAIN1_1D)

  1229  delete (RAIN2_1D)

  1230  delete (NPP1_1D)

  1231  delete (NPP2_1D)

  1232 ;delete (range)

  1233  delete (xvalues)

  1234  delete (yvalues)

  1235  delete (mn_yvalues)

  1236  delete (mx_yvalues)

  1237  delete (good)

  1238  delete (max_bar)

  1239  delete (min_bar)

  1240  delete (max_cap)

  1241  delete (min_cap)

  1242 ;**************************************************************************

  1243 ;(B) histogram 933

  1245 ;--------------------------------------------------------------------

  1246 ;get data

  1248   RAIN1_1D = ndtooned(rain933)

  1249   RAIN2_1D = ndtooned(rainmod933)

  1250   NPP1_1D  = ndtooned(npp933)

  1251   NPP2_1D  = ndtooned(nppmod933)

  1253   nx = 10

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

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

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

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

  1258   dx4          = new((/1/),float)

  1260   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

  1261          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

  1263 ;----------------------------------------

  1264 ;compute correlation coeff and M score

  1266  u = yvalues(0,:)

  1267  v = yvalues(1,:)

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

  1270  uu = u(good)

  1271  vv = v(good)

  1273  ccr933h = esccr(uu,vv,0)

  1274 ;print (ccr933h)

  1276 ;new eq

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

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

  1279  print (M933h)

  1281  delete (u)

  1282  delete (v)

  1283  delete (uu)

  1284  delete (vv)

  1285 ;----------------------------------------------------------------------

  1286 ; histogram res

  1288   delete (resh)

  1289   resh                = True

  1290   resh@gsnMaximize    = True

  1291   resh@gsnDraw        = False

  1292   resh@gsnFrame       = False

  1293   resh@xyMarkLineMode = "Markers"

  1294   resh@xyMarkerSizeF  = 0.014

  1295   resh@xyMarker       = 16

  1296   resh@xyMarkerColors = (/"Brown","Blue"/)

  1297   resh@trYMinF        = min(mn_yvalues) - 10.

  1298   resh@trYMaxF        = max(mx_yvalues) + 10.

  1300   resh@tiYAxisString  = "NPP (g C/m2/year)"

  1301   resh@tiXAxisString  = "Precipitation (m/year)"

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

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

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

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

  1308   lnres = True

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

  1310 ;=================================================================

  1311 ;(B)-3 histogram ob 933 site only

  1313   plot_name = "histogram_ob_933"

  1314   title     = "Observed 933 site"

  1315   resh@tiMainString  = title

  1317   wks   = gsn_open_wks (plot_type,plot_name)

  1319   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

  1321 ;-------------------------------

  1322 ;Attach the vertical bar and the horizontal cap line

  1324   do nd=0,0

  1325     lnres@gsLineColor = line_colors(nd)

  1326     do i=0,nx-1

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

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

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

  1333         x1 = xvalues(nd,i)

  1334         y1 = yvalues(nd,i)

  1335         y2 = mn_yvalues(nd,i)

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

  1338         y2 = mx_yvalues(nd,i)

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

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

  1343         x1 = xvalues(nd,i) - dx4

  1344         x2 = xvalues(nd,i) + dx4

  1345         y1 = mn_yvalues(nd,i)

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

  1348         y1 = mx_yvalues(nd,i)

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

  1350       end if

  1351     end do

  1352   end do

  1354   draw(xy)

  1355   frame(wks)

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

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

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

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

  1362   delete (xy)

  1363   clear (wks)

  1365 ;===========================================================================

  1366 ;(B)-4 histogram model vs ob 933 site

  1368   plot_name = "histogram_mod-ob_933"

  1369   title     = model_name+ " vs Observed 933 site"

  1370   resh@tiMainString  = title

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

  1374 ;-----------------------------

  1375 ; Add a boxed legend using the more simple method, which won't have

  1376 ; vertical lines going through the markers.

  1378   resh@pmLegendDisplayMode    = "Always"

  1379 ; resh@pmLegendWidthF         = 0.1

  1380   resh@pmLegendWidthF         = 0.08

  1381   resh@pmLegendHeightF        = 0.05

  1382   resh@pmLegendOrthogonalPosF = -1.17

  1383 ; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

  1384 ; resh@pmLegendParallelPosF   =  0.18

  1385   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1387 ; resh@lgPerimOn              = False

  1388   resh@lgLabelFontHeightF     = 0.015

  1389   resh@xyExplicitLegendLabels = (/"observed",model_name/)

  1390 ;-----------------------------

  1391   tRes  = True

  1392   tRes@txFontHeightF = 0.025

  1394   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"

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

  1398   xy = gsn_csm_xy(wks,xvalues,yvalues,resh)

  1399 ;-------------------------------

  1400 ;Attach the vertical bar and the horizontal cap line

  1402   do nd=0,1

  1403     lnres@gsLineColor = line_colors(nd)

  1404     do i=0,nx-1

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

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

  1408 ;

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

  1410 ;

  1411         x1 = xvalues(nd,i)

  1412         y1 = yvalues(nd,i)

  1413         y2 = mn_yvalues(nd,i)

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

  1416         y2 = mx_yvalues(nd,i)

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

  1418 ;

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

  1420 ;

  1421         x1 = xvalues(nd,i) - dx4

  1422         x2 = xvalues(nd,i) + dx4

  1423         y1 = mn_yvalues(nd,i)

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

  1426         y1 = mx_yvalues(nd,i)

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

  1428       end if

  1429     end do

  1430   end do

  1431   draw(xy)

  1432   frame(wks)

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

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

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

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

  1439   delete (xy)

  1440   clear (wks)

  1441 ;***********************************************************************-

  1442 ; global contour

  1444 ;res

  1445   resg                     = True             ; Use plot options

  1446   resg@cnFillOn            = True             ; Turn on color fill

  1447   resg@gsnSpreadColors     = True             ; use full colormap

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

  1449 ; resg@lbLabelAutoStride   = True

  1450   resg@cnLinesOn           = False            ; Turn off contourn lines

  1451   resg@mpFillOn            = False            ; Turn off map fill

  1453   resg@gsnSpreadColors      = True            ; use full colormap

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

  1455   resg@cnMinLevelValF       = 0.              ; Min level

  1456   resg@cnMaxLevelValF       = 2200.           ; Max level

  1457   resg@cnLevelSpacingF      = 200.            ; interval

  1458 ;------------------------------------------------------------------------

  1459 ;(C)-1 global contour ob

  1461   delta = 0.00000000001

  1462   nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)

  1464   plot_name = "global_ob"

  1465   title     = "Observed MODIS MOD 17"

  1466   resg@tiMainString  = title

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

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

  1471   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1472   frame(wks)

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

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

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

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

  1479   clear (wks)

  1480 ;------------------------------------------------------------------------

  1481 ;(C)-2 global contour model

  1483   plot_name = "global_model"

  1484   title     = "Model "+ model_name

  1485   resg@tiMainString  = title

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

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

  1490   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1491   frame(wks)

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

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

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

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

  1498   clear (wks)

  1499 ;------------------------------------------------------------------------

  1500 ;(C)-3 global contour model vs ob

  1502   plot_name = "global_model_vs_ob"

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

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

  1507   delete (plot)

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

  1510   resg@gsnFrame             = False          ; Do not draw plot

  1511   resg@gsnDraw              = False          ; Do not advance frame

  1513 ;(d) compute correlation coef and M score

  1515   uu1 = ndtooned(nppmod)

  1516   vv1 = ndtooned(nppglobe)

  1518   delete (good)

  1519   good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))

  1521   ug = uu1(good)

  1522   vg = vv1(good)

  1524   ccrG = esccr(ug,vg,0)

  1525 ; print (ccrG)

  1527   MG   = (ccrG*ccrG)* 5.0

  1528   print (MG)

  1530 ; plot correlation coef

  1532   gRes  = True

  1533   gRes@txFontHeightF = 0.02

  1534   gRes@txAngleF = 90

  1536   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"

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

  1539 ;--------------------------------------------------------------------

  1541 ;(a) ob

  1543   title     = "Observed MODIS MOD 17"

  1544   resg@tiMainString  = title

  1546   plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1548 ;(b) model

  1550   title     = "Model "+ model_name

  1551   resg@tiMainString  = title

  1553   plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1555 ;(c) model-ob

  1557   zz = nppmod

  1558   zz = nppmod - nppglobe

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

  1561   resg@cnMinLevelValF  = -500           ; Min level

  1562   resg@cnMaxLevelValF  =  500.          ; Max level

  1563   resg@cnLevelSpacingF =  50.           ; interval

  1564   resg@tiMainString    = title

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

  1568   pres                            = True        ; panel plot mods desired

  1569   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1570                                                 ; indiv. plots in panel

  1571   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

  1580   frame (wks)

  1581   clear (wks)

  1583   delete (plot)

  1584 ;**********************************************************************

  1585 ;(D)-1 zonal line plot ob

  1587   resz = True

  1589   vv     = zonalAve(nppglobe)

  1590   vv@long_name = nppglobe@long_name

  1592   plot_name = "zonal_ob"

  1593   title     = "Observed MODIS MOD 17"

  1594   resz@tiMainString  = title

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

  1598   plot  = gsn_csm_xy (wks,ym,vv,resz)

  1599   frame(wks)

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

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

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

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

  1606   clear (wks)

  1607 ;---------------------------------------------------------------------

  1608 ;(D)-2 zonal line plot model vs ob

  1610   s  = new ((/2,dimsizes(ym)/), float)

  1612   s(0,:) = zonalAve(nppglobe)

  1613   s(1,:) = zonalAve(nppmod)

  1615   s@long_name = nppglobe@long_name

  1616 ;-------------------------------------------

  1617 ;(d) compute correlation coef and M score

  1619   ccrZ = esccr(s(1,:), s(0,:),0)

  1620 ; print (ccrZ)

  1622   MZ   = (ccrZ*ccrZ)* 5.0

  1623   print (MZ)

  1624 ;-------------------------------------------

  1625   plot_name = "zonal_model_vs_ob"

  1626   title     = "Zonal Average"

  1627   resz@tiMainString  = title

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

  1631 ; resz@vpHeightF          = 0.4               ; change aspect ratio of plot

  1632 ; resz@vpWidthF           = 0.7

  1634   resz@xyMonoLineColor    = "False"           ; want colored lines

  1635   resz@xyLineColors       = (/"black","red"/) ; colors chosen

  1636 ; resz@xyLineThicknesses  = (/3.,3./)      ; line thicknesses

  1637   resz@xyLineThicknesses  = (/2.,2./)      ; line thicknesses

  1638   resz@xyDashPatterns     = (/0.,0./)      ; make all lines solid

  1640   resz@tiYAxisString      = s@long_name      ; add a axis title

  1641   resz@txFontHeightF      = 0.0195            ; change title font heights

  1643 ; Legent

  1644   resz@pmLegendDisplayMode    = "Always"            ; turn on legend

  1645   resz@pmLegendSide           = "Top"               ; Change location of

  1646 ; resz@pmLegendParallelPosF   = .45                 ; move units right

  1647   resz@pmLegendParallelPosF   = .82                 ; move units right

  1648   resz@pmLegendOrthogonalPosF = -0.4                ; move units down

  1650   resz@pmLegendWidthF         = 0.10                ; Change width and

  1651   resz@pmLegendHeightF        = 0.10                ; height of legend.

  1652   resz@lgLabelFontHeightF     = .02                 ; change font height

  1653 ; resz@lgTitleOn              = True                ; turn on legend title

  1654 ; resz@lgTitleString          = "Example"           ; create legend title

  1655 ; resz@lgTitleFontHeightF     = .025                ; font of legend title

  1656   resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels

  1657 ;--------------------------------------------------------------------

  1658   zRes  = True

  1659   zRes@txFontHeightF = 0.025

  1661   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"

  1663   gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)

  1664 ;--------------------------------------------------------------------

  1666   plot  = gsn_csm_xy (wks,ym,s,resz)       ; create plot

  1668   frame(wks)                                            ; advance frame

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

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

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

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

  1675   clear (wks)

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

  1677   temp_name = "temp." + model_name

  1678   system("mkdir -p " + temp_name)

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

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

  1681 ;------------------------------------------------------------------------

  1682 end

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