forrest@0: ; ***********************************************
forrest@0: ; add another model to plot
forrest@0: ; add panel plot to 22.lines.ncl
forrest@0: ; add zone plot to 21.lines.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_98.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:  sta   = chartostring(g->STATION) 
forrest@0:  delete (g)
forrest@0:  
forrest@0: ;print (sta(0))
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 (dimsizes(indexes))
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.061n_1995-2004_MONS_climo_atm.nc"
forrest@0:   fili3 = "b30.061m_401_425_MONS_climo_atm.nc"
forrest@0: ;--------------------------------------------
forrest@0:   g    = addfile(diri2+fili2,"r")
forrest@0:   x1   = g->CO2
forrest@0:   xi   = g->lon
forrest@0:   yi   = g->lat
forrest@0:   delete (g)
forrest@0: 
forrest@0:   xdim  = dimsizes(x1)
forrest@0:   nlev  = xdim(1)
forrest@0: ; y1     = x1(:,0,:,:)
forrest@0:   y1     = x1
forrest@0: ; printVarSummary (y1)
forrest@0:   
forrest@0: ; get the co2 at the lowest level
forrest@0: ; y1     = x1(:,nlev-1,:,:)
forrest@0:   delete (x1)
forrest@0: ;---------------------------------------------
forrest@0: ; g     = addfile(diri2+fili3,"r")
forrest@0: ; x2    = g->CO2
forrest@0: ; delete (g)
forrest@0: ; y2     = x2(:,0,:,:)
forrest@0: ; y2     = x2(:,nlev-1,:,:)
forrest@0: ; delete (x2)
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: ;---------------------------------------------
forrest@0:   y1 = y1 * factor
forrest@0:   y1@_FillValue = 1.e36
forrest@0:   y1@units      = "u mol/mol"
forrest@0: ; y1 = where(y0 .lt. 287.,y1@_FillValue,y1)
forrest@0: ; printVarSummary (y1)
forrest@0: ; print (min(y1)+"/"+max(y1))
forrest@0: ;---------------------------------------------
forrest@0: ; y2 = y2 * factor
forrest@0: ; y2@_FillValue = 1.e36
forrest@0: ; y2@units      = "u mol/mol"
forrest@0: ;---------------------------------------------
forrest@0: ; interpolate into observed station
forrest@0: ; note: model is 0-360E,   90S-90N
forrest@0: ;       ob    is -180-180, 90S-90N
forrest@0: 
forrest@0: ; to be able to handle observation at (-89.98,-24.80)
forrest@0: ; print (yi(0))
forrest@0:   yi(0) = -90.
forrest@0: 
forrest@0:   i = ind(lon_ob .lt. 0.)
forrest@0:   lon_ob(i) = lon_ob(i) + 360.  
forrest@0: ;----------------------------------------------------------------
forrest@0:   yo = linint2_points_Wrap(xi,yi,y1,True,lon_ob,lat_ob,0)
forrest@0:   printVarSummary (yo)
forrest@0: ; yo:[time | 12] x [lev | 26] x [pts | 98]
forrest@0: 
forrest@0:   val_model1 = yo(pts|:,lev|:,time|:)
forrest@0:   delete (yo)
forrest@0:   val_model1_0 = val_model1
forrest@0: ; printVarSummary (val_model1)
forrest@0: ; print (min(val_model1)+"/"+max(val_model1))
forrest@0: 
forrest@0: ; remove annual mean
forrest@0:   val_model1 = val_model1 - conform(val_model1,dim_avg(val_model1),(/0,1/))
forrest@0: ; print (min(val_model1)+"/"+max(val_model1))
forrest@0: 
forrest@0: ;-----------------------------------------------------------------
forrest@0: ;    index of station Barrow, Alaska (71.32,-156.60)
forrest@0:      ind_z = ind(lat_ob .eq. 71.32)
forrest@0:      print (ind_z)
forrest@0:      print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0:      print ("observation at Barrow, Alaska (71.32,-156.60)")
forrest@0:      print (val_ob(ind_z,:))
forrest@0:      print ("model top atm  at Barrow, Alaska (71.32,-156.60)")
forrest@0:      print (val_model1_0(ind_z,0,:))
forrest@0:      print ("model surface at Barrow, Alaska (71.32,-156.60)")
forrest@0:      print (val_model1_0(ind_z,nlev-1,:))
forrest@0:      print ("model top atm  at Barrow, Alaska (71.32,-156.60)")
forrest@0:      print (val_model1(ind_z,0,:))
forrest@0:      print ("model surface at Barrow, Alaska (71.32,-156.60)")
forrest@0:      print (val_model1(ind_z,nlev-1,:))
forrest@0: end