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

     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 set_line(lines:string,nline:integer,newlines:string)

    11 begin

    12 ; add line to ascci/html file

    14   nnewlines = dimsizes(newlines)

    15   if(nline+nnewlines-1.ge.dimsizes(lines))

    16     print("set_line: bad index, not setting anything.")

    17     return

    18   end if

    19   lines(nline:nline+nnewlines-1) = newlines

    20 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    21   nline = nline + nnewlines

    22   return

    23 end

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

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

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

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

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

    29                  ,dx4[1]:numeric \

    30                   )

    31 begin

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

    33 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D

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

    36   nbins       = 15     ; Number of bins to use.

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

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

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

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

    43 ; particular range, and then take an average.

    45   nr           = dimsizes(range)

    46   nx           = nr-1

    47 ; print (nx)

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

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

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

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

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

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

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

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

    58   do nd=0,1

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

    62     if(nd.eq.0) then

    63       data     = RAIN1_1D

    64       npp_data = NPP1_1D

    65     else

    66       data     = RAIN2_1D

    67       npp_data = NPP2_1D

    68     end if

    70 ; Loop through each range and check for values.

    72     do i=0,nr-2

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

    74 ;        print("")

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

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

    77       else

    78 ;        print("")

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

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

    81       end if

    83 ; Calculate average, and get min and max.

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

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

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

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

    89         count = dimsizes(idx)

    90       else

    91         count            = 0

    92         yvalues(nd,i)    = yvalues@_FillValue

    93         mn_yvalues(nd,i) = yvalues@_FillValue

    94         mx_yvalues(nd,i) = yvalues@_FillValue

    95       end if

    97 ; Print out information.

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

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

   101 ; Clean up for next time in loop.

   103       delete(idx)

   104     end do

   105     delete(data)

   106     delete(npp_data)

   107   end do

   108 end

   109 ;****************************************************************************

   110 ; Main code.

   111 ;****************************************************************************

   112 begin

   114  plot_type     = "ps"

   115  plot_type_new = "png"

   117 ;************************************************

   118 ; read model data

   119 ;************************************************

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

   122 ;model_name = "06cn"

   123 ;model_grid = "T42"

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

   126 ;model_name = "06casa"

   127 ;model_grid = "T42"

   129 ;film = "i01.10casa_1948-2004_ANN_climo.nc"

   130 ;model_name = "10casa"

   131 ;model_grid = "T42"

   133  film = "i01.10cn_1948-2004_ANN_climo.nc"

   134  model_name = "10cn"

   135  model_grid = "T42"

   137  html_name = "table.html." + model_name

   138  html_new  = html_name +".new"

   139  system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \

   140          "mv -f "+html_new+" "+html_name)

   141 ;--------------------------------------------------

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

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

   145  nppmod0  = fm->NPP

   146  rainmod0 = fm->RAIN

   147  xm       = fm->lon

   148  ym       = fm->lat

   150 ;************************************************

   151 ; read ob data

   152 ;************************************************

   154 ;(1) data at 81 sites

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

   156  fil81 = "data.81.nc"

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

   159  id81   = f81->SITE_ID

   160  npp81  = f81->TNPP_C

   161  rain81 = tofloat(f81->PREC_ANN)

   162  x81    = f81->LONG_DD

   163  y81    = f81->LAT_DD

   165  id81@long_name  = "SITE_ID"

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

   168 ;for model data interpolation

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

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

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

   173     end if

   174  end do

   175 ;print (x81)

   176 ;-------------------------------------

   177 ;(2) data at 933 sites

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

   179  fil933 = "data.933.nc"

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

   182  id933   = f933->SITE_ID

   183  npp933  = f933->TNPP_C

   184  rain933 = f933->PREC

   185  x933    = f933->LONG_DD

   186  y933    = f933->LAT_DD

   188  id933@long_name  = "SITE_ID"

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

   191 ;for model data interpolation

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

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

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

   196     end if

   197  end do

   198 ;print (x933)

   199 ;----------------------------------------

   200 ;(3) global data, interpolated into model grid

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

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

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

   205  nppglobe0 = fglobe->NPP

   206  nppglobe  = nppglobe0

   208 ;***********************************************************************

   209 ; interpolate model data into ob sites

   210 ;***********************************************************************

   212  nppmod   = nppmod0(0,:,:)

   213  rainmod  = rainmod0(0,:,:)

   214  delete (nppmod0)

   215  delete (rainmod0)

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

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

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

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

   225 ;**********************************************************

   226 ; unified units

   227 ;**********************************************************

   228 ; Units for these variables are:

   229 ;

   230 ; rain81  : mm/year

   231 ; rainmod : mm/s

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

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

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

   235 ;

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

   238   nsec_per_year = 60*60*24*365

   240   rain81    = rain81 / 1000.

   241   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   242   nppmod81  = nppmod81 * nsec_per_year

   244   rain933    = rain933 / 1000.

   245   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   246   nppmod933  = nppmod933 * nsec_per_year

   248   nppmod  = nppmod * nsec_per_year

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

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

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

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

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

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

   256   rain81@units     = "m/yr"

   257   rainmod81@units  = "m/yr"

   258   rain933@units    = "m/yr"

   259   rainmod933@units = "m/yr"

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

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

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

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

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

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

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

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

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

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

   272 ;*******************************************************************

   273 ;(A)-1 html table of site81 -- observed

   274 ;*******************************************************************

   275   output_html = "table_site81_ob.html"

   277   header = (/"<HTML>" \

   278             ,"<HEAD>" \

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

   280             ,"</HEAD>" \

   281             ,"<H1>Observation at 81 sites</H1>" \

   282             /)

   283   footer = "</HTML>"

   285   table_header = (/ \

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

   287        ,"<tr>" \

   288        ,"   <th bgcolor=DDDDDD >Site ID</th>" \

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

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

   291        ,"   <th bgcolor=DDDDDD >NPP(gC/m2.year)</th>" \

   292        ,"   <th bgcolor=DDDDDD >RAIN(m/year)</th>" \

   293        ,"</tr>" \

   294        /)

   295   table_footer = "</table>"

   296   row_header = "<tr>"

   297   row_footer = "</tr>"

   299   lines = new(50000,string)

   300   nline = 0

   302   set_line(lines,nline,header)

   303   set_line(lines,nline,table_header)

   304 ;-----------------------------------------------

   305 ;row of table

   307   nrow = dimsizes(id81)

   308   do n = 0,nrow-1

   309      set_line(lines,nline,row_header)

   311      txt1 = sprintf("%5.0f", id81(n))

   312      txt2 = sprintf("%5.2f", y81(n))

   313      txt3 = sprintf("%5.2f", x81(n))

   314      txt4 = sprintf("%5.2f", npp81(n))

   315      txt5 = sprintf("%5.2f", rain81(n))

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

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

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

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

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

   323      set_line(lines,nline,row_footer)

   324   end do

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

   326   set_line(lines,nline,table_footer)

   327   set_line(lines,nline,footer)

   329 ; Now write to an HTML file.

   330   idx = ind(.not.ismissing(lines))

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

   332     asciiwrite(output_html,lines(idx))

   333   else

   334    print ("error?")

   335   end if

   336 ;*******************************************************************

   337 ;(A)-2 html table of site933 -- observed

   338 ;*******************************************************************

   339   output_html = "table_site933_ob.html"

   341   header = (/"<HTML>" \

   342             ,"<HEAD>" \

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

   344             ,"</HEAD>" \

   345             ,"<H1>Observation at 933 sites</H1>" \

   346             /)

   347   footer = "</HTML>"

   349   table_header = (/ \

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

   351        ,"<tr>" \

   352        ,"   <th bgcolor=DDDDDD >Site ID</th>" \

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

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

   355        ,"   <th bgcolor=DDDDDD >NPP(gC/m2.year)</th>" \

   356        ,"   <th bgcolor=DDDDDD >RAIN(m/year)</th>" \

   357        ,"</tr>" \

   358        /)

   359   table_footer = "</table>"

   360   row_header = "<tr>"

   361   row_footer = "</tr>"

   363   delete (lines)

   364   lines = new(50000,string)

   365   nline = 0

   367   set_line(lines,nline,header)

   368   set_line(lines,nline,table_header)

   369 ;-----------------------------------------------

   370 ;row of table

   372   nrow = dimsizes(id933)

   373   do n = 0,nrow-1

   374      set_line(lines,nline,row_header)

   376      txt1 = sprintf("%5.0f", id933(n))

   377      txt2 = sprintf("%5.2f", y933(n))

   378      txt3 = sprintf("%5.2f", x933(n))

   379      txt4 = sprintf("%5.2f", npp933(n))

   380      txt5 = sprintf("%5.2f", rain933(n))

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

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

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

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

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

   388      set_line(lines,nline,row_footer)

   389   end do

   390 ;-----------------------------------------------

   391   set_line(lines,nline,table_footer)

   392   set_line(lines,nline,footer)

   394 ; Now write to an HTML file.

   395   delete (idx)

   396   idx = ind(.not.ismissing(lines))

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

   398     asciiwrite(output_html,lines(idx))

   399   else

   400    print ("error?")

   401   end if

   402 ;*******************************************************************

   403 ;(A)-3 html table of site81 -- model vs ob

   404 ;*******************************************************************

   405   output_html = "table_site81_model_vs_ob.html"

   407   header_text = "<H1>Model "+model_name+" vs Observed at 81 sites</H1>"

   409   header = (/"<HTML>" \

   410             ,"<HEAD>" \

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

   412             ,"</HEAD>" \

   413             ,header_text \

   414             /)

   415   footer = "</HTML>"

   417   delete (table_header)

   418   table_header_text = "   <th bgcolor=DDDDDD >"+model_name+"</th>"

   419   table_header = (/ \

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

   421        ,"<tr>" \

   422        ,"   <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \

   423        ,"   <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \

   424        ,"   <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \

   425        ,"   <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \

   426        ,"   <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \

   427        ,"</tr>" \

   428        ,"<tr>" \

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

   430        ,     table_header_text \

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

   432        ,     table_header_text \

   433        ,"</tr>" \

   434        /)

   435   table_footer = "</table>"

   436   row_header = "<tr>"

   437   row_footer = "</tr>"

   439   delete (lines)

   440   lines = new(50000,string)

   441   nline = 0

   443   set_line(lines,nline,header)

   444   set_line(lines,nline,table_header)

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

   446 ;row of table

   448   nrow = dimsizes(id81)

   449   do n = 0,nrow-1

   450      set_line(lines,nline,row_header)

   452      txt1 = sprintf("%5.0f", id81(n))

   453      txt2 = sprintf("%5.2f", y81(n))

   454      txt3 = sprintf("%5.2f", x81(n))

   455      txt4 = sprintf("%5.2f", npp81(n))

   456      txt5 = sprintf("%5.2f", nppmod81(n))

   457      txt6 = sprintf("%5.2f", rain81(n))

   458      txt7 = sprintf("%5.2f", rainmod81(n))

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

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

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

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

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

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

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

   468      set_line(lines,nline,row_footer)

   469   end do

   470 ;-----------------------------------------------

   471   set_line(lines,nline,table_footer)

   472   set_line(lines,nline,footer)

   474 ; Now write to an HTML file.

   475   delete (idx)

   476   idx = ind(.not.ismissing(lines))

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

   478     asciiwrite(output_html,lines(idx))

   479   else

   480    print ("error?")

   481   end if

   483 ;*******************************************************************

   484 ;(A)-4 html table of site933 -- model vs ob

   485 ;*******************************************************************

   486   output_html = "table_site933_model_vs_ob.html"

   488   header_text = "<H1>Model "+model_name+" vs Observed at 933 sites</H1>"

   490   header = (/"<HTML>" \

   491             ,"<HEAD>" \

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

   493             ,"</HEAD>" \

   494             ,header_text \

   495             /)

   496   footer = "</HTML>"

   498   delete (table_header)

   499   table_header_text = "   <th bgcolor=DDDDDD >"+model_name+"</th>"

   500   table_header = (/ \

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

   502        ,"<tr>" \

   503        ,"   <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \

   504        ,"   <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \

   505        ,"   <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \

   506        ,"   <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \

   507        ,"   <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \

   508        ,"</tr>" \

   509        ,"<tr>" \

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

   511        ,     table_header_text \

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

   513        ,     table_header_text \

   514        ,"</tr>" \

   515        /)

   516   table_footer = "</table>"

   517   row_header = "<tr>"

   518   row_footer = "</tr>"

   520   delete (lines)

   521   lines = new(50000,string)

   522   nline = 0

   524   set_line(lines,nline,header)

   525   set_line(lines,nline,table_header)

   526 ;-----------------------------------------------

   527 ;row of table

   529   nrow = dimsizes(id933)

   530   do n = 0,nrow-1

   531      set_line(lines,nline,row_header)

   533      txt1 = sprintf("%5.0f", id933(n))

   534      txt2 = sprintf("%5.2f", y933(n))

   535      txt3 = sprintf("%5.2f", x933(n))

   536      txt4 = sprintf("%5.2f", npp933(n))

   537      txt5 = sprintf("%5.2f", nppmod933(n))

   538      txt6 = sprintf("%5.2f", rain933(n))

   539      txt7 = sprintf("%5.2f", rainmod933(n))

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

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

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

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

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

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

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

   549      set_line(lines,nline,row_footer)

   550   end do

   551 ;-----------------------------------------------

   552   set_line(lines,nline,table_footer)

   553   set_line(lines,nline,footer)

   555 ; Now write to an HTML file.

   556   delete (idx)

   557   idx = ind(.not.ismissing(lines))

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

   559     asciiwrite(output_html,lines(idx))

   560   else

   561    print ("error?")

   562   end if

   563 ;***************************************************************************

   564 ;(A)-5 scatter plot, model vs ob 81

   565 ;***************************************************************************

   567  plot_name = "scatter_model_vs_ob_81"

   568  title     = model_name +" vs ob 81 sites"

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

   571  res                   = True                  ; plot mods desired

   572  res@tiMainString      = title                 ; add title

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

   574  res@xyMarkers         =  16                   ; choose type of marker

   575  res@xyMarkerColor     = "red"                 ; Marker color

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

   577  res@tmLabelAutoStride = True                  ; nice tick mark labels

   579  res@gsnDraw           = False

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

   581 ;-------------------------------

   582 ;compute correlation coef. and M

   583  ccr81 = esccr(nppmod81,npp81,0)

   584 ;print (ccr81)

   586 ;new eq

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

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

   590  M_npp_S81  = sprintf("%.2f", M81)

   591  system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \

   592         "mv -f "+html_new+" "+html_name)

   593  print (M_npp_S81)

   595  tRes  = True

   596  tRes@txFontHeightF = 0.025

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

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

   601 ;--------------------------------

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

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

   604 ; add polyline

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

   607 ;-------------------------------

   608  draw (plot)

   609  frame(wks)

   610  clear (wks)

   611  delete (dum)

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

   614         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

   615         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

   617 ;***************************************************************************

   618 ;(A)-6 scatter plot, model vs ob 933

   619 ;***************************************************************************

   621  plot_name = "scatter_model_vs_ob_933"

   622  title     = model_name +" vs ob 933 sites"

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

   625  res                   = True                  ; plot mods desired

   626  res@tiMainString      = title                 ; add title

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

   628  res@xyMarkers         =  16                   ; choose type of marker

   629  res@xyMarkerColor     = "red"                 ; Marker color

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

   631  res@tmLabelAutoStride = True                  ; nice tick mark labels

   633  res@gsnDraw           = False

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

   635 ;-------------------------------

   636 ;compute correlation coef. and M

   637  ccr933 = esccr(nppmod933,npp933,0)

   638 ;print (ccr933)

   640 ;new eq

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

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

   644  M_npp_S933  = sprintf("%.2f", M933)

   645  system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \

   646         "mv -f "+html_new+" "+html_name)

   647  print (M_npp_S933)

   649  tRes  = True

   650  tRes@txFontHeightF = 0.025

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

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

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

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

   657 ;-------------------------------

   658 ; add polyline

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

   661 ;-------------------------------

   662  draw (plot)

   663  frame(wks)

   664  clear (wks)

   665  delete (dum)

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

   668         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

   669         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

   671 ;**************************************************************************

   672 ;(B) histogram 81

   673 ;**************************************************************************

   674 ;get data

   675   RAIN1_1D = ndtooned(rain81)

   676   RAIN2_1D = ndtooned(rainmod81)

   677   NPP1_1D  = ndtooned(npp81)

   678   NPP2_1D  = ndtooned(nppmod81)

   680 ; number of bin

   681   nx = 8

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

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

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

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

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

   689   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

   690          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

   692 ;----------------------------------------

   693 ;compute correlation coeff and M score

   695  u = yvalues(0,:)

   696  v = yvalues(1,:)

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

   699  uu = u(good)

   700  vv = v(good)

   702  ccr81h = esccr(uu,vv,0)

   703 ;print (ccr81h)

   705 ;new eq

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

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

   709  M_npp_H81  = sprintf("%.2f", M81h)

   710  system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \

   711         "mv -f "+html_new+" "+html_name)

   712  print (M_npp_H81)

   714  delete (u)

   715  delete (v)

   716  delete (uu)

   717  delete (vv)

   718 ;----------------------------------------------------------------------

   719 ; histogram res

   720   resh                = True

   721   resh@gsnMaximize    = True

   722   resh@gsnDraw        = False

   723   resh@gsnFrame       = False

   724   resh@xyMarkLineMode = "Markers"

   725   resh@xyMarkerSizeF  = 0.014

   726   resh@xyMarker       = 16

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

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

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

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

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

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

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

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

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

   739   lnres = True

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

   742 ;**************************************************************************

   743 ;(B)-1 histogram plot, ob 81 site

   744 ;**************************************************************************

   746   plot_name = "histogram_ob_81"

   747   title     = "Observed 81 site"

   748   resh@tiMainString  = title

   750   wks   = gsn_open_wks (plot_type,plot_name)

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

   754 ;-------------------------------

   755 ;Attach the vertical bar and the horizontal cap line

   757   do nd=0,0

   758     lnres@gsLineColor = line_colors(nd)

   759     do i=0,nx-1

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

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

   763 ;

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

   765 ;

   766         x1 = xvalues(nd,i)

   767         y1 = yvalues(nd,i)

   768         y2 = mn_yvalues(nd,i)

   769         plx = (/x1,x1/)

   770         ply = (/y1,y2/)

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

   773         y2 = mx_yvalues(nd,i)

   774         plx = (/x1,x1/)

   775         ply = (/y1,y2/)

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

   777 ;

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

   779 ;

   780         x1 = xvalues(nd,i) - dx4

   781         x2 = xvalues(nd,i) + dx4

   782         y1 = mn_yvalues(nd,i)

   783         plx = (/x1,x2/)

   784         ply = (/y1,y1/)

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

   787         y1 = mx_yvalues(nd,i)

   788         plx = (/x1,x2/)

   789         ply = (/y1,y1/)

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

   791       end if

   792     end do

   793   end do

   795   draw(xy)

   796   frame(wks)

   797   clear (wks)

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

   800         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

   801         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

   803 ;****************************************************************************

   804 ;(B)-2 histogram plot, model vs ob 81 site

   805 ;****************************************************************************

   807   plot_name = "histogram_model_vs_ob_81"

   808   title     = model_name+ " vs Observed 81 site"

   809   resh@tiMainString  = title

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

   813 ;-----------------------------

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

   815 ; vertical lines going through the markers.

   817   resh@pmLegendDisplayMode    = "Always"

   818 ; resh@pmLegendWidthF         = 0.1

   819   resh@pmLegendWidthF         = 0.08

   820   resh@pmLegendHeightF        = 0.05

   821   resh@pmLegendOrthogonalPosF = -1.17

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

   823 ; resh@pmLegendParallelPosF   =  0.18

   824   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

   826 ; resh@lgPerimOn              = False

   827   resh@lgLabelFontHeightF     = 0.015

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

   829 ;-----------------------------

   830   tRes  = True

   831   tRes@txFontHeightF = 0.025

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

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

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

   838 ;-------------------------------

   839 ;Attach the vertical bar and the horizontal cap line

   841   do nd=0,1

   842     lnres@gsLineColor = line_colors(nd)

   843     do i=0,nx-1

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

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

   847 ;

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

   849 ;

   850         x1 = xvalues(nd,i)

   851         y1 = yvalues(nd,i)

   852         y2 = mn_yvalues(nd,i)

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

   855         y2 = mx_yvalues(nd,i)

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

   857 ;

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

   859 ;

   860         x1 = xvalues(nd,i) - dx4

   861         x2 = xvalues(nd,i) + dx4

   862         y1 = mn_yvalues(nd,i)

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

   865         y1 = mx_yvalues(nd,i)

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

   867       end if

   868     end do

   869   end do

   871   draw(xy)

   872   frame(wks)

   873   clear(wks)

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

   876         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

   877         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

   880  delete (RAIN1_1D)

   881  delete (RAIN2_1D)

   882  delete (NPP1_1D)

   883  delete (NPP2_1D)

   884 ;delete (range)

   885  delete (xvalues)

   886  delete (yvalues)

   887  delete (mn_yvalues)

   888  delete (mx_yvalues)

   889  delete (good)

   890  delete (max_bar)

   891  delete (min_bar)

   892  delete (max_cap)

   893  delete (min_cap)

   895 ;**************************************************************************

   896 ;(B) histogram 933

   897 ;**************************************************************************

   899 ; get data

   900   RAIN1_1D = ndtooned(rain933)

   901   RAIN2_1D = ndtooned(rainmod933)

   902   NPP1_1D  = ndtooned(npp933)

   903   NPP2_1D  = ndtooned(nppmod933)

   905 ; number of bin

   906   nx = 10

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

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

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

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

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

   914   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

   915          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

   917 ;----------------------------------------

   918 ;compute correlation coeff and M score

   920  u = yvalues(0,:)

   921  v = yvalues(1,:)

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

   924  uu = u(good)

   925  vv = v(good)

   927  ccr933h = esccr(uu,vv,0)

   928 ;print (ccr933h)

   930 ;new eq

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

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

   934  M_npp_H933 = sprintf("%.2f", M933h)

   935  system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \

   936         "mv -f "+html_new+" "+html_name)

   937  print (M_npp_H933)

   939  delete (u)

   940  delete (v)

   941  delete (uu)

   942  delete (vv)

   943 ;----------------------------------------------------------------------

   944 ; histogram res

   945   delete (resh)

   946   resh                = True

   947   resh@gsnMaximize    = True

   948   resh@gsnDraw        = False

   949   resh@gsnFrame       = False

   950   resh@xyMarkLineMode = "Markers"

   951   resh@xyMarkerSizeF  = 0.014

   952   resh@xyMarker       = 16

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

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

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

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

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

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

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

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

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

   965   lnres = True

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

   968 ;**************************************************************************

   969 ;(B)-3 histogram plot, ob 933 site

   970 ;**************************************************************************

   972   plot_name = "histogram_ob_933"

   973   title     = "Observed 933 site"

   974   resh@tiMainString  = title

   976   wks   = gsn_open_wks (plot_type,plot_name)

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

   980 ;-------------------------------

   981 ;Attach the vertical bar and the horizontal cap line

   983   do nd=0,0

   984     lnres@gsLineColor = line_colors(nd)

   985     do i=0,nx-1

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

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

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

   992         x1 = xvalues(nd,i)

   993         y1 = yvalues(nd,i)

   994         y2 = mn_yvalues(nd,i)

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

   997         y2 = mx_yvalues(nd,i)

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

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

  1002         x1 = xvalues(nd,i) - dx4

  1003         x2 = xvalues(nd,i) + dx4

  1004         y1 = mn_yvalues(nd,i)

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

  1007         y1 = mx_yvalues(nd,i)

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

  1009       end if

  1010     end do

  1011   end do

  1013   draw(xy)

  1014   frame(wks)

  1015   delete (xy)

  1016   clear (wks)

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

  1019         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1020         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1022 ;****************************************************************************

  1023 ;(B)-4 histogram plot, model vs ob 933 site

  1024 ;****************************************************************************

  1026   plot_name = "histogram_model_vs_ob_933"

  1027   title     = model_name+ " vs Observed 933 site"

  1028   resh@tiMainString  = title

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

  1032 ;-----------------------------

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

  1034 ; vertical lines going through the markers.

  1036   resh@pmLegendDisplayMode    = "Always"

  1037 ; resh@pmLegendWidthF         = 0.1

  1038   resh@pmLegendWidthF         = 0.08

  1039   resh@pmLegendHeightF        = 0.05

  1040   resh@pmLegendOrthogonalPosF = -1.17

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

  1042 ; resh@pmLegendParallelPosF   =  0.18

  1043   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1045 ; resh@lgPerimOn              = False

  1046   resh@lgLabelFontHeightF     = 0.015

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

  1048 ;-----------------------------

  1049   tRes  = True

  1050   tRes@txFontHeightF = 0.025

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

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

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

  1057 ;-------------------------------

  1058 ;Attach the vertical bar and the horizontal cap line

  1060   do nd=0,1

  1061     lnres@gsLineColor = line_colors(nd)

  1062     do i=0,nx-1

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

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

  1066 ;

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

  1068 ;

  1069         x1 = xvalues(nd,i)

  1070         y1 = yvalues(nd,i)

  1071         y2 = mn_yvalues(nd,i)

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

  1074         y2 = mx_yvalues(nd,i)

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

  1076 ;

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

  1078 ;

  1079         x1 = xvalues(nd,i) - dx4

  1080         x2 = xvalues(nd,i) + dx4

  1081         y1 = mn_yvalues(nd,i)

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

  1084         y1 = mx_yvalues(nd,i)

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

  1086       end if

  1087     end do

  1088   end do

  1090   draw(xy)

  1091   frame(wks)

  1092   delete(xy)

  1093   clear(wks)

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

  1096         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1097         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1099 ;***************************************************************************

  1100 ;(C) global contour

  1101 ;***************************************************************************

  1103 ;res

  1104   resg                     = True             ; Use plot options

  1105   resg@cnFillOn            = True             ; Turn on color fill

  1106   resg@gsnSpreadColors     = True             ; use full colormap

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

  1108 ; resg@lbLabelAutoStride   = True

  1109   resg@cnLinesOn           = False            ; Turn off contourn lines

  1110   resg@mpFillOn            = False            ; Turn off map fill

  1112   resg@gsnSpreadColors      = True            ; use full colormap

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

  1114   resg@cnMinLevelValF       = 0.              ; Min level

  1115   resg@cnMaxLevelValF       = 2200.           ; Max level

  1116   resg@cnLevelSpacingF      = 200.            ; interval

  1118 ;****************************************************************************

  1119 ;(C)-1 global contour plot, ob

  1120 ;****************************************************************************

  1122   delta = 0.00000000001

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

  1125   plot_name = "global_ob"

  1126   title     = "Observed MODIS MOD 17"

  1127   resg@tiMainString  = title

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

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

  1132   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1134   frame(wks)

  1135   clear (wks)

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

  1138         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1139         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1141 ;****************************************************************************

  1142 ;(C)-2 global contour plot, model

  1143 ;****************************************************************************

  1145   plot_name = "global_model"

  1146   title     = "Model "+ model_name

  1147   resg@tiMainString  = title

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

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

  1152   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1154   frame(wks)

  1155   clear (wks)

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

  1158         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1159         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1161 ;****************************************************************************

  1162 ;(C)-3 global contour plot, model vs ob

  1163 ;****************************************************************************

  1165   plot_name = "global_model_vs_ob"

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

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

  1170   delete (plot)

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

  1173   resg@gsnFrame             = False          ; Do not draw plot

  1174   resg@gsnDraw              = False          ; Do not advance frame

  1176 ; compute correlation coef and M score

  1178   uu1 = ndtooned(nppmod)

  1179   vv1 = ndtooned(nppglobe)

  1181   delete (good)

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

  1184   ug = uu1(good)

  1185   vg = vv1(good)

  1187   ccrG = esccr(ug,vg,0)

  1188 ; print (ccrG)

  1190   MG   = (ccrG*ccrG)* 5.0

  1192   M_npp_G = sprintf("%.2f", MG)

  1193   system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \

  1194          "mv -f "+html_new+" "+html_name)

  1195   print (M_npp_G)

  1197 ; plot correlation coef

  1199   gRes  = True

  1200   gRes@txFontHeightF = 0.02

  1201   gRes@txAngleF = 90

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

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

  1206 ;--------------------------------------------------------------------

  1207 ;(a) ob

  1209   title     = "Observed MODIS MOD 17"

  1210   resg@tiMainString  = title

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

  1214 ;(b) model

  1216   title     = "Model "+ model_name

  1217   resg@tiMainString  = title

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

  1221 ;(c) model-ob

  1223   zz = nppmod

  1224   zz = nppmod - nppglobe

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

  1227   resg@cnMinLevelValF  = -500           ; Min level

  1228   resg@cnMaxLevelValF  =  500.          ; Max level

  1229   resg@cnLevelSpacingF =  50.           ; interval

  1230   resg@tiMainString    = title

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

  1234   pres                            = True        ; panel plot mods desired

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

  1236                                                 ; indiv. plots in panel

  1237   pres@gsnMaximize                = True        ; fill the page

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

  1241   frame (wks)

  1242   clear (wks)

  1243   delete (plot)

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

  1246         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1247         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1249 ;***************************************************************************

  1250 ;(D)-1 zonal line plot, ob

  1251 ;***************************************************************************

  1253   vv     = zonalAve(nppglobe)

  1254   vv@long_name = nppglobe@long_name

  1256   plot_name = "zonal_ob"

  1257   title     = "Observed MODIS MOD 17"

  1259   resz = True

  1260   resz@tiMainString  = title

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

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

  1266   frame(wks)

  1267   clear (wks)

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

  1270         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1271         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1273 ;****************************************************************************

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

  1275 ;****************************************************************************

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

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

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

  1282   s@long_name = nppglobe@long_name

  1283 ;-------------------------------------------

  1284 ; compute correlation coef and M score

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

  1287 ; print (ccrZ)

  1289   MZ   = (ccrZ*ccrZ)* 5.0

  1291   M_npp_Z = sprintf("%.2f", MZ)

  1292   system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \

  1293          "mv -f "+html_new+" "+html_name)

  1294   print (M_npp_Z)

  1295 ;-------------------------------------------

  1296   plot_name = "zonal_model_vs_ob"

  1297   title     = "Zonal Average"

  1298   resz@tiMainString  = title

  1300   wks = gsn_open_wks (plot_type,plot_name)

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

  1303 ; resz@vpWidthF           = 0.7

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

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

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

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

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

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

  1312   resz@txFontHeightF      = 0.0195            ; change title font heights

  1314 ; Legent

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

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

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

  1318   resz@pmLegendParallelPosF   = .82           ; move units right

  1319   resz@pmLegendOrthogonalPosF = -0.4          ; move units down

  1321   resz@pmLegendWidthF         = 0.10          ; Change width and

  1322   resz@pmLegendHeightF        = 0.10          ; height of legend.

  1323   resz@lgLabelFontHeightF     = .02           ; change font height

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

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

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

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

  1328 ;--------------------------------------------------------------------

  1329   zRes  = True

  1330   zRes@txFontHeightF = 0.025

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

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

  1335 ;--------------------------------------------------------------------

  1337   plot  = gsn_csm_xy (wks,ym,s,resz)

  1339   frame(wks)

  1340   clear (wks)

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

  1343         "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

  1344         "rm "+plot_name+"-*."+plot_type_new+";"+ \

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

  1346 ;***************************************************************************

  1347 ; tar up output plots

  1348 ;***************************************************************************

  1350   temp_name = "npp." + model_name

  1351   system("mkdir -p " + temp_name+";"+ \

  1352          "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \

  1353          "mv table_*.html " + temp_name +";"+ \

  1354          "mv *.png " + temp_name +";"+ \

  1355          "tar cf "+ temp_name +".tar " + temp_name)

  1356 ;***************************************************************************

  1357 end

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