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:  
forrest@0: ;nclass = 18
forrest@0:  nclass = 20
forrest@0:  
forrest@0: ;************************************************
forrest@0: ; read in data: observed
forrest@0: ;************************************************
forrest@0:  diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0: ;fili1  = "land_class_T42.nc"
forrest@0:  fili1  = "land_class_T42_new.nc"
forrest@0:  fili2  = "LAI_2000-2005_mean_T42.nc"
forrest@0:  data_file_ob1 = addfile(diri1+fili1,"r")
forrest@0:  data_file_ob2 = addfile(diri1+fili2,"r")
forrest@0:  
forrest@0:  RAIN1 = tofloat(data_file_ob1->LAND_CLASS)               
forrest@0:  NPP1  = data_file_ob2->LAI      
forrest@0: ;************************************************
forrest@0: ; read in data: model       
forrest@0: ;************************************************
forrest@0:  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0: ;fili3  = "i01.03cn_1545-1569_ANN_climo.nc"
forrest@0:  fili3  = "i01.04casa_1605-1629_ANN_climo.nc"
forrest@0:  data_file_model = addfile(diri2+fili3,"r")
forrest@0:       
forrest@0:  NPP2  = data_file_model->TLAI      
forrest@0: ;************************************************
forrest@0: ; print min/max and unit
forrest@0: ;************************************************
forrest@0:   pminmax(RAIN1,"RAIN1")
forrest@0:   pminmax(NPP1,"NPP1")
forrest@0:   pminmax(NPP2,"NPP2")
forrest@0: 
forrest@0:   RAIN1_1D = ndtooned(RAIN1)
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: 
forrest@0: ; nbins       = nclass + 1         ; Number of bins to use.
forrest@0: ; nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False)
forrest@0: ; nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
forrest@0: ; range       = fspan(nicevals(0),nicevals(1),nvals)
forrest@0: 
forrest@0:   nclassn     = nclass + 1
forrest@0:   range       = fspan(0,nclassn-1,nclassn)
forrest@0: 
forrest@0: ; print (nicevals)
forrest@0: ; print (nvals)
forrest@0:   print (range)
forrest@0: ; exit
forrest@0: 
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(RAIN1_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(RAIN1_1D))
forrest@0:   mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))
forrest@0:   mx_yvalues   = new((/2,nx/),typeof(RAIN1_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     = RAIN1_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(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: 
forrest@0:  u = yvalues(0,:)
forrest@0:  v = yvalues(1,:)
forrest@0:  print (u)
forrest@0:  print (v)
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: ;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: