|
1 ; *********************************************** |
|
2 ; add another model to plot |
|
3 ; add panel plot to 22.lines.ncl |
|
4 ; add zone plot to 21.lines.ncl |
|
5 ; *********************************************** |
|
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
9 ;************************************************ |
|
10 begin |
|
11 ;************************************************ |
|
12 ; read in data: observed |
|
13 ;************************************************ |
|
14 diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/" |
|
15 fili = "co2_globalView_98.nc" |
|
16 g = addfile (diri+fili,"r") |
|
17 val = g->CO2_SEAS |
|
18 lon = g->LON |
|
19 lat = g->LAT |
|
20 sta = chartostring(g->STATION) |
|
21 delete (g) |
|
22 |
|
23 ;print (sta(0)) |
|
24 |
|
25 ncase = dimsizes(lat) |
|
26 ;print (ncase) |
|
27 |
|
28 ;************************************************************** |
|
29 ; get only the lowest level at each station |
|
30 ;************************************************************** |
|
31 lat_tmp = lat |
|
32 lat_tmp@_FillValue = 1.e+36 |
|
33 |
|
34 do n = 0,ncase-1 |
|
35 if (.not. ismissing(lat_tmp(n))) then |
|
36 indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon) |
|
37 if (dimsizes(indexes) .gt. 1) then |
|
38 lat_tmp(indexes(1:)) = lat_tmp@_FillValue |
|
39 end if |
|
40 delete (indexes) |
|
41 end if |
|
42 end do |
|
43 |
|
44 indexes = ind(.not. ismissing(lat_tmp)) |
|
45 ;print (dimsizes(indexes)) |
|
46 ;print (indexes) |
|
47 |
|
48 lat_ob = lat(indexes) |
|
49 lon_ob = lon(indexes) |
|
50 val_ob = val(indexes,:) |
|
51 ;printVarSummary (val_ob) |
|
52 ;print (lat_ob +"/"+lon_ob) |
|
53 |
|
54 ;************************************************ |
|
55 ; read in model data |
|
56 ;************************************************ |
|
57 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
58 fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc" |
|
59 fili3 = "b30.061m_401_425_MONS_climo_atm.nc" |
|
60 ;-------------------------------------------- |
|
61 g = addfile(diri2+fili2,"r") |
|
62 x1 = g->CO2 |
|
63 xi = g->lon |
|
64 yi = g->lat |
|
65 delete (g) |
|
66 |
|
67 xdim = dimsizes(x1) |
|
68 nlev = xdim(1) |
|
69 ; y1 = x1(:,0,:,:) |
|
70 y1 = x1 |
|
71 ; printVarSummary (y1) |
|
72 |
|
73 ; get the co2 at the lowest level |
|
74 ; y1 = x1(:,nlev-1,:,:) |
|
75 delete (x1) |
|
76 ;--------------------------------------------- |
|
77 ; g = addfile(diri2+fili3,"r") |
|
78 ; x2 = g->CO2 |
|
79 ; delete (g) |
|
80 ; y2 = x2(:,0,:,:) |
|
81 ; y2 = x2(:,nlev-1,:,:) |
|
82 ; delete (x2) |
|
83 ;--------------------------------------------- |
|
84 ; change to unit of observed (u mol/mol) |
|
85 ; Model_units [=] kgCO2 / kgDryAir |
|
86 ; 28.966 = molecular weight of dry air |
|
87 ; 44. = molecular weight of CO2 |
|
88 ; u mol = 1e-6 mol |
|
89 |
|
90 factor = (28.966/44.) * 1e6 |
|
91 ;--------------------------------------------- |
|
92 y1 = y1 * factor |
|
93 y1@_FillValue = 1.e36 |
|
94 y1@units = "u mol/mol" |
|
95 ; y1 = where(y0 .lt. 287.,y1@_FillValue,y1) |
|
96 ; printVarSummary (y1) |
|
97 ; print (min(y1)+"/"+max(y1)) |
|
98 ;--------------------------------------------- |
|
99 ; y2 = y2 * factor |
|
100 ; y2@_FillValue = 1.e36 |
|
101 ; y2@units = "u mol/mol" |
|
102 ;--------------------------------------------- |
|
103 ; interpolate into observed station |
|
104 ; note: model is 0-360E, 90S-90N |
|
105 ; ob is -180-180, 90S-90N |
|
106 |
|
107 ; to be able to handle observation at (-89.98,-24.80) |
|
108 ; print (yi(0)) |
|
109 yi(0) = -90. |
|
110 |
|
111 i = ind(lon_ob .lt. 0.) |
|
112 lon_ob(i) = lon_ob(i) + 360. |
|
113 ;---------------------------------------------------------------- |
|
114 yo = linint2_points_Wrap(xi,yi,y1,True,lon_ob,lat_ob,0) |
|
115 printVarSummary (yo) |
|
116 ; yo:[time | 12] x [lev | 26] x [pts | 98] |
|
117 |
|
118 val_model1 = yo(pts|:,lev|:,time|:) |
|
119 delete (yo) |
|
120 val_model1_0 = val_model1 |
|
121 ; printVarSummary (val_model1) |
|
122 ; print (min(val_model1)+"/"+max(val_model1)) |
|
123 |
|
124 ; remove annual mean |
|
125 val_model1 = val_model1 - conform(val_model1,dim_avg(val_model1),(/0,1/)) |
|
126 ; print (min(val_model1)+"/"+max(val_model1)) |
|
127 |
|
128 ;----------------------------------------------------------------- |
|
129 ; index of station Barrow, Alaska (71.32,-156.60) |
|
130 ind_z = ind(lat_ob .eq. 71.32) |
|
131 print (ind_z) |
|
132 print (lat_ob(ind_z)+"/"+lon_ob(ind_z)) |
|
133 print ("observation at Barrow, Alaska (71.32,-156.60)") |
|
134 print (val_ob(ind_z,:)) |
|
135 print ("model top atm at Barrow, Alaska (71.32,-156.60)") |
|
136 print (val_model1_0(ind_z,0,:)) |
|
137 print ("model surface at Barrow, Alaska (71.32,-156.60)") |
|
138 print (val_model1_0(ind_z,nlev-1,:)) |
|
139 print ("model top atm at Barrow, Alaska (71.32,-156.60)") |
|
140 print (val_model1(ind_z,0,:)) |
|
141 print ("model surface at Barrow, Alaska (71.32,-156.60)") |
|
142 print (val_model1(ind_z,nlev-1,:)) |
|
143 end |