clamp: all/01.npp.ncl@4be95183fbcd

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

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

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

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

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

     6 procedure set_line(lines:string,nline:integer,newlines:string)

     7 begin

     8 ; add line to ascci/html file

    10   nnewlines = dimsizes(newlines)

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

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

    13     return

    14   end if

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

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

    17   nline = nline + nnewlines

    18   return

    19 end

    20 ;****************************************************************************

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

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

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

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

    25                  ,dx4[1]:numeric \

    26                   )

    27 begin

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

    29 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D

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

    32   nbins       = 15     ; Number of bins to use.

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

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

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

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

    39 ; particular range, and then take an average.

    41   nr           = dimsizes(range)

    42   nx           = nr-1

    43 ; print (nx)

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

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

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

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

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

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

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

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

    54   do nd=0,1

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

    58     if(nd.eq.0) then

    59       data     = RAIN1_1D

    60       npp_data = NPP1_1D

    61     else

    62       data     = RAIN2_1D

    63       npp_data = NPP2_1D

    64     end if

    66 ; Loop through each range and check for values.

    68     do i=0,nr-2

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

    70 ;        print("")

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

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

    73       else

    74 ;        print("")

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

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

    77       end if

    79 ; Calculate average, and get min and max.

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

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

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

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

    85         count = dimsizes(idx)

    86       else

    87         count            = 0

    88         yvalues(nd,i)    = yvalues@_FillValue

    89         mn_yvalues(nd,i) = yvalues@_FillValue

    90         mx_yvalues(nd,i) = yvalues@_FillValue

    91       end if

    93 ; Print out information.

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

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

    97 ; Clean up for next time in loop.

    99       delete(idx)

   100     end do

   101     delete(data)

   102     delete(npp_data)

   103   end do

   104 end

   105 ;****************************************************************************

   106 ; Main code.

   107 ;****************************************************************************

   108 begin

   110  plot_type     = "ps"

   111  plot_type_new = "png"

   113 ;-----------------------------------------------------

   114 ; edit table.html of current model for movel1_vs_model2

   116  if (isvar("compare")) then

   117     html_name2 = compare+"/table.html"

   118     html_new2  = html_name2 +".new"

   119  end if

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

   121 ; edit table.html for current model

   123  html_name = model_name+"/table.html"

   124  html_new  = html_name +".new"

   126 ;------------------------------------------------------

   127 ; read model data

   129  fm   = addfile (dirm+film1,"r")

   131  nppmod0  = fm->NPP

   132  rainmod0 = fm->RAIN

   133  xm       = fm->lon

   134  ym       = fm->lat

   136  delete (fm)

   138 ;----------------------------------------------------

   139 ; read ob data

   141 ;(1) data at 81 sites

   143  dir81 = diro + "npp/"

   144  fil81 = "data.81.nc"

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

   147  id81   = f81->SITE_ID

   148  npp81  = f81->TNPP_C

   149  rain81 = tofloat(f81->PREC_ANN)

   150  x81    = f81->LONG_DD

   151  y81    = f81->LAT_DD

   153  delete (f81)

   155  id81@long_name  = "SITE_ID"

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

   158 ;for model data interpolation

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

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

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

   163     end if

   164  end do

   166 ;-------------------------------------

   167 ;(2) data at 933 sites

   169  dir933 = diro + "npp/"

   170  fil933 = "data.933.nc"

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

   173  id933   = f933->SITE_ID

   174  npp933  = f933->TNPP_C

   175  rain933 = f933->PREC

   176  x933    = f933->LONG_DD

   177  y933    = f933->LAT_DD

   179  delete (f933)

   181  id933@long_name  = "SITE_ID"

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

   184 ;for model data interpolation

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

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

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

   189     end if

   190  end do

   192 ;----------------------------------------

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

   195  dirg = diro + "npp/"

   196  ;filg = "Npp_"+model_grid+"_mean.nc"

   197  filg = "npp_"+model_grid+"_mean_2000-2004.nc"

   198  fglobe   = addfile (dirg+filg,"r")

   200  nppglobe0 = fglobe->NPP

   201  nppglobe  = nppglobe0

   203  delete (fglobe)

   205 ;***********************************************************************

   206 ; interpolate model data into ob sites

   207 ;***********************************************************************

   209  nppmod   = nppmod0(0,:,:)

   210  rainmod  = rainmod0(0,:,:)

   211  delete (nppmod0)

   212  delete (rainmod0)

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

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

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

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

   222 ;**********************************************************

   223 ; unified units

   224 ;**********************************************************

   225 ; Units for these variables are:

   226 ;

   227 ; rain81  : mm/year

   228 ; rainmod : mm/s

   229 ; npp81   : gC/m^2/year

   230 ; nppmod81: gC/m^2/s

   231 ; nppglobe: gC/m^2/year

   232 ;

   233 ; We want to convert these to "m/year" and "gC/m^2/year".

   235   nsec_per_year = 60*60*24*365

   237   rain81    = rain81 / 1000.

   238   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   239   nppmod81  = nppmod81 * nsec_per_year

   241   rain933    = rain933 / 1000.

   242   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   243   nppmod933  = nppmod933 * nsec_per_year

   245   nppmod  = nppmod * nsec_per_year

   247   npp81@units      = "gC m~S~-2~N~ y~S~-1~N~"

   248   nppmod81@units   = "gC m~S~-2~N~ y~S~-1~N~"

   249   npp933@units     = "gC m~S~-2~N~ y~S~-1~N~"

   250   nppmod933@units  = "gC m~S~-2~N~ y~S~-1~N~"

   251   nppmod@units     = "gC m~S~-2~N~ y~S~-1~N~"

   252   nppglobe@units   = "gC m~S~-2~N~ y~S~-1~N~"

   253   rain81@units     = "m y~S~-1~N~"

   254   rainmod81@units  = "m y~S~-1~N~"

   255   rain933@units    = "m y~S~-1~N~"

   256   rainmod933@units = "m y~S~-1~N~"

   258   npp81@long_name      = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"

   259   npp933@long_name     = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"

   260   nppmod81@long_name   = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"

   261   nppmod933@long_name  = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"

   262   nppmod@long_name     = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"

   263   nppglobe@long_name   = "NPP (gC m~S~-2~N~ y~S~-1~N~)"

   264   rain81@long_name     = "PREC (m y~S~-1~N~)"

   265   rain933@long_name    = "PREC (m y~S~-1~N~)"

   266   rainmod81@long_name  = "PREC (m y~S~-1~N~)"

   267   rainmod933@long_name = "PREC (m y~S~-1~N~)"

   269 ; change longitude from 0-360 to -180-180

   270   x81  = where(x81  .gt. 180., x81 -360., x81 )

   271   x933 = where(x933 .gt. 180., x933-360., x933)

   273 ;*******************************************************************

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

   275 ;*******************************************************************

   277   output_html = "table_site81_ob.html"

   279   header = (/"<HTML>" \

   280             ,"<HEAD>" \

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

   282             ,"</HEAD>" \

   283             ,"<H1>EMDI Observations Class A (81 sites)</H1>" \

   284             /)

   285   footer = "</HTML>"

   287   table_header = (/ \

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

   289        ,"<tr>" \

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

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

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

   293        ,"   <th bgcolor=DDDDDD >NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \

   294        ,"   <th bgcolor=DDDDDD >PPT (m y<small><sup>-1</sup></small>)</th>" \

   295        ,"</tr>" \

   296        /)

   297   table_footer = "</table>"

   298   row_header = "<tr>"

   299   row_footer = "</tr>"

   301   lines = new(50000,string)

   302   nline = 0

   304   set_line(lines,nline,header)

   305   set_line(lines,nline,table_header)

   306 ;-----------------------------------------------

   307 ;row of table

   309   nrow = dimsizes(id81)

   310   do n = 0,nrow-1

   311      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

   325      set_line(lines,nline,row_footer)

   326   end do

   327 ;-----------------------------------------------

   328   set_line(lines,nline,table_footer)

   329   set_line(lines,nline,footer)

   331 ; Now write to an HTML file.

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

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

   334     asciiwrite(output_html,lines(idx))

   335   else

   336    print ("error?")

   337   end if

   339 ;*******************************************************************

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

   341 ;*******************************************************************

   342   output_html = "table_site933_ob.html"

   344   header = (/"<HTML>" \

   345             ,"<HEAD>" \

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

   347             ,"</HEAD>" \

   348             ,"<H1>EMDI Observations Class B (933 sites)</H1>" \

   349             /)

   350   footer = "</HTML>"

   352   table_header = (/ \

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

   354        ,"<tr>" \

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

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

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

   358        ,"   <th bgcolor=DDDDDD >NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \

   359        ,"   <th bgcolor=DDDDDD >PPT (m y<small><sup>-1</sup></small>)</th>" \

   360        ,"</tr>" \

   361        /)

   362   table_footer = "</table>"

   363   row_header = "<tr>"

   364   row_footer = "</tr>"

   366   delete (lines)

   367   lines = new(50000,string)

   368   nline = 0

   370   set_line(lines,nline,header)

   371   set_line(lines,nline,table_header)

   372 ;-----------------------------------------------

   373 ;row of table

   375   nrow = dimsizes(id933)

   376   do n = 0,nrow-1

   377      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

   391      set_line(lines,nline,row_footer)

   392   end do

   393 ;-----------------------------------------------

   394   set_line(lines,nline,table_footer)

   395   set_line(lines,nline,footer)

   397 ; Now write to an HTML file.

   398   delete (idx)

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

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

   401     asciiwrite(output_html,lines(idx))

   402   else

   403    print ("error?")

   404   end if

   406 ;*******************************************************************

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

   408 ;*******************************************************************

   409   output_html = "table_site81_model_vs_ob.html"

   411   header_text = "<H1>Model "+model_name+" vs Class A Observations (81 sites)</H1>"

   413   header = (/"<HTML>" \

   414             ,"<HEAD>" \

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

   416             ,"</HEAD>" \

   417             ,header_text \

   418             /)

   419   footer = "</HTML>"

   421   delete (table_header)

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

   423   table_header = (/ \

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

   425        ,"<tr>" \

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

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

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

   429        ,"   <th bgcolor=DDDDDD colspan=2>NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \

   430        ,"   <th bgcolor=DDDDDD colspan=2>RAIN (m y<small><sup>-1</sup></small>)</th>" \

   431        ,"</tr>" \

   432        ,"<tr>" \

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

   434        ,     table_header_text \

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

   436        ,     table_header_text \

   437        ,"</tr>" \

   438        /)

   439   table_footer = "</table>"

   440   row_header = "<tr>"

   441   row_footer = "</tr>"

   443   delete (lines)

   444   lines = new(50000,string)

   445   nline = 0

   447   set_line(lines,nline,header)

   448   set_line(lines,nline,table_header)

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

   450 ;row of table

   452   nrow = dimsizes(id81)

   453   do n = 0,nrow-1

   454      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   472      set_line(lines,nline,row_footer)

   473   end do

   474 ;-----------------------------------------------

   475   set_line(lines,nline,table_footer)

   476   set_line(lines,nline,footer)

   478 ; Now write to an HTML file.

   479   delete (idx)

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

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

   482     asciiwrite(output_html,lines(idx))

   483   else

   484    print ("error?")

   485   end if

   487 ;*******************************************************************

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

   489 ;*******************************************************************

   490   output_html = "table_site933_model_vs_ob.html"

   492   header_text = "<H1>Model "+model_name+" vs Class B Observations (933 sites)</H1>"

   494   header = (/"<HTML>" \

   495             ,"<HEAD>" \

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

   497             ,"</HEAD>" \

   498             ,header_text \

   499             /)

   500   footer = "</HTML>"

   502   delete (table_header)

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

   504   table_header = (/ \

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

   506        ,"<tr>" \

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

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

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

   510        ,"   <th bgcolor=DDDDDD colspan=2>NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \

   511        ,"   <th bgcolor=DDDDDD colspan=2>RAIN (m y<small><sup>-1</sup></small>)</th>" \

   512        ,"</tr>" \

   513        ,"<tr>" \

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

   515        ,     table_header_text \

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

   517        ,     table_header_text \

   518        ,"</tr>" \

   519        /)

   520   table_footer = "</table>"

   521   row_header = "<tr>"

   522   row_footer = "</tr>"

   524   delete (lines)

   525   lines = new(50000,string)

   526   nline = 0

   528   set_line(lines,nline,header)

   529   set_line(lines,nline,table_header)

   530 ;-----------------------------------------------

   531 ;row of table

   533   nrow = dimsizes(id933)

   534   do n = 0,nrow-1

   535      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   553      set_line(lines,nline,row_footer)

   554   end do

   555 ;-----------------------------------------------

   556   set_line(lines,nline,table_footer)

   557   set_line(lines,nline,footer)

   559 ; Now write to an HTML file.

   560   delete (idx)

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

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

   563     asciiwrite(output_html,lines(idx))

   564   else

   565    print ("error?")

   566   end if

   568 ;***************************************************************************

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

   570 ;***************************************************************************

   572  plot_name = "scatter_model_vs_ob_81"

   573  ;title     = model_name +" vs Class A observations (81 sites)"

   574  title     = model_name +" (1975-2000) vs Class A observations (81 sites)"

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

   577  res                   = True                  ; plot mods desired

   578  res@tiMainString      = title                 ; add title

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

   580  res@xyMarkers         =  16                   ; choose type of marker

   581  res@xyMarkerColor     = "red"                 ; Marker color

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

   583  res@tmLabelAutoStride = True                  ; nice tick mark labels

   585  res@gsnDraw           = False

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

   587 ;-------------------------------

   588 ;compute correlation coef. and M

   589  ccr81 = esccr(nppmod81,npp81,0)

   590 ;print (ccr81)

   592  score_max = 1.0

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

   595  M81s = (1. - (bias/dimsizes(y81)))*score_max

   596  M_npp_S81  = sprintf("%.2f", M81s)

   598  if (isvar("compare")) then

   599     system("sed -e '1,/M_npp_S81/s/M_npp_S81/"+M_npp_S81+"/' "+html_name2+" > "+html_new2+";"+ \

   600            "mv -f "+html_new2+" "+html_name2)

   601  end if

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

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

   606  tRes  = True

   607  tRes@txFontHeightF = 0.025

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

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

   612 ;--------------------------------

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

   614 ;-------------------------------

   615 ; add polyline

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

   618 ;-------------------------------

   619  draw (plot)

   620  delete (wks)

   621  delete (dum)

   623  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

   625  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   627 ;***************************************************************************

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

   629 ;***************************************************************************

   631  plot_name = "scatter_model_vs_ob_933"

   632  title     = model_name +" (1975-2000) vs Class B Observations (933 sites)"

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

   635  res                   = True                  ; plot mods desired

   636  res@tiMainString      = title                 ; add title

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

   638  res@xyMarkers         =  16                   ; choose type of marker

   639  res@xyMarkerColor     = "red"                 ; Marker color

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

   641  res@tmLabelAutoStride = True                  ; nice tick mark labels

   643  res@gsnDraw           = False

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

   645 ;-------------------------------

   646 ;compute correlation coef. and M

   648  ccr933 = esccr(nppmod933,npp933,0)

   650  score_max = 1.0

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

   653  M933s = (1. - (bias/dimsizes(y933)))*score_max

   654  M_npp_S933  = sprintf("%.2f", M933s)

   656  if (isvar("compare")) then

   657     system("sed -e '1,/M_npp_S933/s/M_npp_S933/"+M_npp_S933+"/' "+html_name2+" > "+html_new2+";"+ \

   658            "mv -f "+html_new2+" "+html_name2)

   659  end if

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

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

   664  tRes  = True

   665  tRes@txFontHeightF = 0.025

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

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

   670 ;--------------------------------

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

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

   673 ; add polyline

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

   676 ;-------------------------------

   677  draw (plot)

   678  delete (wks)

   679  delete (dum)

   681  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

   683  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   685 ;**************************************************************************

   686 ;(B) histogram 81

   687 ;**************************************************************************

   688 ;get data

   689   RAIN1_1D = ndtooned(rain81)

   690   RAIN2_1D = ndtooned(rainmod81)

   691   NPP1_1D  = ndtooned(npp81)

   692   NPP2_1D  = ndtooned(nppmod81)

   694 ; number of bin

   695   nx = 8

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

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

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

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

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

   703   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

   704          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

   706 ;----------------------------------------

   707 ;compute correlation coeff and M score

   709  u = yvalues(0,:)

   710  v = yvalues(1,:)

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

   713  uu = u(good)

   714  vv = v(good)

   716  ccr81h = esccr(uu,vv,0)

   718  score_max = 2.0

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

   721  M81h = (1.- (bias/dimsizes(uu)))*score_max

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

   724  if (isvar("compare")) then

   725     system("sed -e '1,/M_npp_H81/s/M_npp_H81/"+M_npp_H81+"/' "+html_name2+" > "+html_new2+";"+ \

   726            "mv -f "+html_new2+" "+html_name2)

   727  end if

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

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

   732  delete (u)

   733  delete (v)

   734  delete (uu)

   735  delete (vv)

   736 ;----------------------------------------------------------------------

   737 ; histogram res

   738   resh                = True

   739   resh@gsnMaximize    = True

   740   resh@gsnDraw        = False

   741   resh@gsnFrame       = False

   742   resh@xyMarkLineMode = "Markers"

   743   resh@xyMarkerSizeF  = 0.014

   744   resh@xyMarker       = 16

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

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

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

   749   resh@tiYAxisString  = "NPP (gC m~S~-2~N~ y~S~-1~N~)"

   750   resh@tiXAxisString  = "Precipitation (m y~S~-1~N~)"

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

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

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

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

   757   lnres = True

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

   760 ;**************************************************************************

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

   762 ;**************************************************************************

   764   plot_name = "histogram_ob_81"

   765   title     = "Class A Observations (81 sites)"

   766   resh@tiMainString  = title

   768   wks   = gsn_open_wks (plot_type,plot_name)

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

   772   ; Added by Forrest Hoffman to print out values on Wed Feb  4 14:36:00 EST 2009

   773   print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))

   775 ;-------------------------------

   776 ;Attach the vertical bar and the horizontal cap line

   778   do nd=0,0

   779     lnres@gsLineColor = line_colors(nd)

   780     do i=0,nx-1

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

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

   784 ;

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

   786 ;

   787         x1 = xvalues(nd,i)

   788         y1 = yvalues(nd,i)

   789         y2 = mn_yvalues(nd,i)

   790         plx = (/x1,x1/)

   791         ply = (/y1,y2/)

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

   794         y2 = mx_yvalues(nd,i)

   795         plx = (/x1,x1/)

   796         ply = (/y1,y2/)

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

   798 ;

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

   800 ;

   801         x1 = xvalues(nd,i) - dx4

   802         x2 = xvalues(nd,i) + dx4

   803         y1 = mn_yvalues(nd,i)

   804         plx = (/x1,x2/)

   805         ply = (/y1,y1/)

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

   808         y1 = mx_yvalues(nd,i)

   809         plx = (/x1,x2/)

   810         ply = (/y1,y1/)

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

   812       end if

   813     end do

   814   end do

   816   draw(xy)

   817   delete (wks)

   819  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

   821  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   823 ;****************************************************************************

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

   825 ;****************************************************************************

   827   plot_name = "histogram_model_vs_ob_81"

   828   title     = model_name+ " (1975-2000) vs Class A Observations (81 sites)"

   829   resh@tiMainString  = title

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

   833 ;-----------------------------

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

   835 ; vertical lines going through the markers.

   837   resh@pmLegendDisplayMode    = "Always"

   838 ; resh@pmLegendWidthF         = 0.1

   839   resh@pmLegendWidthF         = 0.08

   840   resh@pmLegendHeightF        = 0.05

   841   resh@pmLegendOrthogonalPosF = -1.17

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

   843 ; resh@pmLegendParallelPosF   =  0.18

   844   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

   846 ; resh@lgPerimOn              = False

   847   resh@lgLabelFontHeightF     = 0.015

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

   849 ;-----------------------------

   850   tRes  = True

   851   tRes@txFontHeightF = 0.025

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

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

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

   859   ; Added by Forrest Hoffman to print out values on Wed Feb  4 14:36:00 EST 2009

   860   print(model_name + ": " + xvalues(1,:) + " " + yvalues(1,:))

   861   ;print("All: " + xvalues(:,:) + " " + yvalues(:,:))

   862 ;-------------------------------

   863 ;Attach the vertical bar and the horizontal cap line

   865   do nd=0,1

   866     lnres@gsLineColor = line_colors(nd)

   867     do i=0,nx-1

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

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

   871 ;

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

   873 ;

   874         x1 = xvalues(nd,i)

   875         y1 = yvalues(nd,i)

   876         y2 = mn_yvalues(nd,i)

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

   879         y2 = mx_yvalues(nd,i)

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

   881 ;

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

   883 ;

   884         x1 = xvalues(nd,i) - dx4

   885         x2 = xvalues(nd,i) + dx4

   886         y1 = mn_yvalues(nd,i)

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

   889         y1 = mx_yvalues(nd,i)

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

   891       end if

   892     end do

   893   end do

   895   draw(xy)

   896   delete(wks)

   898  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

   900  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   902  delete (RAIN1_1D)

   903  delete (RAIN2_1D)

   904  delete (NPP1_1D)

   905  delete (NPP2_1D)

   906 ;delete (range)

   907  delete (xvalues)

   908  delete (yvalues)

   909  delete (mn_yvalues)

   910  delete (mx_yvalues)

   911  delete (good)

   912  delete (max_bar)

   913  delete (min_bar)

   914  delete (max_cap)

   915  delete (min_cap)

   917 ;**************************************************************************

   918 ;(B) histogram 933

   919 ;**************************************************************************

   921 ; get data

   922   RAIN1_1D = ndtooned(rain933)

   923   RAIN2_1D = ndtooned(rainmod933)

   924   NPP1_1D  = ndtooned(npp933)

   925   NPP2_1D  = ndtooned(nppmod933)

   927 ; number of bin

   928   nx = 10

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

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

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

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

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

   936   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

   937          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

   939 ;----------------------------------------

   940 ;compute correlation coeff and M score

   942  u = yvalues(0,:)

   943  v = yvalues(1,:)

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

   946  uu = u(good)

   947  vv = v(good)

   949  ccr933h = esccr(uu,vv,0)

   951  score_max = 2.0

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

   954  M933h = (1.- (bias/dimsizes(uu)))*score_max

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

   957  if (isvar("compare")) then

   958     system("sed -e '1,/M_npp_H933/s/M_npp_H933/"+M_npp_H933+"/' "+html_name2+" > "+html_new2+";"+ \

   959            "mv -f "+html_new2+" "+html_name2)

   960  end if

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

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

   965  delete (u)

   966  delete (v)

   967  delete (uu)

   968  delete (vv)

   969 ;----------------------------------------------------------------------

   970 ; histogram res

   971   delete (resh)

   972   resh                = True

   973   resh@gsnMaximize    = True

   974   resh@gsnDraw        = False

   975   resh@gsnFrame       = False

   976   resh@xyMarkLineMode = "Markers"

   977   resh@xyMarkerSizeF  = 0.014

   978   resh@xyMarker       = 16

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

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

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

   983   resh@tiYAxisString  = "NPP (gC m~S~-2~N~ y~S~-1~N~)"

   984   resh@tiXAxisString  = "Precipitation (m y~S~-1~N~)"

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

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

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

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

   991   lnres = True

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

   994 ;**************************************************************************

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

   996 ;**************************************************************************

   998   plot_name = "histogram_ob_933"

   999   title     = "Class B Observations (933 sites)"

  1000   resh@tiMainString  = title

  1002   wks   = gsn_open_wks (plot_type,plot_name)

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

  1006   ; Added by Forrest Hoffman to print out values on Wed Feb  4 14:36:00 EST 2009

  1007   print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))

  1009 ;-------------------------------

  1010 ;Attach the vertical bar and the horizontal cap line

  1012   do nd=0,0

  1013     lnres@gsLineColor = line_colors(nd)

  1014     do i=0,nx-1

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

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

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

  1021         x1 = xvalues(nd,i)

  1022         y1 = yvalues(nd,i)

  1023         y2 = mn_yvalues(nd,i)

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

  1026         y2 = mx_yvalues(nd,i)

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

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

  1031         x1 = xvalues(nd,i) - dx4

  1032         x2 = xvalues(nd,i) + dx4

  1033         y1 = mn_yvalues(nd,i)

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

  1036         y1 = mx_yvalues(nd,i)

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

  1038       end if

  1039     end do

  1040   end do

  1042   draw(xy)

  1043   delete (xy)

  1044   delete (wks)

  1046  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1048  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1050 ;****************************************************************************

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

  1052 ;****************************************************************************

  1054   plot_name = "histogram_model_vs_ob_933"

  1055   title     = model_name+ " (1975-2000) vs Class B Observations (933 sites)"

  1056   resh@tiMainString  = title

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

  1060 ;-----------------------------

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

  1062 ; vertical lines going through the markers.

  1064   resh@pmLegendDisplayMode    = "Always"

  1065 ; resh@pmLegendWidthF         = 0.1

  1066   resh@pmLegendWidthF         = 0.08

  1067   resh@pmLegendHeightF        = 0.05

  1068   resh@pmLegendOrthogonalPosF = -1.17

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

  1070 ; resh@pmLegendParallelPosF   =  0.18

  1071   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1073 ; resh@lgPerimOn              = False

  1074   resh@lgLabelFontHeightF     = 0.015

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

  1076 ;-----------------------------

  1077   tRes  = True

  1078   tRes@txFontHeightF = 0.025

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

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

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

  1086   ; Added by Forrest Hoffman to print out values on Wed Feb  4 14:36:00 EST 2009

  1087   print("Observations: " + xvalues(1,:) + " " + yvalues(1,:))

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

  1089 ;Attach the vertical bar and the horizontal cap line

  1091   do nd=0,1

  1092     lnres@gsLineColor = line_colors(nd)

  1093     do i=0,nx-1

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

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

  1097 ;

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

  1099 ;

  1100         x1 = xvalues(nd,i)

  1101         y1 = yvalues(nd,i)

  1102         y2 = mn_yvalues(nd,i)

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

  1105         y2 = mx_yvalues(nd,i)

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

  1107 ;

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

  1109 ;

  1110         x1 = xvalues(nd,i) - dx4

  1111         x2 = xvalues(nd,i) + dx4

  1112         y1 = mn_yvalues(nd,i)

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

  1115         y1 = mx_yvalues(nd,i)

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

  1117       end if

  1118     end do

  1119   end do

  1121   draw(xy)

  1122   delete(xy)

  1123   delete(wks)

  1125  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1127  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1129 ;***************************************************************************

  1130 ; Read 2000-2004 dataset for MODIS comparison

  1131 ;***************************************************************************

  1133 ;------------------------------------------------------

  1134 ; read model data

  1136  fm   = addfile (dirm+film9,"r")

  1138  nppmod0  = fm->NPP

  1139  xm       = fm->lon

  1140  ym       = fm->lat

  1142  delete (fm)

  1144  nppmod   = nppmod0(0,:,:)

  1145  delete (nppmod0)

  1146  nppmod  = nppmod * nsec_per_year

  1147  nppmod@units     = "gC m~S~-2~N~ y~S~-1~N~"

  1148  nppmod@long_name     = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"

  1150 ;***************************************************************************

  1151 ;(C) global contour

  1152 ;***************************************************************************

  1154 ;res

  1155   resg                     = True             ; Use plot options

  1156   resg@cnFillOn            = True             ; Turn on color fill

  1157   resg@gsnSpreadColors     = True             ; use full colormap

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

  1159 ; resg@lbLabelAutoStride   = True

  1160   resg@cnLinesOn           = False            ; Turn off contourn lines

  1161   resg@mpFillOn            = False            ; Turn off map fill

  1163   resg@gsnSpreadColors      = True            ; use full colormap

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

  1165   resg@cnMinLevelValF       = 0.              ; Min level

  1166   resg@cnMaxLevelValF       = 2200.           ; Max level

  1167   resg@cnLevelSpacingF      = 200.            ; interval

  1169 ;****************************************************************************

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

  1171 ;****************************************************************************

  1173   delta = 0.00000000001

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

  1176   plot_name = "global_ob"

  1177   title     = "MODIS MOD17A3 (2000-2004)"

  1178   resg@tiMainString  = title

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

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

  1183   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1185   delete (wks)

  1187  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1189  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1190 ;****************************************************************************

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

  1192 ;****************************************************************************

  1194   plot_name = "global_model"

  1195   title     = "Model "+ model_name + " (2000-2004)"

  1196   resg@tiMainString  = title

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

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

  1201   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1203   delete (wks)

  1205  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1207  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1209 ;****************************************************************************

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

  1211 ;****************************************************************************

  1213   plot_name = "global_model_vs_ob"

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

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

  1218   delete (plot)

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

  1221   resg@gsnFrame             = False          ; Do not draw plot

  1222   resg@gsnDraw              = False          ; Do not advance frame

  1224 ; compute correlation coef and M score

  1226   uu1 = ndtooned(nppmod)

  1227   vv1 = ndtooned(nppglobe)

  1229   delete (good)

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

  1232   ug = uu1(good)

  1233   vg = vv1(good)

  1235   ccrG = esccr(ug,vg,0)

  1237   score_max = 2.0

  1239   MG   = (ccrG*ccrG)* score_max

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

  1242  if (isvar("compare")) then

  1243     system("sed -e '1,/M_npp_G/s/M_npp_G/"+M_npp_G+"/' "+html_name2+" > "+html_new2+";"+ \

  1244            "mv -f "+html_new2+" "+html_name2)

  1245  end if

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

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

  1250 ; plot correlation coef

  1252   gRes  = True

  1253   gRes@txFontHeightF = 0.02

  1254   gRes@txAngleF = 90

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

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

  1259 ;--------------------------------------------------------------------

  1260 ;(a) ob

  1262   title     = "MODIS MOD17A3 (2000-2004)"

  1263   resg@tiMainString  = title

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

  1267 ;(b) model

  1269   title     = "Model "+ model_name + " (2000-2004)"

  1270   resg@tiMainString  = title

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

  1274 ;(c) model-ob

  1276   zz = nppmod

  1277   zz = nppmod - nppglobe

  1278   title = "Model "+model_name+" - MODIS MOD17A3"

  1280   resg@cnMinLevelValF  = -500           ; Min level

  1281   resg@cnMaxLevelValF  =  500.          ; Max level

  1282   resg@cnLevelSpacingF =  50.           ; interval

  1283   resg@tiMainString    = title

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

  1287   pres                            = True        ; panel plot mods desired

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

  1289                                                 ; indiv. plots in panel

  1290   pres@gsnMaximize                = True        ; fill the page

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

  1294   delete (wks)

  1295   delete (plot)

  1297  system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1299  ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1301 ;***************************************************************************

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

  1303 ;***************************************************************************

  1305   vv     = zonalAve(nppglobe)

  1306   vv@long_name = nppglobe@long_name

  1308   plot_name = "zonal_ob"

  1309   title     = "MODIS MOD17A3 (2000-2004)"

  1311   resz = True

  1312   resz@tiMainString  = title

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

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

  1318   delete (wks)

  1320  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1322  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1324 ;****************************************************************************

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

  1326 ;****************************************************************************

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

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

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

  1333   s@long_name = nppglobe@long_name

  1334 ;-------------------------------------------

  1335 ; compute correlation coef and M score

  1337   score_max = 2.0

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

  1340   MZ   = (ccrZ*ccrZ)* score_max

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

  1343   if (isvar("compare")) then

  1344      system("sed -e '1,/M_npp_Z/s/M_npp_Z/"+M_npp_Z+"/' "+html_name2+" > "+html_new2+";"+ \

  1345             "mv -f "+html_new2+" "+html_name2)

  1346   end if

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

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

  1350 ;-------------------------------------------

  1351   plot_name = "zonal_model_vs_ob"

  1352   title     = "Zonal Average (2000-2004)"

  1353   resz@tiMainString  = title

  1355   wks = gsn_open_wks (plot_type,plot_name)

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

  1358 ; resz@vpWidthF           = 0.7

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

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

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

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

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

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

  1367   resz@txFontHeightF      = 0.0195            ; change title font heights

  1369 ; Legent

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

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

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

  1373   resz@pmLegendParallelPosF   = .75           ; move units right

  1374   resz@pmLegendOrthogonalPosF = -0.4          ; move units down

  1376   resz@pmLegendWidthF         = 0.10          ; Change width and

  1377   resz@pmLegendHeightF        = 0.10          ; height of legend.

  1378   resz@lgLabelFontHeightF     = .015          ; change font height

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

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

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

  1382   resz@xyExplicitLegendLabels = (/"MODIS MOD17A3",model_name/) ; explicit labels

  1383 ;--------------------------------------------------------------------

  1384   zRes  = True

  1385   zRes@txFontHeightF = 0.025

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

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

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

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

  1394   delete (wks)

  1396  system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

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

  1398  ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1400 ;***************************************************************************

  1401 ; add total score and write to file

  1402 ;***************************************************************************

  1403   M_total = M81s + M81h + M933s + M933h + MG + MZ

  1405   asciiwrite("M_save.npp", M_total)

  1407 ;***************************************************************************

  1408 ; output plots

  1409 ;***************************************************************************

  1410   output_dir = model_name+"/npp"

  1412   system("mv table_*.html " + output_dir +";"+ \

  1413          "mv *.png " + output_dir)

  1415 ;***************************************************************************

  1416 end

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