diff -r 000000000000 -r 0c6405ab2ff4 co2/14.metric.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/co2/14.metric.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,197 @@
+; ***********************************************
+; 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"
+
+  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
+
+  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 = val_ob
+  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
+
+ u = ndtooned(val_ob(ind_z,:))
+ v = ndtooned(val_model(ind_z,:))
+
+ good = ind(.not.ismissing(u) .and. .not.ismissing(v))
+
+ uu = u(good)
+ vv = v(good)
+ st = sta(ind_z)
+
+ npts = dimsizes(st)
+ print (npts)
+
+ ccr = esccr(uu,vv,0)
+;print (ccr)
+
+ un =  onedtond(uu,(/npts,12/))
+ vn =  onedtond(vv,(/npts,12/))
+;print (un)
+;print (vn)
+
+ score = new((/npts/),float)
+
+ do n=0,npts-1
+    amp_ob    = max(un(n,:)) - min(un(n,:)) 
+    amp_model = max(vn(n,:)) - min(vn(n,:))
+    score(n) = 1.-abs((amp_model/amp_ob)-1.)
+;   print (amp_ob)
+;   print (amp_model)  
+ end do
+;print (score)
+
+ M = avg(score)
+ print (M)
+
+ M_total = (ccr*ccr + M)*0.5 * score_max
+ print (M_total)
+
+ delete (ind_z)
+ delete (good)
+ delete (u)
+ delete (v)
+ delete (uu)
+ delete (vv)
+ delete (un)
+ delete (vn)
+ delete (st) 
+ delete (score)
+ end do
+end