1 | PRO read_seaice, Files, NumObs=NumObs, $ |
---|
2 | Latitudes=Latitudes, Longitudes=Longitudes, $ |
---|
3 | Obs=Obs, modarr=modarr, qc=qcarr, $ |
---|
4 | Dates=Dates, rmdi=rmdi, iobs=outiobs, jobs=outjobs, $ |
---|
5 | nodates=nodates, quiet=quiet |
---|
6 | ;------------------------------------------------------------ |
---|
7 | ;IDL program to read in netcdf files of sea ice data. |
---|
8 | ; |
---|
9 | ;Author: D. J. Lea Feb 2008 |
---|
10 | ; |
---|
11 | ;------------------------------------------------------------ |
---|
12 | rmdi = -99999. |
---|
13 | NumFiles=n_elements(Files) |
---|
14 | ;RefDate = '1950-01-01' |
---|
15 | ;!DATE_SEPARATOR='-' |
---|
16 | RefDate=JULDAY(1,1,1950,0,0) ; should be at 0 UTC |
---|
17 | |
---|
18 | ; could read in from file |
---|
19 | |
---|
20 | ifile2=0 |
---|
21 | for ifile = 0, NumFiles-1 do begin |
---|
22 | ;------------------------------------------------------------ |
---|
23 | ;1. Open the file containing the data |
---|
24 | ncid = ncdf_open(Files(ifile), /nowrite) |
---|
25 | if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+Files(ifile) |
---|
26 | if ncid ge 0 then begin |
---|
27 | if not keyword_set(quiet) then print, 'Opened file '+Files(ifile)+' successfully' |
---|
28 | |
---|
29 | ;------------------------------------------------------------ |
---|
30 | ;2. Read in the dimensions in the file |
---|
31 | ncinfo = ncdf_inquire(ncid) |
---|
32 | ncdf_diminq, ncid, 0, name, NData |
---|
33 | if Ndata gt 0 then begin |
---|
34 | |
---|
35 | ;------------------------------------------------------------ |
---|
36 | ;3. Allocate the data arrays and read in the data |
---|
37 | lons = dblarr(NData) |
---|
38 | lats = dblarr(NData) |
---|
39 | obsval = fltarr(NData) |
---|
40 | modval = fltarr(NData) |
---|
41 | bytQC = bytarr(NData) |
---|
42 | QC = fltarr(NData) |
---|
43 | Dats = dblarr(NData) |
---|
44 | dts = replicate(!dt_base, NData) |
---|
45 | iobs = intarr(NData) |
---|
46 | jobs = intarr(NData) |
---|
47 | |
---|
48 | ; output attribute and variable info |
---|
49 | ; ncattinfo, id |
---|
50 | |
---|
51 | varid = ncdf_varid(ncid, 'lon') |
---|
52 | ; if not keyword_set(quiet) then print,varid |
---|
53 | if (varid eq -1) then varid = ncdf_varid(ncid, 'LONGITUDE') |
---|
54 | ncdf_varget, ncid, varid, lons |
---|
55 | |
---|
56 | varid = ncdf_varid(ncid, 'lat') |
---|
57 | ; if not keyword_set(quiet) then print,varid |
---|
58 | if (varid eq -1) then varid = ncdf_varid(ncid, 'LATITUDE') |
---|
59 | ncdf_varget, ncid, varid, lats |
---|
60 | |
---|
61 | if (keyword_set(nodates) eq 0) then begin |
---|
62 | |
---|
63 | varid = ncdf_varid(ncid, 'JULD') |
---|
64 | if varid ne -1 then begin |
---|
65 | ncdf_varget, ncid, varid, Dats |
---|
66 | dts = Dats+RefDate |
---|
67 | endif else begin |
---|
68 | varid = ncdf_varid(ncid, 'time') |
---|
69 | ncdf_varget, ncid, varid, secs_from_base |
---|
70 | varid = ncdf_varid(ncid, 'SeaIce_dtime') |
---|
71 | ncdf_varget, ncid, varid, dtime |
---|
72 | ncdf_attget, ncid, varid, 'scale_factor', scale_factor |
---|
73 | if not keyword_set(quiet) then print,'dtime(0): ',dtime(0), scale_factor |
---|
74 | dtime=dtime*scale_factor |
---|
75 | ; RefDate = '1981-01-01' |
---|
76 | RefDate = JULDAY(1,1,1981,0,0) ; should be ref from 0 UTC |
---|
77 | dtime = dtime + secs_from_base |
---|
78 | dts = RefDate + dtime/86400. |
---|
79 | endelse |
---|
80 | |
---|
81 | endif |
---|
82 | |
---|
83 | if not keyword_set(quiet) then print,'reading sea ice data' |
---|
84 | varid = ncdf_varid(ncid, 'sea_ice_concentration') |
---|
85 | if not keyword_set(quiet) then print,varid |
---|
86 | if (varid eq -1) then begin |
---|
87 | varid2 = ncdf_varid(ncid, 'SEAICE') |
---|
88 | if (varid2 ne -1) then ncdf_varget, ncid, varid2, obsval |
---|
89 | endif |
---|
90 | if (varid ne -1) then ncdf_varget, ncid, varid, obsval |
---|
91 | if not keyword_set(quiet) then print,'scale_factor' |
---|
92 | scale_factor=1. |
---|
93 | if (varid ne -1) then ncdf_attget, ncid, varid, 'scale_factor', scale_factor |
---|
94 | |
---|
95 | if not keyword_set(quiet) then print,'_FillValue' |
---|
96 | FillValue=99999 |
---|
97 | |
---|
98 | ; ncdf_attget, ncid, varid, '_FillValue', FillValue |
---|
99 | if not keyword_set(quiet) then print,FillValue |
---|
100 | |
---|
101 | pts = where(obsval eq FillValue, count) |
---|
102 | |
---|
103 | if not keyword_set(quiet) then print,'scale_factor ',scale_factor |
---|
104 | obsval=obsval*scale_factor |
---|
105 | |
---|
106 | if count gt 0 then obsval(pts) = rmdi |
---|
107 | |
---|
108 | |
---|
109 | if not keyword_set(quiet) then print,'reading sea ice model values' |
---|
110 | varid = ncdf_varid(ncid, 'SEAICE_Hx') |
---|
111 | if (varid ne -1) then ncdf_varget, ncid, varid, modval |
---|
112 | |
---|
113 | varid = ncdf_varid(ncid, 'IOBS') |
---|
114 | if (varid ne -1) then ncdf_varget, ncid, varid, iobs |
---|
115 | |
---|
116 | varid = ncdf_varid(ncid, 'JOBS') |
---|
117 | if (varid ne -1) then ncdf_varget, ncid, varid, jobs |
---|
118 | |
---|
119 | |
---|
120 | scale_factor=1 |
---|
121 | pts = where(modval eq FillValue or modval eq -9999, count) |
---|
122 | |
---|
123 | modval=modval*scale_factor |
---|
124 | |
---|
125 | if count gt 0 then modval(pts) = rmdi |
---|
126 | |
---|
127 | if not keyword_set(quiet) then print,'sea ice qc' |
---|
128 | |
---|
129 | varid = ncdf_varid(ncid, 'SEAICE_QC') |
---|
130 | if (varid ne -1) then begin |
---|
131 | if not keyword_set(quiet) then print,'bytQC' |
---|
132 | ncdf_varget, ncid, varid, bytQC |
---|
133 | ncdf_attget, ncid, varid, '_FillValue', FillValue |
---|
134 | QC(*) = 0. |
---|
135 | pts = where(bytQC eq FillValue, count) |
---|
136 | if count gt 0 then QC(pts) = rmdi |
---|
137 | pts = where(bytQC ne 48, count) |
---|
138 | if count gt 0 then QC(pts) = 1. |
---|
139 | endif else begin |
---|
140 | varid = ncdf_varid(ncid, 'confidence_flag') |
---|
141 | ncdf_varget, ncid, varid, QC |
---|
142 | endelse |
---|
143 | |
---|
144 | if ifile2 eq 0 then begin |
---|
145 | Latitudes = [float(lats)] |
---|
146 | Longitudes = [float(lons)] |
---|
147 | Obs = [obsval] |
---|
148 | Modarr = [modval] |
---|
149 | QCarr = [QC] |
---|
150 | Dates = [dts] |
---|
151 | if (n_elements(iobs) gt 0) then outiobs = [iobs] |
---|
152 | if (n_elements(jobs) gt 0) then outjobs = [jobs] |
---|
153 | endif else begin |
---|
154 | Latitudes = [Latitudes, float(lats)] |
---|
155 | Longitudes = [Longitudes, float(lons)] |
---|
156 | Obs = [Obs, obsval] |
---|
157 | Modarr = [Modarr, modval] |
---|
158 | QCarr = [QCarr, QC] |
---|
159 | Dates = [Dates, dts] |
---|
160 | if (n_elements(iobs) gt 0) then outiobs = [outiobs, iobs] |
---|
161 | if (n_elements(jobs) gt 0) then outjobs = [outjobs, jobs] |
---|
162 | endelse |
---|
163 | |
---|
164 | ifile2=ifile2+1 |
---|
165 | endif |
---|
166 | |
---|
167 | if not keyword_set(quiet) then print,'closing file' |
---|
168 | |
---|
169 | ncdf_close, ncid |
---|
170 | |
---|
171 | endif |
---|
172 | endfor ; ifile |
---|
173 | |
---|
174 | NumObs = n_elements(Latitudes) |
---|
175 | |
---|
176 | if (n_elements(Modarr) ne NumObs) then Modarr=replicate(rmdi,NumObs) |
---|
177 | |
---|
178 | END |
---|