; ***********************************************

; 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