1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/99.all.ncl.1 Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1667 @@
1.4 +;*****************************************************
1.5 +; combine scatter, histogram, global and zonal plots
1.6 +; compute all correlation coef and M score
1.7 +; add 1-to-1 line in scatter plots
1.8 +;*****************************************************
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.12 +;*****************************************************
1.13 +procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
1.14 + ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
1.15 + ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
1.16 + ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
1.17 + ,dx4[1]:numeric \
1.18 + )
1.19 +begin
1.20 +; Calculaee "nice" bins for binning the data in equally spaced ranges.
1.21 +; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
1.22 +; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
1.23 +
1.24 + nbins = 15 ; Number of bins to use.
1.25 +
1.26 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.27 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.28 + range = fspan(nicevals(0),nicevals(1),nvals)
1.29 +
1.30 +; Use this range information to grab all the values in a
1.31 +; particular range, and then take an average.
1.32 +
1.33 + nr = dimsizes(range)
1.34 + nx = nr-1
1.35 +; print (nx)
1.36 +
1.37 +; xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.38 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.39 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.40 + dx4 = dx/4 ; 1/4 of the range
1.41 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.42 +; yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.43 +; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.44 +; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.45 +
1.46 + do nd=0,1
1.47 +
1.48 +; See if we are doing model or observational data.
1.49 +
1.50 + if(nd.eq.0) then
1.51 + data = RAIN1_1D
1.52 + npp_data = NPP1_1D
1.53 + else
1.54 + data = RAIN2_1D
1.55 + npp_data = NPP2_1D
1.56 + end if
1.57 +
1.58 +; Loop through each range and check for values.
1.59 +
1.60 + do i=0,nr-2
1.61 + if (i.ne.(nr-2)) then
1.62 +; print("")
1.63 +; print("In range ["+range(i)+","+range(i+1)+")")
1.64 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.65 + else
1.66 +; print("")
1.67 +; print("In range ["+range(i)+",)")
1.68 + idx = ind(range(i).le.data)
1.69 + end if
1.70 +
1.71 +; Calculate average, and get min and max.
1.72 +
1.73 + if(.not.any(ismissing(idx))) then
1.74 + yvalues(nd,i) = avg(npp_data(idx))
1.75 + mn_yvalues(nd,i) = min(npp_data(idx))
1.76 + mx_yvalues(nd,i) = max(npp_data(idx))
1.77 + count = dimsizes(idx)
1.78 + else
1.79 + count = 0
1.80 + yvalues(nd,i) = yvalues@_FillValue
1.81 + mn_yvalues(nd,i) = yvalues@_FillValue
1.82 + mx_yvalues(nd,i) = yvalues@_FillValue
1.83 + end if
1.84 +
1.85 +; Print out information.
1.86 +
1.87 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.88 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.89 +
1.90 +
1.91 +; Clean up for next time in loop.
1.92 +
1.93 + delete(idx)
1.94 + end do
1.95 + delete(data)
1.96 + delete(npp_data)
1.97 + end do
1.98 +end
1.99 +;****************************************************************************
1.100 +; Main code.
1.101 +;****************************************************************************
1.102 +begin
1.103 +
1.104 + plot_type = "ps"
1.105 + plot_type_new = "png"
1.106 +
1.107 +;************************************************
1.108 +; read model data
1.109 +;************************************************
1.110 +
1.111 +;film = "i01.06cn_1798-2004_ANN_climo.nc"
1.112 +;model_name = "06cn"
1.113 +;model_grid = "T42"
1.114 +
1.115 +;film = "i01.06casa_1798-2004_ANN_climo.nc"
1.116 +;model_name = "06casa"
1.117 +;model_grid = "T42"
1.118 +
1.119 + film = "i01.10casa_1948-2004_ANN_climo.nc"
1.120 + model_name = "10casa"
1.121 + model_grid = "T42"
1.122 +
1.123 +;film = "i01.10cn_1948-2004_ANN_climo.nc"
1.124 +;model_name = "10cn"
1.125 +;model_grid = "T42"
1.126 +;--------------------------------------------------
1.127 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.128 + fm = addfile (dirm+film,"r")
1.129 +
1.130 + nppmod0 = fm->NPP
1.131 + rainmod0 = fm->RAIN
1.132 + xm = fm->lon
1.133 + ym = fm->lat
1.134 +
1.135 +;************************************************
1.136 +; read ob data
1.137 +;************************************************
1.138 +
1.139 +;(1) data at 81 sites
1.140 + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.141 + fil81 = "data.81.nc"
1.142 + f81 = addfile (dir81+fil81,"r")
1.143 +
1.144 + id81 = f81->SITE_ID
1.145 + npp81 = f81->TNPP_C
1.146 + rain81 = tofloat(f81->PREC_ANN)
1.147 + x81 = f81->LONG_DD
1.148 + y81 = f81->LAT_DD
1.149 +
1.150 + id81@long_name = "SITE_ID"
1.151 +
1.152 +;change longitude from (-180,180) to (0,360)
1.153 +;for model data interpolation
1.154 +
1.155 + do i= 0,dimsizes(x81)-1
1.156 + if (x81(i) .lt. 0.) then
1.157 + x81(i) = x81(i)+ 360.
1.158 + end if
1.159 + end do
1.160 +;print (x81)
1.161 +;-------------------------------------
1.162 +;(2) data at 933 sites
1.163 + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.164 + fil933 = "data.933.nc"
1.165 + f933 = addfile (dir933+fil933,"r")
1.166 +
1.167 + id933 = f933->SITE_ID
1.168 + npp933 = f933->TNPP_C
1.169 + rain933 = f933->PREC
1.170 + x933 = f933->LONG_DD
1.171 + y933 = f933->LAT_DD
1.172 +
1.173 + id933@long_name = "SITE_ID"
1.174 +
1.175 +;change longitude from (-180,180) to (0,360)
1.176 +;for model data interpolation
1.177 +
1.178 + do i= 0,dimsizes(x933)-1
1.179 + if (x933(i) .lt. 0.) then
1.180 + x933(i) = x933(i)+ 360.
1.181 + end if
1.182 + end do
1.183 +;print (x933)
1.184 +;----------------------------------------
1.185 +;(3) global data, interpolated into model grid
1.186 + dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.187 + filglobe = "Npp_"+model_grid+"_mean.nc"
1.188 + fglobe = addfile (dirglobe+filglobe,"r")
1.189 +
1.190 + nppglobe0 = fglobe->NPP
1.191 + nppglobe = nppglobe0
1.192 +
1.193 +;***********************************************************************
1.194 +; interpolate model data into ob sites
1.195 +;***********************************************************************
1.196 +
1.197 + nppmod = nppmod0(0,:,:)
1.198 + rainmod = rainmod0(0,:,:)
1.199 + delete (nppmod0)
1.200 + delete (rainmod0)
1.201 +
1.202 + nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
1.203 +
1.204 + nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
1.205 +
1.206 + rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
1.207 +
1.208 + rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
1.209 +
1.210 +;**********************************************************
1.211 +; unified units
1.212 +;**********************************************************
1.213 +; Units for these variables are:
1.214 +;
1.215 +; rain81 : mm/year
1.216 +; rainmod : mm/s
1.217 +; npp81 : g C/m^2/year
1.218 +; nppmod81: g C/m^2/s
1.219 +; nppglobe: g C/m^2/year
1.220 +;
1.221 +; We want to convert these to "m/year" and "g C/m^2/year".
1.222 +
1.223 + nsec_per_year = 60*60*24*365
1.224 +
1.225 + rain81 = rain81 / 1000.
1.226 + rainmod81 = (rainmod81/ 1000.) * nsec_per_year
1.227 + nppmod81 = nppmod81 * nsec_per_year
1.228 +
1.229 + rain933 = rain933 / 1000.
1.230 + rainmod933 = (rainmod933/ 1000.) * nsec_per_year
1.231 + nppmod933 = nppmod933 * nsec_per_year
1.232 +
1.233 + nppmod = nppmod * nsec_per_year
1.234 +
1.235 + npp81@units = "gC/m^2/yr"
1.236 + nppmod81@units = "gC/m^2/yr"
1.237 + npp933@units = "gC/m^2/yr"
1.238 + nppmod933@units = "gC/m^2/yr"
1.239 + nppmod@units = "gC/m^2/yr"
1.240 + nppglobe@units = "gC/m^2/yr"
1.241 + rain81@units = "m/yr"
1.242 + rainmod81@units = "m/yr"
1.243 + rain933@units = "m/yr"
1.244 + rainmod933@units = "m/yr"
1.245 +
1.246 + npp81@long_name = "NPP (gC/m2/year)"
1.247 + npp933@long_name = "NPP (gC/m2/year)"
1.248 + nppmod81@long_name = "NPP (gC/m2/year)"
1.249 + nppmod933@long_name = "NPP (gC/m2/year)"
1.250 + nppmod@long_name = "NPP (gC/m2/year)"
1.251 + nppglobe@long_name = "NPP (gC/m2/year)"
1.252 + rain81@long_name = "PREC (m/year)"
1.253 + rain933@long_name = "PREC (m/year)"
1.254 + rainmod81@long_name = "PREC (m/year)"
1.255 + rainmod933@long_name = "PREC (m/year)"
1.256 +
1.257 +;*******************************************************************
1.258 +;(A)-1 table of site81 (1) (the first 41 sites)
1.259 +;*******************************************************************
1.260 +
1.261 + table_length = 0.95
1.262 +
1.263 +; table header name
1.264 + table_header_name = "Site ID"
1.265 +
1.266 +; column (not including header column)
1.267 + col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
1.268 + ncol = dimsizes(col_header)
1.269 +
1.270 +; row (not including header row)
1.271 + nrow = 41
1.272 + row_header = new ((/nrow/),string )
1.273 + row_header(0:nrow-1) = id81(0:nrow-1)
1.274 +
1.275 +; Table header
1.276 + ncr1 = (/1,1/) ; 1 row, 1 column
1.277 + x1 = (/0.005,0.15/) ; Start and end X
1.278 + y1 = (/0.900,0.95/) ; Start and end Y
1.279 + text1 = table_header_name
1.280 + res1 = True
1.281 + res1@txFontHeightF = 0.015
1.282 + res1@gsFillColor = "CornFlowerBlue"
1.283 +
1.284 +; Column header (equally space in x2)
1.285 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.286 + x2 = (/x1(1),0.995/) ; start from end of x1
1.287 + y2 = y1 ; same as y1
1.288 + text2 = col_header
1.289 + res2 = True
1.290 + res2@txFontHeightF = 0.015
1.291 + res2@gsFillColor = "Gray"
1.292 +
1.293 +; Row header (equally space in y2)
1.294 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.295 + x3 = x1 ; same as x1
1.296 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.297 + text3 = row_header
1.298 + res3 = True
1.299 + res3@txFontHeightF = 0.010
1.300 + res3@gsFillColor = "Gray"
1.301 +
1.302 +; Main table body
1.303 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.304 + x4 = x2 ; Start and end x
1.305 + y4 = y3 ; Start and end Y
1.306 + text4 = new((/nrow,ncol/),string)
1.307 +
1.308 + color_fill4 = new((/nrow,ncol/),string)
1.309 + color_fill4 = "white"
1.310 +; color_fill4(:,ncol-1) = "grey"
1.311 +; color_fill4(nrow-1,:) = "green"
1.312 +
1.313 + res4 = True ; Set up resource list
1.314 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.315 + res4@txFontHeightF = 0.015
1.316 + res4@gsFillColor = color_fill4
1.317 +
1.318 + delete (color_fill4)
1.319 +;-------------------------------------------------------------------
1.320 +; for table value
1.321 +
1.322 + do n = 0,nrow-1
1.323 + text4(n,0) = sprintf("%5.2f", y81(n))
1.324 + text4(n,1) = sprintf("%5.2f", x81(n))
1.325 + text4(n,2) = sprintf("%5.2f", npp81(n))
1.326 + text4(n,3) = sprintf("%5.2f", rain81(n))
1.327 + end do
1.328 +;-------------------------------------------------------------------
1.329 +
1.330 + plot_name = "table_site81_ob1"
1.331 +
1.332 + wks = gsn_open_wks (plot_type,plot_name)
1.333 +;------------------------------------------
1.334 +; for table title
1.335 +
1.336 + gRes = True
1.337 + gRes@txFontHeightF = 0.02
1.338 +; gRes@txAngleF = 90
1.339 +
1.340 + title_text = "Observation at 81 Sites (1)"
1.341 +
1.342 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.343 +;------------------------------------------
1.344 +
1.345 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.346 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.347 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.348 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.349 +
1.350 + frame(wks)
1.351 + clear (wks)
1.352 +
1.353 + delete (row_header)
1.354 + delete (res3)
1.355 + delete (ncr3)
1.356 + delete (text3)
1.357 + delete (res4)
1.358 + delete (ncr4)
1.359 + delete (text4)
1.360 +
1.361 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.362 + system("rm "+plot_name+"."+plot_type)
1.363 + system("rm "+plot_name+"-1."+plot_type_new)
1.364 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.365 +
1.366 +;*******************************************************************
1.367 +;(A)-2 table of site81 (2) (the last 40 sites)
1.368 +;*******************************************************************
1.369 +
1.370 + table_length = 0.95
1.371 +
1.372 +; table header name
1.373 + table_header_name = "Site ID"
1.374 +
1.375 +; column (not including header column)
1.376 + col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
1.377 + ncol = dimsizes(col_header)
1.378 +
1.379 +; row (not including header row)
1.380 + nrow = 40
1.381 + row_header = new ((/nrow/),string )
1.382 + row_header(0:nrow-1) = id81(41:41+nrow-1)
1.383 +
1.384 +; Table header
1.385 + ncr1 = (/1,1/) ; 1 row, 1 column
1.386 + x1 = (/0.005,0.15/) ; Start and end X
1.387 + y1 = (/0.900,0.95/) ; Start and end Y
1.388 + text1 = table_header_name
1.389 + res1 = True
1.390 + res1@txFontHeightF = 0.015
1.391 + res1@gsFillColor = "CornFlowerBlue"
1.392 +
1.393 +; Column header (equally space in x2)
1.394 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.395 + x2 = (/x1(1),0.995/) ; start from end of x1
1.396 + y2 = y1 ; same as y1
1.397 + text2 = col_header
1.398 + res2 = True
1.399 + res2@txFontHeightF = 0.015
1.400 + res2@gsFillColor = "Gray"
1.401 +
1.402 +; Row header (equally space in y2)
1.403 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.404 + x3 = x1 ; same as x1
1.405 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.406 + text3 = row_header
1.407 + res3 = True
1.408 + res3@txFontHeightF = 0.010
1.409 + res3@gsFillColor = "Gray"
1.410 +
1.411 +; Main table body
1.412 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.413 + x4 = x2 ; Start and end x
1.414 + y4 = y3 ; Start and end Y
1.415 + text4 = new((/nrow,ncol/),string)
1.416 +
1.417 + color_fill4 = new((/nrow,ncol/),string)
1.418 + color_fill4 = "white"
1.419 +; color_fill4(:,ncol-1) = "grey"
1.420 +; color_fill4(nrow-1,:) = "green"
1.421 +
1.422 + res4 = True ; Set up resource list
1.423 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.424 + res4@txFontHeightF = 0.015
1.425 + res4@gsFillColor = color_fill4
1.426 +
1.427 + delete (color_fill4)
1.428 +;-------------------------------------------------------------------
1.429 +; for table value
1.430 +
1.431 + do n = 0,nrow-1
1.432 + text4(n,0) = sprintf("%5.2f", y81(n+41))
1.433 + text4(n,1) = sprintf("%5.2f", x81(n+41))
1.434 + text4(n,2) = sprintf("%5.2f", npp81(n+41))
1.435 + text4(n,3) = sprintf("%5.2f", rain81(n+41))
1.436 + end do
1.437 +;-------------------------------------------------------------------
1.438 +
1.439 + plot_name = "table_site81_ob2"
1.440 +
1.441 + wks = gsn_open_wks (plot_type,plot_name)
1.442 +;------------------------------------------
1.443 +; for table title
1.444 +
1.445 + gRes = True
1.446 + gRes@txFontHeightF = 0.02
1.447 +; gRes@txAngleF = 90
1.448 +
1.449 + title_text = "Observation at 81 Sites (2)"
1.450 +
1.451 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.452 +;------------------------------------------
1.453 +
1.454 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.455 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.456 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.457 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.458 +
1.459 + frame(wks)
1.460 + clear (wks)
1.461 +
1.462 + delete (row_header)
1.463 + delete (res3)
1.464 + delete (ncr3)
1.465 + delete (text3)
1.466 + delete (res4)
1.467 + delete (ncr4)
1.468 + delete (text4)
1.469 +
1.470 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.471 + system("rm "+plot_name+"."+plot_type)
1.472 + system("rm "+plot_name+"-1."+plot_type_new)
1.473 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.474 +
1.475 +;**************************************************************************
1.476 +;(A)-3 scatter plot, ob 933
1.477 +;**************************************************************************
1.478 +
1.479 + plot_name = "scatter_ob_933"
1.480 + title = "Observed 933 sites"
1.481 +
1.482 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.483 + res = True ; plot mods desired
1.484 + res@tiMainString = title ; add title
1.485 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.486 + res@xyMarkers = 16 ; choose type of marker
1.487 + res@xyMarkerColor = "red" ; Marker color
1.488 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.489 + res@tmLabelAutoStride = True ; nice tick mark labels
1.490 +
1.491 + plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
1.492 + frame(wks)
1.493 + clear (wks)
1.494 +
1.495 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.496 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.497 + system("rm "+plot_name+"-*."+plot_type_new)
1.498 + system("rm "+plot_name+"."+plot_type)
1.499 +
1.500 +;**************************************************************************
1.501 +;(A)-4 scatter plot, model 933
1.502 +;**************************************************************************
1.503 +
1.504 + plot_name = "scatter_model_933"
1.505 + title = "Model "+ model_name +" 933 sites"
1.506 +
1.507 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.508 + res = True ; plot mods desired
1.509 + res@tiMainString = title ; add title
1.510 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.511 + res@xyMarkers = 16 ; choose type of marker
1.512 + res@xyMarkerColor = "red" ; Marker color
1.513 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.514 + res@tmLabelAutoStride = True ; nice tick mark labels
1.515 +
1.516 + plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
1.517 + frame(wks)
1.518 + clear (wks)
1.519 +
1.520 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.521 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.522 + system("rm "+plot_name+"-*."+plot_type_new)
1.523 + system("rm "+plot_name+"."+plot_type)
1.524 +
1.525 +;***************************************************************************
1.526 +;(A)-5 scatter plot, model vs ob 81
1.527 +;***************************************************************************
1.528 +
1.529 + plot_name = "scatter_model_vs_ob_81"
1.530 + title = model_name +" vs ob 81 sites"
1.531 +
1.532 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.533 + res = True ; plot mods desired
1.534 + res@tiMainString = title ; add title
1.535 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.536 + res@xyMarkers = 16 ; choose type of marker
1.537 + res@xyMarkerColor = "red" ; Marker color
1.538 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.539 + res@tmLabelAutoStride = True ; nice tick mark labels
1.540 +
1.541 + res@gsnDraw = False
1.542 + res@gsnFrame = False ; don't advance frame yet
1.543 +;-------------------------------
1.544 +;compute correlation coef. and M
1.545 + ccr81 = esccr(nppmod81,npp81,0)
1.546 +;print (ccr81)
1.547 +
1.548 +;new eq
1.549 + bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
1.550 + M81 = (1. - (bias/dimsizes(y81)))*5.
1.551 + print (M81)
1.552 +
1.553 + tRes = True
1.554 + tRes@txFontHeightF = 0.025
1.555 +
1.556 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.557 +
1.558 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.559 +;--------------------------------
1.560 + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
1.561 +;-------------------------------
1.562 +; add polyline
1.563 +
1.564 + dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
1.565 +;-------------------------------
1.566 + draw (plot)
1.567 + frame(wks)
1.568 + clear (wks)
1.569 + delete (dum)
1.570 +
1.571 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.572 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.573 + system("rm "+plot_name+"-*."+plot_type_new)
1.574 + system("rm "+plot_name+"."+plot_type)
1.575 +
1.576 +;***************************************************************************
1.577 +;(A)-6 scatter plot, model vs ob 933
1.578 +;***************************************************************************
1.579 +
1.580 + plot_name = "scatter_model_vs_ob_933"
1.581 + title = model_name +" vs ob 933 sites"
1.582 +
1.583 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.584 + res = True ; plot mods desired
1.585 + res@tiMainString = title ; add title
1.586 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.587 + res@xyMarkers = 16 ; choose type of marker
1.588 + res@xyMarkerColor = "red" ; Marker color
1.589 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.590 + res@tmLabelAutoStride = True ; nice tick mark labels
1.591 +
1.592 + res@gsnDraw = False
1.593 + res@gsnFrame = False ; don't advance frame yet
1.594 +;-------------------------------
1.595 +;compute correlation coef. and M
1.596 + ccr933 = esccr(nppmod933,npp933,0)
1.597 +;print (ccr933)
1.598 +
1.599 +;new eq
1.600 + bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
1.601 + M933 = (1. - (bias/dimsizes(y933)))*5.
1.602 + print (M933)
1.603 +
1.604 + tRes = True
1.605 + tRes@txFontHeightF = 0.025
1.606 +
1.607 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.608 +
1.609 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.610 +;--------------------------------
1.611 + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
1.612 +;-------------------------------
1.613 +; add polyline
1.614 +
1.615 + dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
1.616 +;-------------------------------
1.617 + draw (plot)
1.618 + frame(wks)
1.619 + clear (wks)
1.620 + delete (dum)
1.621 +
1.622 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.623 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.624 + system("rm "+plot_name+"-*."+plot_type_new)
1.625 + system("rm "+plot_name+"."+plot_type)
1.626 +
1.627 +;*******************************************************************
1.628 +;(A)-7 table of site81 (1) (the first 41 sites), model vs ob
1.629 +;*******************************************************************
1.630 +
1.631 + table_length = 0.95
1.632 +
1.633 +; table header name
1.634 + table_header_name = "Site ID"
1.635 +
1.636 +; row (not including header row)
1.637 + nrow = 41
1.638 + row_header = new ((/nrow/),string )
1.639 + row_header(0:nrow-1) = id81(0:nrow-1)
1.640 +
1.641 +; Table header
1.642 + ncr1 = (/1,1/) ; 1 row, 1 column
1.643 + x1 = (/0.005,0.15/) ; Start and end X
1.644 + y1 = (/0.85 ,0.95/) ; Start and end Y
1.645 + text1 = table_header_name
1.646 + res1 = True
1.647 + res1@txFontHeightF = 0.015
1.648 + res1@gsFillColor = "CornFlowerBlue"
1.649 +
1.650 +; Column header1 (equally space in x2)
1.651 +
1.652 + delete (col_header)
1.653 + delete (ncr2)
1.654 + delete (text2)
1.655 +
1.656 + col_header = (/"Latitude","Longitude"/)
1.657 + ncol = dimsizes(col_header)
1.658 +
1.659 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.660 + x2 = (/x1(1),0.35/) ; start from end of x1
1.661 + y2 = y1 ; same as y1
1.662 + text2 = col_header
1.663 + res2 = True
1.664 + res2@txFontHeightF = 0.015
1.665 + res2@gsFillColor = "Gray"
1.666 +
1.667 +; Column header2 (equally space in x2)
1.668 +
1.669 + col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
1.670 + col_header2 = (/"ob",model_name \
1.671 + ,"ob",model_name \
1.672 + /)
1.673 +
1.674 + ncol1 = dimsizes(col_header1)
1.675 + ncol2 = dimsizes(col_header2)
1.676 +
1.677 + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
1.678 + x21 = (/x2(1),0.995/) ; start from end of x2
1.679 + yhalf = (y1(0)+y1(1))*0.5
1.680 + y21 = (/yhalf,y1(1)/) ; top half of y1
1.681 + text21 = col_header1
1.682 + res21 = True
1.683 + res21@txFontHeightF = 0.015
1.684 + res21@gsFillColor = "Gray"
1.685 +
1.686 + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
1.687 + x22 = x21 ; start from end of x1
1.688 + y22 = (/y1(0),yhalf/) ; bottom half of y1
1.689 + text22 = col_header2
1.690 + res22 = True
1.691 + res22@txFontHeightF = 0.015
1.692 + res22@gsFillColor = "Gray"
1.693 +
1.694 +; Row header (equally space in y2)
1.695 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.696 + x3 = x1 ; same as x1
1.697 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.698 + text3 = row_header
1.699 + res3 = True
1.700 + res3@txFontHeightF = 0.010
1.701 + res3@gsFillColor = "Gray"
1.702 +
1.703 +; Main table body-1
1.704 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.705 + x4 = x2 ; Start and end x
1.706 + y4 = y3 ; Start and end Y
1.707 + text4 = new((/nrow,ncol/),string)
1.708 +
1.709 + color_fill4 = new((/nrow,ncol/),string)
1.710 + color_fill4 = "white"
1.711 +; color_fill4(:,ncol-1) = "grey"
1.712 +; color_fill4(nrow-1,:) = "green"
1.713 +
1.714 + res4 = True ; Set up resource list
1.715 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.716 + res4@txFontHeightF = 0.015
1.717 + res4@gsFillColor = color_fill4
1.718 +
1.719 + delete (color_fill4)
1.720 +
1.721 +; Main table body-2
1.722 + ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
1.723 + x42 = x21 ; Start and end x
1.724 + y42 = y3 ; Start and end Y
1.725 + text42 = new((/nrow,ncol2/),string)
1.726 +
1.727 + color_fill42 = new((/nrow,ncol2/),string)
1.728 + color_fill42 = "white"
1.729 +; color_fill42(:,ncol-1) = "grey"
1.730 +; color_fill42(nrow-1,:) = "green"
1.731 +
1.732 + res42 = True ; Set up resource list
1.733 +; res42@gsnDebug = True ; Useful to print NDC row,col values used.
1.734 + res42@txFontHeightF = 0.015
1.735 + res42@gsFillColor = color_fill42
1.736 +
1.737 + delete (color_fill42)
1.738 +;-------------------------------------------------------------------
1.739 +; for table value
1.740 +
1.741 + do n = 0,nrow-1
1.742 + text4(n,0) = sprintf("%5.2f", y81(n))
1.743 + text4(n,1) = sprintf("%5.2f", x81(n))
1.744 + end do
1.745 +
1.746 + do n = 0,nrow-1
1.747 + text42(n,0) = sprintf("%5.2f", npp81(n))
1.748 + text42(n,1) = sprintf("%5.2f", nppmod81(n))
1.749 + text42(n,2) = sprintf("%5.2f", rain81(n))
1.750 + text42(n,3) = sprintf("%5.2f", rainmod81(n))
1.751 + end do
1.752 +;---------------------------------------------------------------------------
1.753 +
1.754 + plot_name = "table_site81_model_vs_ob1"
1.755 +
1.756 + wks = gsn_open_wks (plot_type,plot_name)
1.757 +;------------------------------------------
1.758 +; for table title
1.759 +
1.760 + gRes = True
1.761 + gRes@txFontHeightF = 0.02
1.762 +; gRes@txAngleF = 90
1.763 +
1.764 + title_text = "Model vs Observation at 81 Sites (1)"
1.765 +
1.766 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.767 +;------------------------------------------
1.768 +
1.769 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.770 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.771 + gsn_table(wks,ncr21,x21,y21,text21,res21)
1.772 + gsn_table(wks,ncr22,x22,y22,text22,res22)
1.773 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.774 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.775 + gsn_table(wks,ncr42,x42,y42,text42,res42)
1.776 +
1.777 + frame(wks)
1.778 + clear (wks)
1.779 +
1.780 + delete (col_header)
1.781 + delete (row_header)
1.782 + delete (res1)
1.783 + delete (ncr1)
1.784 + delete (text1)
1.785 + delete (res2)
1.786 + delete (ncr2)
1.787 + delete (text2)
1.788 + delete (res21)
1.789 + delete (ncr21)
1.790 + delete (text21)
1.791 + delete (res22)
1.792 + delete (ncr22)
1.793 + delete (text22)
1.794 + delete (res3)
1.795 + delete (ncr3)
1.796 + delete (text3)
1.797 + delete (res4)
1.798 + delete (ncr4)
1.799 + delete (text4)
1.800 + delete (res42)
1.801 + delete (ncr42)
1.802 + delete (text42)
1.803 +
1.804 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.805 + system("rm "+plot_name+"."+plot_type)
1.806 + system("rm "+plot_name+"-1."+plot_type_new)
1.807 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.808 +
1.809 +;*******************************************************************
1.810 +;(A)-8 table of site81 (2) (the last 40 sites), model vs ob
1.811 +;*******************************************************************
1.812 +
1.813 + table_length = 0.95
1.814 +
1.815 +; table header name
1.816 + table_header_name = "Site ID"
1.817 +
1.818 +; row (not including header row)
1.819 + nrow = 40
1.820 + row_header = new ((/nrow/),string )
1.821 + row_header(0:nrow-1) = id81(41:41+nrow-1)
1.822 +
1.823 +; Table header
1.824 + ncr1 = (/1,1/) ; 1 row, 1 column
1.825 + x1 = (/0.005,0.15/) ; Start and end X
1.826 + y1 = (/0.85 ,0.95/) ; Start and end Y
1.827 + text1 = table_header_name
1.828 + res1 = True
1.829 + res1@txFontHeightF = 0.015
1.830 + res1@gsFillColor = "CornFlowerBlue"
1.831 +
1.832 +; Column header1 (equally space in x2)
1.833 +
1.834 + col_header = (/"Latitude","Longitude"/)
1.835 + ncol = dimsizes(col_header)
1.836 +
1.837 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.838 + x2 = (/x1(1),0.35/) ; start from end of x1
1.839 + y2 = y1 ; same as y1
1.840 + text2 = col_header
1.841 + res2 = True
1.842 + res2@txFontHeightF = 0.015
1.843 + res2@gsFillColor = "Gray"
1.844 +
1.845 +; Column header2 (equally space in x2)
1.846 +
1.847 +; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
1.848 + col_header2 = (/"ob",model_name \
1.849 + ,"ob",model_name \
1.850 + /)
1.851 +
1.852 + ncol1 = dimsizes(col_header1)
1.853 + ncol2 = dimsizes(col_header2)
1.854 +
1.855 + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
1.856 + x21 = (/x2(1),0.995/) ; start from end of x2
1.857 + yhalf = (y1(0)+y1(1))*0.5
1.858 + y21 = (/yhalf,y1(1)/) ; top half of y1
1.859 + text21 = col_header1
1.860 + res21 = True
1.861 + res21@txFontHeightF = 0.015
1.862 + res21@gsFillColor = "Gray"
1.863 +
1.864 + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
1.865 + x22 = x21 ; start from end of x1
1.866 + y22 = (/y1(0),yhalf/) ; bottom half of y1
1.867 + text22 = col_header2
1.868 + res22 = True
1.869 + res22@txFontHeightF = 0.015
1.870 + res22@gsFillColor = "Gray"
1.871 +
1.872 +; Row header (equally space in y2)
1.873 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.874 + x3 = x1 ; same as x1
1.875 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.876 + text3 = row_header
1.877 + res3 = True
1.878 + res3@txFontHeightF = 0.010
1.879 + res3@gsFillColor = "Gray"
1.880 +
1.881 +; Main table body-1
1.882 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.883 + x4 = x2 ; Start and end x
1.884 + y4 = y3 ; Start and end Y
1.885 + text4 = new((/nrow,ncol/),string)
1.886 +
1.887 + color_fill4 = new((/nrow,ncol/),string)
1.888 + color_fill4 = "white"
1.889 +; color_fill4(:,ncol-1) = "grey"
1.890 +; color_fill4(nrow-1,:) = "green"
1.891 +
1.892 + res4 = True ; Set up resource list
1.893 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.894 + res4@txFontHeightF = 0.015
1.895 + res4@gsFillColor = color_fill4
1.896 +
1.897 + delete (color_fill4)
1.898 +
1.899 +; Main table body-2
1.900 + ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
1.901 + x42 = x21 ; Start and end x
1.902 + y42 = y3 ; Start and end Y
1.903 + text42 = new((/nrow,ncol2/),string)
1.904 +
1.905 + color_fill42 = new((/nrow,ncol2/),string)
1.906 + color_fill42 = "white"
1.907 +; color_fill42(:,ncol-1) = "grey"
1.908 +; color_fill42(nrow-1,:) = "green"
1.909 +
1.910 + res42 = True ; Set up resource list
1.911 +; res42@gsnDebug = True ; Useful to print NDC row,col values used.
1.912 + res42@txFontHeightF = 0.015
1.913 + res42@gsFillColor = color_fill42
1.914 +
1.915 + delete (color_fill42)
1.916 +;-------------------------------------------------------------------
1.917 +; for table value
1.918 +
1.919 + do n = 0,nrow-1
1.920 + text4(n,0) = sprintf("%5.2f", y81(n+41))
1.921 + text4(n,1) = sprintf("%5.2f", x81(n+41))
1.922 + end do
1.923 +
1.924 + do n = 0,nrow-1
1.925 + text42(n,0) = sprintf("%5.2f", npp81(n+41))
1.926 + text42(n,1) = sprintf("%5.2f", nppmod81(n+41))
1.927 + text42(n,2) = sprintf("%5.2f", rain81(n+41))
1.928 + text42(n,3) = sprintf("%5.2f", rainmod81(n+41))
1.929 + end do
1.930 +;-------------------------------------------------------------------
1.931 +
1.932 + plot_name = "table_site81_model_vs_ob2"
1.933 +
1.934 + wks = gsn_open_wks (plot_type,plot_name)
1.935 +;------------------------------------------
1.936 +; for table title
1.937 +
1.938 + gRes = True
1.939 + gRes@txFontHeightF = 0.02
1.940 +; gRes@txAngleF = 90
1.941 +
1.942 + title_text = "Model vs Observation at 81 Sites (2)"
1.943 +
1.944 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.945 +;------------------------------------------
1.946 +
1.947 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.948 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.949 + gsn_table(wks,ncr21,x21,y21,text21,res21)
1.950 + gsn_table(wks,ncr22,x22,y22,text22,res22)
1.951 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.952 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.953 + gsn_table(wks,ncr42,x42,y42,text42,res42)
1.954 +
1.955 + frame(wks)
1.956 + clear (wks)
1.957 +
1.958 + delete (col_header)
1.959 + delete (row_header)
1.960 + delete (res1)
1.961 + delete (ncr1)
1.962 + delete (text1)
1.963 + delete (res2)
1.964 + delete (ncr2)
1.965 + delete (text2)
1.966 + delete (res21)
1.967 + delete (ncr21)
1.968 + delete (text21)
1.969 + delete (res22)
1.970 + delete (ncr22)
1.971 + delete (text22)
1.972 + delete (res3)
1.973 + delete (ncr3)
1.974 + delete (text3)
1.975 + delete (res4)
1.976 + delete (ncr4)
1.977 + delete (text4)
1.978 + delete (res42)
1.979 + delete (ncr42)
1.980 + delete (text42)
1.981 +
1.982 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.983 + system("rm "+plot_name+"."+plot_type)
1.984 + system("rm "+plot_name+"-1."+plot_type_new)
1.985 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.986 +
1.987 +;**************************************************************************
1.988 +;(B) histogram 81
1.989 +;**************************************************************************
1.990 +;get data
1.991 + RAIN1_1D = ndtooned(rain81)
1.992 + RAIN2_1D = ndtooned(rainmod81)
1.993 + NPP1_1D = ndtooned(npp81)
1.994 + NPP2_1D = ndtooned(nppmod81)
1.995 +
1.996 +; number of bin
1.997 + nx = 8
1.998 +
1.999 + xvalues = new((/2,nx/),float)
1.1000 + yvalues = new((/2,nx/),float)
1.1001 + mn_yvalues = new((/2,nx/),float)
1.1002 + mx_yvalues = new((/2,nx/),float)
1.1003 + dx4 = new((/1/),float)
1.1004 +
1.1005 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.1006 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.1007 +
1.1008 +;----------------------------------------
1.1009 +;compute correlation coeff and M score
1.1010 +
1.1011 + u = yvalues(0,:)
1.1012 + v = yvalues(1,:)
1.1013 +
1.1014 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.1015 + uu = u(good)
1.1016 + vv = v(good)
1.1017 +
1.1018 + ccr81h = esccr(uu,vv,0)
1.1019 +;print (ccr81h)
1.1020 +
1.1021 +;new eq
1.1022 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.1023 + M81h = (1.- (bias/dimsizes(uu)))*5.
1.1024 + print (M81h)
1.1025 +
1.1026 + delete (u)
1.1027 + delete (v)
1.1028 + delete (uu)
1.1029 + delete (vv)
1.1030 +;----------------------------------------------------------------------
1.1031 +; histogram res
1.1032 + resh = True
1.1033 + resh@gsnMaximize = True
1.1034 + resh@gsnDraw = False
1.1035 + resh@gsnFrame = False
1.1036 + resh@xyMarkLineMode = "Markers"
1.1037 + resh@xyMarkerSizeF = 0.014
1.1038 + resh@xyMarker = 16
1.1039 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.1040 + resh@trYMinF = min(mn_yvalues) - 10.
1.1041 + resh@trYMaxF = max(mx_yvalues) + 10.
1.1042 +
1.1043 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.1044 + resh@tiXAxisString = "Precipitation (m/year)"
1.1045 +
1.1046 + max_bar = new((/2,nx/),graphic)
1.1047 + min_bar = new((/2,nx/),graphic)
1.1048 + max_cap = new((/2,nx/),graphic)
1.1049 + min_cap = new((/2,nx/),graphic)
1.1050 +
1.1051 + lnres = True
1.1052 + line_colors = (/"brown","blue"/)
1.1053 +
1.1054 +;**************************************************************************
1.1055 +;(B)-1 histogram plot, ob 81 site
1.1056 +;**************************************************************************
1.1057 +
1.1058 + plot_name = "histogram_ob_81"
1.1059 + title = "Observed 81 site"
1.1060 + resh@tiMainString = title
1.1061 +
1.1062 + wks = gsn_open_wks (plot_type,plot_name)
1.1063 +
1.1064 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.1065 +
1.1066 +;-------------------------------
1.1067 +;Attach the vertical bar and the horizontal cap line
1.1068 +
1.1069 + delete (x1)
1.1070 + delete (x2)
1.1071 + delete (y1)
1.1072 + delete (y2)
1.1073 +
1.1074 + do nd=0,0
1.1075 + lnres@gsLineColor = line_colors(nd)
1.1076 + do i=0,nx-1
1.1077 +
1.1078 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1079 + .not.ismissing(mx_yvalues(nd,i))) then
1.1080 +;
1.1081 +; Attach the vertical bar, both above and below the marker.
1.1082 +;
1.1083 + x1 = xvalues(nd,i)
1.1084 + y1 = yvalues(nd,i)
1.1085 + y2 = mn_yvalues(nd,i)
1.1086 + plx = (/x1,x1/)
1.1087 + ply = (/y1,y2/)
1.1088 + min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1089 +
1.1090 + y2 = mx_yvalues(nd,i)
1.1091 + plx = (/x1,x1/)
1.1092 + ply = (/y1,y2/)
1.1093 + max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1094 +;
1.1095 +; Attach the horizontal cap line, both above and below the marker.
1.1096 +;
1.1097 + x1 = xvalues(nd,i) - dx4
1.1098 + x2 = xvalues(nd,i) + dx4
1.1099 + y1 = mn_yvalues(nd,i)
1.1100 + plx = (/x1,x2/)
1.1101 + ply = (/y1,y1/)
1.1102 + min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1103 +
1.1104 + y1 = mx_yvalues(nd,i)
1.1105 + plx = (/x1,x2/)
1.1106 + ply = (/y1,y1/)
1.1107 + max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1108 + end if
1.1109 + end do
1.1110 + end do
1.1111 +
1.1112 + draw(xy)
1.1113 + frame(wks)
1.1114 + clear (wks)
1.1115 +
1.1116 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1117 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1118 + system("rm "+plot_name+"-*."+plot_type_new)
1.1119 + system("rm "+plot_name+"."+plot_type)
1.1120 +
1.1121 +;****************************************************************************
1.1122 +;(B)-2 histogram plot, model vs ob 81 site
1.1123 +;****************************************************************************
1.1124 +
1.1125 + plot_name = "histogram_model_vs_ob_81"
1.1126 + title = model_name+ " vs Observed 81 site"
1.1127 + resh@tiMainString = title
1.1128 +
1.1129 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1130 +
1.1131 +;-----------------------------
1.1132 +; Add a boxed legend using the more simple method, which won't have
1.1133 +; vertical lines going through the markers.
1.1134 +
1.1135 + resh@pmLegendDisplayMode = "Always"
1.1136 +; resh@pmLegendWidthF = 0.1
1.1137 + resh@pmLegendWidthF = 0.08
1.1138 + resh@pmLegendHeightF = 0.05
1.1139 + resh@pmLegendOrthogonalPosF = -1.17
1.1140 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1141 +; resh@pmLegendParallelPosF = 0.18
1.1142 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.1143 +
1.1144 +; resh@lgPerimOn = False
1.1145 + resh@lgLabelFontHeightF = 0.015
1.1146 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.1147 +;-----------------------------
1.1148 + tRes = True
1.1149 + tRes@txFontHeightF = 0.025
1.1150 +
1.1151 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1.1152 +
1.1153 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1154 +
1.1155 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.1156 +;-------------------------------
1.1157 +;Attach the vertical bar and the horizontal cap line
1.1158 +
1.1159 + do nd=0,1
1.1160 + lnres@gsLineColor = line_colors(nd)
1.1161 + do i=0,nx-1
1.1162 +
1.1163 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1164 + .not.ismissing(mx_yvalues(nd,i))) then
1.1165 +;
1.1166 +; Attach the vertical bar, both above and below the marker.
1.1167 +;
1.1168 + x1 = xvalues(nd,i)
1.1169 + y1 = yvalues(nd,i)
1.1170 + y2 = mn_yvalues(nd,i)
1.1171 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1172 +
1.1173 + y2 = mx_yvalues(nd,i)
1.1174 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1175 +;
1.1176 +; Attach the horizontal cap line, both above and below the marker.
1.1177 +;
1.1178 + x1 = xvalues(nd,i) - dx4
1.1179 + x2 = xvalues(nd,i) + dx4
1.1180 + y1 = mn_yvalues(nd,i)
1.1181 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1182 +
1.1183 + y1 = mx_yvalues(nd,i)
1.1184 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1185 + end if
1.1186 + end do
1.1187 + end do
1.1188 +
1.1189 + draw(xy)
1.1190 + frame(wks)
1.1191 + clear(wks)
1.1192 +
1.1193 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1194 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1195 + system("rm "+plot_name+"-*."+plot_type_new)
1.1196 + system("rm "+plot_name+"."+plot_type)
1.1197 +
1.1198 + delete (RAIN1_1D)
1.1199 + delete (RAIN2_1D)
1.1200 + delete (NPP1_1D)
1.1201 + delete (NPP2_1D)
1.1202 +;delete (range)
1.1203 + delete (xvalues)
1.1204 + delete (yvalues)
1.1205 + delete (mn_yvalues)
1.1206 + delete (mx_yvalues)
1.1207 + delete (good)
1.1208 + delete (max_bar)
1.1209 + delete (min_bar)
1.1210 + delete (max_cap)
1.1211 + delete (min_cap)
1.1212 +
1.1213 +;**************************************************************************
1.1214 +;(B) histogram 933
1.1215 +;**************************************************************************
1.1216 +
1.1217 +; get data
1.1218 + RAIN1_1D = ndtooned(rain933)
1.1219 + RAIN2_1D = ndtooned(rainmod933)
1.1220 + NPP1_1D = ndtooned(npp933)
1.1221 + NPP2_1D = ndtooned(nppmod933)
1.1222 +
1.1223 +; number of bin
1.1224 + nx = 10
1.1225 +
1.1226 + xvalues = new((/2,nx/),float)
1.1227 + yvalues = new((/2,nx/),float)
1.1228 + mn_yvalues = new((/2,nx/),float)
1.1229 + mx_yvalues = new((/2,nx/),float)
1.1230 + dx4 = new((/1/),float)
1.1231 +
1.1232 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.1233 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.1234 +
1.1235 +;----------------------------------------
1.1236 +;compute correlation coeff and M score
1.1237 +
1.1238 + u = yvalues(0,:)
1.1239 + v = yvalues(1,:)
1.1240 +
1.1241 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.1242 + uu = u(good)
1.1243 + vv = v(good)
1.1244 +
1.1245 + ccr933h = esccr(uu,vv,0)
1.1246 +;print (ccr933h)
1.1247 +
1.1248 +;new eq
1.1249 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.1250 + M933h = (1.- (bias/dimsizes(uu)))*5.
1.1251 + print (M933h)
1.1252 +
1.1253 + delete (u)
1.1254 + delete (v)
1.1255 + delete (uu)
1.1256 + delete (vv)
1.1257 +;----------------------------------------------------------------------
1.1258 +; histogram res
1.1259 + delete (resh)
1.1260 + resh = True
1.1261 + resh@gsnMaximize = True
1.1262 + resh@gsnDraw = False
1.1263 + resh@gsnFrame = False
1.1264 + resh@xyMarkLineMode = "Markers"
1.1265 + resh@xyMarkerSizeF = 0.014
1.1266 + resh@xyMarker = 16
1.1267 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.1268 + resh@trYMinF = min(mn_yvalues) - 10.
1.1269 + resh@trYMaxF = max(mx_yvalues) + 10.
1.1270 +
1.1271 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.1272 + resh@tiXAxisString = "Precipitation (m/year)"
1.1273 +
1.1274 + max_bar = new((/2,nx/),graphic)
1.1275 + min_bar = new((/2,nx/),graphic)
1.1276 + max_cap = new((/2,nx/),graphic)
1.1277 + min_cap = new((/2,nx/),graphic)
1.1278 +
1.1279 + lnres = True
1.1280 + line_colors = (/"brown","blue"/)
1.1281 +
1.1282 +;**************************************************************************
1.1283 +;(B)-3 histogram plot, ob 933 site
1.1284 +;**************************************************************************
1.1285 +
1.1286 + plot_name = "histogram_ob_933"
1.1287 + title = "Observed 933 site"
1.1288 + resh@tiMainString = title
1.1289 +
1.1290 + wks = gsn_open_wks (plot_type,plot_name)
1.1291 +
1.1292 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.1293 +
1.1294 +;-------------------------------
1.1295 +;Attach the vertical bar and the horizontal cap line
1.1296 +
1.1297 + do nd=0,0
1.1298 + lnres@gsLineColor = line_colors(nd)
1.1299 + do i=0,nx-1
1.1300 +
1.1301 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1302 + .not.ismissing(mx_yvalues(nd,i))) then
1.1303 +
1.1304 +; Attach the vertical bar, both above and below the marker.
1.1305 +
1.1306 + x1 = xvalues(nd,i)
1.1307 + y1 = yvalues(nd,i)
1.1308 + y2 = mn_yvalues(nd,i)
1.1309 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1310 +
1.1311 + y2 = mx_yvalues(nd,i)
1.1312 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1313 +
1.1314 +; Attach the horizontal cap line, both above and below the marker.
1.1315 +
1.1316 + x1 = xvalues(nd,i) - dx4
1.1317 + x2 = xvalues(nd,i) + dx4
1.1318 + y1 = mn_yvalues(nd,i)
1.1319 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1320 +
1.1321 + y1 = mx_yvalues(nd,i)
1.1322 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1323 + end if
1.1324 + end do
1.1325 + end do
1.1326 +
1.1327 + draw(xy)
1.1328 + frame(wks)
1.1329 + delete (xy)
1.1330 + clear (wks)
1.1331 +
1.1332 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1333 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1334 + system("rm "+plot_name+"-*."+plot_type_new)
1.1335 + system("rm "+plot_name+"."+plot_type)
1.1336 +
1.1337 +;****************************************************************************
1.1338 +;(B)-4 histogram plot, model vs ob 933 site
1.1339 +;****************************************************************************
1.1340 +
1.1341 + plot_name = "histogram_model_vs_ob_933"
1.1342 + title = model_name+ " vs Observed 933 site"
1.1343 + resh@tiMainString = title
1.1344 +
1.1345 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1346 +
1.1347 +;-----------------------------
1.1348 +; Add a boxed legend using the more simple method, which won't have
1.1349 +; vertical lines going through the markers.
1.1350 +
1.1351 + resh@pmLegendDisplayMode = "Always"
1.1352 +; resh@pmLegendWidthF = 0.1
1.1353 + resh@pmLegendWidthF = 0.08
1.1354 + resh@pmLegendHeightF = 0.05
1.1355 + resh@pmLegendOrthogonalPosF = -1.17
1.1356 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1357 +; resh@pmLegendParallelPosF = 0.18
1.1358 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.1359 +
1.1360 +; resh@lgPerimOn = False
1.1361 + resh@lgLabelFontHeightF = 0.015
1.1362 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.1363 +;-----------------------------
1.1364 + tRes = True
1.1365 + tRes@txFontHeightF = 0.025
1.1366 +
1.1367 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1.1368 +
1.1369 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1370 +
1.1371 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.1372 +;-------------------------------
1.1373 +;Attach the vertical bar and the horizontal cap line
1.1374 +
1.1375 + do nd=0,1
1.1376 + lnres@gsLineColor = line_colors(nd)
1.1377 + do i=0,nx-1
1.1378 +
1.1379 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1380 + .not.ismissing(mx_yvalues(nd,i))) then
1.1381 +;
1.1382 +; Attach the vertical bar, both above and below the marker.
1.1383 +;
1.1384 + x1 = xvalues(nd,i)
1.1385 + y1 = yvalues(nd,i)
1.1386 + y2 = mn_yvalues(nd,i)
1.1387 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1388 +
1.1389 + y2 = mx_yvalues(nd,i)
1.1390 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1391 +;
1.1392 +; Attach the horizontal cap line, both above and below the marker.
1.1393 +;
1.1394 + x1 = xvalues(nd,i) - dx4
1.1395 + x2 = xvalues(nd,i) + dx4
1.1396 + y1 = mn_yvalues(nd,i)
1.1397 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1398 +
1.1399 + y1 = mx_yvalues(nd,i)
1.1400 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1401 + end if
1.1402 + end do
1.1403 + end do
1.1404 +
1.1405 + draw(xy)
1.1406 + frame(wks)
1.1407 + delete(xy)
1.1408 + clear(wks)
1.1409 +
1.1410 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1411 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1412 + system("rm "+plot_name+"-*."+plot_type_new)
1.1413 + system("rm "+plot_name+"."+plot_type)
1.1414 +
1.1415 +;***************************************************************************
1.1416 +;(C) global contour
1.1417 +;***************************************************************************
1.1418 +
1.1419 +;res
1.1420 + resg = True ; Use plot options
1.1421 + resg@cnFillOn = True ; Turn on color fill
1.1422 + resg@gsnSpreadColors = True ; use full colormap
1.1423 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1424 +; resg@lbLabelAutoStride = True
1.1425 + resg@cnLinesOn = False ; Turn off contourn lines
1.1426 + resg@mpFillOn = False ; Turn off map fill
1.1427 +
1.1428 + resg@gsnSpreadColors = True ; use full colormap
1.1429 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1430 + resg@cnMinLevelValF = 0. ; Min level
1.1431 + resg@cnMaxLevelValF = 2200. ; Max level
1.1432 + resg@cnLevelSpacingF = 200. ; interval
1.1433 +
1.1434 +;****************************************************************************
1.1435 +;(C)-1 global contour plot, ob
1.1436 +;****************************************************************************
1.1437 +
1.1438 + delta = 0.00000000001
1.1439 + nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1.1440 +
1.1441 + plot_name = "global_ob"
1.1442 + title = "Observed MODIS MOD 17"
1.1443 + resg@tiMainString = title
1.1444 +
1.1445 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1446 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1447 +
1.1448 + plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1449 +
1.1450 + frame(wks)
1.1451 + clear (wks)
1.1452 +
1.1453 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1454 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1455 + system("rm "+plot_name+"-*."+plot_type_new)
1.1456 + system("rm "+plot_name+"."+plot_type)
1.1457 +
1.1458 +;****************************************************************************
1.1459 +;(C)-2 global contour plot, model
1.1460 +;****************************************************************************
1.1461 +
1.1462 + plot_name = "global_model"
1.1463 + title = "Model "+ model_name
1.1464 + resg@tiMainString = title
1.1465 +
1.1466 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1467 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1468 +
1.1469 + plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1470 +
1.1471 + frame(wks)
1.1472 + clear (wks)
1.1473 +
1.1474 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1475 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1476 + system("rm "+plot_name+"-*."+plot_type_new)
1.1477 + system("rm "+plot_name+"."+plot_type)
1.1478 +
1.1479 +;****************************************************************************
1.1480 +;(C)-3 global contour plot, model vs ob
1.1481 +;****************************************************************************
1.1482 +
1.1483 + plot_name = "global_model_vs_ob"
1.1484 +
1.1485 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1486 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1487 +
1.1488 + delete (plot)
1.1489 + plot=new(3,graphic) ; create graphic array
1.1490 +
1.1491 + resg@gsnFrame = False ; Do not draw plot
1.1492 + resg@gsnDraw = False ; Do not advance frame
1.1493 +
1.1494 +; compute correlation coef and M score
1.1495 +
1.1496 + uu1 = ndtooned(nppmod)
1.1497 + vv1 = ndtooned(nppglobe)
1.1498 +
1.1499 + delete (good)
1.1500 + good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1.1501 +
1.1502 + ug = uu1(good)
1.1503 + vg = vv1(good)
1.1504 +
1.1505 + ccrG = esccr(ug,vg,0)
1.1506 +; print (ccrG)
1.1507 +
1.1508 + MG = (ccrG*ccrG)* 5.0
1.1509 + print (MG)
1.1510 +
1.1511 +; plot correlation coef
1.1512 +
1.1513 + gRes = True
1.1514 + gRes@txFontHeightF = 0.02
1.1515 + gRes@txAngleF = 90
1.1516 +
1.1517 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.1518 +
1.1519 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1520 +;--------------------------------------------------------------------
1.1521 +;(a) ob
1.1522 +
1.1523 + title = "Observed MODIS MOD 17"
1.1524 + resg@tiMainString = title
1.1525 +
1.1526 + plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1527 +
1.1528 +;(b) model
1.1529 +
1.1530 + title = "Model "+ model_name
1.1531 + resg@tiMainString = title
1.1532 +
1.1533 + plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1534 +
1.1535 +;(c) model-ob
1.1536 +
1.1537 + zz = nppmod
1.1538 + zz = nppmod - nppglobe
1.1539 + title = "Model_"+model_name+" - Observed"
1.1540 +
1.1541 + resg@cnMinLevelValF = -500 ; Min level
1.1542 + resg@cnMaxLevelValF = 500. ; Max level
1.1543 + resg@cnLevelSpacingF = 50. ; interval
1.1544 + resg@tiMainString = title
1.1545 +
1.1546 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1547 +
1.1548 + pres = True ; panel plot mods desired
1.1549 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1550 + ; indiv. plots in panel
1.1551 + pres@gsnMaximize = True ; fill the page
1.1552 +
1.1553 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1554 +
1.1555 + frame (wks)
1.1556 + clear (wks)
1.1557 + delete (plot)
1.1558 +
1.1559 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1560 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1561 + system("rm "+plot_name+"-*."+plot_type_new)
1.1562 + system("rm "+plot_name+"."+plot_type)
1.1563 +
1.1564 +;***************************************************************************
1.1565 +;(D)-1 zonal line plot, ob
1.1566 +;***************************************************************************
1.1567 +
1.1568 + vv = zonalAve(nppglobe)
1.1569 + vv@long_name = nppglobe@long_name
1.1570 +
1.1571 + plot_name = "zonal_ob"
1.1572 + title = "Observed MODIS MOD 17"
1.1573 +
1.1574 + resz = True
1.1575 + resz@tiMainString = title
1.1576 +
1.1577 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1578 +
1.1579 + plot = gsn_csm_xy (wks,ym,vv,resz)
1.1580 +
1.1581 + frame(wks)
1.1582 + clear (wks)
1.1583 +
1.1584 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1585 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1586 + system("rm "+plot_name+"-*."+plot_type_new)
1.1587 + system("rm "+plot_name+"."+plot_type)
1.1588 +
1.1589 +;****************************************************************************
1.1590 +;(D)-2 zonal line plot, model vs ob
1.1591 +;****************************************************************************
1.1592 +
1.1593 + s = new ((/2,dimsizes(ym)/), float)
1.1594 +
1.1595 + s(0,:) = zonalAve(nppglobe)
1.1596 + s(1,:) = zonalAve(nppmod)
1.1597 +
1.1598 + s@long_name = nppglobe@long_name
1.1599 +;-------------------------------------------
1.1600 +; compute correlation coef and M score
1.1601 +
1.1602 + ccrZ = esccr(s(1,:), s(0,:),0)
1.1603 +; print (ccrZ)
1.1604 +
1.1605 + MZ = (ccrZ*ccrZ)* 5.0
1.1606 + print (MZ)
1.1607 +;-------------------------------------------
1.1608 + plot_name = "zonal_model_vs_ob"
1.1609 + title = "Zonal Average"
1.1610 + resz@tiMainString = title
1.1611 +
1.1612 + wks = gsn_open_wks (plot_type,plot_name)
1.1613 +
1.1614 +; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1.1615 +; resz@vpWidthF = 0.7
1.1616 +
1.1617 + resz@xyMonoLineColor = "False" ; want colored lines
1.1618 + resz@xyLineColors = (/"black","red"/) ; colors chosen
1.1619 +; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1.1620 + resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1.1621 + resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1.1622 +
1.1623 + resz@tiYAxisString = s@long_name ; add a axis title
1.1624 + resz@txFontHeightF = 0.0195 ; change title font heights
1.1625 +
1.1626 +; Legent
1.1627 + resz@pmLegendDisplayMode = "Always" ; turn on legend
1.1628 + resz@pmLegendSide = "Top" ; Change location of
1.1629 +; resz@pmLegendParallelPosF = .45 ; move units right
1.1630 + resz@pmLegendParallelPosF = .82 ; move units right
1.1631 + resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1.1632 +
1.1633 + resz@pmLegendWidthF = 0.10 ; Change width and
1.1634 + resz@pmLegendHeightF = 0.10 ; height of legend.
1.1635 + resz@lgLabelFontHeightF = .02 ; change font height
1.1636 +; resz@lgTitleOn = True ; turn on legend title
1.1637 +; resz@lgTitleString = "Example" ; create legend title
1.1638 +; resz@lgTitleFontHeightF = .025 ; font of legend title
1.1639 + resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
1.1640 +;--------------------------------------------------------------------
1.1641 + zRes = True
1.1642 + zRes@txFontHeightF = 0.025
1.1643 +
1.1644 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1.1645 +
1.1646 + gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1.1647 +;--------------------------------------------------------------------
1.1648 +
1.1649 + plot = gsn_csm_xy (wks,ym,s,resz)
1.1650 +
1.1651 + frame(wks)
1.1652 + clear (wks)
1.1653 +
1.1654 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1655 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1656 + system("rm "+plot_name+"-*."+plot_type_new)
1.1657 + system("rm "+plot_name+"."+plot_type)
1.1658 +
1.1659 +;***************************************************************************
1.1660 +; tar up output plots
1.1661 +;***************************************************************************
1.1662 +
1.1663 + temp_name = "temp." + model_name
1.1664 + system("mkdir -p " + temp_name)
1.1665 + system("mv *.png " + temp_name)
1.1666 + system("tar cf "+ temp_name +".tar " + temp_name)
1.1667 +
1.1668 +;***************************************************************************
1.1669 +end
1.1670 +