npp/99.all.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
forrest@0
     1
;*****************************************************
forrest@0
     2
; combine scatter, histogram, global and zonal plots
forrest@0
     3
; compute all correlation coef and M score
forrest@0
     4
; add 1-to-1 line in scatter plots
forrest@0
     5
;*****************************************************
forrest@0
     6
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
forrest@0
     7
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
forrest@0
     8
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0
     9
;****************************************************************************
forrest@0
    10
procedure set_line(lines:string,nline:integer,newlines:string) 
forrest@0
    11
begin
forrest@0
    12
; add line to ascci/html file
forrest@0
    13
    
forrest@0
    14
  nnewlines = dimsizes(newlines)
forrest@0
    15
  if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0
    16
    print("set_line: bad index, not setting anything.") 
forrest@0
    17
    return
forrest@0
    18
  end if 
forrest@0
    19
  lines(nline:nline+nnewlines-1) = newlines
forrest@0
    20
;  print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0
    21
  nline = nline + nnewlines
forrest@0
    22
  return 
forrest@0
    23
end
forrest@0
    24
;****************************************************************************
forrest@0
    25
procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
forrest@0
    26
                 ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
forrest@0
    27
                 ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
forrest@0
    28
                 ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
forrest@0
    29
                 ,dx4[1]:numeric \
forrest@0
    30
                  )
forrest@0
    31
begin
forrest@0
    32
; Calculaee "nice" bins for binning the data in equally spaced ranges.
forrest@0
    33
; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
forrest@0
    34
; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
forrest@0
    35
forrest@0
    36
  nbins       = 15     ; Number of bins to use.
forrest@0
    37
forrest@0
    38
  nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
forrest@0
    39
  nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
forrest@0
    40
  range       = fspan(nicevals(0),nicevals(1),nvals)
forrest@0
    41
forrest@0
    42
; Use this range information to grab all the values in a
forrest@0
    43
; particular range, and then take an average.
forrest@0
    44
forrest@0
    45
  nr           = dimsizes(range)
forrest@0
    46
  nx           = nr-1
forrest@0
    47
; print (nx)
forrest@0
    48
forrest@0
    49
; xvalues      = new((/2,nx/),typeof(RAIN1_1D))
forrest@0
    50
  xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
forrest@0
    51
  dx           = xvalues(0,1) - xvalues(0,0)       ; range width
forrest@0
    52
  dx4          = dx/4                              ; 1/4 of the range
forrest@0
    53
  xvalues(1,:) = xvalues(0,:) - dx/5.
forrest@0
    54
; yvalues      = new((/2,nx/),typeof(RAIN1_1D))
forrest@0
    55
; mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))
forrest@0
    56
; mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))
forrest@0
    57
forrest@0
    58
  do nd=0,1
forrest@0
    59
forrest@0
    60
; See if we are doing model or observational data.
forrest@0
    61
forrest@0
    62
    if(nd.eq.0) then
forrest@0
    63
      data     = RAIN1_1D
forrest@0
    64
      npp_data = NPP1_1D
forrest@0
    65
    else
forrest@0
    66
      data     = RAIN2_1D
forrest@0
    67
      npp_data = NPP2_1D
forrest@0
    68
    end if
forrest@0
    69
forrest@0
    70
; Loop through each range and check for values.
forrest@0
    71
forrest@0
    72
    do i=0,nr-2
forrest@0
    73
      if (i.ne.(nr-2)) then
forrest@0
    74
;        print("")
forrest@0
    75
;        print("In range ["+range(i)+","+range(i+1)+")")
forrest@0
    76
        idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
forrest@0
    77
      else
forrest@0
    78
;        print("")
forrest@0
    79
;        print("In range ["+range(i)+",)")
forrest@0
    80
        idx = ind(range(i).le.data)
forrest@0
    81
      end if
forrest@0
    82
forrest@0
    83
; Calculate average, and get min and max.
forrest@0
    84
forrest@0
    85
      if(.not.any(ismissing(idx))) then
forrest@0
    86
        yvalues(nd,i)    = avg(npp_data(idx))
forrest@0
    87
        mn_yvalues(nd,i) = min(npp_data(idx))
forrest@0
    88
        mx_yvalues(nd,i) = max(npp_data(idx))
forrest@0
    89
        count = dimsizes(idx)
forrest@0
    90
      else
forrest@0
    91
        count            = 0
forrest@0
    92
        yvalues(nd,i)    = yvalues@_FillValue
forrest@0
    93
        mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0
    94
        mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0
    95
      end if
forrest@0
    96
forrest@0
    97
; Print out information.
forrest@0
    98
;      print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0
    99
;      print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0
   100
forrest@0
   101
; Clean up for next time in loop.
forrest@0
   102
forrest@0
   103
      delete(idx)
forrest@0
   104
    end do
forrest@0
   105
    delete(data)
forrest@0
   106
    delete(npp_data)
forrest@0
   107
  end do 
forrest@0
   108
end
forrest@0
   109
;****************************************************************************
forrest@0
   110
; Main code.
forrest@0
   111
;****************************************************************************
forrest@0
   112
begin
forrest@0
   113
forrest@0
   114
 plot_type     = "ps"
forrest@0
   115
 plot_type_new = "png"
forrest@0
   116
forrest@0
   117
;************************************************
forrest@0
   118
; read model data
forrest@0
   119
;************************************************
forrest@0
   120
forrest@0
   121
;film = "i01.06cn_1798-2004_ANN_climo.nc"
forrest@0
   122
;model_name = "06cn"
forrest@0
   123
;model_grid = "T42"
forrest@0
   124
forrest@0
   125
;film = "i01.06casa_1798-2004_ANN_climo.nc"
forrest@0
   126
;model_name = "06casa"
forrest@0
   127
;model_grid = "T42"
forrest@0
   128
forrest@0
   129
;film = "i01.10casa_1948-2004_ANN_climo.nc"
forrest@0
   130
;model_name = "10casa"
forrest@0
   131
;model_grid = "T42"
forrest@0
   132
forrest@0
   133
 film = "i01.10cn_1948-2004_ANN_climo.nc"
forrest@0
   134
 model_name = "10cn"
forrest@0
   135
 model_grid = "T42"
forrest@0
   136
forrest@0
   137
 html_name = "table.html." + model_name
forrest@0
   138
 html_new  = html_name +".new"
forrest@0
   139
 system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   140
         "mv -f "+html_new+" "+html_name)
forrest@0
   141
;--------------------------------------------------
forrest@0
   142
 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0
   143
 fm   = addfile (dirm+film,"r")
forrest@0
   144
  
forrest@0
   145
 nppmod0  = fm->NPP
forrest@0
   146
 rainmod0 = fm->RAIN
forrest@0
   147
 xm       = fm->lon  
forrest@0
   148
 ym       = fm->lat
forrest@0
   149
forrest@0
   150
;************************************************
forrest@0
   151
; read ob data
forrest@0
   152
;************************************************
forrest@0
   153
forrest@0
   154
;(1) data at 81 sites
forrest@0
   155
 dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
forrest@0
   156
 fil81 = "data.81.nc"
forrest@0
   157
 f81   = addfile (dir81+fil81,"r")
forrest@0
   158
forrest@0
   159
 id81   = f81->SITE_ID  
forrest@0
   160
 npp81  = f81->TNPP_C
forrest@0
   161
 rain81 = tofloat(f81->PREC_ANN)
forrest@0
   162
 x81    = f81->LONG_DD  
forrest@0
   163
 y81    = f81->LAT_DD
forrest@0
   164
forrest@0
   165
 id81@long_name  = "SITE_ID"
forrest@0
   166
forrest@0
   167
;change longitude from (-180,180) to (0,360)
forrest@0
   168
;for model data interpolation
forrest@0
   169
forrest@0
   170
 do i= 0,dimsizes(x81)-1
forrest@0
   171
    if (x81(i) .lt. 0.) then
forrest@0
   172
        x81(i) = x81(i)+ 360.
forrest@0
   173
    end if
forrest@0
   174
 end do
forrest@0
   175
;print (x81)
forrest@0
   176
;-------------------------------------
forrest@0
   177
;(2) data at 933 sites
forrest@0
   178
 dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
forrest@0
   179
 fil933 = "data.933.nc"
forrest@0
   180
 f933   = addfile (dir933+fil933,"r")
forrest@0
   181
forrest@0
   182
 id933   = f933->SITE_ID  
forrest@0
   183
 npp933  = f933->TNPP_C
forrest@0
   184
 rain933 = f933->PREC
forrest@0
   185
 x933    = f933->LONG_DD  
forrest@0
   186
 y933    = f933->LAT_DD 
forrest@0
   187
forrest@0
   188
 id933@long_name  = "SITE_ID"
forrest@0
   189
forrest@0
   190
;change longitude from (-180,180) to (0,360)
forrest@0
   191
;for model data interpolation
forrest@0
   192
forrest@0
   193
 do i= 0,dimsizes(x933)-1
forrest@0
   194
    if (x933(i) .lt. 0.) then
forrest@0
   195
        x933(i) = x933(i)+ 360.
forrest@0
   196
    end if
forrest@0
   197
 end do
forrest@0
   198
;print (x933)
forrest@0
   199
;----------------------------------------
forrest@0
   200
;(3) global data, interpolated into model grid
forrest@0
   201
 dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
forrest@0
   202
 filglobe = "Npp_"+model_grid+"_mean.nc"
forrest@0
   203
 fglobe   = addfile (dirglobe+filglobe,"r")
forrest@0
   204
 
forrest@0
   205
 nppglobe0 = fglobe->NPP
forrest@0
   206
 nppglobe  = nppglobe0
forrest@0
   207
forrest@0
   208
;***********************************************************************
forrest@0
   209
; interpolate model data into ob sites
forrest@0
   210
;***********************************************************************
forrest@0
   211
forrest@0
   212
 nppmod   = nppmod0(0,:,:)
forrest@0
   213
 rainmod  = rainmod0(0,:,:)
forrest@0
   214
 delete (nppmod0)
forrest@0
   215
 delete (rainmod0)
forrest@0
   216
forrest@0
   217
 nppmod81  =linint2_points(xm,ym,nppmod,True,x81,y81,0)
forrest@0
   218
forrest@0
   219
 nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
forrest@0
   220
forrest@0
   221
 rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
forrest@0
   222
forrest@0
   223
 rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
forrest@0
   224
forrest@0
   225
;**********************************************************
forrest@0
   226
; unified units
forrest@0
   227
;**********************************************************
forrest@0
   228
; Units for these variables are:
forrest@0
   229
;
forrest@0
   230
; rain81  : mm/year
forrest@0
   231
; rainmod : mm/s
forrest@0
   232
; npp81   : g C/m^2/year
forrest@0
   233
; nppmod81: g C/m^2/s
forrest@0
   234
; nppglobe: g C/m^2/year
forrest@0
   235
;
forrest@0
   236
; We want to convert these to "m/year" and "g C/m^2/year".
forrest@0
   237
forrest@0
   238
  nsec_per_year = 60*60*24*365                 
forrest@0
   239
forrest@0
   240
  rain81    = rain81 / 1000.
forrest@0
   241
  rainmod81 = (rainmod81/ 1000.) * nsec_per_year
forrest@0
   242
  nppmod81  = nppmod81 * nsec_per_year
forrest@0
   243
forrest@0
   244
  rain933    = rain933 / 1000.
forrest@0
   245
  rainmod933 = (rainmod933/ 1000.) * nsec_per_year
forrest@0
   246
  nppmod933  = nppmod933 * nsec_per_year
forrest@0
   247
forrest@0
   248
  nppmod  = nppmod * nsec_per_year
forrest@0
   249
forrest@0
   250
  npp81@units      = "gC/m^2/yr"
forrest@0
   251
  nppmod81@units   = "gC/m^2/yr"
forrest@0
   252
  npp933@units     = "gC/m^2/yr"
forrest@0
   253
  nppmod933@units  = "gC/m^2/yr"
forrest@0
   254
  nppmod@units     = "gC/m^2/yr"
forrest@0
   255
  nppglobe@units   = "gC/m^2/yr"
forrest@0
   256
  rain81@units     = "m/yr"
forrest@0
   257
  rainmod81@units  = "m/yr"
forrest@0
   258
  rain933@units    = "m/yr"
forrest@0
   259
  rainmod933@units = "m/yr"
forrest@0
   260
forrest@0
   261
  npp81@long_name      = "NPP (gC/m2/year)"
forrest@0
   262
  npp933@long_name     = "NPP (gC/m2/year)"
forrest@0
   263
  nppmod81@long_name   = "NPP (gC/m2/year)"
forrest@0
   264
  nppmod933@long_name  = "NPP (gC/m2/year)"
forrest@0
   265
  nppmod@long_name     = "NPP (gC/m2/year)"
forrest@0
   266
  nppglobe@long_name   = "NPP (gC/m2/year)"
forrest@0
   267
  rain81@long_name     = "PREC (m/year)"
forrest@0
   268
  rain933@long_name    = "PREC (m/year)"
forrest@0
   269
  rainmod81@long_name  = "PREC (m/year)"
forrest@0
   270
  rainmod933@long_name = "PREC (m/year)"
forrest@0
   271
forrest@0
   272
;*******************************************************************
forrest@0
   273
;(A)-1 html table of site81 -- observed
forrest@0
   274
;*******************************************************************
forrest@0
   275
  output_html = "table_site81_ob.html"
forrest@0
   276
forrest@0
   277
  header = (/"<HTML>" \
forrest@0
   278
            ,"<HEAD>" \
forrest@0
   279
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   280
            ,"</HEAD>" \
forrest@0
   281
            ,"<H1>Observation at 81 sites</H1>" \
forrest@0
   282
            /) 
forrest@0
   283
  footer = "</HTML>"
forrest@0
   284
forrest@0
   285
  table_header = (/ \
forrest@0
   286
        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
forrest@0
   287
       ,"<tr>" \
forrest@0
   288
       ,"   <th bgcolor=DDDDDD >Site ID</th>" \
forrest@0
   289
       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
forrest@0
   290
       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
forrest@0
   291
       ,"   <th bgcolor=DDDDDD >NPP(gC/m2.year)</th>" \
forrest@0
   292
       ,"   <th bgcolor=DDDDDD >RAIN(m/year)</th>" \
forrest@0
   293
       ,"</tr>" \
forrest@0
   294
       /)
forrest@0
   295
  table_footer = "</table>"
forrest@0
   296
  row_header = "<tr>"
forrest@0
   297
  row_footer = "</tr>"
forrest@0
   298
forrest@0
   299
  lines = new(50000,string)
forrest@0
   300
  nline = 0
forrest@0
   301
forrest@0
   302
  set_line(lines,nline,header)
forrest@0
   303
  set_line(lines,nline,table_header)
forrest@0
   304
;-----------------------------------------------
forrest@0
   305
;row of table
forrest@0
   306
  
forrest@0
   307
  nrow = dimsizes(id81)
forrest@0
   308
  do n = 0,nrow-1
forrest@0
   309
     set_line(lines,nline,row_header)
forrest@0
   310
forrest@0
   311
     txt1 = sprintf("%5.0f", id81(n))
forrest@0
   312
     txt2 = sprintf("%5.2f", y81(n))  
forrest@0
   313
     txt3 = sprintf("%5.2f", x81(n))    
forrest@0
   314
     txt4 = sprintf("%5.2f", npp81(n))     
forrest@0
   315
     txt5 = sprintf("%5.2f", rain81(n)) 
forrest@0
   316
forrest@0
   317
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   318
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   319
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   320
     set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0
   321
     set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0
   322
forrest@0
   323
     set_line(lines,nline,row_footer)
forrest@0
   324
  end do
forrest@0
   325
;-----------------------------------------------
forrest@0
   326
  set_line(lines,nline,table_footer)
forrest@0
   327
  set_line(lines,nline,footer) 
forrest@0
   328
forrest@0
   329
; Now write to an HTML file.
forrest@0
   330
  idx = ind(.not.ismissing(lines))
forrest@0
   331
  if(.not.any(ismissing(idx))) then
forrest@0
   332
    asciiwrite(output_html,lines(idx))
forrest@0
   333
  else
forrest@0
   334
   print ("error?")
forrest@0
   335
  end if
forrest@0
   336
;*******************************************************************
forrest@0
   337
;(A)-2 html table of site933 -- observed
forrest@0
   338
;*******************************************************************
forrest@0
   339
  output_html = "table_site933_ob.html"
forrest@0
   340
forrest@0
   341
  header = (/"<HTML>" \
forrest@0
   342
            ,"<HEAD>" \
forrest@0
   343
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   344
            ,"</HEAD>" \
forrest@0
   345
            ,"<H1>Observation at 933 sites</H1>" \
forrest@0
   346
            /) 
forrest@0
   347
  footer = "</HTML>"
forrest@0
   348
forrest@0
   349
  table_header = (/ \
forrest@0
   350
        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
forrest@0
   351
       ,"<tr>" \
forrest@0
   352
       ,"   <th bgcolor=DDDDDD >Site ID</th>" \
forrest@0
   353
       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
forrest@0
   354
       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
forrest@0
   355
       ,"   <th bgcolor=DDDDDD >NPP(gC/m2.year)</th>" \
forrest@0
   356
       ,"   <th bgcolor=DDDDDD >RAIN(m/year)</th>" \
forrest@0
   357
       ,"</tr>" \
forrest@0
   358
       /)
forrest@0
   359
  table_footer = "</table>"
forrest@0
   360
  row_header = "<tr>"
forrest@0
   361
  row_footer = "</tr>"
forrest@0
   362
forrest@0
   363
  delete (lines)  
forrest@0
   364
  lines = new(50000,string)
forrest@0
   365
  nline = 0
forrest@0
   366
forrest@0
   367
  set_line(lines,nline,header)
forrest@0
   368
  set_line(lines,nline,table_header)
forrest@0
   369
;-----------------------------------------------
forrest@0
   370
;row of table
forrest@0
   371
  
forrest@0
   372
  nrow = dimsizes(id933)
forrest@0
   373
  do n = 0,nrow-1
forrest@0
   374
     set_line(lines,nline,row_header)
forrest@0
   375
forrest@0
   376
     txt1 = sprintf("%5.0f", id933(n))
forrest@0
   377
     txt2 = sprintf("%5.2f", y933(n))  
forrest@0
   378
     txt3 = sprintf("%5.2f", x933(n))    
forrest@0
   379
     txt4 = sprintf("%5.2f", npp933(n))     
forrest@0
   380
     txt5 = sprintf("%5.2f", rain933(n))
forrest@0
   381
forrest@0
   382
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   383
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   384
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   385
     set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0
   386
     set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0
   387
forrest@0
   388
     set_line(lines,nline,row_footer)
forrest@0
   389
  end do
forrest@0
   390
;-----------------------------------------------
forrest@0
   391
  set_line(lines,nline,table_footer)
forrest@0
   392
  set_line(lines,nline,footer) 
forrest@0
   393
forrest@0
   394
; Now write to an HTML file.
forrest@0
   395
  delete (idx)
forrest@0
   396
  idx = ind(.not.ismissing(lines))
forrest@0
   397
  if(.not.any(ismissing(idx))) then
forrest@0
   398
    asciiwrite(output_html,lines(idx))
forrest@0
   399
  else
forrest@0
   400
   print ("error?")
forrest@0
   401
  end if
forrest@0
   402
;*******************************************************************
forrest@0
   403
;(A)-3 html table of site81 -- model vs ob
forrest@0
   404
;*******************************************************************
forrest@0
   405
  output_html = "table_site81_model_vs_ob.html"
forrest@0
   406
forrest@0
   407
  header_text = "<H1>Model "+model_name+" vs Observed at 81 sites</H1>" 
forrest@0
   408
forrest@0
   409
  header = (/"<HTML>" \
forrest@0
   410
            ,"<HEAD>" \
forrest@0
   411
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   412
            ,"</HEAD>" \
forrest@0
   413
            ,header_text \
forrest@0
   414
            /) 
forrest@0
   415
  footer = "</HTML>"
forrest@0
   416
forrest@0
   417
  delete (table_header)
forrest@0
   418
  table_header_text = "   <th bgcolor=DDDDDD >"+model_name+"</th>"
forrest@0
   419
  table_header = (/ \
forrest@0
   420
        "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
forrest@0
   421
       ,"<tr>" \
forrest@0
   422
       ,"   <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
forrest@0
   423
       ,"   <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
forrest@0
   424
       ,"   <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
forrest@0
   425
       ,"   <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \
forrest@0
   426
       ,"   <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \
forrest@0
   427
       ,"</tr>" \
forrest@0
   428
       ,"<tr>" \
forrest@0
   429
       ,"   <th bgcolor=DDDDDD >observed</th>" \
forrest@0
   430
       ,     table_header_text \
forrest@0
   431
       ,"   <th bgcolor=DDDDDD >observed</th>" \
forrest@0
   432
       ,     table_header_text \
forrest@0
   433
       ,"</tr>" \
forrest@0
   434
       /)
forrest@0
   435
  table_footer = "</table>"
forrest@0
   436
  row_header = "<tr>"
forrest@0
   437
  row_footer = "</tr>"
forrest@0
   438
forrest@0
   439
  delete (lines)
forrest@0
   440
  lines = new(50000,string)
forrest@0
   441
  nline = 0
forrest@0
   442
forrest@0
   443
  set_line(lines,nline,header)
forrest@0
   444
  set_line(lines,nline,table_header)
forrest@0
   445
;-----------------------------------------------
forrest@0
   446
;row of table
forrest@0
   447
  
forrest@0
   448
  nrow = dimsizes(id81)
forrest@0
   449
  do n = 0,nrow-1
forrest@0
   450
     set_line(lines,nline,row_header)
forrest@0
   451
forrest@0
   452
     txt1 = sprintf("%5.0f", id81(n))
forrest@0
   453
     txt2 = sprintf("%5.2f", y81(n))  
forrest@0
   454
     txt3 = sprintf("%5.2f", x81(n))    
forrest@0
   455
     txt4 = sprintf("%5.2f", npp81(n))
forrest@0
   456
     txt5 = sprintf("%5.2f", nppmod81(n))     
forrest@0
   457
     txt6 = sprintf("%5.2f", rain81(n))
forrest@0
   458
     txt7 = sprintf("%5.2f", rainmod81(n)) 
forrest@0
   459
forrest@0
   460
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   461
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   462
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   463
     set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0
   464
     set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0
   465
     set_line(lines,nline,"<th>"+txt6+"</th>")
forrest@0
   466
     set_line(lines,nline,"<th>"+txt7+"</th>")
forrest@0
   467
forrest@0
   468
     set_line(lines,nline,row_footer)
forrest@0
   469
  end do
forrest@0
   470
;-----------------------------------------------
forrest@0
   471
  set_line(lines,nline,table_footer)
forrest@0
   472
  set_line(lines,nline,footer) 
forrest@0
   473
forrest@0
   474
; Now write to an HTML file.
forrest@0
   475
  delete (idx)
forrest@0
   476
  idx = ind(.not.ismissing(lines))
forrest@0
   477
  if(.not.any(ismissing(idx))) then
forrest@0
   478
    asciiwrite(output_html,lines(idx))
forrest@0
   479
  else
forrest@0
   480
   print ("error?")
forrest@0
   481
  end if
forrest@0
   482
forrest@0
   483
;*******************************************************************
forrest@0
   484
;(A)-4 html table of site933 -- model vs ob
forrest@0
   485
;*******************************************************************
forrest@0
   486
  output_html = "table_site933_model_vs_ob.html"
forrest@0
   487
forrest@0
   488
  header_text = "<H1>Model "+model_name+" vs Observed at 933 sites</H1>" 
forrest@0
   489
forrest@0
   490
  header = (/"<HTML>" \
forrest@0
   491
            ,"<HEAD>" \
forrest@0
   492
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   493
            ,"</HEAD>" \
forrest@0
   494
            ,header_text \
forrest@0
   495
            /) 
forrest@0
   496
  footer = "</HTML>"
forrest@0
   497
forrest@0
   498
  delete (table_header)
forrest@0
   499
  table_header_text = "   <th bgcolor=DDDDDD >"+model_name+"</th>"
forrest@0
   500
  table_header = (/ \
forrest@0
   501
        "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
forrest@0
   502
       ,"<tr>" \
forrest@0
   503
       ,"   <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
forrest@0
   504
       ,"   <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
forrest@0
   505
       ,"   <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
forrest@0
   506
       ,"   <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \
forrest@0
   507
       ,"   <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \
forrest@0
   508
       ,"</tr>" \
forrest@0
   509
       ,"<tr>" \
forrest@0
   510
       ,"   <th bgcolor=DDDDDD >observed</th>" \
forrest@0
   511
       ,     table_header_text \
forrest@0
   512
       ,"   <th bgcolor=DDDDDD >observed</th>" \
forrest@0
   513
       ,     table_header_text \
forrest@0
   514
       ,"</tr>" \
forrest@0
   515
       /)
forrest@0
   516
  table_footer = "</table>"
forrest@0
   517
  row_header = "<tr>"
forrest@0
   518
  row_footer = "</tr>"
forrest@0
   519
forrest@0
   520
  delete (lines)
forrest@0
   521
  lines = new(50000,string)
forrest@0
   522
  nline = 0
forrest@0
   523
forrest@0
   524
  set_line(lines,nline,header)
forrest@0
   525
  set_line(lines,nline,table_header)
forrest@0
   526
;-----------------------------------------------
forrest@0
   527
;row of table
forrest@0
   528
  
forrest@0
   529
  nrow = dimsizes(id933)
forrest@0
   530
  do n = 0,nrow-1
forrest@0
   531
     set_line(lines,nline,row_header)
forrest@0
   532
forrest@0
   533
     txt1 = sprintf("%5.0f", id933(n))
forrest@0
   534
     txt2 = sprintf("%5.2f", y933(n))  
forrest@0
   535
     txt3 = sprintf("%5.2f", x933(n))    
forrest@0
   536
     txt4 = sprintf("%5.2f", npp933(n))
forrest@0
   537
     txt5 = sprintf("%5.2f", nppmod933(n))     
forrest@0
   538
     txt6 = sprintf("%5.2f", rain933(n))
forrest@0
   539
     txt7 = sprintf("%5.2f", rainmod933(n)) 
forrest@0
   540
forrest@0
   541
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   542
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   543
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   544
     set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0
   545
     set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0
   546
     set_line(lines,nline,"<th>"+txt6+"</th>")
forrest@0
   547
     set_line(lines,nline,"<th>"+txt7+"</th>")
forrest@0
   548
forrest@0
   549
     set_line(lines,nline,row_footer)
forrest@0
   550
  end do
forrest@0
   551
;-----------------------------------------------
forrest@0
   552
  set_line(lines,nline,table_footer)
forrest@0
   553
  set_line(lines,nline,footer) 
forrest@0
   554
forrest@0
   555
; Now write to an HTML file.
forrest@0
   556
  delete (idx)
forrest@0
   557
  idx = ind(.not.ismissing(lines))
forrest@0
   558
  if(.not.any(ismissing(idx))) then
forrest@0
   559
    asciiwrite(output_html,lines(idx))
forrest@0
   560
  else
forrest@0
   561
   print ("error?")
forrest@0
   562
  end if
forrest@0
   563
;***************************************************************************
forrest@0
   564
;(A)-5 scatter plot, model vs ob 81
forrest@0
   565
;***************************************************************************
forrest@0
   566
  
forrest@0
   567
 plot_name = "scatter_model_vs_ob_81"
forrest@0
   568
 title     = model_name +" vs ob 81 sites"
forrest@0
   569
forrest@0
   570
 wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
forrest@0
   571
 res                   = True                  ; plot mods desired
forrest@0
   572
 res@tiMainString      = title                 ; add title
forrest@0
   573
 res@xyMarkLineModes   = "Markers"             ; choose which have markers
forrest@0
   574
 res@xyMarkers         =  16                   ; choose type of marker  
forrest@0
   575
 res@xyMarkerColor     = "red"                 ; Marker color
forrest@0
   576
 res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
forrest@0
   577
 res@tmLabelAutoStride = True                  ; nice tick mark labels
forrest@0
   578
forrest@0
   579
 res@gsnDraw           = False
forrest@0
   580
 res@gsnFrame          = False                 ; don't advance frame yet
forrest@0
   581
;-------------------------------
forrest@0
   582
;compute correlation coef. and M
forrest@0
   583
 ccr81 = esccr(nppmod81,npp81,0)
forrest@0
   584
;print (ccr81)
forrest@0
   585
forrest@0
   586
;new eq
forrest@0
   587
 bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81))) 
forrest@0
   588
 M81  = (1. - (bias/dimsizes(y81)))*5.
forrest@0
   589
forrest@0
   590
 M_npp_S81  = sprintf("%.2f", M81)
forrest@0
   591
 system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   592
        "mv -f "+html_new+" "+html_name)
forrest@0
   593
 print (M_npp_S81)
forrest@0
   594
 
forrest@0
   595
 tRes  = True
forrest@0
   596
 tRes@txFontHeightF = 0.025
forrest@0
   597
forrest@0
   598
 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
forrest@0
   599
forrest@0
   600
 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
forrest@0
   601
;--------------------------------
forrest@0
   602
 plot  = gsn_csm_xy (wks,npp81,nppmod81,res)       ; create plot
forrest@0
   603
;-------------------------------
forrest@0
   604
; add polyline
forrest@0
   605
forrest@0
   606
 dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
forrest@0
   607
;-------------------------------
forrest@0
   608
 draw (plot)
forrest@0
   609
 frame(wks)
forrest@0
   610
 clear (wks)
forrest@0
   611
 delete (dum)
forrest@0
   612
forrest@0
   613
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   614
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   615
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
   616
        "rm "+plot_name+"."+plot_type)
forrest@0
   617
;***************************************************************************
forrest@0
   618
;(A)-6 scatter plot, model vs ob 933
forrest@0
   619
;***************************************************************************
forrest@0
   620
  
forrest@0
   621
 plot_name = "scatter_model_vs_ob_933"
forrest@0
   622
 title     = model_name +" vs ob 933 sites"
forrest@0
   623
forrest@0
   624
 wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
forrest@0
   625
 res                   = True                  ; plot mods desired
forrest@0
   626
 res@tiMainString      = title                 ; add title
forrest@0
   627
 res@xyMarkLineModes   = "Markers"             ; choose which have markers
forrest@0
   628
 res@xyMarkers         =  16                   ; choose type of marker  
forrest@0
   629
 res@xyMarkerColor     = "red"                 ; Marker color
forrest@0
   630
 res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
forrest@0
   631
 res@tmLabelAutoStride = True                  ; nice tick mark labels
forrest@0
   632
forrest@0
   633
 res@gsnDraw           = False
forrest@0
   634
 res@gsnFrame          = False                 ; don't advance frame yet
forrest@0
   635
;-------------------------------
forrest@0
   636
;compute correlation coef. and M
forrest@0
   637
 ccr933 = esccr(nppmod933,npp933,0)
forrest@0
   638
;print (ccr933)
forrest@0
   639
forrest@0
   640
;new eq
forrest@0
   641
 bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
forrest@0
   642
 M933 = (1. - (bias/dimsizes(y933)))*5.
forrest@0
   643
forrest@0
   644
 M_npp_S933  = sprintf("%.2f", M933)
forrest@0
   645
 system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   646
        "mv -f "+html_new+" "+html_name)
forrest@0
   647
 print (M_npp_S933)
forrest@0
   648
forrest@0
   649
 tRes  = True
forrest@0
   650
 tRes@txFontHeightF = 0.025
forrest@0
   651
forrest@0
   652
 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
forrest@0
   653
forrest@0
   654
 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
forrest@0
   655
;--------------------------------
forrest@0
   656
 plot  = gsn_csm_xy (wks,npp933,nppmod933,res)    ; create plot
forrest@0
   657
;-------------------------------
forrest@0
   658
; add polyline
forrest@0
   659
forrest@0
   660
 dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
forrest@0
   661
;-------------------------------
forrest@0
   662
 draw (plot)
forrest@0
   663
 frame(wks)
forrest@0
   664
 clear (wks)
forrest@0
   665
 delete (dum)
forrest@0
   666
forrest@0
   667
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   668
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   669
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
   670
        "rm "+plot_name+"."+plot_type)
forrest@0
   671
;**************************************************************************
forrest@0
   672
;(B) histogram 81
forrest@0
   673
;**************************************************************************
forrest@0
   674
;get data
forrest@0
   675
  RAIN1_1D = ndtooned(rain81)
forrest@0
   676
  RAIN2_1D = ndtooned(rainmod81)
forrest@0
   677
  NPP1_1D  = ndtooned(npp81)
forrest@0
   678
  NPP2_1D  = ndtooned(nppmod81)
forrest@0
   679
forrest@0
   680
; number of bin
forrest@0
   681
  nx = 8
forrest@0
   682
forrest@0
   683
  xvalues      = new((/2,nx/),float)
forrest@0
   684
  yvalues      = new((/2,nx/),float)
forrest@0
   685
  mn_yvalues   = new((/2,nx/),float)
forrest@0
   686
  mx_yvalues   = new((/2,nx/),float)
forrest@0
   687
  dx4          = new((/1/),float)
forrest@0
   688
forrest@0
   689
  get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
forrest@0
   690
         ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
forrest@0
   691
forrest@0
   692
;----------------------------------------
forrest@0
   693
;compute correlation coeff and M score
forrest@0
   694
 
forrest@0
   695
 u = yvalues(0,:)
forrest@0
   696
 v = yvalues(1,:)
forrest@0
   697
forrest@0
   698
 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0
   699
 uu = u(good)
forrest@0
   700
 vv = v(good)
forrest@0
   701
forrest@0
   702
 ccr81h = esccr(uu,vv,0)
forrest@0
   703
;print (ccr81h)
forrest@0
   704
forrest@0
   705
;new eq
forrest@0
   706
 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
forrest@0
   707
 M81h = (1.- (bias/dimsizes(uu)))*5.
forrest@0
   708
forrest@0
   709
 M_npp_H81  = sprintf("%.2f", M81h)
forrest@0
   710
 system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   711
        "mv -f "+html_new+" "+html_name)
forrest@0
   712
 print (M_npp_H81)
forrest@0
   713
 
forrest@0
   714
 delete (u)
forrest@0
   715
 delete (v)
forrest@0
   716
 delete (uu)
forrest@0
   717
 delete (vv)
forrest@0
   718
;----------------------------------------------------------------------
forrest@0
   719
; histogram res
forrest@0
   720
  resh                = True
forrest@0
   721
  resh@gsnMaximize    = True
forrest@0
   722
  resh@gsnDraw        = False
forrest@0
   723
  resh@gsnFrame       = False
forrest@0
   724
  resh@xyMarkLineMode = "Markers"
forrest@0
   725
  resh@xyMarkerSizeF  = 0.014
forrest@0
   726
  resh@xyMarker       = 16
forrest@0
   727
  resh@xyMarkerColors = (/"Brown","Blue"/)
forrest@0
   728
  resh@trYMinF        = min(mn_yvalues) - 10.
forrest@0
   729
  resh@trYMaxF        = max(mx_yvalues) + 10.
forrest@0
   730
forrest@0
   731
  resh@tiYAxisString  = "NPP (g C/m2/year)"
forrest@0
   732
  resh@tiXAxisString  = "Precipitation (m/year)"
forrest@0
   733
forrest@0
   734
  max_bar = new((/2,nx/),graphic)
forrest@0
   735
  min_bar = new((/2,nx/),graphic)
forrest@0
   736
  max_cap = new((/2,nx/),graphic)
forrest@0
   737
  min_cap = new((/2,nx/),graphic)
forrest@0
   738
forrest@0
   739
  lnres = True
forrest@0
   740
  line_colors = (/"brown","blue"/)
forrest@0
   741
forrest@0
   742
;**************************************************************************
forrest@0
   743
;(B)-1 histogram plot, ob 81 site 
forrest@0
   744
;**************************************************************************
forrest@0
   745
forrest@0
   746
  plot_name = "histogram_ob_81"
forrest@0
   747
  title     = "Observed 81 site"
forrest@0
   748
  resh@tiMainString  = title
forrest@0
   749
forrest@0
   750
  wks   = gsn_open_wks (plot_type,plot_name)    
forrest@0
   751
forrest@0
   752
  xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
forrest@0
   753
forrest@0
   754
;-------------------------------
forrest@0
   755
;Attach the vertical bar and the horizontal cap line 
forrest@0
   756
forrest@0
   757
  do nd=0,0
forrest@0
   758
    lnres@gsLineColor = line_colors(nd)
forrest@0
   759
    do i=0,nx-1
forrest@0
   760
     
forrest@0
   761
      if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0
   762
         .not.ismissing(mx_yvalues(nd,i))) then
forrest@0
   763
;
forrest@0
   764
; Attach the vertical bar, both above and below the marker.
forrest@0
   765
;
forrest@0
   766
        x1 = xvalues(nd,i)
forrest@0
   767
        y1 = yvalues(nd,i)
forrest@0
   768
        y2 = mn_yvalues(nd,i)
forrest@0
   769
        plx = (/x1,x1/)
forrest@0
   770
        ply = (/y1,y2/)
forrest@0
   771
        min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0
   772
forrest@0
   773
        y2 = mx_yvalues(nd,i)
forrest@0
   774
        plx = (/x1,x1/)
forrest@0
   775
        ply = (/y1,y2/)
forrest@0
   776
        max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0
   777
;
forrest@0
   778
; Attach the horizontal cap line, both above and below the marker.
forrest@0
   779
;
forrest@0
   780
        x1 = xvalues(nd,i) - dx4
forrest@0
   781
        x2 = xvalues(nd,i) + dx4
forrest@0
   782
        y1 = mn_yvalues(nd,i)
forrest@0
   783
        plx = (/x1,x2/)
forrest@0
   784
        ply = (/y1,y1/)
forrest@0
   785
        min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0
   786
forrest@0
   787
        y1 = mx_yvalues(nd,i)
forrest@0
   788
        plx = (/x1,x2/)
forrest@0
   789
        ply = (/y1,y1/)
forrest@0
   790
        max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0
   791
      end if
forrest@0
   792
    end do
forrest@0
   793
  end do
forrest@0
   794
forrest@0
   795
  draw(xy)
forrest@0
   796
  frame(wks)
forrest@0
   797
  clear (wks)
forrest@0
   798
forrest@0
   799
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   800
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   801
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
   802
        "rm "+plot_name+"."+plot_type)
forrest@0
   803
;****************************************************************************
forrest@0
   804
;(B)-2 histogram plot, model vs ob 81 site 
forrest@0
   805
;****************************************************************************
forrest@0
   806
forrest@0
   807
  plot_name = "histogram_model_vs_ob_81"
forrest@0
   808
  title     = model_name+ " vs Observed 81 site"
forrest@0
   809
  resh@tiMainString  = title
forrest@0
   810
forrest@0
   811
  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
forrest@0
   812
forrest@0
   813
;-----------------------------
forrest@0
   814
; Add a boxed legend using the more simple method, which won't have
forrest@0
   815
; vertical lines going through the markers.
forrest@0
   816
forrest@0
   817
  resh@pmLegendDisplayMode    = "Always"
forrest@0
   818
; resh@pmLegendWidthF         = 0.1
forrest@0
   819
  resh@pmLegendWidthF         = 0.08
forrest@0
   820
  resh@pmLegendHeightF        = 0.05
forrest@0
   821
  resh@pmLegendOrthogonalPosF = -1.17
forrest@0
   822
; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0
   823
; resh@pmLegendParallelPosF   =  0.18
forrest@0
   824
  resh@pmLegendParallelPosF   =  0.88  ;(rightward)
forrest@0
   825
forrest@0
   826
; resh@lgPerimOn              = False
forrest@0
   827
  resh@lgLabelFontHeightF     = 0.015
forrest@0
   828
  resh@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0
   829
;-----------------------------
forrest@0
   830
  tRes  = True
forrest@0
   831
  tRes@txFontHeightF = 0.025
forrest@0
   832
forrest@0
   833
  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
forrest@0
   834
forrest@0
   835
  gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0
   836
forrest@0
   837
  xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
forrest@0
   838
;-------------------------------
forrest@0
   839
;Attach the vertical bar and the horizontal cap line 
forrest@0
   840
forrest@0
   841
  do nd=0,1
forrest@0
   842
    lnres@gsLineColor = line_colors(nd)
forrest@0
   843
    do i=0,nx-1
forrest@0
   844
     
forrest@0
   845
      if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0
   846
         .not.ismissing(mx_yvalues(nd,i))) then
forrest@0
   847
;
forrest@0
   848
; Attach the vertical bar, both above and below the marker.
forrest@0
   849
;
forrest@0
   850
        x1 = xvalues(nd,i)
forrest@0
   851
        y1 = yvalues(nd,i)
forrest@0
   852
        y2 = mn_yvalues(nd,i)
forrest@0
   853
        min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0
   854
forrest@0
   855
        y2 = mx_yvalues(nd,i)
forrest@0
   856
        max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0
   857
;
forrest@0
   858
; Attach the horizontal cap line, both above and below the marker.
forrest@0
   859
;
forrest@0
   860
        x1 = xvalues(nd,i) - dx4
forrest@0
   861
        x2 = xvalues(nd,i) + dx4
forrest@0
   862
        y1 = mn_yvalues(nd,i)
forrest@0
   863
        min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0
   864
forrest@0
   865
        y1 = mx_yvalues(nd,i)
forrest@0
   866
        max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0
   867
      end if
forrest@0
   868
    end do
forrest@0
   869
  end do
forrest@0
   870
forrest@0
   871
  draw(xy)
forrest@0
   872
  frame(wks)
forrest@0
   873
  clear(wks)
forrest@0
   874
forrest@0
   875
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   876
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   877
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
   878
        "rm "+plot_name+"."+plot_type)
forrest@0
   879
forrest@0
   880
 delete (RAIN1_1D)
forrest@0
   881
 delete (RAIN2_1D)
forrest@0
   882
 delete (NPP1_1D)
forrest@0
   883
 delete (NPP2_1D)
forrest@0
   884
;delete (range)
forrest@0
   885
 delete (xvalues) 
forrest@0
   886
 delete (yvalues)
forrest@0
   887
 delete (mn_yvalues)
forrest@0
   888
 delete (mx_yvalues)
forrest@0
   889
 delete (good)
forrest@0
   890
 delete (max_bar)
forrest@0
   891
 delete (min_bar)
forrest@0
   892
 delete (max_cap)
forrest@0
   893
 delete (min_cap)
forrest@0
   894
   
forrest@0
   895
;**************************************************************************
forrest@0
   896
;(B) histogram 933
forrest@0
   897
;**************************************************************************
forrest@0
   898
forrest@0
   899
; get data
forrest@0
   900
  RAIN1_1D = ndtooned(rain933)
forrest@0
   901
  RAIN2_1D = ndtooned(rainmod933)
forrest@0
   902
  NPP1_1D  = ndtooned(npp933)
forrest@0
   903
  NPP2_1D  = ndtooned(nppmod933)
forrest@0
   904
forrest@0
   905
; number of bin
forrest@0
   906
  nx = 10
forrest@0
   907
forrest@0
   908
  xvalues      = new((/2,nx/),float)
forrest@0
   909
  yvalues      = new((/2,nx/),float)
forrest@0
   910
  mn_yvalues   = new((/2,nx/),float)
forrest@0
   911
  mx_yvalues   = new((/2,nx/),float)
forrest@0
   912
  dx4          = new((/1/),float)
forrest@0
   913
forrest@0
   914
  get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
forrest@0
   915
         ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
forrest@0
   916
 
forrest@0
   917
;----------------------------------------
forrest@0
   918
;compute correlation coeff and M score
forrest@0
   919
 
forrest@0
   920
 u = yvalues(0,:)
forrest@0
   921
 v = yvalues(1,:)
forrest@0
   922
forrest@0
   923
 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0
   924
 uu = u(good)
forrest@0
   925
 vv = v(good)
forrest@0
   926
forrest@0
   927
 ccr933h = esccr(uu,vv,0)
forrest@0
   928
;print (ccr933h)
forrest@0
   929
forrest@0
   930
;new eq
forrest@0
   931
 bias  = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
forrest@0
   932
 M933h = (1.- (bias/dimsizes(uu)))*5.
forrest@0
   933
forrest@0
   934
 M_npp_H933 = sprintf("%.2f", M933h)
forrest@0
   935
 system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   936
        "mv -f "+html_new+" "+html_name)
forrest@0
   937
 print (M_npp_H933)
forrest@0
   938
forrest@0
   939
 delete (u)
forrest@0
   940
 delete (v)
forrest@0
   941
 delete (uu)
forrest@0
   942
 delete (vv)
forrest@0
   943
;----------------------------------------------------------------------
forrest@0
   944
; histogram res
forrest@0
   945
  delete (resh)
forrest@0
   946
  resh                = True
forrest@0
   947
  resh@gsnMaximize    = True
forrest@0
   948
  resh@gsnDraw        = False
forrest@0
   949
  resh@gsnFrame       = False
forrest@0
   950
  resh@xyMarkLineMode = "Markers"
forrest@0
   951
  resh@xyMarkerSizeF  = 0.014
forrest@0
   952
  resh@xyMarker       = 16
forrest@0
   953
  resh@xyMarkerColors = (/"Brown","Blue"/)
forrest@0
   954
  resh@trYMinF        = min(mn_yvalues) - 10.
forrest@0
   955
  resh@trYMaxF        = max(mx_yvalues) + 10.
forrest@0
   956
forrest@0
   957
  resh@tiYAxisString  = "NPP (g C/m2/year)"
forrest@0
   958
  resh@tiXAxisString  = "Precipitation (m/year)"
forrest@0
   959
forrest@0
   960
  max_bar = new((/2,nx/),graphic)
forrest@0
   961
  min_bar = new((/2,nx/),graphic)
forrest@0
   962
  max_cap = new((/2,nx/),graphic)
forrest@0
   963
  min_cap = new((/2,nx/),graphic)
forrest@0
   964
forrest@0
   965
  lnres = True
forrest@0
   966
  line_colors = (/"brown","blue"/)
forrest@0
   967
forrest@0
   968
;**************************************************************************
forrest@0
   969
;(B)-3 histogram plot, ob 933 site 
forrest@0
   970
;**************************************************************************
forrest@0
   971
forrest@0
   972
  plot_name = "histogram_ob_933"
forrest@0
   973
  title     = "Observed 933 site"
forrest@0
   974
  resh@tiMainString  = title
forrest@0
   975
forrest@0
   976
  wks   = gsn_open_wks (plot_type,plot_name)    
forrest@0
   977
forrest@0
   978
  xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
forrest@0
   979
forrest@0
   980
;-------------------------------
forrest@0
   981
;Attach the vertical bar and the horizontal cap line 
forrest@0
   982
forrest@0
   983
  do nd=0,0
forrest@0
   984
    lnres@gsLineColor = line_colors(nd)
forrest@0
   985
    do i=0,nx-1
forrest@0
   986
     
forrest@0
   987
      if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0
   988
         .not.ismissing(mx_yvalues(nd,i))) then
forrest@0
   989
forrest@0
   990
; Attach the vertical bar, both above and below the marker.
forrest@0
   991
forrest@0
   992
        x1 = xvalues(nd,i)
forrest@0
   993
        y1 = yvalues(nd,i)
forrest@0
   994
        y2 = mn_yvalues(nd,i)
forrest@0
   995
        min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0
   996
forrest@0
   997
        y2 = mx_yvalues(nd,i)
forrest@0
   998
        max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0
   999
forrest@0
  1000
; Attach the horizontal cap line, both above and below the marker.
forrest@0
  1001
forrest@0
  1002
        x1 = xvalues(nd,i) - dx4
forrest@0
  1003
        x2 = xvalues(nd,i) + dx4
forrest@0
  1004
        y1 = mn_yvalues(nd,i)
forrest@0
  1005
        min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0
  1006
forrest@0
  1007
        y1 = mx_yvalues(nd,i)
forrest@0
  1008
        max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0
  1009
      end if
forrest@0
  1010
    end do
forrest@0
  1011
  end do
forrest@0
  1012
forrest@0
  1013
  draw(xy)
forrest@0
  1014
  frame(wks)
forrest@0
  1015
  delete (xy)
forrest@0
  1016
  clear (wks)
forrest@0
  1017
forrest@0
  1018
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1019
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1020
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1021
        "rm "+plot_name+"."+plot_type)
forrest@0
  1022
;****************************************************************************
forrest@0
  1023
;(B)-4 histogram plot, model vs ob 933 site
forrest@0
  1024
;**************************************************************************** 
forrest@0
  1025
forrest@0
  1026
  plot_name = "histogram_model_vs_ob_933"
forrest@0
  1027
  title     = model_name+ " vs Observed 933 site"
forrest@0
  1028
  resh@tiMainString  = title
forrest@0
  1029
forrest@0
  1030
  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
forrest@0
  1031
forrest@0
  1032
;-----------------------------
forrest@0
  1033
; Add a boxed legend using the more simple method, which won't have
forrest@0
  1034
; vertical lines going through the markers.
forrest@0
  1035
forrest@0
  1036
  resh@pmLegendDisplayMode    = "Always"
forrest@0
  1037
; resh@pmLegendWidthF         = 0.1
forrest@0
  1038
  resh@pmLegendWidthF         = 0.08
forrest@0
  1039
  resh@pmLegendHeightF        = 0.05
forrest@0
  1040
  resh@pmLegendOrthogonalPosF = -1.17
forrest@0
  1041
; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0
  1042
; resh@pmLegendParallelPosF   =  0.18
forrest@0
  1043
  resh@pmLegendParallelPosF   =  0.88  ;(rightward)
forrest@0
  1044
forrest@0
  1045
; resh@lgPerimOn              = False
forrest@0
  1046
  resh@lgLabelFontHeightF     = 0.015
forrest@0
  1047
  resh@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0
  1048
;-----------------------------
forrest@0
  1049
  tRes  = True
forrest@0
  1050
  tRes@txFontHeightF = 0.025
forrest@0
  1051
forrest@0
  1052
  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
forrest@0
  1053
forrest@0
  1054
  gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0
  1055
forrest@0
  1056
  xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
forrest@0
  1057
;-------------------------------
forrest@0
  1058
;Attach the vertical bar and the horizontal cap line 
forrest@0
  1059
forrest@0
  1060
  do nd=0,1
forrest@0
  1061
    lnres@gsLineColor = line_colors(nd)
forrest@0
  1062
    do i=0,nx-1
forrest@0
  1063
     
forrest@0
  1064
      if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0
  1065
         .not.ismissing(mx_yvalues(nd,i))) then
forrest@0
  1066
;
forrest@0
  1067
; Attach the vertical bar, both above and below the marker.
forrest@0
  1068
;
forrest@0
  1069
        x1 = xvalues(nd,i)
forrest@0
  1070
        y1 = yvalues(nd,i)
forrest@0
  1071
        y2 = mn_yvalues(nd,i)
forrest@0
  1072
        min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0
  1073
forrest@0
  1074
        y2 = mx_yvalues(nd,i)
forrest@0
  1075
        max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0
  1076
;
forrest@0
  1077
; Attach the horizontal cap line, both above and below the marker.
forrest@0
  1078
;
forrest@0
  1079
        x1 = xvalues(nd,i) - dx4
forrest@0
  1080
        x2 = xvalues(nd,i) + dx4
forrest@0
  1081
        y1 = mn_yvalues(nd,i)
forrest@0
  1082
        min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0
  1083
forrest@0
  1084
        y1 = mx_yvalues(nd,i)
forrest@0
  1085
        max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0
  1086
      end if
forrest@0
  1087
    end do
forrest@0
  1088
  end do
forrest@0
  1089
forrest@0
  1090
  draw(xy)
forrest@0
  1091
  frame(wks)
forrest@0
  1092
  delete(xy)
forrest@0
  1093
  clear(wks)
forrest@0
  1094
forrest@0
  1095
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1096
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1097
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1098
        "rm "+plot_name+"."+plot_type)
forrest@0
  1099
;***************************************************************************
forrest@0
  1100
;(C) global contour 
forrest@0
  1101
;***************************************************************************
forrest@0
  1102
forrest@0
  1103
;res
forrest@0
  1104
  resg                     = True             ; Use plot options
forrest@0
  1105
  resg@cnFillOn            = True             ; Turn on color fill
forrest@0
  1106
  resg@gsnSpreadColors     = True             ; use full colormap
forrest@0
  1107
; resg@cnFillMode          = "RasterFill"     ; Turn on raster color
forrest@0
  1108
; resg@lbLabelAutoStride   = True
forrest@0
  1109
  resg@cnLinesOn           = False            ; Turn off contourn lines
forrest@0
  1110
  resg@mpFillOn            = False            ; Turn off map fill
forrest@0
  1111
forrest@0
  1112
  resg@gsnSpreadColors      = True            ; use full colormap
forrest@0
  1113
  resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
forrest@0
  1114
  resg@cnMinLevelValF       = 0.              ; Min level
forrest@0
  1115
  resg@cnMaxLevelValF       = 2200.           ; Max level
forrest@0
  1116
  resg@cnLevelSpacingF      = 200.            ; interval
forrest@0
  1117
forrest@0
  1118
;****************************************************************************
forrest@0
  1119
;(C)-1 global contour plot, ob
forrest@0
  1120
;****************************************************************************
forrest@0
  1121
forrest@0
  1122
  delta = 0.00000000001
forrest@0
  1123
  nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
forrest@0
  1124
  
forrest@0
  1125
  plot_name = "global_ob"
forrest@0
  1126
  title     = "Observed MODIS MOD 17"
forrest@0
  1127
  resg@tiMainString  = title
forrest@0
  1128
forrest@0
  1129
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
  1130
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
  1131
forrest@0
  1132
  plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
forrest@0
  1133
   
forrest@0
  1134
  frame(wks)
forrest@0
  1135
  clear (wks)
forrest@0
  1136
forrest@0
  1137
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1138
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1139
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1140
        "rm "+plot_name+"."+plot_type)
forrest@0
  1141
;****************************************************************************
forrest@0
  1142
;(C)-2 global contour plot, model
forrest@0
  1143
;****************************************************************************
forrest@0
  1144
forrest@0
  1145
  plot_name = "global_model"
forrest@0
  1146
  title     = "Model "+ model_name
forrest@0
  1147
  resg@tiMainString  = title
forrest@0
  1148
forrest@0
  1149
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
  1150
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
  1151
forrest@0
  1152
  plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
forrest@0
  1153
   
forrest@0
  1154
  frame(wks)
forrest@0
  1155
  clear (wks)
forrest@0
  1156
forrest@0
  1157
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1158
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1159
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1160
        "rm "+plot_name+"."+plot_type)
forrest@0
  1161
;****************************************************************************
forrest@0
  1162
;(C)-3 global contour plot, model vs ob
forrest@0
  1163
;****************************************************************************
forrest@0
  1164
forrest@0
  1165
  plot_name = "global_model_vs_ob"
forrest@0
  1166
forrest@0
  1167
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
  1168
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
  1169
forrest@0
  1170
  delete (plot)
forrest@0
  1171
  plot=new(3,graphic)                        ; create graphic array
forrest@0
  1172
forrest@0
  1173
  resg@gsnFrame             = False          ; Do not draw plot 
forrest@0
  1174
  resg@gsnDraw              = False          ; Do not advance frame
forrest@0
  1175
forrest@0
  1176
; compute correlation coef and M score
forrest@0
  1177
forrest@0
  1178
  uu1 = ndtooned(nppmod)
forrest@0
  1179
  vv1 = ndtooned(nppglobe)
forrest@0
  1180
forrest@0
  1181
  delete (good) 
forrest@0
  1182
  good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
forrest@0
  1183
forrest@0
  1184
  ug = uu1(good)
forrest@0
  1185
  vg = vv1(good)
forrest@0
  1186
forrest@0
  1187
  ccrG = esccr(ug,vg,0)
forrest@0
  1188
; print (ccrG)
forrest@0
  1189
forrest@0
  1190
  MG   = (ccrG*ccrG)* 5.0
forrest@0
  1191
forrest@0
  1192
  M_npp_G = sprintf("%.2f", MG)
forrest@0
  1193
  system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \
forrest@0
  1194
         "mv -f "+html_new+" "+html_name)
forrest@0
  1195
  print (M_npp_G)
forrest@0
  1196
forrest@0
  1197
; plot correlation coef
forrest@0
  1198
forrest@0
  1199
  gRes  = True
forrest@0
  1200
  gRes@txFontHeightF = 0.02
forrest@0
  1201
  gRes@txAngleF = 90
forrest@0
  1202
forrest@0
  1203
  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
forrest@0
  1204
forrest@0
  1205
  gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0
  1206
;--------------------------------------------------------------------  
forrest@0
  1207
;(a) ob
forrest@0
  1208
forrest@0
  1209
  title     = "Observed MODIS MOD 17"
forrest@0
  1210
  resg@tiMainString  = title
forrest@0
  1211
forrest@0
  1212
  plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)       
forrest@0
  1213
forrest@0
  1214
;(b) model
forrest@0
  1215
forrest@0
  1216
  title     = "Model "+ model_name
forrest@0
  1217
  resg@tiMainString  = title
forrest@0
  1218
forrest@0
  1219
  plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg) 
forrest@0
  1220
forrest@0
  1221
;(c) model-ob
forrest@0
  1222
forrest@0
  1223
  zz = nppmod
forrest@0
  1224
  zz = nppmod - nppglobe
forrest@0
  1225
  title = "Model_"+model_name+" - Observed"
forrest@0
  1226
forrest@0
  1227
  resg@cnMinLevelValF  = -500           ; Min level
forrest@0
  1228
  resg@cnMaxLevelValF  =  500.          ; Max level
forrest@0
  1229
  resg@cnLevelSpacingF =  50.           ; interval
forrest@0
  1230
  resg@tiMainString    = title
forrest@0
  1231
forrest@0
  1232
  plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 
forrest@0
  1233
forrest@0
  1234
  pres                            = True        ; panel plot mods desired
forrest@0
  1235
; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around
forrest@0
  1236
                                                ; indiv. plots in panel
forrest@0
  1237
  pres@gsnMaximize                = True        ; fill the page
forrest@0
  1238
forrest@0
  1239
  gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot
forrest@0
  1240
forrest@0
  1241
  frame (wks)
forrest@0
  1242
  clear (wks)
forrest@0
  1243
  delete (plot)
forrest@0
  1244
forrest@0
  1245
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1246
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1247
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1248
        "rm "+plot_name+"."+plot_type)
forrest@0
  1249
;***************************************************************************
forrest@0
  1250
;(D)-1 zonal line plot, ob
forrest@0
  1251
;***************************************************************************
forrest@0
  1252
 
forrest@0
  1253
  vv     = zonalAve(nppglobe)
forrest@0
  1254
  vv@long_name = nppglobe@long_name
forrest@0
  1255
forrest@0
  1256
  plot_name = "zonal_ob"
forrest@0
  1257
  title     = "Observed MODIS MOD 17"
forrest@0
  1258
forrest@0
  1259
  resz = True
forrest@0
  1260
  resz@tiMainString  = title
forrest@0
  1261
forrest@0
  1262
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
  1263
forrest@0
  1264
  plot  = gsn_csm_xy (wks,ym,vv,resz)
forrest@0
  1265
   
forrest@0
  1266
  frame(wks)
forrest@0
  1267
  clear (wks)
forrest@0
  1268
forrest@0
  1269
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1270
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1271
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1272
        "rm "+plot_name+"."+plot_type)
forrest@0
  1273
;****************************************************************************
forrest@0
  1274
;(D)-2 zonal line plot, model vs ob
forrest@0
  1275
;****************************************************************************
forrest@0
  1276
forrest@0
  1277
  s  = new ((/2,dimsizes(ym)/), float)  
forrest@0
  1278
forrest@0
  1279
  s(0,:) = zonalAve(nppglobe) 
forrest@0
  1280
  s(1,:) = zonalAve(nppmod)
forrest@0
  1281
forrest@0
  1282
  s@long_name = nppglobe@long_name
forrest@0
  1283
;-------------------------------------------
forrest@0
  1284
; compute correlation coef and M score
forrest@0
  1285
forrest@0
  1286
  ccrZ = esccr(s(1,:), s(0,:),0)
forrest@0
  1287
; print (ccrZ)
forrest@0
  1288
forrest@0
  1289
  MZ   = (ccrZ*ccrZ)* 5.0
forrest@0
  1290
forrest@0
  1291
  M_npp_Z = sprintf("%.2f", MZ)
forrest@0
  1292
  system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \
forrest@0
  1293
         "mv -f "+html_new+" "+html_name)
forrest@0
  1294
  print (M_npp_Z)
forrest@0
  1295
;-------------------------------------------
forrest@0
  1296
  plot_name = "zonal_model_vs_ob"
forrest@0
  1297
  title     = "Zonal Average"
forrest@0
  1298
  resz@tiMainString  = title
forrest@0
  1299
forrest@0
  1300
  wks = gsn_open_wks (plot_type,plot_name)   
forrest@0
  1301
forrest@0
  1302
; resz@vpHeightF          = 0.4               ; change aspect ratio of plot
forrest@0
  1303
; resz@vpWidthF           = 0.7
forrest@0
  1304
forrest@0
  1305
  resz@xyMonoLineColor    = "False"           ; want colored lines
forrest@0
  1306
  resz@xyLineColors       = (/"black","red"/) ; colors chosen
forrest@0
  1307
; resz@xyLineThicknesses  = (/3.,3./)         ; line thicknesses
forrest@0
  1308
  resz@xyLineThicknesses  = (/2.,2./)         ; line thicknesses
forrest@0
  1309
  resz@xyDashPatterns     = (/0.,0./)         ; make all lines solid
forrest@0
  1310
forrest@0
  1311
  resz@tiYAxisString      = s@long_name       ; add a axis title    
forrest@0
  1312
  resz@txFontHeightF      = 0.0195            ; change title font heights
forrest@0
  1313
forrest@0
  1314
; Legent
forrest@0
  1315
  resz@pmLegendDisplayMode    = "Always"      ; turn on legend
forrest@0
  1316
  resz@pmLegendSide           = "Top"         ; Change location of 
forrest@0
  1317
; resz@pmLegendParallelPosF   = .45           ; move units right
forrest@0
  1318
  resz@pmLegendParallelPosF   = .82           ; move units right
forrest@0
  1319
  resz@pmLegendOrthogonalPosF = -0.4          ; move units down
forrest@0
  1320
forrest@0
  1321
  resz@pmLegendWidthF         = 0.10          ; Change width and
forrest@0
  1322
  resz@pmLegendHeightF        = 0.10          ; height of legend.
forrest@0
  1323
  resz@lgLabelFontHeightF     = .02           ; change font height
forrest@0
  1324
; resz@lgTitleOn              = True          ; turn on legend title
forrest@0
  1325
; resz@lgTitleString          = "Example"     ; create legend title
forrest@0
  1326
; resz@lgTitleFontHeightF     = .025          ; font of legend title
forrest@0
  1327
  resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
forrest@0
  1328
;--------------------------------------------------------------------
forrest@0
  1329
  zRes  = True
forrest@0
  1330
  zRes@txFontHeightF = 0.025
forrest@0
  1331
forrest@0
  1332
  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
forrest@0
  1333
forrest@0
  1334
  gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
forrest@0
  1335
;--------------------------------------------------------------------
forrest@0
  1336
  
forrest@0
  1337
  plot  = gsn_csm_xy (wks,ym,s,resz)      
forrest@0
  1338
forrest@0
  1339
  frame(wks)                                            
forrest@0
  1340
  clear (wks)
forrest@0
  1341
forrest@0
  1342
 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1343
        "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
  1344
        "rm "+plot_name+"-*."+plot_type_new+";"+ \
forrest@0
  1345
        "rm "+plot_name+"."+plot_type)
forrest@0
  1346
;***************************************************************************
forrest@0
  1347
; tar up output plots
forrest@0
  1348
;***************************************************************************
forrest@0
  1349
forrest@0
  1350
  temp_name = "npp." + model_name
forrest@0
  1351
  system("mkdir -p " + temp_name+";"+ \
forrest@0
  1352
         "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
forrest@0
  1353
         "mv table_*.html " + temp_name +";"+ \ 
forrest@0
  1354
         "mv *.png " + temp_name +";"+ \ 
forrest@0
  1355
         "tar cf "+ temp_name +".tar " + temp_name)
forrest@0
  1356
;***************************************************************************
forrest@0
  1357
end
forrest@0
  1358