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