1 ; ***********************************************
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 ;************************************************
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")
18 sta = chartostring(g->STATION)
26 ;**************************************************************
27 ; get only the lowest level at each station
28 ;**************************************************************
30 lat_tmp@_FillValue = 1.e+36
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
42 indexes = ind(.not. ismissing(lat_tmp))
43 ;print (dimsizes(indexes))
48 val_ob = val(indexes,:)
49 ;printVarSummary (val_ob)
50 ;print (lat_ob +"/"+lon_ob)
52 ;************************************************
54 ;************************************************
55 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
56 ; fili2 = "b30.061m_401_425_MONS_climo.nc"
57 fili2 = "b30.061n_1995-2004_MONS_climo.nc"
59 g = addfile(diri2+fili2,"r")
68 ; get the co2 at the lowest level
71 ; change to unit of observed (u mol/mol)
72 ; Model_units [=] kgCO2 / kgDryAir
73 ; 28.966 = molecular weight of dry air
74 ; 44. = molecular weight of CO2
77 factor = (28.966/44.) * 1e6
82 ; y = where(y0 .lt. 287.,y@_FillValue,y)
84 ; print (min(y)+"/"+max(y))
86 ; interpolate into observed station
87 ; note: model is 0-360E, 90S-90N
89 ; to be able to handle observation at (-89.98,-24.80)
93 i = ind(lon_ob .lt. 0.)
94 lon_ob(i) = lon_ob(i) + 360.
96 yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
98 val_model = yo(pts|:,time|:)
99 ; printVarSummary (val_model)
100 ; print (min(val_model)+"/"+max(val_model))
103 val_model = val_model - conform(val_model,dim_avg(val_model),0)
104 ; print (min(val_model)+"/"+max(val_model))
110 ; maximum score for the zone, 60N-90N
112 ; index of stations in this zone
113 ind_z = ind(lat_ob .ge. 60.)
115 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
116 ; print (val_ob(ind_z,:))
117 ; print (val_model(ind_z,:))
121 ; maximum score for the zone, 30N-60N
123 ; index of stations in this zone
124 ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
126 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
127 ; print (val_ob(ind_z,:))
128 ; print (val_model(ind_z,:))
132 ; maximum score for the zone, EQ-30N
134 ; index of stations in this zone
135 ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
137 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
138 ; print (val_ob(ind_z,:))
139 ; print (val_model(ind_z,:))
143 ; maximum score for the zone, 90S-EQ
145 ; index of stations in this zone
146 ind_z = ind(lat_ob .lt. 0. )
148 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
149 ; print (val_ob(ind_z,:))
150 ; print (val_model(ind_z,:))
153 npts = dimsizes(ind_z)
156 amp_ob = new((/npts/),float)
157 amp_model = new((/npts/),float)
159 amp_ratio_sta = new((/npts/),float)
160 ccr_sta = new((/npts/),float)
161 M_sta = new((/npts/),float)
162 score_sta = new((/npts/),float)
165 amp_ob(n) = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))
166 amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
168 amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
169 ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
170 M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
171 score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
173 print (sta(ind_z(n))+"/"+lat(ind_z(n))+"/"+lon(ind_z(n))+"/"+amp_ratio_sta(n)+"/"+ccr_sta(n)+"/"+M_sta(n)+"/"+score_sta(n))
176 amp_ratio_zone = avg(amp_ratio_sta)
177 ccr_zone = avg(ccr_sta)
178 M_zone = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)
179 score_zone = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
181 print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)
186 delete (amp_ratio_sta)