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_feedback.pro in branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/TOOLS/OBSTOOLS/dataplot – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/TOOLS/OBSTOOLS/dataplot/read_feedback.pro @ 5985

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

Reinstate keywords before upgrading to head of trunk

  • Property svn:keywords set to Id
File size: 9.2 KB
Line 
1PRO read_feedback, Files, DepRange=DepRange, VarName=VarNames, NumData=NumData, $
2                OutLats=OutLats, OutLons=OutLons, $
3                OutInstruments=OutInst, OutPlatform=OutPlatform, $
4                OutDeps=OutDeps,  $
5                OutObs1=OutObs1, OutObs2=OutObs2, $
6                OutMod1=OutMod1, OutMod2=OutMod2, $
7                OutQC1=OutQC1, OutQC2=OutQC2, $
8                MDT=MDT, OutDates=OutDates, $
9                rmdi=rmdi, quiet=quiet, $
10                OutProfileNum=OutProfileNum
11
12;------------------------------------------------------------------------------------------
13;Program to read in data from a feedback format file
14;
15; Inputs:
16;        File        => File name to be read in
17;        DepRange    => (optional) Range of depths for which the data is to be extracted.
18;
19; Outputs:
20
21;
22;Author: Matt Martin. 29/09/09
23;------------------------------------------------------------------------------------------
24rmdi = 99999.
25if NOT keyword_set(DepRange) then DepRange=[0.,1.e1]   ;Default to top 10 metres.
26NumFiles=n_elements(Files)
27
28ifile2=0
29numdata=0
30for ifile = 0, NumFiles-1 do begin
31  File = Files(ifile)
32 
33;1. Open the netcdf file
34  ncid = ncdf_open(File, /nowrite)
35  if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+File
36  if ncid ge 0 then begin
37  if not keyword_set(quiet) then print, 'Opened file '+File+' successfully'
38 
39  ;2. Get info about file
40  ncinfo = ncdf_inquire(ncid)
41
42  ;3. Read in the relevant dimensions
43  NumExtra = 0
44  for idim = 0, ncinfo.ndims-1 do begin
45    ncdf_diminq, ncid, idim, name, dimlen
46    if name eq 'N_OBS' then NumObs = dimlen $
47    else if name eq 'N_LEVELS' then NumLevels = dimlen $
48    else if name eq 'N_VARS' then NumVars = dimlen $
49    else if name eq 'N_QCF' then NumFlags = dimlen $
50    else if name eq 'N_ENTRIES' then NumEntries = dimlen $
51    else if name eq 'N_EXTRA' then NumExtra = dimlen $
52    else if name eq 'STRINGNAM' then strnam = dimlen $
53    else if name eq 'STRINGGRID' then strgrid = dimlen $
54    else if name eq 'STRINGWMO' then strwmo = dimlen $
55    else if name eq 'STRINGTYP' then strtyp = dimlen $
56    else if name eq 'STRINGJULD' then strjuld = dimlen
57  endfor
58 
59  ;4. Set up arrays
60
61  if (NumObs gt 0) then begin 
62    ByteNam = intarr(strnam, NumVars)
63    if (n_elements(NumEntries) eq 1) then ByteEntries = intarr(strnam, NumEntries)
64    ByteId = intarr(strwmo, NumObs)
65    ByteType = intarr(strtyp, NumObs)
66    ByteJulDRef = intarr(strjuld)
67    VarNames = strarr(NumVars)
68    if (n_elements(NumEntries) eq 1) then Entries = strarr(NumEntries)
69    Id = strarr(NumObs)
70    Type = strarr(NumObs)
71    if NumExtra gt 0 then begin
72      ByteExtra = intarr(strnam, NumExtra)   
73      Extra = strarr(NumExtra)   
74    endif
75    Longitude = fltarr(NumObs)
76    Latitude  = fltarr(NumObs)
77    Depth     = fltarr(NumLevels, NumObs)
78    DepQC     = intarr(NumLevels, NumObs)
79    JulD      = dblarr(NumObs)
80    ObsQC     = intarr(NumObs)
81    PosQC     = intarr(NumObs)
82    JulDQC    = intarr(NumObs)
83    Obs       = fltarr(NumLevels, NumObs, NumVars)
84    Hx        = fltarr(NumLevels, NumObs, NumVars)
85    VarQC     = intarr(NumObs, NumVars)
86    LevelQC   = intarr(NumLevels, NumObs, NumVars)
87   
88  ;5. Read in the data 
89    varid = ncdf_varid(ncid, 'VARIABLES')
90    ncdf_varget, ncid, varid, ByteNam
91;    info, VarNames
92;    info, ByteNam
93    for ivar= 0, NumVars-1 do begin
94      VarNames(ivar) = string(ByteNam(*,ivar))
95    endfor
96   
97    if ifile2 eq 0 then VarName = VarNames(0)
98    if VarName ne VarNames(0) then message, 'Can only read in from files containing the same variables'
99   
100    varid = ncdf_varid(ncid, 'JULD_REFERENCE')
101    ncdf_varget, ncid, varid, ByteJulDRef
102    JulDRef = string(ByteJulDRef)
103    RefDate = JULDAY(fix(strmid(JulDRef,4,2)), fix(strmid(JulDRef,6,2)), fix(strmid(JulDRef,0,4)), $
104                     fix(strmid(JulDRef,8,2)), fix(strmid(JulDRef,10,2)), fix(strmid(JulDRef,12,2)))
105   print, 'RefDate: ',RefDate
106    varid = ncdf_varid(ncid, 'JULD')
107    ncdf_varget, ncid, varid, JulD
108;    print, JulD
109    Dates = dblarr(NumLevels, NumObs)
110    for iob = 0L, long(NumObs)-1L do begin
111      dt = RefDate + JulD(iob)
112;      print,dt, JulD(iob)
113      Dates(0:NumLevels-1, iob) = replicate(dt, NumLevels)
114    endfor
115   
116    varid = ncdf_varid(ncid, 'STATION_IDENTIFIER')
117    ncdf_varget, ncid, varid, ByteId
118    Identifier = string(ByteId)
119   
120    varid = ncdf_varid(ncid, 'STATION_TYPE')
121    ncdf_varget, ncid, varid, ByteType
122    Type = string(ByteType)
123   
124    varid = ncdf_varid(ncid, 'LONGITUDE')
125    ncdf_varget, ncid, varid, Longitude
126 
127    varid = ncdf_varid(ncid, 'LATITUDE')
128    ncdf_varget, ncid, varid, Latitude
129 
130    varid = ncdf_varid(ncid, 'DEPTH')
131    ncdf_varget, ncid, varid, Depth
132    ncdf_attget, ncid, varid, '_Fillvalue', FillVal
133    pts = where(Depth eq FillVal, count)
134    if count gt 0 then Depth(pts) = rmdi
135   
136    varid = ncdf_varid(ncid, 'OBSERVATION_QC')
137    ncdf_varget, ncid, varid, ObsQC
138    ncdf_attget, ncid, varid, '_Fillvalue', FillVal
139    pts = where(ObsQC eq FillVal, count)
140    if count gt 0 then ObsQC(pts) = rmdi     
141   
142    for ivar = 0, NumVars-1 do begin
143   
144      varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_OBS')
145      ncdf_varget, ncid, varid, tmp
146      Obs(*,*,ivar) = tmp
147      ncdf_attget, ncid, varid, '_Fillvalue', FillValObs
148 
149      varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_Hx')
150      if (varid gt -1) then begin
151         ncdf_varget, ncid, varid, tmp
152         Hx(*,*,ivar) = tmp
153         ncdf_attget, ncid, varid, '_Fillvalue', FillValHx
154      endif else begin
155         FillValHx=0
156      endelse
157     
158      varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_QC')
159      ncdf_varget, ncid, varid, tmp
160      VarQC(*,ivar) = tmp
161      ncdf_attget, ncid, varid, '_Fillvalue', FillValQC
162     
163     
164      varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_LEVEL_QC')
165      ncdf_varget, ncid, varid, tmp
166      LevelQC(*,*,ivar) = tmp
167      ncdf_attget, ncid, varid, '_Fillvalue', FillValLevQC
168         
169    endfor       
170 
171;  print,' DJL levelqc(*,218,0) ',levelqc(*,218,0)
172;  print,' DJL levelqc(*,218,1) ',levelqc(*,218,1)
173   
174    pts = where(Obs eq FillValObs, count)
175    if count gt 0 then Obs(pts) = rmdi
176    pts = where(Hx eq FillValHx, count)
177    if count gt 0 then Hx(pts) = rmdi   
178    pts = where(VarQC eq FillValQC, count)
179    if count gt 0 then VarQC(pts) = rmdi 
180    pts = where(LevelQC eq FillValLevQC, count)
181    if count gt 0 then LevelQC(pts) = rmdi       
182
183    Obs1 = Obs(*,*,0)
184    Hx1  = Hx(*,*,0)
185    VarQC1=LevelQC(*,*,0)   
186    if NumVars gt 1 then begin
187      Obs2 = Obs(*,*,1)
188      Hx2  = Hx(*,*,1)
189      VarQC2=LevelQC(*,*,1)   
190    endif
191   
192    if strtrim(VarName,2) eq 'SLA' then begin
193      MDT = fltarr(NumLevels, NumObs)
194      print, NumLevels, NumObs
195      varid = ncdf_varid(ncid, 'MDT')
196      ncdf_varget, ncid, varid, MDT
197      ncdf_attget, ncid, varid, '_Fillvalue', FillVal
198      pts = where(MDT eq FillVal, count)
199      if count gt 0 then MDT(pts) = rmdi   
200      Obs2 = MDT
201    endif
202         
203    ;6. Put the data into the correct form to be output from this routine
204    Platformsa = strarr(NumLevels, NumObs)
205    Instrumentsa = strarr(NumLevels, NumObs)
206    Lons = fltarr(NumLevels,NumObs)
207    Lats = fltarr(NumLevels,NumObs)
208    ProfileNum = lonarr(NumLevels,NumObs)
209       
210    for i=0,NumLevels-1 do begin
211      Platformsa(i,*)=Type
212      Instrumentsa(i,*)=Identifier
213    endfor
214     
215    for i=0L,long(NumObs)-1L do begin
216      Lats(*,i)=Latitude(i)
217      Lons(*,i)=Longitude(i)
218      ProfileNum(*,i) = i+1
219    endfor 
220;   
221   pts = where(Depth ge DepRange(0) and Depth le DepRange(1), NumDataInc)
222   NumData=NumData+NumDataInc
223;   
224   if ifile2 eq 0 then begin 
225     OutObs1 = [Obs1(pts)]
226     OutMod1 = [Hx1(pts)]
227     OutQC1  = [VarQC1(pts)]
228     OutDeps = [Depth(pts)]
229     OutLats = [Lats(pts)]
230     OutLons = [Lons(pts)]
231     OutDates= [Dates(pts)]
232     OutInst= [Instrumentsa(pts)]   
233     OutPlatform= [Platformsa(pts)]
234     OutProfileNum = [ProfileNum(pts)]
235     if NumVars eq 2 then begin
236       OutObs2 = [Obs2(pts)]     
237       OutMod2 = [Hx2(pts)]     
238       OutQC2  = [VarQC2(pts)]
239     endif
240   endif else begin
241     OutObs1 = [OutObs1, Obs1(pts)]
242     OutMod1 = [OutMod1, Hx1(pts)]
243     OutQC1  = [OutQC1, VarQC1(pts)]
244     OutDeps = [OutDeps, Depth(pts)]
245     OutLats = [OutLats, Lats(pts)]
246     OutLons = [OutLons, Lons(pts)]
247     OutDates= [OutDates, Dates(pts)]
248     OutInst= [OutInst, Instrumentsa(pts)]   
249     OutPlatform= [OutPlatform, Platformsa(pts)]
250     OutProfileNum = [OutProfileNum, ProfileNum(pts)]
251     if NumVars eq 2 then begin
252       OutObs2 = [OutObs2, Obs2(pts)]     
253       OutMod2 = [OutMod2, Hx2(pts)]     
254       OutQC2  = [OutQC2, VarQC2(pts)]
255     endif
256   endelse
257   
258   if NumVars eq 1 then begin
259     if VarName ne 'SLA' then OutObs2 = fltarr(n_elements(OutObs1)) + rmdi
260     OutMod2 = fltarr(n_elements(OutMod1)) + rmdi
261     OutQC2  = fltarr(n_elements(OutQC1)) + rmdi
262   endif
263;
264   ifile2=ifile2+1
265   endif
266  ncdf_close, ncid
267 
268  endif
269 
270endfor ;ifile
271
272pts1 = where(OutQC1 eq 1, count1)
273pts2 = where(OutQC1 ne 1, count2)
274if count1 gt 0 then OutQc1(pts1) = 0
275if count2 gt 0 then OutQc1(pts2) = 1
276if NumVars gt 1 then begin
277  pts1 = where(OutQC2 eq 1, count1)
278  pts2 = where(OutQC2 ne 1, count2)
279  if count1 gt 0 then OutQc2(pts1) = 0
280  if count2 gt 0 then OutQc2(pts2) = 1
281endif
282
283END
Note: See TracBrowser for help on using the repository browser.