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: ; note: model is 0-360E, 90S-90N
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,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: 
forrest@0:   nzone = 4
forrest@0:   do z = 0,nzone-1
forrest@0: 
forrest@0:   if (z .eq. 0) then 
forrest@0: ;    maximum score for the zone, 60N-90N 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 60.)
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 1) then 
forrest@0: ;    maximum score for the zone, 30N-60N 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 2) then 
forrest@0: ;    maximum score for the zone, EQ-30N 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 3) then 
forrest@0: ;    maximum score for the zone, 90S-EQ 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .lt. 0. )
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:  u = ndtooned(val_ob(ind_z,:))
forrest@0:  v = ndtooned(val_model(ind_z,:))
forrest@0: 
forrest@0:  good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0: 
forrest@0:  uu = u(good)
forrest@0:  vv = v(good)
forrest@0: 
forrest@0:  npts = dimsizes(uu)/12
forrest@0:  print (npts)
forrest@0: 
forrest@0:  ccr = esccr(uu,vv,0)
forrest@0:  print (ccr)
forrest@0: 
forrest@0:  un =  onedtond(uu,(/npts,12/))
forrest@0:  vn =  onedtond(vv,(/npts,12/))
forrest@0: ;print (un)
forrest@0: ;print (vn)
forrest@0: 
forrest@0:  score = new((/npts/),float)
forrest@0: 
forrest@0:  do n=0,npts-1
forrest@0:     amp_ob    = max(un(n,:)) - min(un(n,:)) 
forrest@0:     amp_model = max(vn(n,:)) - min(vn(n,:))
forrest@0:     score(n) = 1.-abs((amp_model/amp_ob)-1.)
forrest@0: ;   print (amp_ob)
forrest@0: ;   print (amp_model)  
forrest@0:  end do
forrest@0: ;print (score)
forrest@0: 
forrest@0:  M = avg(score)
forrest@0:  print (M)
forrest@0: 
forrest@0:  M_total = (ccr*ccr + M)*0.5 * score_max
forrest@0:  print (M_total)
forrest@0: 
forrest@0:  delete (ind_z)
forrest@0:  delete (good)
forrest@0:  delete (u)
forrest@0:  delete (v)
forrest@0:  delete (uu)
forrest@0:  delete (vv)
forrest@0:  delete (un)
forrest@0:  delete (vn) 
forrest@0:  delete (score)
forrest@0:  end do
forrest@0: end