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

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

Update documentation for obstools and dataplot. Removal of dataplot code not needed. Addition of headers to some dataplot code. Addition of .exe to command example in obstools.

File size: 7.6 KB
Line 
1;+---------------------------------------------------------------------------
2PRO read_cdfobs, Files, NumObs=NumObs, $
3               Latitudes=Latitudes, Longitudes=Longitudes, Depths=Depths, $
4               Obs=Observations, ModelVals=ModelVals, qcs=QCs, $
5               Ob2=Observations2, ModelVal2=ModelVals2, qc2=QCs2, $
6               Ob3=Observations3, $
7               Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, kobs=kobs, $
8               Salinity=Salinity, nodates=nodates, types=types, $
9               filetype=ObsType, quiet=quiet, MDT=MDT, $
10               ProfileNum=ProfileNum, error=error, $
11               notfussy=notfussy,VarName=VarName 
12;+--------------------------------------------------------------------------
13; Read in observation and feedback files
14; detects filetype and calls the appropriate reading routine
15;
16; Author:  D. J. Lea       Feb 2008
17;+--------------------------------------------------------------------------
18
19; Declare error label.
20ON_IOERROR, IOERROR
21         
22error=0
23title=''
24
25; Set types to undefined
26types=-1 & tempvar=size(temporary(types))
27
28; Read in netcdf data file and observation operator
29;2. Work out which type of data is in the files by looking at the first one.
30ncid = ncdf_open(Files(0), /nowrite)
31if ncid lt 0 then message, 'Error opening file '+File
32result = NCDF_ATTINQ( ncid, 'title', /global)
33if result.datatype ne 'UNKNOWN' then ncdf_attget, ncid, 'title', Title, /global
34
35if string(Title) eq "NEMO observation operator output" then ObsType='feedback' $
36else $
37if string(Title) eq "Forecast class 4 file" then ObsType='ForecastClass4' $
38else $
39ObsType = 'none'
40
41if ObsType eq 'none' then begin
42  varid = ncdf_varid(ncid, 'POTM_CORRECTED') 
43  if varid ge 0 then ObsType = 'Prof'
44 
45  varid = ncdf_varid(ncid, 'SLA')
46  if varid ge 0 then ObsType = 'SSH'
47 
48  varid = ncdf_varid(ncid, 'SST')
49  varid2 = ncdf_varid(ncid, 'sea_surface_temperature')
50  varid=max([varid,varid2])
51  if varid ge 0 then ObsType = 'SST'
52 
53  varid = ncdf_varid(ncid, 'SEAICE')
54  varid2 = ncdf_varid(ncid, 'sea_ice_concentration')
55  varid=max([varid,varid2])
56  if varid ge 0 then ObsType = 'SEAICE'   
57 
58  varid = ncdf_varid(ncid, 'LOGCHL')
59  varid2 = ncdf_varid(ncid, 'LogChl')
60  varid3 = ncdf_varid(ncid, 'CHL1_mean')
61  varid=max([varid,varid2,varid3])
62  if varid ge 0 then ObsType = 'LOGCHL'   
63 
64  varid = ncdf_varid(ncid, 'PCO2')
65  if varid ge 0 then ObsType = 'PCO2'
66
67  varid = ncdf_varid(ncid, 'UCRT')
68  varid2 = ncdf_varid(ncid, 'VCRT')
69  varid=max([varid,varid2])
70  if varid ge 0 then ObsType = 'CRT'   
71endif
72 
73ncdf_close, ncid
74
75if not keyword_set(quiet) then print, 'Reading in data as type ',ObsType
76
77  if ObsType eq 'feedback' then begin
78 
79    if n_elements(DepRange) eq 0 then DepRange=[0,5000] 
80 
81    read_feedback, Files, DepRange=DepRange, VarName=VarName, NumData=numobs, $
82                OutLats=Latitudes, OutLons=Longitudes, $
83                OutInstruments=Instruments, OutPlatform=Platform, $
84                OutDeps=Depths,  $
85                OutObs1=Observations, OutObs2=Observations2, $
86                OutMod1=ModelVals, OutMod2=ModelVals2, $
87                OutQC1=QCs, OutQC2=QCs2, $
88                MDT=MDT, OutDates=Dates, $
89                rmdi=rmdi, quiet=quiet, $
90                OutProfileNum=ProfileNum
91               
92;                info, instruments
93;                info, platform
94               
95                types=Platform
96                if n_elements(instruments) gt 0 then types=Platform+' '+Instruments
97               
98  print, 'Varname: ',VarName
99  endif else if ObsType eq 'Prof' then begin
100 
101    if keyword_set(PlotModel) then begin
102      OutTMod = 1
103      OutSMod = 1
104    endif
105
106   if n_elements(DepRange) eq 0 then DepRange=[0,5000]
107   
108  read_enact, Files, DepRange=DepRange, NumData=NumData, $
109                OutLats=Latitudes, OutLons=Longitudes, $
110                Instruments=Instruments, Platform=Platform, $
111                OutDeps=Depths, OutSObs=OutSObs, OutSQC=OutSQC, $
112                OutTObs=OutTObs, OutTQC=OutTQC, OutDates=Dates, $
113                OutTMod=OutTMod, OutSMod=OutSMod, rmdi=rmdi, quiet=quiet, $
114                iobs=iobs, jobs=jobs, kobs=kobs, $
115                ProfileNum=ProfileNum
116
117   types=Instruments+' '+Platform
118
119; if salinity keyword is set then put salinity values first 
120    if keyword_set(Salinity) then begin     
121      Observations = OutSObs
122      ModelVals    = OutSMod   
123      QCs          = OutSQC
124      Observations2 = OutTObs
125      ModelVals2    = OutTMod
126      QCs2          = OutTQC
127    endif else begin
128      Observations = OutTObs
129      ModelVals    = OutTMod
130      QCs          = OutTQC
131      Observations2 = OutSObs
132      ModelVals2    = OutSMod
133      QCs2          = OutSQC
134    endelse
135
136   numobs=numdata
137        numobs=n_elements(latitudes)
138                 
139  endif else if ObsType eq 'SSH' then begin
140   
141    read_sla, Files, NumObs=NumObs, $
142              Latitudes=Latitudes, Longitudes=Longitudes, $
143              ObsSLA=Observations, ModSLA=ModelVals, SLAQC=QCs, $
144              MDT=MDT, Satellites=Satellites, types=types, Dates=Dates, rmdi=rmdi, $
145         iobs=iobs, jobs=jobs, mstp=mstp, quiet=quiet
146    Depths = fltarr(NumObs)
147   
148  endif else if ObsType eq 'SST' then begin
149     
150    read_sst, Files, NumObs=NumObs, $
151              Latitudes=Latitudes, Longitudes=Longitudes, $
152              ObsSST=Observations, ModSST=ModelVals, SSTQC=QCs, $
153              Dates=Dates, rmdi=rmdi, types=types, iobs=iobs, jobs=jobs, $
154              quiet=quiet
155    Depths = fltarr(NumObs)
156
157  endif else if ObsType eq 'SEAICE' then begin
158 
159    read_seaice, Files, NumObs=NumObs, $
160              Latitudes=Latitudes, Longitudes=Longitudes, $
161              Obs=Observations, Modarr=ModelVals, QC=QCs, $
162              Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
163              quiet=quiet
164    Depths = fltarr(NumObs)
165 
166  endif else if ObsType eq 'LOGCHL' then begin
167     
168    read_chl, Files, NumObs=NumObs, $
169              Latitudes=Latitudes, Longitudes=Longitudes, $
170              Obs=Observations, Modarr=ModelVals, QC=QCs, $
171              Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
172              quiet=quiet
173    Depths = fltarr(NumObs)
174             
175   endif else if ObsType eq 'PCO2' then begin
176     
177    read_pco2, Files, NumObs=NumObs, $
178              Latitudes=Latitudes, Longitudes=Longitudes, $
179              Obs=Observations, Modarr=ModelVals, QC=QCs, $
180              Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
181              quiet=quiet
182    Depths = fltarr(NumObs)
183
184   endif else if ObsType eq 'CRT' then begin
185
186    read_crt, Files, NumObs=NumObs, $
187              OutLats=Latitudes, OutLons=Longitudes, $
188              OutTypes=Types, OutDates=Dates, $
189              OutUObs=Observations, OutVObs=Observations2, $
190              OutUMod=ModelVals, OutVMod=ModelVals2, Quiet=Quiet, $
191              FloatNum=FloatNum, QC=QCs, OutSpeed=Observations3
192    rmdi=0.0
193    QCs2=QCs
194    Depths = fltarr(NumObs)
195
196   endif else if ObsType eq 'ForecastClass4' then begin
197   
198    print, 'reading forecast class 4 files'
199    read_forc_class4, Files, NumObs=NumObs, $
200              Latitudes=Latitudes, Longitudes=Longitudes, $
201              Depths=Depths, $
202              Types=types, Dates=Dates, $
203              Obsarr=Observations, Obs2arr=Observations2, $
204              Modarr=ModelVals, Mod2arr=ModelVals2, $
205              QCs1=QCs, QCs2=QCs2, rmdi=rmdi, $
206              notfussy=notfussy
207                 
208  endif else message, 'Error: ObsType is not set correctly'                 
209
210  if (n_elements(types) eq 0) then begin
211     types=replicate(rmdi,NumObs)
212  endif
213
214  goto, NOERROR
215 
216IOERROR: error=1
217    print,'read_cdfobs: an error occurred trying read files: ',files
218
219NOERROR:
220   
221END
Note: See TracBrowser for help on using the repository browser.