clamp: co2/14.metric.ncl@4be95183fbcd

     1 ; ***********************************************

     2 ; xy_4.ncl

     3 ; ***********************************************

     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     7 ;************************************************

     8 begin

     9 ;************************************************

    10 ; read in data: observed

    11 ;************************************************

    12  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"

    13  fili  = "co2_globalView_98.nc"

    14  g     = addfile (diri+fili,"r")

    15  val   = g->CO2_SEAS

    16  lon   = g->LON

    17  lat   = g->LAT

    18  sta   = chartostring(g->STATION)

    19  delete (g)

    21 ;print (sta(0))

    23  ncase = dimsizes(lat)

    24 ;print (ncase)

    26 ;**************************************************************

    27 ; get only the lowest level at each station

    28 ;**************************************************************

    29  lat_tmp = lat

    30  lat_tmp@_FillValue = 1.e+36

    32  do n = 0,ncase-1

    33     if (.not. ismissing(lat_tmp(n))) then

    34        indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)

    35        if (dimsizes(indexes) .gt. 1) then

    36           lat_tmp(indexes(1:)) = lat_tmp@_FillValue

    37        end if

    38        delete (indexes)

    39     end if

    40  end do

    42  indexes = ind(.not. ismissing(lat_tmp))

    43 ;print (dimsizes(indexes))

    44 ;print (indexes)

    46  lat_ob = lat(indexes)

    47  lon_ob = lon(indexes)

    48  val_ob = val(indexes,:)

    49 ;printVarSummary (val_ob)

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

    52 ;************************************************

    53 ; read in model data

    54 ;************************************************

    55   diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    56   fili2 = "b30.061m_401_425_MONS_climo.nc"

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

    59   x     = g->CO2

    60   xi    = g->lon

    61   yi    = g->lat

    62   xdim  = dimsizes(x)

    63   nlev  = xdim(1)

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

    65 ; printVarSummary (y)

    67 ; get the co2 at the lowest level

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

    70 ; change to unit of observed (u mol/mol)

    71 ; Model_units [=] kgCO2 / kgDryAir

    72 ; 28.966 = molecular weight of dry air

    73 ; 44.       = molecular weight of CO2

    74 ; u mol = 1e-6 mol

    76   factor = (28.966/44.) * 1e6

    77   y      = y * factor

    79   y@_FillValue = 1.e36

    80   y@units      = "u mol/mol"

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

    82 ; printVarSummary (y)

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

    85 ; interpolate into observed station

    86 ; note: model is 0-360E, 90S-90N

    88   i = ind(lon_ob .lt. 0.)

    89   lon_ob(i) = lon_ob(i) + 360.

    91   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)

    93   val_model = val_ob

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

    95 ; printVarSummary (val_model)

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

    98 ; remove annual mean

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

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

   102   nzone = 4

   103   do z = 0,nzone-1

   105   if (z .eq. 0) then

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

   107      score_max = 5.0

   108 ;    index of stations in this zone

   109      ind_z = ind(lat_ob .ge. 60.)

   110 ;    print (ind_z)

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

   112 ;    print (val_ob(ind_z,:))

   113 ;    print (val_model(ind_z,:))

   114   end if

   116   if (z .eq. 1) then

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

   118      score_max = 5.0

   119 ;    index of stations in this zone

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

   121 ;    print (ind_z)

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

   123 ;    print (val_ob(ind_z,:))

   124 ;    print (val_model(ind_z,:))

   125   end if

   127   if (z .eq. 2) then

   128 ;    maximum score for the zone, EQ-30N

   129      score_max = 5.0

   130 ;    index of stations in this zone

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

   132 ;    print (ind_z)

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

   134 ;    print (val_ob(ind_z,:))

   135 ;    print (val_model(ind_z,:))

   136   end if

   138   if (z .eq. 3) then

   139 ;    maximum score for the zone, 90S-EQ

   140      score_max = 5.0

   141 ;    index of stations in this zone

   142      ind_z = ind(lat_ob .lt. 0. )

   143 ;    print (ind_z)

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

   145 ;    print (val_ob(ind_z,:))

   146 ;    print (val_model(ind_z,:))

   147   end if

   149  u = ndtooned(val_ob(ind_z,:))

   150  v = ndtooned(val_model(ind_z,:))

   152  good = ind(.not.ismissing(u) .and. .not.ismissing(v))

   154  uu = u(good)

   155  vv = v(good)

   156  st = sta(ind_z)

   158  npts = dimsizes(st)

   159  print (npts)

   161  ccr = esccr(uu,vv,0)

   162 ;print (ccr)

   164  un =  onedtond(uu,(/npts,12/))

   165  vn =  onedtond(vv,(/npts,12/))

   166 ;print (un)

   167 ;print (vn)

   169  score = new((/npts/),float)

   171  do n=0,npts-1

   172     amp_ob    = max(un(n,:)) - min(un(n,:))

   173     amp_model = max(vn(n,:)) - min(vn(n,:))

   174     score(n) = 1.-abs((amp_model/amp_ob)-1.)

   175 ;   print (amp_ob)

   176 ;   print (amp_model)

   177  end do

   178 ;print (score)

   180  M = avg(score)

   181  print (M)

   183  M_total = (ccr*ccr + M)*0.5 * score_max

   184  print (M_total)

   186  delete (ind_z)

   187  delete (good)

   188  delete (u)

   189  delete (v)

   190  delete (uu)

   191  delete (vv)

   192  delete (un)

   193  delete (vn)

   194  delete (st)

   195  delete (score)

   196  end do

   197 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--