1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/co2/24.check.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,143 @@
1.4 +; ***********************************************
1.5 +; add another model to plot
1.6 +; add panel plot to 22.lines.ncl
1.7 +; add zone plot to 21.lines.ncl
1.8 +; ***********************************************
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.12 +;************************************************
1.13 +begin
1.14 +;************************************************
1.15 +; read in data: observed
1.16 +;************************************************
1.17 + diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
1.18 + fili = "co2_globalView_98.nc"
1.19 + g = addfile (diri+fili,"r")
1.20 + val = g->CO2_SEAS
1.21 + lon = g->LON
1.22 + lat = g->LAT
1.23 + sta = chartostring(g->STATION)
1.24 + delete (g)
1.25 +
1.26 +;print (sta(0))
1.27 +
1.28 + ncase = dimsizes(lat)
1.29 +;print (ncase)
1.30 +
1.31 +;**************************************************************
1.32 +; get only the lowest level at each station
1.33 +;**************************************************************
1.34 + lat_tmp = lat
1.35 + lat_tmp@_FillValue = 1.e+36
1.36 +
1.37 + do n = 0,ncase-1
1.38 + if (.not. ismissing(lat_tmp(n))) then
1.39 + indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
1.40 + if (dimsizes(indexes) .gt. 1) then
1.41 + lat_tmp(indexes(1:)) = lat_tmp@_FillValue
1.42 + end if
1.43 + delete (indexes)
1.44 + end if
1.45 + end do
1.46 +
1.47 + indexes = ind(.not. ismissing(lat_tmp))
1.48 +;print (dimsizes(indexes))
1.49 +;print (indexes)
1.50 +
1.51 + lat_ob = lat(indexes)
1.52 + lon_ob = lon(indexes)
1.53 + val_ob = val(indexes,:)
1.54 +;printVarSummary (val_ob)
1.55 +;print (lat_ob +"/"+lon_ob)
1.56 +
1.57 +;************************************************
1.58 +; read in model data
1.59 +;************************************************
1.60 + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.61 + fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"
1.62 + fili3 = "b30.061m_401_425_MONS_climo_atm.nc"
1.63 +;--------------------------------------------
1.64 + g = addfile(diri2+fili2,"r")
1.65 + x1 = g->CO2
1.66 + xi = g->lon
1.67 + yi = g->lat
1.68 + delete (g)
1.69 +
1.70 + xdim = dimsizes(x1)
1.71 + nlev = xdim(1)
1.72 +; y1 = x1(:,0,:,:)
1.73 + y1 = x1
1.74 +; printVarSummary (y1)
1.75 +
1.76 +; get the co2 at the lowest level
1.77 +; y1 = x1(:,nlev-1,:,:)
1.78 + delete (x1)
1.79 +;---------------------------------------------
1.80 +; g = addfile(diri2+fili3,"r")
1.81 +; x2 = g->CO2
1.82 +; delete (g)
1.83 +; y2 = x2(:,0,:,:)
1.84 +; y2 = x2(:,nlev-1,:,:)
1.85 +; delete (x2)
1.86 +;---------------------------------------------
1.87 +; change to unit of observed (u mol/mol)
1.88 +; Model_units [=] kgCO2 / kgDryAir
1.89 +; 28.966 = molecular weight of dry air
1.90 +; 44. = molecular weight of CO2
1.91 +; u mol = 1e-6 mol
1.92 +
1.93 + factor = (28.966/44.) * 1e6
1.94 +;---------------------------------------------
1.95 + y1 = y1 * factor
1.96 + y1@_FillValue = 1.e36
1.97 + y1@units = "u mol/mol"
1.98 +; y1 = where(y0 .lt. 287.,y1@_FillValue,y1)
1.99 +; printVarSummary (y1)
1.100 +; print (min(y1)+"/"+max(y1))
1.101 +;---------------------------------------------
1.102 +; y2 = y2 * factor
1.103 +; y2@_FillValue = 1.e36
1.104 +; y2@units = "u mol/mol"
1.105 +;---------------------------------------------
1.106 +; interpolate into observed station
1.107 +; note: model is 0-360E, 90S-90N
1.108 +; ob is -180-180, 90S-90N
1.109 +
1.110 +; to be able to handle observation at (-89.98,-24.80)
1.111 +; print (yi(0))
1.112 + yi(0) = -90.
1.113 +
1.114 + i = ind(lon_ob .lt. 0.)
1.115 + lon_ob(i) = lon_ob(i) + 360.
1.116 +;----------------------------------------------------------------
1.117 + yo = linint2_points_Wrap(xi,yi,y1,True,lon_ob,lat_ob,0)
1.118 + printVarSummary (yo)
1.119 +; yo:[time | 12] x [lev | 26] x [pts | 98]
1.120 +
1.121 + val_model1 = yo(pts|:,lev|:,time|:)
1.122 + delete (yo)
1.123 + val_model1_0 = val_model1
1.124 +; printVarSummary (val_model1)
1.125 +; print (min(val_model1)+"/"+max(val_model1))
1.126 +
1.127 +; remove annual mean
1.128 + val_model1 = val_model1 - conform(val_model1,dim_avg(val_model1),(/0,1/))
1.129 +; print (min(val_model1)+"/"+max(val_model1))
1.130 +
1.131 +;-----------------------------------------------------------------
1.132 +; index of station Barrow, Alaska (71.32,-156.60)
1.133 + ind_z = ind(lat_ob .eq. 71.32)
1.134 + print (ind_z)
1.135 + print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
1.136 + print ("observation at Barrow, Alaska (71.32,-156.60)")
1.137 + print (val_ob(ind_z,:))
1.138 + print ("model top atm at Barrow, Alaska (71.32,-156.60)")
1.139 + print (val_model1_0(ind_z,0,:))
1.140 + print ("model surface at Barrow, Alaska (71.32,-156.60)")
1.141 + print (val_model1_0(ind_z,nlev-1,:))
1.142 + print ("model top atm at Barrow, Alaska (71.32,-156.60)")
1.143 + print (val_model1(ind_z,0,:))
1.144 + print ("model surface at Barrow, Alaska (71.32,-156.60)")
1.145 + print (val_model1(ind_z,nlev-1,:))
1.146 +end