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_enact.pro in branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/TOOLS/OBSTOOLS/dataplot – NEMO

source: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/TOOLS/OBSTOOLS/dataplot/read_enact.pro @ 7152

Last change on this file since 7152 was 7152, checked in by jcastill, 7 years ago

Initial implementation of wave coupling branch - INGV wave branch + UKMO wave coupling branch
File size: 10.9 KB
Line 
1PRO read_enact, Files, DepRange=DepRange, $
2                NumData=NumData, $
3                OutLats=OutLats, OutLons=OutLons, Instruments=OutInst, $
4      Platform=OutPlatform, $
5                OutDeps=OutDeps, OutSObs=OutSObs, OutSQC=OutSQC, $
6                OutTObs=OutTObs, OutTQC=OutTQC, OutDates=OutDates, $
7                OutTMod=OutTMod, OutSMod=OutSMod, rmdi=rmdi, quiet=quiet, $
8                iobs=outiobs, jobs=outjobs, kobs=outkobs, $
9                ProfileNum=OutProfileNum
10
11;------------------------------------------------------------------------------------------
12;Program to read in data from an ENACT format file containing profile data.
13;
14; Inputs:
15;        File        => File name to be read in
16;        DepRange    => (optional) Range of depths for which the data is to be extracted.
17;
18; Outputs:
19;        NumProfs    => Number of profiles  -Integer
20;        Latitudes   => Latitudes(NumProfs) -Real
21;        Longitudes  => Longitudes(NumProfs) -Real
22;        Dates       => Dates(NumProfs)  -Date/time structure.
23;        Instruments => Instrument type (NumProfs) -String
24;        Depths      => Depths of each ob (NumProfs)  -Real
25;        Salinity    => Salinity (psu) (NumProfs)  -Real
26;        SalQC       => QC flag for salinity (NumProfs) 1=> Good  -Int
27;        Temperature => Temperature (deg C) (NumProfs)  -Real
28;        TempQC      => QC flag for temperature (NumProfs) 1=> Good  -Int
29;        ModelT      => Model temperature (deg C) (NumProfs)  -Real
30;        ModelS      => Model salinity (psu) (NumProfs)  -Real
31;        rmdi        => Missing data indicator  -Real
32;
33;Author: Matt Martin. 5/10/07.
34;------------------------------------------------------------------------------------------
35rmdi = 99999.
36if NOT keyword_set(DepRange) then DepRange=[0.,1.e1]   ;Default to top 10 metres.
37NumFiles=n_elements(Files)
38
39ifile2=0
40for ifile = 0, NumFiles-1 do begin
41  File = Files(ifile)
42 
43;1. Open the netcdf file
44  ncid = ncdf_open(File, /nowrite)
45  if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+File
46  if ncid ge 0 then begin
47  if not keyword_set(quiet) then print, 'Opened file '+File+' successfully'
48 
49  ;2. Get info about file
50  ncinfo = ncdf_inquire(ncid)
51;  if not keyword_set(quiet) then , ncinfo
52  ;3. Read in the relevant dimensions
53  for idim = 0, ncinfo.ndims-1 do begin
54    ncdf_diminq, ncid, idim, name, dimlen
55    if name eq 'N_PROF' then NumProfs = dimlen $
56    else if name eq 'N_LEVELS' then NumLevels = dimlen $
57    else if name eq 'STRING4' then str4 = dimlen
58  endfor
59 
60  ;4. Set up arrays and read in the data
61
62  if (NumProfs gt 0) then begin 
63  JulDays = fltarr(NumProfs)
64  Latitudes = fltarr(NumProfs)
65  Longitudes = fltarr(NumProfs)
66  BytInsts = intarr(4, NumProfs)
67  Instruments = strarr(NumProfs)
68  Platforms = strarr(NumProfs)
69  Instrumentsa = strarr(NumLevels, NumProfs)
70  Platformsa = strarr(NumLevels, NumProfs)
71  Depths = fltarr(NumLevels, NumProfs)
72  Salinity = fltarr(NumLevels, NumProfs)
73  SalQC = intarr(NumLevels, NumProfs)
74  Temperature = fltarr(NumLevels, NumProfs)
75  TempQC = intarr(NumLevels, NumProfs)
76  ModelT = fltarr(NumLevels, NumProfs)
77  ModelS = fltarr(NumLevels, NumProfs)
78  ProfileNum = replicate(1,NumLevels)#indgen(NumProfs)
79   
80  varid = ncdf_varid(ncid, 'JULD')
81  ncdf_varget, ncid, varid, JulDays
82 
83;  print,JulDays
84 
85;  BaseDate = var_to_dt(1950,1,1)
86;  BaseDate=JULDAY(1,1,1950,0,0)
87  BaseDate=double(JULDAY(1,1,1950,0,0))      ;should be at 0 UTC
88;  info,BaseDate
89  Dates = replicate(BaseDate, NumLevels,NumProfs)
90;  info,dates
91  for iprof = 0L, NumProfs-1 do begin
92    secs = JulDays(iprof)* 24. * 60. * 60.
93;    dt = dt_add(BaseDate, second=secs)
94    dt=BaseDate+JulDays(iprof)
95;    CALDAT, dt, mon, day, year, hour, minute, second
96;    print,dt-BaseDate, mon, day, year, hour, minute, second
97    Dates(0:NumLevels-1,iprof) = replicate(dt,NumLevels)
98  endfor
99
100 
101 
102  varid = ncdf_varid(ncid, 'LONGITUDE')
103  ncdf_varget, ncid, varid, Longitudes
104
105  varid = ncdf_varid(ncid, 'LATITUDE')
106  ncdf_varget, ncid, varid, Latitudes
107
108  varid = ncdf_varid(ncid, 'DEPH_CORRECTED')
109  ncdf_varget, ncid, varid, Depths
110  res = ncdf_attinq( ncid, varid , '_FillValue')
111  if res.length gt 0 and res.datatype ne "UNKNOWN" then ncdf_attget, ncid, varid, '_FillValue', FillVal $
112  else begin
113    res = ncdf_attinq( ncid, varid , '_fillvalue')
114    ncdf_attget, ncid, varid, '_fillvalue', FillVal
115  endelse
116  pts = where(Depths eq FillVal, count)
117  if count gt 0 then Depths(pts) = rmdi
118
119;  if keyword_set(Salinity) then begin
120    varid = ncdf_varid(ncid, 'PSAL_CORRECTED')
121    if (varid ne -1) then  begin
122;    if not keyword_set(quiet) then print,'reading psal_corrected'
123    ncdf_varget, ncid, varid, Salinity
124    res = ncdf_attinq( ncid, varid , '_FillValue')
125    if res.length gt 0 and res.datatype ne "UNKNOWN" then ncdf_attget, ncid, varid, '_FillValue', FillVal $
126    else begin
127      res = ncdf_attinq( ncid, varid , '_fillvalue')
128      ncdf_attget, ncid, varid, '_fillvalue', FillVal
129    endelse
130    pts = where(Salinity eq FillVal, count)
131    if count gt 0 then Salinity(pts) = rmdi
132   endif     
133    varid = ncdf_varid(ncid, 'PSAL_CORRECTED_QC')
134    ncdf_varget, ncid, varid, SalQC
135    res = ncdf_attinq( ncid, varid , '_FillValue')
136    if res.length gt 0 and res.datatype ne "UNKNOWN" then ncdf_attget, ncid, varid, '_FillValue', FillVal $
137    else begin
138      res = ncdf_attinq( ncid, varid , '_fillvalue')
139      ncdf_attget, ncid, varid, '_fillvalue', FillVal
140    endelse
141    SalQC = SalQC - 48
142;  endif
143
144;  if keyword_set(Temperature) then begin
145    varid = ncdf_varid(ncid, 'POTM_CORRECTED')
146    if (varid ne -1) then  begin
147;    if not keyword_set(quiet) then print,'reading potm_corrected'
148    ncdf_varget, ncid, varid, Temperature
149    res = ncdf_attinq( ncid, varid , '_FillValue')
150    if res.length gt 0 and res.datatype ne "UNKNOWN" then ncdf_attget, ncid, varid, '_FillValue', FillVal $
151    else begin
152      res = ncdf_attinq( ncid, varid , '_fillvalue')
153      ncdf_attget, ncid, varid, '_fillvalue', FillVal
154    endelse
155;    if not keyword_set(quiet) then print,'Temp fill value ',fillval
156    pts = where(Temperature eq FillVal, count)
157    if count gt 0 then Temperature(pts) = rmdi
158   endif   
159    varid = ncdf_varid(ncid, 'POTM_CORRECTED_QC')
160    ncdf_varget, ncid, varid, TempQC
161    TempQC = TempQC - 48
162;   endif
163
164;  if keyword_set(OutTMod) then begin
165    varid = ncdf_varid(ncid, 'POTM_Hx')
166    if (varid ne -1) then begin
167    ncdf_varget, ncid, varid, ModelT
168    res = ncdf_attinq( ncid, varid , '_FillValue')
169    if res.length gt 0 and res.datatype ne "UNKNOWN" then ncdf_attget, ncid, varid, '_FillValue', FillVal $
170    else begin
171      res = ncdf_attinq( ncid, varid , '_fillvalue')
172      ncdf_attget, ncid, varid, '_fillvalue', FillVal
173    endelse
174    pts = where(ModelT eq FillVal, count)
175    if count gt 0 then ModelT(pts) = rmdi
176    endif else begin
177   ModelT=rmdi
178    endelse
179;   endif
180
181;  if keyword_set(OutSMod) then begin
182    varid = ncdf_varid(ncid, 'PSAL_Hx')
183    if (varid ne -1) then begin
184    ncdf_varget, ncid, varid, ModelS
185    res = ncdf_attinq( ncid, varid , '_FillValue')
186    if res.length gt 0 and res.datatype ne "UNKNOWN" then ncdf_attget, ncid, varid, '_FillValue', FillVal $
187    else begin
188      res = ncdf_attinq( ncid, varid , '_fillvalue')
189      ncdf_attget, ncid, varid, '_fillvalue', FillVal
190    endelse
191    pts = where(ModelS eq FillVal, count)
192    if count gt 0 then ModelS(pts) = rmdi
193    endif else begin
194   ModelS=rmdi
195    endelse
196;  endif
197
198  varid = ncdf_varid(ncid, 'WMO_INST_TYPE')
199  ncdf_varget, ncid, varid, BytInsts
200 
201  pts = where(BytInsts(0,*) eq 32 and $
202              BytInsts(1,*) eq 56 and $
203              BytInsts(2,*) eq 50 and $
204              BytInsts(3,*) eq 48, count)
205  if count gt 0 then Instruments(pts) = 'BUOYS'  ; 820
206 
207  pts = where(BytInsts(0,*) eq 32 and $
208              BytInsts(1,*) eq 52 and $
209              BytInsts(2,*) eq 48 and $
210              BytInsts(3,*) eq 49, count)
211  if count gt 0 then Instruments(pts) = 'XBT'    ; 401
212 
213  pts = where(BytInsts(0,*) eq 32 and $
214              BytInsts(1,*) eq 55 and $
215              BytInsts(2,*) eq 52 and $
216              BytInsts(3,*) eq 49, count)
217  if count gt 0 then Instruments(pts) = 'TESAC'  ; 741
218 
219  pts = where(BytInsts(0,*) eq 32 and $
220              BytInsts(1,*) eq 56 and $
221              BytInsts(2,*) eq 51 and $
222              BytInsts(3,*) eq 49, count)
223  if count gt 0 then Instruments(pts) = 'ARGO'   ; 831
224 
225  pts = where(Instruments eq '', count)
226  if count gt 0 then Instruments(pts) = 'UNKNOWN'
227 
228  varid = ncdf_varid(ncid, 'PLATFORM_NUMBER')
229  ncdf_varget, ncid, varid, platforms
230  platforms=string(platforms)
231 
232  for i=0,numlevels-1 do begin
233    platformsa(i,*)=platforms
234    instrumentsa(i,*)=instruments
235  endfor
236
237  varid = ncdf_varid(ncid, 'IOBS')
238  if (varid ne -1) then ncdf_varget, ncid, varid, iobsval
239
240  varid = ncdf_varid(ncid, 'JOBS')
241  if (varid ne -1) then ncdf_varget, ncid, varid, jobsval
242
243  varid = ncdf_varid(ncid, 'KOBS')
244  if (varid ne -1) then ncdf_varget, ncid, varid, kobsval
245   
246  ;Select those obs in the required depth range
247;  Lats = replv(Latitudes(*), [NumLevels,NumProfs],1)
248;  Lons = replv(Longitudes(*), [NumLevels,NumProfs],1)
249
250   Lons=fltarr(NumLevels,NumProfs)
251        Lats=fltarr(NumLevels,NumProfs)
252   iobs=fltarr(NumLevels,NumProfs)
253        jobs=fltarr(NumLevels,NumProfs)
254       
255   for i=0L,NumProfs-1 do begin
256          Lats(*,i)=Latitudes(i)
257          Lons(*,i)=Longitudes(i)
258          if (n_elements(iobsval) gt 0) then $
259             iobs(*,i)=iobsval(i)
260          if (n_elements(jobsval) gt 0) then $
261              jobs(*,i)=jobsval(i)
262        endfor 
263 
264  pts = where(Depths ge DepRange(0) and Depths le DepRange(1), NumData)
265 
266  if ifile2 eq 0 then begin 
267    OutTObs = [Temperature(pts)]
268    OutSObs = [Salinity(pts)]
269    OutTMod = [ModelT(pts)]
270    OutSMod = [ModelS(pts)]
271    OutTQC  = [TempQC(pts)]
272    OutSQC  = [SalQC(pts)]
273    OutDeps = [Depths(pts)]
274    OutLats = [Lats(pts)]
275    OutLons = [Lons(pts)]
276    OutDates= [Dates(pts)]
277    OutInst= [Instrumentsa(pts)]   
278    OutPlatform= [Platformsa(pts)]
279    OutProfileNum = [ProfileNum(pts)]
280    if (n_elements(iobsval) gt 0) then outiobs = [iobs(pts)]
281    if (n_elements(jobsval) gt 0) then outjobs = [jobs(pts)]
282    if (n_elements(kobsval) gt 0) then outkobs = [kobsval(pts)]
283   
284  endif else begin
285    OutTObs = [OutTObs,Temperature(pts)]
286    OutSObs = [OutSObs,Salinity(pts)]
287    OutTMod = [OutTMod,ModelT(pts)]
288    OutSMod = [OutSMod,ModelS(pts)]
289    OutTQC  = [OutTQC,TempQC(pts)]
290    OutSQC  = [OutSQC,SalQC(pts)]
291    OutDeps = [OutDeps,Depths(pts)]
292    OutLats = [OutLats,Lats(pts)]
293    OutLons = [OutLons,Lons(pts)]
294    OutDates= [OutDates,Dates(pts)]
295    OutInst= [OutInst,Instrumentsa(pts)]
296    OutPlatform = [OutPlatform, Platformsa(pts)]
297    OutProfileNum = [OutProfileNum, ProfileNum(pts)]
298    if (n_elements(iobsval) gt 0) then outiobs = [outiobs, iobs(pts)]
299    if (n_elements(jobsval) gt 0) then outjobs = [outjobs, jobs(pts)]
300    if (n_elements(kobsval) gt 0) then outkobs = [outkobs, kobsval(pts)]
301
302  endelse
303
304   ifile2=ifile2+1
305   endif
306  ncdf_close, ncid
307
308
309   endif
310 
311endfor ;ifile
312
313END
Note: See TracBrowser for help on using the repository browser.