diff -r 000000000000 -r 0c6405ab2ff4 co2/15.metric.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/co2/15.metric.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,191 @@
+; ***********************************************
+; xy_4.ncl
+; ***********************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+;************************************************
+begin
+;************************************************
+; read in data: observed
+;************************************************
+ diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
+ fili  = "co2_globalView_98.nc"
+ g     = addfile (diri+fili,"r")
+ val   = g->CO2_SEAS  
+ lon   = g->LON 
+ lat   = g->LAT
+ sta   = chartostring(g->STATION) 
+ delete (g)
+ 
+;print (sta(0))
+
+ ncase = dimsizes(lat)
+;print (ncase)
+
+;**************************************************************
+; get only the lowest level at each station 
+;**************************************************************
+ lat_tmp = lat
+ lat_tmp@_FillValue = 1.e+36
+ 
+ do n = 0,ncase-1
+    if (.not. ismissing(lat_tmp(n))) then 
+       indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
+       if (dimsizes(indexes) .gt. 1) then
+          lat_tmp(indexes(1:)) = lat_tmp@_FillValue
+       end if
+       delete (indexes)
+    end if
+ end do
+
+ indexes = ind(.not. ismissing(lat_tmp))
+;print (dimsizes(indexes))
+;print (indexes)
+ 
+ lat_ob = lat(indexes)
+ lon_ob = lon(indexes)
+ val_ob = val(indexes,:)
+;printVarSummary (val_ob)
+;print (lat_ob +"/"+lon_ob)
+
+;************************************************
+; read in model data
+;************************************************
+  diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+; fili2 = "b30.061m_401_425_MONS_climo.nc"
+  fili2 = "b30.061n_1995-2004_MONS_climo.nc"
+
+  g     = addfile(diri2+fili2,"r")
+  x     = g->CO2
+  xi    = g->lon
+  yi    = g->lat
+  xdim  = dimsizes(x)
+  nlev  = xdim(1)
+  y     = x(:,0,:,:)
+; printVarSummary (y)
+  
+; get the co2 at the lowest level
+  y     = x(:,nlev-1,:,:)
+
+; change to unit of observed (u mol/mol)
+; Model_units [=] kgCO2 / kgDryAir
+; 28.966 = molecular weight of dry air
+; 44.       = molecular weight of CO2
+; u mol = 1e-6 mol
+
+  factor = (28.966/44.) * 1e6
+  y      = y * factor
+
+  y@_FillValue = 1.e36
+  y@units      = "u mol/mol"
+; y = where(y0 .lt. 287.,y@_FillValue,y)
+; printVarSummary (y)
+; print (min(y)+"/"+max(y))
+
+; interpolate into observed station
+; note: model is 0-360E, 90S-90N
+
+; to be able to handle observation at (-89.98,-24.80)
+  print (yi(0))
+  yi(0) = -90.
+
+  i = ind(lon_ob .lt. 0.)
+  lon_ob(i) = lon_ob(i) + 360.  
+
+  yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
+
+  val_model = yo(pts|:,time|:)
+; printVarSummary (val_model)
+; print (min(val_model)+"/"+max(val_model))
+
+; remove annual mean
+  val_model = val_model - conform(val_model,dim_avg(val_model),0)
+; print (min(val_model)+"/"+max(val_model))
+
+  nzone = 4
+  do z = 0,nzone-1
+
+  if (z .eq. 0) then 
+;    maximum score for the zone, 60N-90N 
+     score_max = 5.0
+;    index of stations in this zone
+     ind_z = ind(lat_ob .ge. 60.)
+;    print (ind_z)
+;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
+;    print (val_ob(ind_z,:))
+;    print (val_model(ind_z,:))
+  end if
+
+  if (z .eq. 1) then 
+;    maximum score for the zone, 30N-60N 
+     score_max = 5.0
+;    index of stations in this zone
+     ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
+;    print (ind_z)
+;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
+;    print (val_ob(ind_z,:))
+;    print (val_model(ind_z,:))
+  end if
+
+  if (z .eq. 2) then 
+;    maximum score for the zone, EQ-30N 
+     score_max = 5.0
+;    index of stations in this zone
+     ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
+;    print (ind_z)
+;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
+;    print (val_ob(ind_z,:))
+;    print (val_model(ind_z,:))
+  end if
+
+  if (z .eq. 3) then 
+;    maximum score for the zone, 90S-EQ 
+     score_max = 5.0
+;    index of stations in this zone
+     ind_z = ind(lat_ob .lt. 0. )
+;    print (ind_z)
+;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
+;    print (val_ob(ind_z,:))
+;    print (val_model(ind_z,:))
+  end if
+
+ npts = dimsizes(ind_z)
+ print (npts)
+
+ amp_ob        = new((/npts/),float)
+ amp_model     = new((/npts/),float)
+
+ amp_ratio_sta = new((/npts/),float)
+ ccr_sta       = new((/npts/),float)
+ M_sta         = new((/npts/),float)
+ score_sta     = new((/npts/),float)
+
+ do n=0,npts-1
+    amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 
+    amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
+
+    amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
+    ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
+    M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
+    score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
+
+    print (sta(ind_z(n))+"/"+lat(ind_z(n))+"/"+lon(ind_z(n))+"/"+amp_ratio_sta(n)+"/"+ccr_sta(n)+"/"+M_sta(n)+"/"+score_sta(n))  
+ end do
+
+ amp_ratio_zone = avg(amp_ratio_sta)
+ ccr_zone       = avg(ccr_sta)
+ M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 
+ score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
+
+ print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)  
+
+ delete (ind_z)
+ delete (amp_model)
+ delete (amp_ob)
+ delete (amp_ratio_sta)
+ delete (ccr_sta)
+ delete (M_sta)
+ delete (score_sta)
+ end do
+end