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

 ; xy_4.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.061m_401_425_MONS_climo.nc"

   g     = addfile(diri2+fili2,"r")

   x     = g->CO2

   xi    = g->lon

   yi    = g->lat

   xdim  = dimsizes(x)

   nlev  = xdim(1)

   y     = x(:,0,:,:)

 ; printVarSummary (y)

 ; get the co2 at the lowest level

   y     = x(:,nlev-1,:,:)

 ; 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

   y      = y * factor

   y@_FillValue = 1.e36

   y@units      = "u mol/mol"

 ; y = where(y0 .lt. 287.,y@_FillValue,y)

 ; printVarSummary (y)

 ; print (min(y)+"/"+max(y))

 ; interpolate into observed station

 ; note: model is 0-360E, 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,y,True,lon_ob,lat_ob,0)

   val_model = yo(pts|:,time|:)

 ; printVarSummary (val_model)

 ; print (min(val_model)+"/"+max(val_model))

 ; remove annual mean

   val_model = val_model - conform(val_model,dim_avg(val_model),0)

 ; print (min(val_model)+"/"+max(val_model))

   plot_sta      = new((/2,12/),float)

   plot_sta@long_name = "Seasonal CO2"

   mon = ispan(1,12,1)

   mon@long_name = "month"

   plot_type = "ps"

   plot_type_new = "png"

   res                   = True                      ; plot mods desired

   res@xyLineThicknesses = (/1.0,2.0/)               ; make 2nd lines thicker

   res@xyLineColors      = (/"blue","red"/)          ; change line color

 ;------------------------------------------------------------------

 ; Add a boxed legend using the more simple method, which won't have

 ; vertical lines going through the markers.

   res@pmLegendDisplayMode    = "Always"

 ; res@pmLegendWidthF         = 0.1

   res@pmLegendWidthF         = 0.08

   res@pmLegendHeightF        = 0.05

 ; res@pmLegendOrthogonalPosF = -1.17

 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

 ; res@pmLegendParallelPosF   =  0.18

   res@pmLegendParallelPosF   =  0.23  ;(rightward)

 ; res@lgPerimOn             = False

   res@lgLabelFontHeightF     = 0.015

   res@xyExplicitLegendLabels = (/"observed","model_b30.061m"/)

 ;-------------------------------------------------------------------

   nzone = 4

   do z = 0,nzone-1

   if (z .eq. 0) then 

 ;    maximum score for the zone, 60N-90N 

      score_max = 5.0

 ;    index of stations in this zone

      ind_z = ind(lat_ob .ge. 60.)

 ;    print (ind_z)

 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

 ;    print (val_ob(ind_z,:))

 ;    print (val_model(ind_z,:))

   end if

   if (z .eq. 1) then 

 ;    maximum score for the zone, 30N-60N 

      score_max = 5.0

 ;    index of stations in this zone

      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)

 ;    print (ind_z)

 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

 ;    print (val_ob(ind_z,:))

 ;    print (val_model(ind_z,:))

   end if

   if (z .eq. 2) then 

 ;    maximum score for the zone, EQ-30N 

      score_max = 5.0

 ;    index of stations in this zone

      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)

 ;    print (ind_z)

 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

 ;    print (val_ob(ind_z,:))

 ;    print (val_model(ind_z,:))

   end if

   if (z .eq. 3) then 

 ;    maximum score for the zone, 90S-EQ 

      score_max = 5.0

 ;    index of stations in this zone

      ind_z = ind(lat_ob .lt. 0. )

 ;    print (ind_z)

 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

 ;    print (val_ob(ind_z,:))

 ;    print (val_model(ind_z,:))

   end if

  npts = dimsizes(ind_z)

  print (npts)

  do n=0,npts-1

     plot_sta(0,:) = (/val_ob(ind_z(n),:)/)

     plot_sta(1,:) = (/val_model(ind_z(n),:)/)

     title = sta(ind_z(n))+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"

     plot_name = sta(ind_z(n))

 ;   print (title)

 ;   print (plot_name)

     wks = gsn_open_wks (plot_type,plot_name)          ; open workstation

     res@tiMainString      = title                     ; add title

     plot  = gsn_csm_xy (wks,mon,plot_sta,res)         ; create plot

     frame(wks)

     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

     system("rm "+plot_name+"."+plot_type)

     system("rm "+plot_name+"-1."+plot_type_new)

     system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

     clear (wks)

  end do

  delete (ind_z)

  end do

 end