clamp: co2/13.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_222.nc"

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

    15  val   = g->CO2_SEAS

    16  lon   = g->LON

    17  lat   = g->LAT

    18  delete (g)

    20  ncase = dimsizes(lat)

    21 ;print (ncase)

    23 ;**************************************************************

    24 ; get only the lowest level at each station

    25 ;**************************************************************

    26  lat_tmp = lat

    27  lat_tmp@_FillValue = 1.e+36

    29  do n = 0,ncase-1

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

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

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

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

    34        end if

    35        delete (indexes)

    36     end if

    37  end do

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

    40 ;print (indexes)

    42  lat_ob = lat(indexes)

    43  lon_ob = lon(indexes)

    44  val_ob = val(indexes,:)

    45 ;printVarSummary (val_ob)

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

    48 ;************************************************

    49 ; read in model data

    50 ;************************************************

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

    52   fili2 = "b30.061m_401_425_MONS_climo.nc"

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

    55   x     = g->CO2

    56   xi    = g->lon

    57   yi    = g->lat

    58   xdim  = dimsizes(x)

    59   nlev  = xdim(1)

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

    61 ; printVarSummary (y)

    63 ; get the co2 at the lowest level

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

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

    67 ; Model_units [=] kgCO2 / kgDryAir

    68 ; 28.966 = molecular weight of dry air

    69 ; 44.       = molecular weight of CO2

    70 ; u mol = 1e-6 mol

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

    73   y      = y * factor

    75   y@_FillValue = 1.e36

    76   y@units      = "u mol/mol"

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

    78 ; printVarSummary (y)

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

    81 ; interpolate into observed station

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

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

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

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

    89   val_model = val_ob

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

    91 ; printVarSummary (val_model)

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

    94 ; remove annual mean

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

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

    98   nzone = 4

    99   do z = 0,nzone-1

   101   if (z .eq. 0) then

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

   103      score_max = 5.0

   104 ;    index of stations in this zone

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

   106 ;    print (ind_z)

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

   108 ;    print (val_ob(ind_z,:))

   109 ;    print (val_model(ind_z,:))

   110   end if

   112   if (z .eq. 1) then

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

   114      score_max = 5.0

   115 ;    index of stations in this zone

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

   117 ;    print (ind_z)

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

   119 ;    print (val_ob(ind_z,:))

   120 ;    print (val_model(ind_z,:))

   121   end if

   123   if (z .eq. 2) then

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

   125      score_max = 5.0

   126 ;    index of stations in this zone

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

   128 ;    print (ind_z)

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

   130 ;    print (val_ob(ind_z,:))

   131 ;    print (val_model(ind_z,:))

   132   end if

   134   if (z .eq. 3) then

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

   136      score_max = 5.0

   137 ;    index of stations in this zone

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

   139 ;    print (ind_z)

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

   141 ;    print (val_ob(ind_z,:))

   142 ;    print (val_model(ind_z,:))

   143   end if

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

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

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

   150  uu = u(good)

   151  vv = v(good)

   153  npts = dimsizes(uu)/12

   154  print (npts)

   156  ccr = esccr(uu,vv,0)

   157  print (ccr)

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

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

   161 ;print (un)

   162 ;print (vn)

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

   166  do n=0,npts-1

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

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

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

   170 ;   print (amp_ob)

   171 ;   print (amp_model)

   172  end do

   173 ;print (score)

   175  M = avg(score)

   176  print (M)

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

   179  print (M_total)

   181  delete (ind_z)

   182  delete (good)

   183  delete (u)

   184  delete (v)

   185  delete (uu)

   186  delete (vv)

   187  delete (un)

   188  delete (vn)

   189  delete (score)

   190  end do

   191 end

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