New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
read_sla.pro in branches/UKMO/icebergs_restart_single_file/NEMOGCM/TOOLS/OBSTOOLS/dataplot – NEMO

source: branches/UKMO/icebergs_restart_single_file/NEMOGCM/TOOLS/OBSTOOLS/dataplot/read_sla.pro @ 6019

Last change on this file since 6019 was 6019, checked in by timgraham, 8 years ago

Reinstated svn keywords before upgrading to head of trunk
  • Property svn:keywords set to Id
File size: 9.6 KB
Line 
1;+
2PRO  read_sla, Files, NumObs=NumObs, $
3               Latitudes=Latitudes, Longitudes=Longitudes, $
4               ObsSLA=ObsSLA, ModSLA=ModSLA, SLAQC=SLAQC, $
5               MDT=MDT, Satellites=Satellites, Types=Types, Dates=Dates, rmdi=rmdi, $
6          iobs=iobs, jobs=jobs, mstp=mstp, quiet=quiet, track=track
7;+-----------------------------------------------------------
8;IDL program to read in netcdf files of altimeter data.
9;
10;Author:  D. J. Lea       - Feb 2008
11;
12;------------------------------------------------------------
13rmdi = -99999.
14NumFiles=n_elements(Files)
15;RefDate = '1950-01-01'
16;!DATE_SEPARATOR='-'
17RefDate=JULDAY(1,1,1950,0,0)     ; should be at 0 UTC
18
19ifile2=0
20for ifile = 0, NumFiles-1 do begin
21;------------------------------------------------------------
22;1. Open the file containing the data
23  ncid = ncdf_open(Files(ifile), /nowrite)
24  if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+Files(ifile)
25  if ncid ge 0 then begin
26  if not keyword_set(quiet) then print, 'Opened file '+Files(ifile)+' successfully'
27
28;------------------------------------------------------------
29;2. Read in the dimensions in the file
30  ncinfo = ncdf_inquire(ncid)
31;  var = ncdf_varinq( ncid, "SLA" )
32  if not keyword_set(quiet) then print,'ncinfo.ndims ',ncinfo.ndims
33  ncdf_diminq, ncid, 0, name, NData
34  NData2=0
35  cycles=1
36  if (ncinfo.ndims gt 2) then ncdf_diminq, ncid, 2, name, NData2
37  if (ncinfo.ndims gt 1) then ncdf_diminq, ncid, 1, name, cycles
38
39
40  if not keyword_set(quiet) then print,NData, NData2
41
42 NData=max([NData*cycles, NData2*cycles])
43
44; print,NData, NData2
45
46  if (NData gt 0) then begin
47;------------------------------------------------------------
48;3. Allocate the data arrays and read in the data
49  lons   = dblarr(NData)
50  lats   = dblarr(NData)
51  obsval = fltarr(NData)
52  modval = fltarr(NData)
53  mdtval = fltarr(NData) 
54  bytQC  = bytarr(NData)
55  QC     = fltarr(NData)
56  Dats   = dblarr(NData)
57  trackval  = intarr(NData)
58  Sats   = intarr(NData)
59  strSats= strarr(NData)
60  StrVal = strarr(8)
61;  dts    = replicate(!dt_base, NData)
62  dts = dblarr(NData)
63 
64;  info,lons
65  varid = ncdf_varid(ncid, 'LONGITUDE')
66  scale_factor=1.
67  if (varid eq -1) then begin
68   varid = ncdf_varid(ncid, 'Longitudes')
69        if (varid eq -1) then begin
70             varid = ncdf_varid(ncid, 'longitude')
71        endif
72        ncdf_attget, ncid, varid, 'scale_factor', scale_factor
73  endif
74  ncdf_varget, ncid, varid, lons
75  lons=lons*scale_factor
76;  info, lons
77  if (cycles gt 1) then begin
78    lonsn=lons
79    for i=1,cycles-1 do begin
80       lonsn=[lonsn,lons]
81    endfor
82    lons=lonsn
83  endif
84
85;  info, lons
86
87;  info,lats     
88  varid = ncdf_varid(ncid, 'LATITUDE')
89  scale_factor=1.
90  if (varid eq -1) then begin
91   varid = ncdf_varid(ncid, 'Latitudes')
92        if (varid eq -1) then begin
93             varid = ncdf_varid(ncid, 'latitude')
94        endif
95        ncdf_attget, ncid, varid, 'scale_factor', scale_factor
96  endif
97  ncdf_varget, ncid, varid, lats
98  lats=lats*scale_factor
99;  info, lats
100  if (cycles gt 1) then begin
101    latsn=lats
102    for i=1,cycles-1 do begin
103       latsn=[latsn,lats]
104    endfor
105    lats=latsn
106  endif
107 
108;  info, lats
109   
110  varid = ncdf_varid(ncid, 'JULD')
111  if (varid eq -1) then begin
112        varid = ncdf_varid(ncid, 'time')
113  endif
114  if (varid ne -1) then begin
115   ncdf_varget, ncid, varid, Dats
116  endif else begin
117        varid = ncdf_varid(ncid, 'BeginDates')
118        ncdf_varget, ncid, varid, Bds
119;        print,bds
120;        info, bds
121        if (cycles gt 1) then begin
122           bds=transpose(bds)
123        endif
124        varid = ncdf_varid(ncid, 'NbPoints')
125        ncdf_varget, ncid, varid, nbpoints
126
127; Read in dataindexes and deltat in cls format file
128; These are used to adjust dates
129        varid = ncdf_varid(ncid, 'DeltaT')
130        if (varid ne -1) then begin
131     ncdf_varget, ncid, varid, deltat
132          ncdf_attget, ncid, varid, 'scale_factor', scale_factor
133          deltat=deltat*scale_factor
134          varid = ncdf_varid(ncid, 'DataIndexes')
135          ncdf_varget, ncid, varid, dataindexes
136         
137        endif
138
139   print, 'deltat ',deltat
140;        print, 'dataindexes ',dataindexes
141
142;        info, nbpoints
143;        print,nbpoints
144;        info, cycles
145        cumbds=[0]
146        nbpointst=nbpoints
147        dataindexest=dataindexes
148        for i=1,cycles-1 do begin
149      nbpointst=[nbpointst,nbpoints]
150                dataindexest=[dataindexest,dataindexes]
151        endfor
152        cumbds=[0,total(nbpointst,/cumulative)]
153;       if not keyword_set(quiet) then print,'cumbds ',cumbds
154        nel=n_elements(Bds)
155;        info,nel
156;       if not keyword_set(quiet) then print,n_elements(dats)
157        for i=0,nel-1 do begin
158;        print,'*',i, cumbds(i+1)-cumbds(i)
159;           print,i,cumbds(i),cumbds(i+1)-1
160           dats(cumbds(i):cumbds(i+1)-1)=bds(i)
161           trackval(cumbds(i):cumbds(i+1)-1)=i
162           dataindexest(cumbds(i):cumbds(i+1)-1)=dataindexest(cumbds(i):cumbds(i+1)-1)-dataindexest(cumbds(i))
163        endfor
164       
165; adjust dats based on dataindex
166   dats=dats+dataindexest*deltat/86400.
167       
168       
169;        print, dats
170  endelse
171
172
173         
174  dts=RefDate + Dats
175   
176  varid = ncdf_varid(ncid, 'MISSION')
177  if (varid ne -1) then begin
178  ncdf_varget, ncid, varid, Sats
179  ncdf_attget, ncid, varid, 'Value_0', StrVal0           
180  ncdf_attget, ncid, varid, 'Value_1', StrVal1 
181  ncdf_attget, ncid, varid, 'Value_2', StrVal2 
182  ncdf_attget, ncid, varid, 'Value_3', StrVal3     
183  ncdf_attget, ncid, varid, 'Value_4', StrVal4 
184  ncdf_attget, ncid, varid, 'Value_5', StrVal5
185  ncdf_attget, ncid, varid, 'Value_6', StrVal6   
186  ncdf_attget, ncid, varid, 'Value_7', StrVal7
187  StrVal=[StrVal0, StrVal1, StrVal2, StrVal3, StrVal4, $
188   StrVal5, StrVal6, StrVal7]   
189  for ival = 0, 7 do begin
190    pts = where(Sats eq ival, count)
191    if count gt 0 then strSats(pts) = StrVal(ival)
192  endfor
193  endif
194
195  varid = ncdf_varid(ncid, 'SLA')
196  ncdf_varget, ncid, varid, obsval
197  ncdf_attget, ncid, varid, '_FillValue', FillValue
198  scale_factor=1.
199  add_offset=0.
200  ;if (ncinfo.ndims gt 2) then ncdf_attget, ncid, varid, 'scale_factor', scale_factor
201  ncviq=ncdf_varinq(ncid, varid)
202  for iatt=0,ncviq.natts-1 do begin
203    ncaiq=ncdf_attname(ncid, varid, iatt)
204    if ( ncaiq eq 'scale_factor') then ncdf_attget, ncid, varid, 'scale_factor', scale_factor
205    if ( ncaiq eq 'add_factor') then ncdf_attget, ncid, varid, 'add_factor', add_factor
206  endfor
207
208  if not keyword_set(quiet) then print,'SLA scale factor: ',scale_factor
209  if not keyword_set(quiet) then print,'SLA add offset: ',add_offset
210  if not keyword_set(quiet) then print,'SLA FillValue: ',FillValue
211  pts = where(obsval eq FillValue, count)
212;  if not keyword_set(quiet) then print,count
213  obsval=obsval*scale_factor+add_offset
214  if count gt 0 then obsval(pts) = rmdi
215 
216  varid = ncdf_varid(ncid, 'SLA_Hx')
217  if (varid ne -1) then begin
218  ncdf_varget, ncid, varid, modval
219  ncdf_attget, ncid, varid, '_FillValue', FillValue
220;  if not keyword_set(quiet) then print,'SLA_Hx FillValue ',FillValue
221  pts = where(modval eq FillValue or modval eq -9999.0, count)
222;  if not keyword_set(quiet) then print,count
223  if count gt 0 then modval(pts) = rmdi
224  endif
225   
226  varid = ncdf_varid(ncid, 'SLA_QC')
227  if (varid ne -1) then begin
228  ncdf_varget, ncid, varid, bytQC
229  ncdf_attget, ncid, varid, '_FillValue', FillValue
230  QC(*) = 0.
231  pts = where(bytQC eq FillValue, count)
232  if count gt 0 then QC(pts) = rmdi
233  pts = where(bytQC ne 48, count)
234  if count gt 0 then QC(pts) = 1.
235  endif
236 
237  varid = ncdf_varid(ncid, 'MDT')
238  if (varid ne -1) then ncdf_varget, ncid, varid, mdtval
239
240  varid = ncdf_varid(ncid, 'data_source')
241  if varid ne -1 then begin
242     ncdf_varget, ncid, varid, obstyp
243  endif else begin
244     varid = ncdf_varid(ncid, 'MISSION')
245     if varid ne -1 then ncdf_varget, ncid, varid, obstyp
246  endelse
247
248; reads in multi cycle data the wrong way round!
249  if (cycles gt 1) then begin
250      obsval=reform(obsval,[cycles,ndata/cycles])
251      obsval=transpose(obsval)
252      modval=reform(modval,[cycles,ndata/cycles])
253      modval=transpose(modval)
254      QC=reform(QC,[cycles,ndata/cycles])
255      QC=transpose(QC)
256  endif
257
258  varid = ncdf_varid(ncid, 'IOBS')
259  if (varid ne -1) then ncdf_varget, ncid, varid, iobsval
260
261  varid = ncdf_varid(ncid, 'JOBS')
262  if (varid ne -1) then ncdf_varget, ncid, varid, jobsval
263
264  varid = ncdf_varid(ncid, 'MSTP')
265  if (varid ne -1) then ncdf_varget, ncid, varid, mstpval
266
267;  info, obssla
268;  info, obsval
269;  info, modsla
270;  info, modval
271 
272  if ifile2 eq 0 then begin
273    Latitudes = [float(lats)]
274    Longitudes = [float(lons)]
275    ObsSLA = [obsval(*)]
276    ModSLA = [modval(*)]
277    MDT = [mdtval]
278    SLAQC = [QC]
279    Dates = [dts]
280    Satellites = [strSats]
281    if (n_elements(obstyp) gt 0) then Types = [obstyp]
282    if (n_elements(iobsval) gt 0) then iobs = [iobsval]
283    if (n_elements(jobsval) gt 0) then jobs = [jobsval]
284   if (n_elements(mstpval) gt 0) then mstp = [mstpval]
285   if (n_elements(trackval) gt 0) then track = [trackval]
286  endif else begin
287    Latitudes = [Latitudes, float(lats)]
288    Longitudes = [Longitudes, float(lons)]
289    ObsSLA = [ObsSLA, obsval(*)]
290    ModSLA = [ModSLA, modval(*)]
291    MDT = [MDT, mdtval]   
292    SLAQC = [SLAQC, QC]
293    Dates = [Dates, dts]   
294    Satellites = [Satellites, strSats]
295    if (n_elements(obstyp) gt 0) then Types = [Types, obstyp]
296    if (n_elements(iobsval) gt 0) then iobs = [iobs, iobsval]
297    if (n_elements(jobsval) gt 0) then jobs = [jobs, jobsval]
298    if (n_elements(mstpval) gt 0) then mstp = [mstp, mstpval]
299    if (n_elements(trackval) gt 0) then tracks = [track, trackval]
300  endelse     
301
302ifile2=ifile2+1
303endif
304  ncdf_close, ncid
305
306endif
307
308endfor ; ifile 
309 
310NumObs = n_elements(Latitudes)
311
312END
Note: See TracBrowser for help on using the repository browser.