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

     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  fglobe   = addfile (dirg+filg,"r")

   199  nppglobe0 = fglobe->NPP

   200  nppglobe  = nppglobe0

   202  delete (fglobe)

   204 ;***********************************************************************

   205 ; interpolate model data into ob sites

   206 ;***********************************************************************

   208  nppmod   = nppmod0(0,:,:)

   209  rainmod  = rainmod0(0,:,:)

   210  delete (nppmod0)

   211  delete (rainmod0)

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

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

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

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

   221 ;**********************************************************

   222 ; unified units

   223 ;**********************************************************

   224 ; Units for these variables are:

   225 ;

   226 ; rain81  : mm/year

   227 ; rainmod : mm/s

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

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

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

   231 ;

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

   234   nsec_per_year = 60*60*24*365

   236   rain81    = rain81 / 1000.

   237   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   238   nppmod81  = nppmod81 * nsec_per_year

   240   rain933    = rain933 / 1000.

   241   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   242   nppmod933  = nppmod933 * nsec_per_year

   244   nppmod  = nppmod * nsec_per_year

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

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

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

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

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

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

   252   rain81@units     = "m/yr"

   253   rainmod81@units  = "m/yr"

   254   rain933@units    = "m/yr"

   255   rainmod933@units = "m/yr"

   257   npp81@long_name      = "Obs. NPP (gC/m2/year)"

   258   npp933@long_name     = "Obs. NPP (gC/m2/year)"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   276   output_html = "table_site81_ob.html"

   278   header = (/"<HTML>" \

   279             ,"<HEAD>" \

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

   281             ,"</HEAD>" \

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

   283             /)

   284   footer = "</HTML>"

   286   table_header = (/ \

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

   288        ,"<tr>" \

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

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

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

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

   293        ,"   <th bgcolor=DDDDDD >PPT(m/year)</th>" \

   294        ,"</tr>" \

   295        /)

   296   table_footer = "</table>"

   297   row_header = "<tr>"

   298   row_footer = "</tr>"

   300   lines = new(50000,string)

   301   nline = 0

   303   set_line(lines,nline,header)

   304   set_line(lines,nline,table_header)

   305 ;-----------------------------------------------

   306 ;row of table

   308   nrow = dimsizes(id81)

   309   do n = 0,nrow-1

   310      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

   324      set_line(lines,nline,row_footer)

   325   end do

   326 ;-----------------------------------------------

   327   set_line(lines,nline,table_footer)

   328   set_line(lines,nline,footer)

   330 ; Now write to an HTML file.

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

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

   333     asciiwrite(output_html,lines(idx))

   334   else

   335    print ("error?")

   336   end if

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

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

   340 ;*******************************************************************

   341   output_html = "table_site933_ob.html"

   343   header = (/"<HTML>" \

   344             ,"<HEAD>" \

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

   346             ,"</HEAD>" \

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

   348             /)

   349   footer = "</HTML>"

   351   table_header = (/ \

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

   353        ,"<tr>" \

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

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

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

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

   358        ,"   <th bgcolor=DDDDDD >PPT(m/year)</th>" \

   359        ,"</tr>" \

   360        /)

   361   table_footer = "</table>"

   362   row_header = "<tr>"

   363   row_footer = "</tr>"

   365   delete (lines)

   366   lines = new(50000,string)

   367   nline = 0

   369   set_line(lines,nline,header)

   370   set_line(lines,nline,table_header)

   371 ;-----------------------------------------------

   372 ;row of table

   374   nrow = dimsizes(id933)

   375   do n = 0,nrow-1

   376      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

   390      set_line(lines,nline,row_footer)

   391   end do

   392 ;-----------------------------------------------

   393   set_line(lines,nline,table_footer)

   394   set_line(lines,nline,footer)

   396 ; Now write to an HTML file.

   397   delete (idx)

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

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

   400     asciiwrite(output_html,lines(idx))

   401   else

   402    print ("error?")

   403   end if

   405 ;*******************************************************************

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

   407 ;*******************************************************************

   408   output_html = "table_site81_model_vs_ob.html"

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

   412   header = (/"<HTML>" \

   413             ,"<HEAD>" \

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

   415             ,"</HEAD>" \

   416             ,header_text \

   417             /)

   418   footer = "</HTML>"

   420   delete (table_header)

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

   422   table_header = (/ \

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

   424        ,"<tr>" \

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

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

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

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

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

   430        ,"</tr>" \

   431        ,"<tr>" \

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

   433        ,     table_header_text \

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

   435        ,     table_header_text \

   436        ,"</tr>" \

   437        /)

   438   table_footer = "</table>"

   439   row_header = "<tr>"

   440   row_footer = "</tr>"

   442   delete (lines)

   443   lines = new(50000,string)

   444   nline = 0

   446   set_line(lines,nline,header)

   447   set_line(lines,nline,table_header)

   448 ;-----------------------------------------------

   449 ;row of table

   451   nrow = dimsizes(id81)

   452   do n = 0,nrow-1

   453      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   471      set_line(lines,nline,row_footer)

   472   end do

   473 ;-----------------------------------------------

   474   set_line(lines,nline,table_footer)

   475   set_line(lines,nline,footer)

   477 ; Now write to an HTML file.

   478   delete (idx)

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

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

   481     asciiwrite(output_html,lines(idx))

   482   else

   483    print ("error?")

   484   end if

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

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

   488 ;*******************************************************************

   489   output_html = "table_site933_model_vs_ob.html"

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

   493   header = (/"<HTML>" \

   494             ,"<HEAD>" \

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

   496             ,"</HEAD>" \

   497             ,header_text \

   498             /)

   499   footer = "</HTML>"

   501   delete (table_header)

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

   503   table_header = (/ \

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

   505        ,"<tr>" \

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

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

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

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

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

   511        ,"</tr>" \

   512        ,"<tr>" \

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

   514        ,     table_header_text \

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

   516        ,     table_header_text \

   517        ,"</tr>" \

   518        /)

   519   table_footer = "</table>"

   520   row_header = "<tr>"

   521   row_footer = "</tr>"

   523   delete (lines)

   524   lines = new(50000,string)

   525   nline = 0

   527   set_line(lines,nline,header)

   528   set_line(lines,nline,table_header)

   529 ;-----------------------------------------------

   530 ;row of table

   532   nrow = dimsizes(id933)

   533   do n = 0,nrow-1

   534      set_line(lines,nline,row_header)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   552      set_line(lines,nline,row_footer)

   553   end do

   554 ;-----------------------------------------------

   555   set_line(lines,nline,table_footer)

   556   set_line(lines,nline,footer)

   558 ; Now write to an HTML file.

   559   delete (idx)

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

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

   562     asciiwrite(output_html,lines(idx))

   563   else

   564    print ("error?")

   565   end if

   567 ;***************************************************************************

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

   569 ;***************************************************************************

   571  plot_name = "scatter_model_vs_ob_81"

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

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

   575  res                   = True                  ; plot mods desired

   576  res@tiMainString      = title                 ; add title

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

   578  res@xyMarkers         =  16                   ; choose type of marker

   579  res@xyMarkerColor     = "red"                 ; Marker color

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

   581  res@tmLabelAutoStride = True                  ; nice tick mark labels

   583  res@gsnDraw           = False

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

   585 ;-------------------------------

   586 ;compute correlation coef. and M

   587  ccr81 = esccr(nppmod81,npp81,0)

   588 ;print (ccr81)

   590  score_max = 2.5

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

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

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

   596  if (isvar("compare")) then

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

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

   599  end if

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

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

   604  tRes  = True

   605  tRes@txFontHeightF = 0.025

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

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

   610 ;--------------------------------

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

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

   613 ; add polyline

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

   616 ;-------------------------------

   617  draw (plot)

   618  delete (wks)

   619  delete (dum)

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

   622         "rm "+plot_name+"."+plot_type)

   624 ;***************************************************************************

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

   626 ;***************************************************************************

   628  plot_name = "scatter_model_vs_ob_933"

   629  title     = model_name +" vs Class B Observations (933 sites)"

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

   632  res                   = True                  ; plot mods desired

   633  res@tiMainString      = title                 ; add title

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

   635  res@xyMarkers         =  16                   ; choose type of marker

   636  res@xyMarkerColor     = "red"                 ; Marker color

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

   638  res@tmLabelAutoStride = True                  ; nice tick mark labels

   640  res@gsnDraw           = False

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

   642 ;-------------------------------

   643 ;compute correlation coef. and M

   645  ccr933 = esccr(nppmod933,npp933,0)

   647  score_max = 2.5

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

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

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

   653  if (isvar("compare")) then

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

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

   656  end if

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

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

   661  tRes  = True

   662  tRes@txFontHeightF = 0.025

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

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

   667 ;--------------------------------

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

   669 ;-------------------------------

   670 ; add polyline

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

   673 ;-------------------------------

   674  draw (plot)

   675  delete (wks)

   676  delete (dum)

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

   679         "rm "+plot_name+"."+plot_type)

   681 ;**************************************************************************

   682 ;(B) histogram 81

   683 ;**************************************************************************

   684 ;get data

   685   RAIN1_1D = ndtooned(rain81)

   686   RAIN2_1D = ndtooned(rainmod81)

   687   NPP1_1D  = ndtooned(npp81)

   688   NPP2_1D  = ndtooned(nppmod81)

   690 ; number of bin

   691   nx = 8

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

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

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

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

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

   699   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

   700          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

   702 ;----------------------------------------

   703 ;compute correlation coeff and M score

   705  u = yvalues(0,:)

   706  v = yvalues(1,:)

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

   709  uu = u(good)

   710  vv = v(good)

   712  ccr81h = esccr(uu,vv,0)

   714  score_max = 2.5

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

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

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

   720  if (isvar("compare")) then

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

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

   723  end if

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

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

   728  delete (u)

   729  delete (v)

   730  delete (uu)

   731  delete (vv)

   732 ;----------------------------------------------------------------------

   733 ; histogram res

   734   resh                = True

   735   resh@gsnMaximize    = True

   736   resh@gsnDraw        = False

   737   resh@gsnFrame       = False

   738   resh@xyMarkLineMode = "Markers"

   739   resh@xyMarkerSizeF  = 0.014

   740   resh@xyMarker       = 16

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

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

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

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

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

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

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

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

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

   753   lnres = True

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

   756 ;**************************************************************************

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

   758 ;**************************************************************************

   760   plot_name = "histogram_ob_81"

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

   762   resh@tiMainString  = title

   764   wks   = gsn_open_wks (plot_type,plot_name)

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

   768 ;-------------------------------

   769 ;Attach the vertical bar and the horizontal cap line

   771   do nd=0,0

   772     lnres@gsLineColor = line_colors(nd)

   773     do i=0,nx-1

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

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

   777 ;

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

   779 ;

   780         x1 = xvalues(nd,i)

   781         y1 = yvalues(nd,i)

   782         y2 = mn_yvalues(nd,i)

   783         plx = (/x1,x1/)

   784         ply = (/y1,y2/)

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

   787         y2 = mx_yvalues(nd,i)

   788         plx = (/x1,x1/)

   789         ply = (/y1,y2/)

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

   791 ;

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

   793 ;

   794         x1 = xvalues(nd,i) - dx4

   795         x2 = xvalues(nd,i) + dx4

   796         y1 = mn_yvalues(nd,i)

   797         plx = (/x1,x2/)

   798         ply = (/y1,y1/)

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

   801         y1 = mx_yvalues(nd,i)

   802         plx = (/x1,x2/)

   803         ply = (/y1,y1/)

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

   805       end if

   806     end do

   807   end do

   809   draw(xy)

   810   delete (wks)

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

   813         "rm "+plot_name+"."+plot_type)

   815 ;****************************************************************************

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

   817 ;****************************************************************************

   819   plot_name = "histogram_model_vs_ob_81"

   820   title     = model_name+ " vs Class A Observations (81 sites)"

   821   resh@tiMainString  = title

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

   825 ;-----------------------------

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

   827 ; vertical lines going through the markers.

   829   resh@pmLegendDisplayMode    = "Always"

   830 ; resh@pmLegendWidthF         = 0.1

   831   resh@pmLegendWidthF         = 0.08

   832   resh@pmLegendHeightF        = 0.05

   833   resh@pmLegendOrthogonalPosF = -1.17

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

   835 ; resh@pmLegendParallelPosF   =  0.18

   836   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

   838 ; resh@lgPerimOn              = False

   839   resh@lgLabelFontHeightF     = 0.015

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

   841 ;-----------------------------

   842   tRes  = True

   843   tRes@txFontHeightF = 0.025

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

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

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

   850 ;-------------------------------

   851 ;Attach the vertical bar and the horizontal cap line

   853   do nd=0,1

   854     lnres@gsLineColor = line_colors(nd)

   855     do i=0,nx-1

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

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

   859 ;

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

   861 ;

   862         x1 = xvalues(nd,i)

   863         y1 = yvalues(nd,i)

   864         y2 = mn_yvalues(nd,i)

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

   867         y2 = mx_yvalues(nd,i)

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

   869 ;

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

   871 ;

   872         x1 = xvalues(nd,i) - dx4

   873         x2 = xvalues(nd,i) + dx4

   874         y1 = mn_yvalues(nd,i)

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

   877         y1 = mx_yvalues(nd,i)

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

   879       end if

   880     end do

   881   end do

   883   draw(xy)

   884   delete(wks)

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

   887         "rm "+plot_name+"."+plot_type)

   889  delete (RAIN1_1D)

   890  delete (RAIN2_1D)

   891  delete (NPP1_1D)

   892  delete (NPP2_1D)

   893 ;delete (range)

   894  delete (xvalues)

   895  delete (yvalues)

   896  delete (mn_yvalues)

   897  delete (mx_yvalues)

   898  delete (good)

   899  delete (max_bar)

   900  delete (min_bar)

   901  delete (max_cap)

   902  delete (min_cap)

   904 ;**************************************************************************

   905 ;(B) histogram 933

   906 ;**************************************************************************

   908 ; get data

   909   RAIN1_1D = ndtooned(rain933)

   910   RAIN2_1D = ndtooned(rainmod933)

   911   NPP1_1D  = ndtooned(npp933)

   912   NPP2_1D  = ndtooned(nppmod933)

   914 ; number of bin

   915   nx = 10

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

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

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

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

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

   923   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

   924          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

   926 ;----------------------------------------

   927 ;compute correlation coeff and M score

   929  u = yvalues(0,:)

   930  v = yvalues(1,:)

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

   933  uu = u(good)

   934  vv = v(good)

   936  ccr933h = esccr(uu,vv,0)

   938  score_max = 2.5

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

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

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

   944  if (isvar("compare")) then

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

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

   947  end if

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

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

   952  delete (u)

   953  delete (v)

   954  delete (uu)

   955  delete (vv)

   956 ;----------------------------------------------------------------------

   957 ; histogram res

   958   delete (resh)

   959   resh                = True

   960   resh@gsnMaximize    = True

   961   resh@gsnDraw        = False

   962   resh@gsnFrame       = False

   963   resh@xyMarkLineMode = "Markers"

   964   resh@xyMarkerSizeF  = 0.014

   965   resh@xyMarker       = 16

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

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

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

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

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

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

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

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

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

   978   lnres = True

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

   981 ;**************************************************************************

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

   983 ;**************************************************************************

   985   plot_name = "histogram_ob_933"

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

   987   resh@tiMainString  = title

   989   wks   = gsn_open_wks (plot_type,plot_name)

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

   993 ;-------------------------------

   994 ;Attach the vertical bar and the horizontal cap line

   996   do nd=0,0

   997     lnres@gsLineColor = line_colors(nd)

   998     do i=0,nx-1

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

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

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

  1005         x1 = xvalues(nd,i)

  1006         y1 = yvalues(nd,i)

  1007         y2 = mn_yvalues(nd,i)

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

  1010         y2 = mx_yvalues(nd,i)

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

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

  1015         x1 = xvalues(nd,i) - dx4

  1016         x2 = xvalues(nd,i) + dx4

  1017         y1 = mn_yvalues(nd,i)

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

  1020         y1 = mx_yvalues(nd,i)

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

  1022       end if

  1023     end do

  1024   end do

  1026   draw(xy)

  1027   delete (xy)

  1028   delete (wks)

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

  1031         "rm "+plot_name+"."+plot_type)

  1033 ;****************************************************************************

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

  1035 ;****************************************************************************

  1037   plot_name = "histogram_model_vs_ob_933"

  1038   title     = model_name+ " vs Class B Observations (933 sites)"

  1039   resh@tiMainString  = title

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

  1043 ;-----------------------------

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

  1045 ; vertical lines going through the markers.

  1047   resh@pmLegendDisplayMode    = "Always"

  1048 ; resh@pmLegendWidthF         = 0.1

  1049   resh@pmLegendWidthF         = 0.08

  1050   resh@pmLegendHeightF        = 0.05

  1051   resh@pmLegendOrthogonalPosF = -1.17

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

  1053 ; resh@pmLegendParallelPosF   =  0.18

  1054   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1056 ; resh@lgPerimOn              = False

  1057   resh@lgLabelFontHeightF     = 0.015

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

  1059 ;-----------------------------

  1060   tRes  = True

  1061   tRes@txFontHeightF = 0.025

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

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

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

  1068 ;-------------------------------

  1069 ;Attach the vertical bar and the horizontal cap line

  1071   do nd=0,1

  1072     lnres@gsLineColor = line_colors(nd)

  1073     do i=0,nx-1

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

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

  1077 ;

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

  1079 ;

  1080         x1 = xvalues(nd,i)

  1081         y1 = yvalues(nd,i)

  1082         y2 = mn_yvalues(nd,i)

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

  1085         y2 = mx_yvalues(nd,i)

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

  1087 ;

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

  1089 ;

  1090         x1 = xvalues(nd,i) - dx4

  1091         x2 = xvalues(nd,i) + dx4

  1092         y1 = mn_yvalues(nd,i)

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

  1095         y1 = mx_yvalues(nd,i)

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

  1097       end if

  1098     end do

  1099   end do

  1101   draw(xy)

  1102   delete(xy)

  1103   delete(wks)

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

  1106         "rm "+plot_name+"."+plot_type)

  1108 ;***************************************************************************

  1109 ;(C) global contour

  1110 ;***************************************************************************

  1112 ;res

  1113   resg                     = True             ; Use plot options

  1114   resg@cnFillOn            = True             ; Turn on color fill

  1115   resg@gsnSpreadColors     = True             ; use full colormap

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

  1117 ; resg@lbLabelAutoStride   = True

  1118   resg@cnLinesOn           = False            ; Turn off contourn lines

  1119   resg@mpFillOn            = False            ; Turn off map fill

  1121   resg@gsnSpreadColors      = True            ; use full colormap

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

  1123   resg@cnMinLevelValF       = 0.              ; Min level

  1124   resg@cnMaxLevelValF       = 2200.           ; Max level

  1125   resg@cnLevelSpacingF      = 200.            ; interval

  1127 ;****************************************************************************

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

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

  1131   delta = 0.00000000001

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

  1134   plot_name = "global_ob"

  1135   title     = "Observed MODIS MOD 17"

  1136   resg@tiMainString  = title

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

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

  1141   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1143   delete (wks)

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

  1146         "rm "+plot_name+"."+plot_type)

  1147 ;****************************************************************************

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

  1149 ;****************************************************************************

  1151   plot_name = "global_model"

  1152   title     = "Model "+ model_name

  1153   resg@tiMainString  = title

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

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

  1158   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1160   delete (wks)

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

  1163         "rm "+plot_name+"."+plot_type)

  1165 ;****************************************************************************

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

  1167 ;****************************************************************************

  1169   plot_name = "global_model_vs_ob"

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

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

  1174   delete (plot)

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

  1177   resg@gsnFrame             = False          ; Do not draw plot

  1178   resg@gsnDraw              = False          ; Do not advance frame

  1180 ; compute correlation coef and M score

  1182   uu1 = ndtooned(nppmod)

  1183   vv1 = ndtooned(nppglobe)

  1185   delete (good)

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

  1188   ug = uu1(good)

  1189   vg = vv1(good)

  1191   ccrG = esccr(ug,vg,0)

  1193   score_max = 5.0

  1195   MG   = (ccrG*ccrG)* score_max

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

  1198  if (isvar("compare")) then

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

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

  1201  end if

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

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

  1206 ; plot correlation coef

  1208   gRes  = True

  1209   gRes@txFontHeightF = 0.02

  1210   gRes@txAngleF = 90

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

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

  1215 ;--------------------------------------------------------------------

  1216 ;(a) ob

  1218   title     = "Observed MODIS MOD 17"

  1219   resg@tiMainString  = title

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

  1223 ;(b) model

  1225   title     = "Model "+ model_name

  1226   resg@tiMainString  = title

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

  1230 ;(c) model-ob

  1232   zz = nppmod

  1233   zz = nppmod - nppglobe

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

  1236   resg@cnMinLevelValF  = -500           ; Min level

  1237   resg@cnMaxLevelValF  =  500.          ; Max level

  1238   resg@cnLevelSpacingF =  50.           ; interval

  1239   resg@tiMainString    = title

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

  1243   pres                            = True        ; panel plot mods desired

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

  1245                                                 ; indiv. plots in panel

  1246   pres@gsnMaximize                = True        ; fill the page

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

  1250   delete (wks)

  1251   delete (plot)

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

  1254         "rm "+plot_name+"."+plot_type)

  1256 ;***************************************************************************

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

  1258 ;***************************************************************************

  1260   vv     = zonalAve(nppglobe)

  1261   vv@long_name = nppglobe@long_name

  1263   plot_name = "zonal_ob"

  1264   title     = "Observed MODIS MOD 17"

  1266   resz = True

  1267   resz@tiMainString  = title

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

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

  1273   delete (wks)

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

  1276         "rm "+plot_name+"."+plot_type)

  1278 ;****************************************************************************

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

  1280 ;****************************************************************************

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

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

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

  1287   s@long_name = nppglobe@long_name

  1288 ;-------------------------------------------

  1289 ; compute correlation coef and M score

  1291   score_max = 5.0

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

  1294   MZ   = (ccrZ*ccrZ)* score_max

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

  1297   if (isvar("compare")) then

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

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

  1300   end if

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

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

  1304 ;-------------------------------------------

  1305   plot_name = "zonal_model_vs_ob"

  1306   title     = "Zonal Average"

  1307   resz@tiMainString  = title

  1309   wks = gsn_open_wks (plot_type,plot_name)

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

  1312 ; resz@vpWidthF           = 0.7

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

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

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

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

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

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

  1321   resz@txFontHeightF      = 0.0195            ; change title font heights

  1323 ; Legent

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

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

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

  1327   resz@pmLegendParallelPosF   = .82           ; move units right

  1328   resz@pmLegendOrthogonalPosF = -0.4          ; move units down

  1330   resz@pmLegendWidthF         = 0.10          ; Change width and

  1331   resz@pmLegendHeightF        = 0.10          ; height of legend.

  1332   resz@lgLabelFontHeightF     = .02           ; change font height

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

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

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

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

  1337 ;--------------------------------------------------------------------

  1338   zRes  = True

  1339   zRes@txFontHeightF = 0.025

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

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

  1344 ;--------------------------------------------------------------------

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

  1348   delete (wks)

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

  1351         "rm "+plot_name+"."+plot_type)

  1353 ;***************************************************************************

  1354 ; add total score and write to file

  1355 ;***************************************************************************

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

  1358   asciiwrite("M_save.npp", M_total)

  1360 ;***************************************************************************

  1361 ; output plots

  1362 ;***************************************************************************

  1363   output_dir = model_name+"/npp"

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

  1366          "mv *.png " + output_dir)

  1368 ;***************************************************************************

  1369 end

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