forrest@0: ;********************************************************
forrest@0: ; histogram normalized by rain and compute correleration
forrest@0: ;********************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: 
forrest@0: procedure pminmax(data:numeric,name:string)
forrest@0: begin
forrest@0:   print ("min/max " + name + " = " + min(data) + "/" + max(data))
forrest@0:   if(isatt(data,"units")) then
forrest@0:     print (name + " units = " + data@units)
forrest@0:   end if
forrest@0: end
forrest@0: 
forrest@0: ;
forrest@0: ; Main code.
forrest@0: ;
forrest@0: begin
forrest@0:   data_types    = (/ "Obs",        "Model"       /)
forrest@0:   data_names    = (/ "data.81.nc" , "i01.03cn_1545-1569_ANN_climo.nc" /)
forrest@0: ; data_names    = (/ "data.933.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)
forrest@0: ; data_names    = (/ "data.81.nc" , "i01.04casa_1605-1629_ANN_climo.nc" /)
forrest@0: ; data_names    = (/ "data.933.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)
forrest@0:   filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)
forrest@0:   ndata_types   = dimsizes(data_types)
forrest@0:       
forrest@0:   data_file_obs = addfile(data_names(0),"r")    ; Open obs file
forrest@0:   data_file_mod = addfile(data_names(1),"r")    ; Open model file
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read in data: observed
forrest@0: ;************************************************
forrest@0:  RAIN1 = tofloat(data_file_obs->PREC_ANN) ; for data.81
forrest@0: ;RAIN1 = data_file_obs->PREC              ; for data.933
forrest@0:  NPP1  = data_file_obs->TNPP_C
forrest@0:  xo    = data_file_obs->LONG_DD  
forrest@0:  yo    = data_file_obs->LAT_DD   
forrest@0: 
forrest@0: ; change longitude from (-180,180) to (0,360)
forrest@0:  print (xo)
forrest@0:  nx = dimsizes(xo)
forrest@0:  do i= 0,nx-1
forrest@0:     if (xo(i) .lt. 0.) then
forrest@0:         xo(i) = xo(i)+ 360.
forrest@0:     end if
forrest@0:  end do
forrest@0:  print (xo)
forrest@0:  
forrest@0: ;************************************************
forrest@0: ; read in data: model       
forrest@0: ;************************************************
forrest@0:  ai    = data_file_mod->RAIN
forrest@0:  bi    = data_file_mod->NPP
forrest@0:  xi    = data_file_mod->lon      
forrest@0:  yi    = data_file_mod->lat      
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; interpolate from model grid to observed grid
forrest@0: ;************************************************
forrest@0:  RAIN2 = linint2_points(xi,yi,ai,True,xo,yo,0) 
forrest@0:  NPP2  = linint2_points(xi,yi,bi,True,xo,yo,0) 
forrest@0:  
forrest@0: ;************************************************
forrest@0: ; convert unit
forrest@0: ;************************************************
forrest@0: ; Units for these four variables are:
forrest@0: ;
forrest@0: ; RAIN1 : mm/year
forrest@0: ; RAIN2 : mm/s
forrest@0: ; NPP1  : g C/m^2/year
forrest@0: ; NPP2  : g C/m^2/s
forrest@0: ;
forrest@0: ; We want to convert these to "m/year" and "g C/m^2/year".
forrest@0: ;
forrest@0:   nsec_per_year = 60*60*24*365                  ; # seconds per year
forrest@0: 
forrest@0: ; Do the necessary conversions.
forrest@0:   RAIN1 = RAIN1 / 1000.
forrest@0:   RAIN2 = (RAIN2/ 1000.) * nsec_per_year
forrest@0:   NPP2     = NPP2 * nsec_per_year
forrest@0: 
forrest@0: ; Redo the units.
forrest@0:   RAIN1@units = "m/yr"
forrest@0:   RAIN2@units = "m/yr"
forrest@0:   NPP1@units  = "gC/m^2/yr"
forrest@0:   NPP2@units  = "gC/m^2/yr"
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; print min/max and unit
forrest@0: ;************************************************
forrest@0:   pminmax(RAIN1,"RAIN1")
forrest@0:   pminmax(RAIN2,"RAIN2")
forrest@0:   pminmax(NPP1,"NPP1")
forrest@0:   pminmax(NPP2,"NPP2")
forrest@0: 
forrest@0:   RAIN1_1D = ndtooned(RAIN1)
forrest@0:   RAIN2_1D = ndtooned(RAIN2)
forrest@0:   NPP1_1D  = ndtooned(NPP1)
forrest@0:   NPP2_1D  = ndtooned(NPP2)
forrest@0: ;
forrest@0: ; Calculate some "nice" bins for binning the data in equally spaced
forrest@0: ; ranges.
forrest@0: ;
forrest@0:   nbins       = 15     ; Number of bins to use.
forrest@0: 
forrest@0:   nicevals    = nice_mnmxintvl(min(RAIN2_1D),max(RAIN2_1D),nbins,True)
forrest@0:   nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
forrest@0:   range       = fspan(nicevals(0),nicevals(1),nvals)
forrest@0: ;
forrest@0: ; Use this range information to grab all the values in a
forrest@0: ; particular range, and then take an average.
forrest@0: ;
forrest@0:   nr      = dimsizes(range)
forrest@0:   nx      = nr-1
forrest@0:   xvalues     = new((/2,nx/),typeof(RAIN2_1D))
forrest@0:   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
forrest@0:   dx           = xvalues(0,1) - xvalues(0,0)       ; range width
forrest@0:   dx4          = dx/4                              ; 1/4 of the range
forrest@0:   xvalues(1,:) = xvalues(0,:) - dx/5.
forrest@0:   yvalues      = new((/2,nx/),typeof(RAIN2_1D))
forrest@0:   mn_yvalues   = new((/2,nx/),typeof(RAIN2_1D))
forrest@0:   mx_yvalues   = new((/2,nx/),typeof(RAIN2_1D))
forrest@0: 
forrest@0:   do nd=0,1
forrest@0: ;
forrest@0: ; See if we are doing model or observational data.
forrest@0: ;
forrest@0:     if(nd.eq.0) then
forrest@0:       data     = RAIN1_1D
forrest@0:       npp_data = NPP1_1D
forrest@0:     else
forrest@0:       data     = RAIN2_1D
forrest@0:       npp_data = NPP2_1D
forrest@0:     end if
forrest@0: ;
forrest@0: ; Loop through each range and check for values.
forrest@0: ;
forrest@0:     do i=0,nr-2
forrest@0:       if (i.ne.(nr-2)) then
forrest@0: ;        print("")
forrest@0: ;        print("In range ["+range(i)+","+range(i+1)+")")
forrest@0:         idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
forrest@0:       else
forrest@0: ;        print("")
forrest@0: ;        print("In range ["+range(i)+",)")
forrest@0:         idx = ind(range(i).le.data)
forrest@0:       end if
forrest@0: ;
forrest@0: ; Calculate average, and get min and max.
forrest@0: ;
forrest@0:       if(.not.any(ismissing(idx))) then
forrest@0:         yvalues(nd,i)    = avg(npp_data(idx))
forrest@0:         mn_yvalues(nd,i) = min(npp_data(idx))
forrest@0:         mx_yvalues(nd,i) = max(npp_data(idx))
forrest@0:         count = dimsizes(idx)
forrest@0:       else
forrest@0:         count            = 0
forrest@0:         yvalues(nd,i)    = yvalues@_FillValue
forrest@0:         mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0:         mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0:       end if
forrest@0: ;
forrest@0: ; Print out information.
forrest@0: ;
forrest@0: ;      print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0: ;      print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0: 
forrest@0: ;
forrest@0: ; Clean up for next time in loop.
forrest@0: ;
forrest@0:       delete(idx)
forrest@0:     end do
forrest@0:     delete(data)
forrest@0:     delete(npp_data)
forrest@0:   end do
forrest@0: 
forrest@0: ;
forrest@0: ; Start the graphics.
forrest@0: ;
forrest@0:   wks = gsn_open_wks("png","xy")
forrest@0: 
forrest@0:   res             = True
forrest@0:   res@gsnMaximize = True
forrest@0:   res@gsnDraw     = False
forrest@0:   res@gsnFrame    = False
forrest@0:   res@xyMarkLineMode = "Markers"
forrest@0:   res@xyMarkerSizeF   = 0.014
forrest@0:   res@xyMarker       = 16
forrest@0:   res@xyMarkerColors = (/"Brown","Blue"/)
forrest@0:   res@trYMinF        = min(mn_yvalues) - 10.
forrest@0:   res@trYMaxF        = max(mx_yvalues) + 10.
forrest@0: 
forrest@0: ; res@tiMainString  = "Observed vs i01.03cn 81 site"
forrest@0: ; res@tiMainString  = "Observed vs i01.03cn 933 site"
forrest@0:   res@tiMainString  = "Observed vs i01.04casa 81 site"
forrest@0: ; res@tiMainString  = "Observed vs i01.04casa 933 site"
forrest@0:   res@tiYAxisString = "NPP (g C/m2/year)"
forrest@0:   res@tiXAxisString = "Precipitation (m/year)"
forrest@0: 
forrest@0: ;
forrest@0: ; Add a boxed legend using the more simple method, which won't have
forrest@0: ; vertical lines going through the markers.
forrest@0: ;
forrest@0:   res@pmLegendDisplayMode    = "Always"
forrest@0: ; res@pmLegendWidthF         = 0.1
forrest@0:   res@pmLegendWidthF         = 0.08
forrest@0:   res@pmLegendHeightF        = 0.05
forrest@0:   res@pmLegendOrthogonalPosF = -1.17
forrest@0: ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0: ; res@pmLegendParallelPosF   =  0.18
forrest@0:   res@pmLegendParallelPosF   =  0.23  ;(rightward)
forrest@0: 
forrest@0: ;  res@lgPerimOn             = False
forrest@0:   res@lgLabelFontHeightF     = 0.015
forrest@0: ; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)
forrest@0:   res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)
forrest@0: 
forrest@0:   xy = gsn_csm_xy(wks,xvalues,yvalues,res)
forrest@0: 
forrest@0:   max_bar = new((/2,nx/),graphic)
forrest@0:   min_bar = new((/2,nx/),graphic)
forrest@0:   max_cap = new((/2,nx/),graphic)
forrest@0:   min_cap = new((/2,nx/),graphic)
forrest@0: 
forrest@0:   lnres = True
forrest@0: 
forrest@0:   line_colors = (/"brown","blue"/)
forrest@0:   do nd=0,1
forrest@0:     lnres@gsLineColor = line_colors(nd)
forrest@0:     do i=0,nx-1
forrest@0:      
forrest@0:       if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0:          .not.ismissing(mx_yvalues(nd,i))) then
forrest@0: ;
forrest@0: ; Attach the vertical bar, both above and below the marker.
forrest@0: ;
forrest@0:         x1 = xvalues(nd,i)
forrest@0:         y1 = yvalues(nd,i)
forrest@0:         y2 = mn_yvalues(nd,i)
forrest@0:         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0: 
forrest@0:         y2 = mx_yvalues(nd,i)
forrest@0:         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0: ;
forrest@0: ; Attach the horizontal cap line, both above and below the marker.
forrest@0: ;
forrest@0:         x1 = xvalues(nd,i) - dx4
forrest@0:         x2 = xvalues(nd,i) + dx4
forrest@0:         y1 = mn_yvalues(nd,i)
forrest@0:         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0: 
forrest@0:         y1 = mx_yvalues(nd,i)
forrest@0:         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0:       end if
forrest@0:     end do
forrest@0:   end do
forrest@0: 
forrest@0: ;
forrest@0: ; Here's how to do the legend by hand.
forrest@0: ;
forrest@0: ; mkres               = True    ; Marker resources
forrest@0: ; txres               = True    ; Text resources
forrest@0: ; mkres@gsMarkerIndex = 16
forrest@0: ; mkres@gsMarkerSizeF = 0.02
forrest@0: ; txres@txFontHeightF = 0.02
forrest@0: ; txres@txJust        = "CenterLeft"
forrest@0: ;
forrest@0: ; Change these values if you want to move the marker legend location.
forrest@0: ; These values are in the same data space as the plot. 
forrest@0: ;
forrest@0: ; xlg1_cen = 0.2
forrest@0: ; ylg1_cen = 900.
forrest@0: 
forrest@0: ; xlg2_cen = 0.2
forrest@0: ; ylg2_cen = 760.
forrest@0: 
forrest@0: ; mkres@gsMarkerColor = "brown"
forrest@0: ; lnres@gsLineColor   = "brown"
forrest@0: 
forrest@0: ; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)
forrest@0: ; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \
forrest@0: ;                                           (/ylg1_cen-60,ylg1_cen+60/),lnres)
forrest@0: ; lg_cap_legend11  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
forrest@0: ;                                           (/ylg1_cen-60,ylg1_cen-60/),lnres)
forrest@0: ; lg_cap_legend12  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
forrest@0: ;                                           (/ylg1_cen+60,ylg1_cen+60/),lnres)
forrest@0: 
forrest@0: ; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)
forrest@0: 
forrest@0: ; mkres@gsMarkerColor = "blue"
forrest@0: ; lnres@gsLineColor   = "blue"
forrest@0: 
forrest@0: ; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)
forrest@0: ; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \
forrest@0: ;                                           (/ylg2_cen-60,ylg2_cen+60/),lnres)
forrest@0: ; lg_cap_legend21  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
forrest@0: ;                                           (/ylg2_cen-60,ylg2_cen-60/),lnres)
forrest@0: ; lg_cap_legend22  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
forrest@0: ;                                           (/ylg2_cen+60,ylg2_cen+60/),lnres)
forrest@0: ; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)
forrest@0: 
forrest@0:   draw(xy)
forrest@0:   frame(wks)
forrest@0: 
forrest@0:  u = yvalues(0,:)
forrest@0:  v = yvalues(1,:)
forrest@0: 
forrest@0:  good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0:  uu = u(good)
forrest@0:  vv = v(good)
forrest@0:  nz = dimsizes(uu)
forrest@0:  print (nz)
forrest@0: 
forrest@0:  ccr = esccr(uu,vv,0)
forrest@0:  print (ccr)
forrest@0: 
forrest@0: ;old eq
forrest@0: ;bias = sum(((vv-uu)/uu)^2)
forrest@0: ;M    = (1.- sqrt(bias/nz))*5.
forrest@0: 
forrest@0: ;new eq
forrest@0:  bias = sum(abs(vv-uu)/(vv+uu))
forrest@0:  M    = (1.- (bias/nz))*5.
forrest@0:  print (bias)
forrest@0:  print (M)
forrest@0: 
forrest@0: end
forrest@0: