forrest@0: ; ***********************************************
forrest@0: ; xy_4.ncl
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: begin
forrest@0: ;************************************************
forrest@0: ; read in data: observed
forrest@0: ;************************************************
forrest@0:  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
forrest@0:  fili  = "co2_globalView_222.nc"
forrest@0:  g     = addfile (diri+fili,"r")
forrest@0:  val   = g->CO2_SEAS  
forrest@0:  lon   = g->LON 
forrest@0:  lat   = g->LAT 
forrest@0:  delete (g)
forrest@0:  
forrest@0:  ncase = dimsizes(lat)
forrest@0: ;print (ncase)
forrest@0: 
forrest@0: ;**************************************************************
forrest@0: ; get only the lowest level at each station 
forrest@0: ;**************************************************************
forrest@0:  lat_tmp = lat
forrest@0:  lat_tmp@_FillValue = 1.e+36
forrest@0:  
forrest@0:  do n = 0,ncase-1
forrest@0:     if (.not. ismissing(lat_tmp(n))) then 
forrest@0:        indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
forrest@0:        if (dimsizes(indexes) .gt. 1) then
forrest@0:           lat_tmp(indexes(1:)) = lat_tmp@_FillValue
forrest@0:        end if
forrest@0:        delete (indexes)
forrest@0:     end if
forrest@0:  end do
forrest@0: 
forrest@0:  indexes = ind(.not. ismissing(lat_tmp))
forrest@0: ;print (indexes)
forrest@0:  
forrest@0:  lat_ob = lat(indexes)
forrest@0:  lon_ob = lon(indexes)
forrest@0:  val_ob = val(indexes,:)
forrest@0: ;printVarSummary (val_ob)
forrest@0: ;print (lat_ob +"/"+lon_ob)
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read in model data
forrest@0: ;************************************************
forrest@0:   diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:   fili2 = "b30.061m_401_425_MONS_climo.nc"
forrest@0: 
forrest@0:   g     = addfile(diri2+fili2,"r")
forrest@0:   x     = g->CO2
forrest@0:   xi    = g->lon
forrest@0:   yi    = g->lat
forrest@0:   xdim  = dimsizes(x)
forrest@0:   nlev  = xdim(1)
forrest@0:   y     = x(:,0,:,:)
forrest@0: ; printVarSummary (y)
forrest@0:   
forrest@0: ; get the co2 at the lowest level
forrest@0:   y     = x(:,nlev-1,:,:)
forrest@0: 
forrest@0: ; change to unit of observed (u mol/mol)
forrest@0: ; Model_units [=] kgCO2 / kgDryAir
forrest@0: ; 28.966 = molecular weight of dry air
forrest@0: ; 44.       = molecular weight of CO2
forrest@0: ; u mol = 1e-6 mol
forrest@0: 
forrest@0:   factor = (28.966/44.) * 1e6
forrest@0:   y      = y * factor
forrest@0: 
forrest@0:   y@_FillValue = 1.e36
forrest@0:   y@units      = "u mol/mol"
forrest@0: ; y = where(y0 .lt. 287.,y@_FillValue,y)
forrest@0: ; printVarSummary (y)
forrest@0: ; print (min(y)+"/"+max(y))
forrest@0: 
forrest@0: ; interpolate into observed station
forrest@0: 
forrest@0:   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
forrest@0: 
forrest@0:   val_model = val_ob
forrest@0:   val_model = yo(pts|:,time|:)
forrest@0:   printVarSummary (val_model)
forrest@0:   print (min(val_model)+"/"+max(val_model))
forrest@0: 
forrest@0: ; remove annual mean
forrest@0:   val_model = val_model - conform(val_model,dim_avg(val_model),0)
forrest@0:   print (min(val_model)+"/"+max(val_model))
forrest@0: exit
forrest@0: ;**************************************************************
forrest@0: ; get stations of 60N-90N
forrest@0: ;**************************************************************
forrest@0: 
forrest@0:  ind_1 = ind(lat_new .gt. 60.)
forrest@0:  print (ind_1)
forrest@0:  print (lat_new(ind_1)+"/"+lon_new(ind_1))
forrest@0:  print (val_new(ind_1(0),:))
forrest@0: exit
forrest@0: 
forrest@0: end