1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/co2/12.test.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,104 @@
     1.4 +; ***********************************************
     1.5 +; xy_4.ncl
     1.6 +; ***********************************************
     1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
     1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
     1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    1.10 +;************************************************
    1.11 +begin
    1.12 +;************************************************
    1.13 +; read in data: observed
    1.14 +;************************************************
    1.15 + diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
    1.16 + fili  = "co2_globalView_222.nc"
    1.17 + g     = addfile (diri+fili,"r")
    1.18 + val   = g->CO2_SEAS  
    1.19 + lon   = g->LON 
    1.20 + lat   = g->LAT 
    1.21 + delete (g)
    1.22 + 
    1.23 + ncase = dimsizes(lat)
    1.24 +;print (ncase)
    1.25 +
    1.26 +;**************************************************************
    1.27 +; get only the lowest level at each station 
    1.28 +;**************************************************************
    1.29 + lat_tmp = lat
    1.30 + lat_tmp@_FillValue = 1.e+36
    1.31 + 
    1.32 + do n = 0,ncase-1
    1.33 +    if (.not. ismissing(lat_tmp(n))) then 
    1.34 +       indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
    1.35 +       if (dimsizes(indexes) .gt. 1) then
    1.36 +          lat_tmp(indexes(1:)) = lat_tmp@_FillValue
    1.37 +       end if
    1.38 +       delete (indexes)
    1.39 +    end if
    1.40 + end do
    1.41 +
    1.42 + indexes = ind(.not. ismissing(lat_tmp))
    1.43 +;print (indexes)
    1.44 + 
    1.45 + lat_ob = lat(indexes)
    1.46 + lon_ob = lon(indexes)
    1.47 + val_ob = val(indexes,:)
    1.48 +;printVarSummary (val_ob)
    1.49 +;print (lat_ob +"/"+lon_ob)
    1.50 +
    1.51 +;************************************************
    1.52 +; read in model data
    1.53 +;************************************************
    1.54 +  diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
    1.55 +  fili2 = "b30.061m_401_425_MONS_climo.nc"
    1.56 +
    1.57 +  g     = addfile(diri2+fili2,"r")
    1.58 +  x     = g->CO2
    1.59 +  xi    = g->lon
    1.60 +  yi    = g->lat
    1.61 +  xdim  = dimsizes(x)
    1.62 +  nlev  = xdim(1)
    1.63 +  y     = x(:,0,:,:)
    1.64 +; printVarSummary (y)
    1.65 +  
    1.66 +; get the co2 at the lowest level
    1.67 +  y     = x(:,nlev-1,:,:)
    1.68 +
    1.69 +; change to unit of observed (u mol/mol)
    1.70 +; Model_units [=] kgCO2 / kgDryAir
    1.71 +; 28.966 = molecular weight of dry air
    1.72 +; 44.       = molecular weight of CO2
    1.73 +; u mol = 1e-6 mol
    1.74 +
    1.75 +  factor = (28.966/44.) * 1e6
    1.76 +  y      = y * factor
    1.77 +
    1.78 +  y@_FillValue = 1.e36
    1.79 +  y@units      = "u mol/mol"
    1.80 +; y = where(y0 .lt. 287.,y@_FillValue,y)
    1.81 +; printVarSummary (y)
    1.82 +; print (min(y)+"/"+max(y))
    1.83 +
    1.84 +; interpolate into observed station
    1.85 +
    1.86 +  yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
    1.87 +
    1.88 +  val_model = val_ob
    1.89 +  val_model = yo(pts|:,time|:)
    1.90 +  printVarSummary (val_model)
    1.91 +  print (min(val_model)+"/"+max(val_model))
    1.92 +
    1.93 +; remove annual mean
    1.94 +  val_model = val_model - conform(val_model,dim_avg(val_model),0)
    1.95 +  print (min(val_model)+"/"+max(val_model))
    1.96 +exit
    1.97 +;**************************************************************
    1.98 +; get stations of 60N-90N
    1.99 +;**************************************************************
   1.100 +
   1.101 + ind_1 = ind(lat_new .gt. 60.)
   1.102 + print (ind_1)
   1.103 + print (lat_new(ind_1)+"/"+lon_new(ind_1))
   1.104 + print (val_new(ind_1(0),:))
   1.105 +exit
   1.106 +
   1.107 +end