source: utils/tools_CPL_OASIS_ticket2379_2020/OBSTOOLS/dataplot/read_seaice.pro

Last change on this file was 3026, checked in by djlea, 9 years ago

Remove unneeded sec_to_dt from dataplot. Also add obstools build instructions to documentation.

File size: 5.4 KB
Line 
1PRO  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;------------------------------------------------------------
12rmdi = -99999.
13NumFiles=n_elements(Files)
14;RefDate = '1950-01-01'
15;!DATE_SEPARATOR='-'
16RefDate=JULDAY(1,1,1950,0,0)     ; should be at 0 UTC
17
18; could read in from file
19
20ifile2=0
21for 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
164ifile2=ifile2+1
165endif
166
167  if not keyword_set(quiet) then print,'closing file'
168
169  ncdf_close, ncid
170
171endif
172endfor ; ifile 
173 
174NumObs = n_elements(Latitudes)
175
176if (n_elements(Modarr) ne NumObs) then Modarr=replicate(rmdi,NumObs)
177
178END
Note: See TracBrowser for help on using the repository browser.