diff -r 000000000000 -r 0c6405ab2ff4 co2/24.check.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/co2/24.check.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,143 @@
+; ***********************************************
+; add another model to plot
+; add panel plot to 22.lines.ncl
+; add zone plot to 21.lines.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.061n_1995-2004_MONS_climo_atm.nc"
+  fili3 = "b30.061m_401_425_MONS_climo_atm.nc"
+;--------------------------------------------
+  g    = addfile(diri2+fili2,"r")
+  x1   = g->CO2
+  xi   = g->lon
+  yi   = g->lat
+  delete (g)
+
+  xdim  = dimsizes(x1)
+  nlev  = xdim(1)
+; y1     = x1(:,0,:,:)
+  y1     = x1
+; printVarSummary (y1)
+  
+; get the co2 at the lowest level
+; y1     = x1(:,nlev-1,:,:)
+  delete (x1)
+;---------------------------------------------
+; g     = addfile(diri2+fili3,"r")
+; x2    = g->CO2
+; delete (g)
+; y2     = x2(:,0,:,:)
+; y2     = x2(:,nlev-1,:,:)
+; delete (x2)
+;---------------------------------------------
+; 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
+;---------------------------------------------
+  y1 = y1 * factor
+  y1@_FillValue = 1.e36
+  y1@units      = "u mol/mol"
+; y1 = where(y0 .lt. 287.,y1@_FillValue,y1)
+; printVarSummary (y1)
+; print (min(y1)+"/"+max(y1))
+;---------------------------------------------
+; y2 = y2 * factor
+; y2@_FillValue = 1.e36
+; y2@units      = "u mol/mol"
+;---------------------------------------------
+; interpolate into observed station
+; note: model is 0-360E,   90S-90N
+;       ob    is -180-180, 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,y1,True,lon_ob,lat_ob,0)
+  printVarSummary (yo)
+; yo:[time | 12] x [lev | 26] x [pts | 98]
+
+  val_model1 = yo(pts|:,lev|:,time|:)
+  delete (yo)
+  val_model1_0 = val_model1
+; printVarSummary (val_model1)
+; print (min(val_model1)+"/"+max(val_model1))
+
+; remove annual mean
+  val_model1 = val_model1 - conform(val_model1,dim_avg(val_model1),(/0,1/))
+; print (min(val_model1)+"/"+max(val_model1))
+
+;-----------------------------------------------------------------
+;    index of station Barrow, Alaska (71.32,-156.60)
+     ind_z = ind(lat_ob .eq. 71.32)
+     print (ind_z)
+     print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
+     print ("observation at Barrow, Alaska (71.32,-156.60)")
+     print (val_ob(ind_z,:))
+     print ("model top atm  at Barrow, Alaska (71.32,-156.60)")
+     print (val_model1_0(ind_z,0,:))
+     print ("model surface at Barrow, Alaska (71.32,-156.60)")
+     print (val_model1_0(ind_z,nlev-1,:))
+     print ("model top atm  at Barrow, Alaska (71.32,-156.60)")
+     print (val_model1(ind_z,0,:))
+     print ("model surface at Barrow, Alaska (71.32,-156.60)")
+     print (val_model1(ind_z,nlev-1,:))
+end