Context navigation

Back to Ticket #431

Ticket #431: ropp_io_read_ncdf_get.f90.30062016

File ropp_io_read_ncdf_get.f90.30062016, 124.3 KB (added by Ian Culverwell, 8 years ago)
ropp_io_read_ncdf_get.f90.30062016

Line
1	! $Id: ropp_io_read_ncdf_get.f90 4461 2015-02-04 15:11:06Z idculv $
2
3	!***is Reading/ropp_io_read_ncdf_get *
4	!
5	! NAME
6	! ropp_io_read_ncdf_get - Get data from an (already defined) netCDF.
7	!
8	! SYNOPSIS
9	! call ropp_io_read_ncdf_get(data, rec)
10	!
11	! DESCRIPTION
12	! This subroutine gets variables from an already openend / netCDF data
13	! file.
14	! Reads core parameters from V1.0 format version and new parameters
15	! added in V1.1. Additional variables contained in the netCDF file are
16	! read in and define elements of the Vlist structures in ROprof.
17	!
18	! NOTES
19	! A netCDF file must have been opened before; this subroutine only works
20	! on the current netcdf file. The netCDF data file is left open.
21	!
22	! AUTHOR
23	! Met Office, Exeter, UK.
24	! Any comments on this software should be given via the ROM SAF
25	! Helpdesk at http://www.romsaf.org
26	!
27	! COPYRIGHT
28	! (c) EUMETSAT. All rights reserved.
29	! For further details please refer to the file COPYRIGHT
30	! which you should have received as part of this distribution.
31	!
32	!****
33
34	!-------------------------------------------------------------------------------
35	! 1. Core RO data
36	!-------------------------------------------------------------------------------
37
38	SUBROUTINE ropp_io_read_ncdf_get_rodata(DATA, rec)
39
40	! 1.1 Declarations
41	! ----------------
42
43	USE ropp_utils
44	USE ncdf
45	USE ropp_io, not_this => ropp_io_read_ncdf_get_rodata
46	USE ropp_io_types, ONLY: ROprof, &
47	ThisFmtVer
48
49	IMPLICIT NONE
50
51	TYPE(ROprof), INTENT(inout) :: DATA
52	INTEGER, OPTIONAL :: rec
53
54	INTEGER :: n, ierr
55	INTEGER :: irec
56	INTEGER, DIMENSION(8) :: DT8 ! Date/time array
57
58	CHARACTER(len = 23) :: proc_date
59	CHARACTER(LEN=15) :: cval
60
61	REAL(dp) :: time, start_time, dtocc_time, now_time !dp defined in DateTimeTypes
62	REAL :: fmtver
63
64	INTEGER :: status, varid, ndim, TYPE
65	CHARACTER(len = 256) :: routine
66
67	! 1.2 Error handling
68	! ------------------
69
70	CALL message_get_routine(routine)
71	CALL message_set_routine('ropp_io_read_ncdf_get')
72
73	! 1.3 Default parameters
74	! ----------------------
75
76	IF (PRESENT(rec)) THEN
77	irec = rec
78	ELSE
79	irec = 1
80	ENDIF
81
82	! 1.4 (Global) Attributes
83	! ------------------------
84
85	data%FmtVersion = ' ' ; CALL ncdf_getatt('format_version', data%FmtVersion)
86	READ ( data%FmtVersion(11:), fmt=*, iostat=ierr ) fmtver
87	IF ( ierr /= 0 ) data%FmtVersion = ThisFmtVer
88	data%processing_centre = ' ' ; CALL ncdf_getatt('processing_centre', data%processing_centre)
89	IF (ncdf_isatt('processing_software')) THEN ! added at V8.0
90	data%processing_software = ' ' ; CALL ncdf_getatt('processing_software', data%processing_software)
91	ENDIF
92	proc_date = ' ' ; CALL ncdf_getatt('processing_date', proc_date)
93	data%pod_method = ' ' ; CALL ncdf_getatt('pod_method', data%pod_method)
94	data%phase_method = ' ' ; CALL ncdf_getatt('phase_method', data%phase_method)
95	data%bangle_method = ' ' ; CALL ncdf_getatt('bangle_method', data%bangle_method)
96	data%refrac_method = ' ' ; CALL ncdf_getatt('refrac_method', data%refrac_method)
97	data%meteo_method = ' ' ; CALL ncdf_getatt('meteo_method', data%meteo_method)
98	IF(ncdf_isatt('thin_method'))THEN ! added at V1.1
99	data%thin_method = ' ' ; CALL ncdf_getatt('thin_method', data%thin_method)
100	ENDIF
101	data%software_version = ' ' ; CALL ncdf_getatt('software_version', data%software_version)
102
103	IF (proc_date( 1: 4) /= ' ') READ(proc_date( 1: 4), *) data%DTpro%Year
104	IF (proc_date( 6: 7) /= ' ') READ(proc_date( 6: 7), *) data%DTpro%Month
105	IF (proc_date( 9:10) /= ' ') READ(proc_date( 9:10), *) data%DTpro%Day
106	IF (proc_date(12:13) /= ' ') READ(proc_date(12:13), *) data%DTpro%Hour
107	IF (proc_date(15:16) /= ' ') READ(proc_date(15:16), *) data%DTpro%Minute
108	IF (proc_date(18:19) /= ' ') READ(proc_date(18:19), *) data%DTpro%Second
109	IF (proc_date(21:23) /= ' ') READ(proc_date(21:23), *) data%DTpro%Msec
110
111	! 1.5 Header variables
112	! --------------------
113
114	CALL ncdf_getvar('occ_id', data%occ_id, rec = irec)
115	CALL ncdf_getvar('gns_id', data%gns_id, rec = irec)
116	CALL ncdf_getvar('leo_id', data%leo_id, rec = irec)
117	CALL ncdf_getvar('stn_id', data%stn_id, rec = irec)
118
119	! 1.6 Date and time
120	! -----------------
121
122	CALL ncdf_getvar('start_time', start_time, rec=irec)
123	CALL ncdf_getvar('year', data%DTocc%Year, &
124	units = data%DTocc%units%Year, &
125	range = data%DTocc%range%Year, &
126	rec = irec)
127	CALL ncdf_getvar('month', data%DTocc%Month, &
128	units = data%DTocc%units%Month, &
129	range = data%DTocc%range%Month, &
130	rec = irec)
131	CALL ncdf_getvar('day', data%DTocc%Day, &
132	units = data%DTocc%units%Day, &
133	range = data%DTocc%range%Day, &
134	rec = irec)
135	CALL ncdf_getvar('hour', data%DTocc%Hour, &
136	units = data%DTocc%units%Hour, &
137	range = data%DTocc%range%Hour, &
138	rec = irec)
139	CALL ncdf_getvar('minute', data%DTocc%Minute, &
140	units = data%DTocc%units%Minute, &
141	range = data%DTocc%range%Minute, &
142	rec = irec)
143	CALL ncdf_getvar('second', data%DTocc%Second, &
144	units = data%DTocc%units%Second, &
145	range = data%DTocc%range%Second, &
146	rec = irec)
147	CALL ncdf_getvar('msec', data%DTocc%Msec, &
148	units = data%DTocc%units%Msec, &
149	range = data%DTocc%range%Msec, &
150	rec = irec)
151
152	! 1.6.1 Check consistency: start_time and DTocc (if both are valid)
153	! should refer to the same epoch within 1ms. Issue a warning if not.
154
155	CALL Date_and_Time_UTC(Values=DT8)
156
157	CALL TimeSince(DT8, now_time, 1, Base="JS2000")
158
159	IF (isroppinrange(data%DTocc)) THEN
160	IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN
161	DT8 = (/ data%DTocc%Year, data%DTocc%Month, data%DTocc%Day, 0, &
162	data%DTocc%Hour, data%DTocc%Minute, data%DTocc%Second, &
163	data%DTocc%Msec /)
164	CALL TimeSince(DT8, dtocc_time, 1, Base="JS2000")
165	time = ABS(start_time - dtocc_time)
166	IF ( time > 0.0005_dp ) THEN
167	WRITE ( cval, FMT="(F15.3)" ) time
168	CALL message(msg_warn, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // &
169	"timestamps differ by " // &
170	TRIM(ADJUSTL(cval))//" seconds - using yr/../ms timestamp")
171	END IF
172	END IF
173
174	! If any DTocc element is invalid, substitute converted start_time
175	! (if that is valid) and issue a warning.
176	ELSE
177
178	IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN
179	CALL message(msg_warn, "One or more of yr/mo/dy/hr/mn/sc/ms times " // &
180	"are invalid - using 'start_time' instead")
181	CALL TimeSince(DT8, start_time, -1, Base="JS2000")
182	data%DTocc%Year = DT8(1)
183	data%DTocc%Month = DT8(2)
184	data%DTocc%Day = DT8(3)
185	data%DTocc%Hour = DT8(5)
186	data%DTocc%Minute = DT8(6)
187	data%DTocc%Second = DT8(7)
188	data%DTocc%Msec = DT8(8)
189	END IF
190
191	! NB if neither DTocc nor start_time are valid, do nothing here; missing
192	! coordinates should be dealt with as part of generic Q/C.
193	END IF
194
195	! 1.7 Overall quality
196	! -------------------
197
198	CALL ncdf_getvar('pcd', data%pcd, &
199	units = data%units%pcd, &
200	range = data%range%pcd, &
201	rec = irec)
202	CALL ncdf_getvar('overall_qual', data%overall_qual, &
203	units = data%units%overall_qual, &
204	range = data%range%overall_qual, &
205	rec = irec)
206
207	! 1.8 Georeferencing
208	! ------------------
209
210	CALL ncdf_getvar('time', time, rec=irec)
211	CALL ncdf_getvar('lat', data%georef%lat, &
212	units = data%georef%units%lat, &
213	range = data%georef%range%lat, &
214	rec = irec)
215	CALL ncdf_getvar('lon', data%georef%lon, &
216	units = data%georef%units%lon, &
217	range = data%georef%range%lon, &
218	rec = irec)
219	CALL ncdf_getvar('time_offset', data%georef%time_offset, &
220	units = data%georef%units%time_offset, &
221	range = data%georef%range%time_offset, &
222	rec = irec)
223	CALL ncdf_getvar('undulation', data%georef%Undulation, &
224	units = data%georef%units%Undulation, &
225	range = data%georef%range%undulation, &
226	rec = irec)
227	CALL ncdf_getvar('roc', data%georef%roc, &
228	units = data%georef%units%roc, &
229	range = data%georef%range%roc, &
230	rec = irec)
231	CALL ncdf_getvar('r_coc', data%georef%r_coc, &
232	units = data%georef%units%r_coc, &
233	range = data%georef%range%r_coc, &
234	rec = irec)
235	CALL ncdf_getvar('azimuth', data%georef%azimuth, &
236	units = data%georef%units%azimuth, &
237	range = data%georef%range%azimuth, &
238	rec = irec)
239
240	! 1.8.1 Other attributes
241
242	CALL ncdf_getatt('reference_frame', data%georef%reference_frame%r_coc, varname= 'r_coc')
243
244	! 1.9 Background characterisation (if any)
245	! ----------------------------------------
246
247	IF (ncdf_isvar('bg_source')) THEN
248	data%BG%Source = 'TBD'
249	ELSE
250	data%BG%Source = 'NONE'
251	ENDIF
252
253	IF (data%BG%Source /= 'NONE') THEN
254	CALL ncdf_getvar('bg_source', data%BG%Source, rec = irec)
255	CALL ncdf_getvar('bg_year', data%BG%Year, &
256	units = data%BG%units%Year, &
257	range = data%BG%range%Year, &
258	rec = irec)
259	CALL ncdf_getvar('bg_month', data%BG%Month, &
260	units = data%BG%units%Month, &
261	range = data%BG%range%Month, &
262	rec = irec)
263	CALL ncdf_getvar('bg_day', data%BG%Day, &
264	units = data%BG%units%Day, &
265	range = data%BG%range%Day, &
266	rec = irec)
267	CALL ncdf_getvar('bg_hour', data%BG%Hour, &
268	units = data%BG%units%Hour, &
269	range = data%BG%range%Hour, &
270	rec = irec)
271	CALL ncdf_getvar('bg_minute', data%BG%Minute, &
272	units = data%BG%units%Minute, &
273	range = data%BG%range%Minute, &
274	rec = irec)
275	CALL ncdf_getvar('bg_fcperiod', data%BG%fcPeriod, &
276	units = data%BG%units%fcPeriod, &
277	range = data%BG%range%fcPeriod, &
278	rec = irec)
279	ENDIF
280
281	! 1.10 Level1a variables (if any)
282	! ------------------------------
283
284	IF (ncdf_isvar('dtime')) THEN
285	CALL ncdf_getsize('dtime', n, dim = 1)
286	CALL ropp_io_init(data%Lev1a, n)
287	ELSE
288	data%Lev1a%Npoints = 0
289	ENDIF
290
291	IF (data%Lev1a%Npoints > 0) THEN
292	CALL ncdf_getvar('dtime', data%Lev1a%dtime, &
293	units = data%Lev1a%units%dtime, &
294	range = data%Lev1a%range%dtime, &
295	rec = irec)
296	CALL ncdf_getvar('snr_L1ca', data%Lev1a%snr_L1ca, &
297	units = data%Lev1a%units%snr, &
298	range = data%Lev1a%range%snr, &
299	rec = irec)
300	CALL ncdf_getvar('snr_L1p', data%Lev1a%snr_L1p, &
301	units = data%Lev1a%units%snr, &
302	range = data%Lev1a%range%snr, &
303	rec = irec)
304	CALL ncdf_getvar('snr_L2p', data%Lev1a%snr_L2p, &
305	units = data%Lev1a%units%snr, &
306	range = data%Lev1a%range%snr, &
307	rec = irec)
308	CALL ncdf_getvar('phase_L1', data%Lev1a%phase_L1, &
309	units = data%Lev1a%units%phase, &
310	range = data%Lev1a%range%phase, &
311	rec = irec)
312	CALL ncdf_getvar('phase_L2', data%Lev1a%phase_L2, &
313	units = data%Lev1a%units%phase, &
314	range = data%Lev1a%range%phase, &
315	rec = irec)
316	CALL ncdf_getvar('r_gns', data%Lev1a%r_gns, &
317	units = data%Lev1a%units%r_gns, &
318	range = data%Lev1a%range%r_gns, &
319	rec = irec)
320	CALL ncdf_getvar('v_gns', data%Lev1a%v_gns, &
321	units = data%Lev1a%units%v_gns, &
322	range = data%Lev1a%range%v_gns, &
323	rec = irec)
324	CALL ncdf_getvar('r_leo', data%Lev1a%r_leo, &
325	units = data%Lev1a%units%r_leo, &
326	range = data%Lev1a%range%r_leo, &
327	rec = irec)
328	CALL ncdf_getvar('v_leo', data%Lev1a%v_leo, &
329	units = data%Lev1a%units%v_leo, &
330	range = data%Lev1a%range%v_leo, &
331	rec = irec)
332	CALL ncdf_getvar('phase_qual', data%Lev1a%phase_qual, &
333	units = data%Lev1a%units%phase_qual, &
334	range = data%Lev1a%range%phase_qual, &
335	rec = irec)
336	! 1.10.1 Other attributes
337
338	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_gns, varname = 'r_gns')
339	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_gns, varname = 'v_gns')
340	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_leo, varname = 'r_leo')
341	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_leo, varname = 'v_leo')
342
343	data%Lev1a%Missing = .FALSE.
344	ENDIF
345
346	! 1.11 Level1b variables (if any)
347	! -------------------------------
348
349	IF (ncdf_isvar('lat_tp')) THEN
350	CALL ncdf_getsize('lat_tp', n, dim = 1)
351	CALL ropp_io_init(data%Lev1b, n)
352	ELSE
353	data%Lev1b%Npoints = 0
354	ENDIF
355
356	IF (data%Lev1b%Npoints > 0) THEN
357
358	CALL ncdf_getvar('lat_tp', data%Lev1b%lat_tp, &
359	units = data%Lev1b%units%lat_tp, &
360	range = data%Lev1b%range%lat_tp, &
361	rec = irec)
362	CALL ncdf_getvar('lon_tp', data%Lev1b%lon_tp, &
363	units = data%Lev1b%units%lon_tp, &
364	range = data%Lev1b%range%lon_tp, &
365	rec = irec)
366	CALL ncdf_getvar('azimuth_tp', data%Lev1b%azimuth_tp, &
367	units = data%Lev1b%units%azimuth_tp, &
368	range = data%Lev1b%range%azimuth_tp, &
369	rec = irec)
370
371	CALL ncdf_getvar('impact_L1', data%Lev1b%impact_L1, &
372	units = data%Lev1b%units%impact, &
373	range = data%Lev1b%range%impact, &
374	rec = irec)
375	CALL ncdf_getvar('impact_L2', data%Lev1b%impact_L2, &
376	units = data%Lev1b%units%impact, &
377	range = data%Lev1b%range%impact, &
378	rec = irec)
379	CALL ncdf_getvar('impact', data%Lev1b%impact, &
380	units = data%Lev1b%units%impact, &
381	range = data%Lev1b%range%impact, &
382	rec = irec)
383	IF (ncdf_isvar('impact_opt')) & ! added at v1.1
384	CALL ncdf_getvar('impact_opt', data%Lev1b%impact_opt, &
385	units = data%Lev1b%units%impact, &
386	range = data%Lev1b%range%impact, &
387	rec = irec)
388
389	CALL ncdf_getvar('bangle_L1', data%Lev1b%bangle_L1, &
390	units = data%Lev1b%units%bangle, &
391	range = data%Lev1b%range%bangle, &
392	rec = irec)
393	CALL ncdf_getvar('bangle_L2', data%Lev1b%bangle_L2, &
394	units = data%Lev1b%units%bangle, &
395	range = data%Lev1b%range%bangle, &
396	rec = irec)
397	CALL ncdf_getvar('bangle', data%Lev1b%bangle, &
398	units = data%Lev1b%units%bangle, &
399	range = data%Lev1b%range%bangle, &
400	rec = irec)
401	IF (ncdf_isvar('bangle_opt')) & ! added at v1.1
402	CALL ncdf_getvar('bangle_opt', data%Lev1b%bangle_opt, &
403	units = data%Lev1b%units%bangle, &
404	range = data%Lev1b%range%bangle, &
405	rec = irec)
406
407	CALL ncdf_getvar('bangle_L1_sigma', data%Lev1b%bangle_L1_sigma, &
408	units = data%Lev1b%units%bangle_sigma, &
409	range = data%Lev1b%range%bangle_sigma, &
410	rec = irec)
411	CALL ncdf_getvar('bangle_L2_sigma', data%Lev1b%bangle_L2_sigma, &
412	units = data%Lev1b%units%bangle_sigma, &
413	range = data%Lev1b%range%bangle_sigma, &
414	rec = irec)
415	CALL ncdf_getvar('bangle_sigma', data%Lev1b%bangle_sigma, &
416	units = data%Lev1b%units%bangle_sigma, &
417	range = data%Lev1b%range%bangle_sigma, &
418	rec = irec)
419	IF (ncdf_isvar('bangle_opt_sigma')) & ! added at v1.1
420	CALL ncdf_getvar('bangle_opt_sigma', data%Lev1b%bangle_opt_sigma, &
421	units = data%Lev1b%units%bangle_sigma, &
422	range = data%Lev1b%range%bangle_sigma, &
423	rec = irec)
424
425	CALL ncdf_getvar('bangle_L1_qual', data%Lev1b%bangle_L1_qual, &
426	units = data%Lev1b%units%bangle_qual, &
427	range = data%Lev1b%range%bangle_qual, &
428	rec = irec)
429	CALL ncdf_getvar('bangle_L2_qual', data%Lev1b%bangle_L2_qual, &
430	units = data%Lev1b%units%bangle_qual, &
431	range = data%Lev1b%range%bangle_qual, &
432	rec = irec)
433	CALL ncdf_getvar('bangle_qual', data%Lev1b%bangle_qual, &
434	units = data%Lev1b%units%bangle_qual, &
435	range = data%Lev1b%range%bangle_qual, &
436	rec = irec)
437	IF (ncdf_isvar('bangle_opt_qual')) & ! added at v1.1
438	CALL ncdf_getvar('bangle_opt_qual', data%Lev1b%bangle_opt_qual, &
439	units = data%Lev1b%units%bangle_qual, &
440	range = data%Lev1b%range%bangle_qual, &
441	rec = irec)
442	data%Lev1b%Missing = .FALSE.
443	ENDIF
444
445	! 1.12 Level2a variables (if any)
446	! -------------------------------
447
448	IF (ncdf_isvar('alt_refrac')) THEN
449	CALL ncdf_getsize('alt_refrac', n, dim = 1)
450	CALL ropp_io_init(data%Lev2a, n)
451	ELSE
452	data%Lev2a%Npoints = 0
453	ENDIF
454
455	IF (data%Lev2a%Npoints > 0) THEN
456
457	CALL ncdf_getvar('alt_refrac', data%Lev2a%alt_refrac, &
458	units = data%Lev2a%units%alt_refrac, &
459	range = data%Lev2a%range%alt_refrac, &
460	rec = irec)
461	CALL ncdf_getvar('geop_refrac', data%Lev2a%geop_refrac, &
462	units = data%Lev2a%units%geop_refrac, &
463	range = data%Lev2a%range%geop_refrac, &
464	rec = irec)
465	CALL ncdf_getvar('refrac', data%Lev2a%refrac, &
466	units = data%Lev2a%units%refrac, &
467	range = data%Lev2a%range%refrac, &
468	rec = irec)
469	CALL ncdf_getvar('refrac_sigma', data%Lev2a%refrac_sigma, &
470	units = data%Lev2a%units%refrac_sigma, &
471	range = data%Lev2a%range%refrac_sigma, &
472	rec = irec)
473	CALL ncdf_getvar('refrac_qual', data%Lev2a%refrac_qual, &
474	units = data%Lev2a%units%refrac_qual, &
475	range = data%Lev2a%range%refrac_qual, &
476	rec = irec)
477	IF (ncdf_isvar('dry_temp')) THEN !For backward compatibility
478	CALL ncdf_getvar('dry_temp', data%Lev2a%dry_temp, &
479	units = data%Lev2a%units%dry_temp, &
480	range = data%Lev2a%range%dry_temp, &
481	rec = irec)
482	CALL ncdf_getvar('dry_temp_sigma', data%Lev2a%dry_temp_sigma, &
483	units = data%Lev2a%units%dry_temp_sigma, &
484	range = data%Lev2a%range%dry_temp_sigma, &
485	rec = irec)
486	CALL ncdf_getvar('dry_temp_qual', data%Lev2a%dry_temp_qual, &
487	units = data%Lev2a%units%dry_temp_qual, &
488	range = data%Lev2a%range%dry_temp_qual, &
489	rec = irec)
490	ENDIF
491	data%Lev2a%Missing = .FALSE.
492	ENDIF
493
494	! 1.13 Level2b variables (if any)
495	! -------------------------------
496
497	IF (ncdf_isvar('geop')) THEN
498	CALL ncdf_getsize('geop', n, dim = 1)
499	CALL ropp_io_init(data%Lev2b, n)
500	ELSE
501	data%Lev2b%Npoints = 0
502	ENDIF
503
504	IF (data%Lev2b%Npoints > 0) THEN
505
506	CALL ncdf_getvar('geop', data%Lev2b%geop, &
507	units = data%Lev2b%units%geop, &
508	range = data%Lev2b%range%geop, &
509	rec = irec)
510	CALL ncdf_getvar('geop_sigma', data%Lev2b%geop_sigma, &
511	units = data%Lev2b%units%geop_sigma, &
512	range = data%Lev2b%range%geop_sigma, &
513	rec = irec)
514	CALL ncdf_getvar('press', data%Lev2b%press, &
515	units = data%Lev2b%units%press, &
516	range = data%Lev2b%range%press, &
517	rec = irec)
518	CALL ncdf_getvar('press_sigma', data%Lev2b%press_sigma, &
519	units = data%Lev2b%units%press_sigma, &
520	range = data%Lev2b%range%press_sigma, &
521	rec = irec)
522	CALL ncdf_getvar('temp', data%Lev2b%temp, &
523	units = data%Lev2b%units%temp, &
524	range = data%Lev2b%range%temp, &
525	rec = irec)
526	CALL ncdf_getvar('temp_sigma', data%Lev2b%temp_sigma, &
527	units = data%Lev2b%units%temp_sigma, &
528	range = data%Lev2b%range%temp_sigma, &
529	rec = irec)
530	CALL ncdf_getvar('shum', data%Lev2b%shum, &
531	units = data%Lev2b%units%shum, &
532	range = data%Lev2b%range%shum, &
533	rec = irec)
534	CALL ncdf_getvar('shum_sigma', data%Lev2b%shum_sigma, &
535	units = data%Lev2b%units%shum_sigma, &
536	range = data%Lev2b%range%shum_sigma, &
537	rec = irec)
538	CALL ncdf_getvar('meteo_qual', data%Lev2b%meteo_qual, &
539	units = data%Lev2b%units%meteo_qual, &
540	range = data%Lev2b%range%meteo_qual, &
541	rec = irec)
542
543	data%Lev2b%Missing = .FALSE.
544	ENDIF
545
546	! 1.14 Level2c variables (if any)
547	! -------------------------------
548
549	IF (ncdf_isvar('geop_sfc')) THEN
550	data%Lev2c%Npoints = 1
551	ELSE
552	data%Lev2c%Npoints = 0
553	ENDIF
554
555	IF (data%Lev2c%Npoints > 0) THEN
556
557	CALL ncdf_getvar('geop_sfc', data%Lev2c%geop_sfc, &
558	units = data%Lev2c%units%geop_sfc, &
559	range = data%Lev2c%range%geop_sfc, &
560	rec = irec)
561	CALL ncdf_getvar('press_sfc', data%Lev2c%press_sfc, &
562	units = data%Lev2c%units%press_sfc, &
563	range = data%Lev2c%range%press_sfc, &
564	rec = irec)
565	CALL ncdf_getvar('press_sfc_sigma', data%Lev2c%press_sfc_sigma, &
566	units = data%Lev2c%units%press_sfc_sigma, &
567	range = data%Lev2c%range%press_sfc_sigma, &
568	rec = irec)
569	CALL ncdf_getvar('press_sfc_qual', data%Lev2c%press_sfc_qual, &
570	units = data%Lev2c%units%press_sfc_qual, &
571	range = data%Lev2c%range%press_sfc_qual, &
572	rec = irec)
573
574	IF (ncdf_isvar('Ne_max')) THEN !For backward compatibility
575	CALL ncdf_getvar('Ne_max', data%Lev2c%Ne_max, &
576	units = data%Lev2c%units%Ne_max, &
577	range = data%Lev2c%range%Ne_max, &
578	rec = irec)
579	ENDIF
580	IF (ncdf_isvar('Ne_max_sigma')) THEN !For backward compatibility
581	CALL ncdf_getvar('Ne_max_sigma', data%Lev2c%Ne_max_sigma, &
582	units = data%Lev2c%units%Ne_max_sigma, &
583	range = data%Lev2c%range%Ne_max_sigma, &
584	rec = irec)
585	ENDIF
586
587	IF (ncdf_isvar('H_peak')) THEN !For backward compatibility
588	CALL ncdf_getvar('H_peak', data%Lev2c%H_peak, &
589	units = data%Lev2c%units%H_peak, &
590	range = data%Lev2c%range%H_peak, &
591	rec = irec)
592	ENDIF
593	IF (ncdf_isvar('H_peak_sigma')) THEN !For backward compatibility
594	CALL ncdf_getvar('H_peak_sigma', data%Lev2c%H_peak_sigma, &
595	units = data%Lev2c%units%H_peak_sigma, &
596	range = data%Lev2c%range%H_peak_sigma, &
597	rec = irec)
598	ENDIF
599
600	IF (ncdf_isvar('H_width')) THEN !For backward compatibility
601	CALL ncdf_getvar('H_width', data%Lev2c%H_width, &
602	units = data%Lev2c%units%H_width, &
603	range = data%Lev2c%range%H_width, &
604	rec = irec)
605	ENDIF
606	IF (ncdf_isvar('H_width_sigma')) THEN !For backward compatibility
607	CALL ncdf_getvar('H_width_sigma', data%Lev2c%H_width_sigma, &
608	units = data%Lev2c%units%H_width_sigma, &
609	range = data%Lev2c%range%H_width_sigma, &
610	rec = irec)
611	ENDIF
612
613	IF (ncdf_isvar('tph_bangle')) THEN !For backward compatibility
614	CALL ncdf_getvar('tph_bangle', data%Lev2c%tph_bangle, &
615	units = data%Lev2c%units%tph_bangle, &
616	range = data%Lev2c%range%tph_bangle, &
617	rec = irec)
618	ENDIF
619	IF (ncdf_isvar('tpa_bangle')) THEN !For backward compatibility
620	CALL ncdf_getvar('tpa_bangle', data%Lev2c%tpa_bangle, &
621	units = data%Lev2c%units%tpa_bangle, &
622	range = data%Lev2c%range%tpa_bangle, &
623	rec = irec)
624	ENDIF
625	IF (ncdf_isvar('tph_bangle_flag')) THEN !For backward compatibility
626	CALL ncdf_getvar('tph_bangle_flag', data%Lev2c%tph_bangle_flag, &
627	units = data%Lev2c%units%tph_bangle_flag, &
628	range = data%Lev2c%range%tph_bangle_flag, &
629	rec = irec)
630	ENDIF
631
632	IF (ncdf_isvar('tph_refrac')) THEN !For backward compatibility
633	CALL ncdf_getvar('tph_refrac', data%Lev2c%tph_refrac, &
634	units = data%Lev2c%units%tph_refrac, &
635	range = data%Lev2c%range%tph_refrac, &
636	rec = irec)
637	ENDIF
638	IF (ncdf_isvar('tpn_refrac')) THEN !For backward compatibility
639	CALL ncdf_getvar('tpn_refrac', data%Lev2c%tpn_refrac, &
640	units = data%Lev2c%units%tpn_refrac, &
641	range = data%Lev2c%range%tpn_refrac, &
642	rec = irec)
643	ENDIF
644	IF (ncdf_isvar('tph_refrac_flag')) THEN !For backward compatibility
645	CALL ncdf_getvar('tph_refrac_flag', data%Lev2c%tph_refrac_flag, &
646	units = data%Lev2c%units%tph_refrac_flag, &
647	range = data%Lev2c%range%tph_refrac_flag, &
648	rec = irec)
649	ENDIF
650
651	IF (ncdf_isvar('tph_tdry_lrt')) THEN !For backward compatibility
652	CALL ncdf_getvar('tph_tdry_lrt', data%Lev2c%tph_tdry_lrt, &
653	units = data%Lev2c%units%tph_tdry_lrt, &
654	range = data%Lev2c%range%tph_tdry_lrt, &
655	rec = irec)
656	ENDIF
657	IF (ncdf_isvar('tpt_tdry_lrt')) THEN !For backward compatibility
658	CALL ncdf_getvar('tpt_tdry_lrt', data%Lev2c%tpt_tdry_lrt, &
659	units = data%Lev2c%units%tpt_tdry_lrt, &
660	range = data%Lev2c%range%tpt_tdry_lrt, &
661	rec = irec)
662	ENDIF
663	IF (ncdf_isvar('tph_tdry_lrt_flag')) THEN !For backward compatibility
664	CALL ncdf_getvar('tph_tdry_lrt_flag', data%Lev2c%tph_tdry_lrt_flag, &
665	units = data%Lev2c%units%tph_tdry_lrt_flag, &
666	range = data%Lev2c%range%tph_tdry_lrt_flag, &
667	rec = irec)
668	ENDIF
669
670	IF (ncdf_isvar('tph_tdry_cpt')) THEN !For backward compatibility
671	CALL ncdf_getvar('tph_tdry_cpt', data%Lev2c%tph_tdry_cpt, &
672	units = data%Lev2c%units%tph_tdry_cpt, &
673	range = data%Lev2c%range%tph_tdry_cpt, &
674	rec = irec)
675	ENDIF
676	IF (ncdf_isvar('tpt_tdry_cpt')) THEN !For backward compatibility
677	CALL ncdf_getvar('tpt_tdry_cpt', data%Lev2c%tpt_tdry_cpt, &
678	units = data%Lev2c%units%tpt_tdry_cpt, &
679	range = data%Lev2c%range%tpt_tdry_cpt, &
680	rec = irec)
681	ENDIF
682	IF (ncdf_isvar('tph_tdry_cpt_flag')) THEN !For backward compatibility
683	CALL ncdf_getvar('tph_tdry_cpt_flag', data%Lev2c%tph_tdry_cpt_flag, &
684	units = data%Lev2c%units%tph_tdry_cpt_flag, &
685	range = data%Lev2c%range%tph_tdry_cpt_flag, &
686	rec = irec)
687	ENDIF
688
689	IF (ncdf_isvar('prh_tdry_cpt')) THEN !For backward compatibility
690	CALL ncdf_getvar('prh_tdry_cpt', data%Lev2c%prh_tdry_cpt, &
691	units = data%Lev2c%units%prh_tdry_cpt, &
692	range = data%Lev2c%range%prh_tdry_cpt, &
693	rec = irec)
694	ENDIF
695	IF (ncdf_isvar('prt_tdry_cpt')) THEN !For backward compatibility
696	CALL ncdf_getvar('prt_tdry_cpt', data%Lev2c%prt_tdry_cpt, &
697	units = data%Lev2c%units%prt_tdry_cpt, &
698	range = data%Lev2c%range%prt_tdry_cpt, &
699	rec = irec)
700	ENDIF
701	IF (ncdf_isvar('prh_tdry_cpt_flag')) THEN !For backward compatibility
702	CALL ncdf_getvar('prh_tdry_cpt_flag', data%Lev2c%prh_tdry_cpt_flag, &
703	units = data%Lev2c%units%prh_tdry_cpt_flag, &
704	range = data%Lev2c%range%prh_tdry_cpt_flag, &
705	rec = irec)
706	ENDIF
707
708	IF (ncdf_isvar('tph_temp_lrt')) THEN !For backward compatibility
709	CALL ncdf_getvar('tph_temp_lrt', data%Lev2c%tph_temp_lrt, &
710	units = data%Lev2c%units%tph_temp_lrt, &
711	range = data%Lev2c%range%tph_temp_lrt, &
712	rec = irec)
713	ENDIF
714	IF (ncdf_isvar('tpt_temp_lrt')) THEN !For backward compatibility
715	CALL ncdf_getvar('tpt_temp_lrt', data%Lev2c%tpt_temp_lrt, &
716	units = data%Lev2c%units%tpt_temp_lrt, &
717	range = data%Lev2c%range%tpt_temp_lrt, &
718	rec = irec)
719	ENDIF
720	IF (ncdf_isvar('tph_temp_lrt_flag')) THEN !For backward compatibility
721	CALL ncdf_getvar('tph_temp_lrt_flag', data%Lev2c%tph_temp_lrt_flag, &
722	units = data%Lev2c%units%tph_temp_lrt_flag, &
723	range = data%Lev2c%range%tph_temp_lrt_flag, &
724	rec = irec)
725	ENDIF
726
727	IF (ncdf_isvar('tph_temp_cpt')) THEN !For backward compatibility
728	CALL ncdf_getvar('tph_temp_cpt', data%Lev2c%tph_temp_cpt, &
729	units = data%Lev2c%units%tph_temp_cpt, &
730	range = data%Lev2c%range%tph_temp_cpt, &
731	rec = irec)
732	ENDIF
733	IF (ncdf_isvar('tpt_temp_cpt')) THEN !For backward compatibility
734	CALL ncdf_getvar('tpt_temp_cpt', data%Lev2c%tpt_temp_cpt, &
735	units = data%Lev2c%units%tpt_temp_cpt, &
736	range = data%Lev2c%range%tpt_temp_cpt, &
737	rec = irec)
738	ENDIF
739	IF (ncdf_isvar('tph_temp_cpt_flag')) THEN !For backward compatibility
740	CALL ncdf_getvar('tph_temp_cpt_flag', data%Lev2c%tph_temp_cpt_flag, &
741	units = data%Lev2c%units%tph_temp_cpt_flag, &
742	range = data%Lev2c%range%tph_temp_cpt_flag, &
743	rec = irec)
744	ENDIF
745
746	IF (ncdf_isvar('prh_temp_cpt')) THEN !For backward compatibility
747	CALL ncdf_getvar('prh_temp_cpt', data%Lev2c%prh_temp_cpt, &
748	units = data%Lev2c%units%prh_temp_cpt, &
749	range = data%Lev2c%range%prh_temp_cpt, &
750	rec = irec)
751	ENDIF
752	IF (ncdf_isvar('prt_temp_cpt')) THEN !For backward compatibility
753	CALL ncdf_getvar('prt_temp_cpt', data%Lev2c%prt_temp_cpt, &
754	units = data%Lev2c%units%prt_temp_cpt, &
755	range = data%Lev2c%range%prt_temp_cpt, &
756	rec = irec)
757	ENDIF
758	IF (ncdf_isvar('prh_temp_cpt_flag')) THEN !For backward compatibility
759	CALL ncdf_getvar('prh_temp_cpt_flag', data%Lev2c%prh_temp_cpt_flag, &
760	units = data%Lev2c%units%prh_temp_cpt_flag, &
761	range = data%Lev2c%range%prh_temp_cpt_flag, &
762	rec = irec)
763	ENDIF
764
765	data%Lev2c%Missing = .FALSE.
766
767	ENDIF
768
769	! 1.15 Level2d variables (if any)
770	! -------------------------------
771
772	IF (ncdf_isvar('level_coeff_a')) THEN
773	CALL ncdf_getsize('level_coeff_a', n, dim = 1)
774	CALL ropp_io_init(data%Lev2d, n)
775	ELSE
776	data%Lev2d%Npoints = 0
777	ENDIF
778
779	IF (data%Lev2d%Npoints > 0) THEN
780
781	CALL ncdf_getvar('level_type', data%Lev2d%level_type, &
782	rec = irec)
783	CALL ncdf_getvar('level_coeff_a', data%Lev2d%level_coeff_a, &
784	units = data%Lev2d%units%level_coeff_a, &
785	range = data%Lev2d%range%level_coeff_a, &
786	rec = irec)
787	CALL ncdf_getvar('level_coeff_b', data%Lev2d%level_coeff_b, &
788	units = data%Lev2d%units%level_coeff_b, &
789	range = data%Lev2d%range%level_coeff_b, &
790	rec = irec)
791
792	data%Lev2d%Missing = .FALSE.
793	ENDIF
794
795	! 1.16 Additional variables (if any)
796	! ----------------------------------
797
798	CALL ropp_io_init(data%vlist)
799
800	DO varid=1,ncdf_nvars
801
802	IF(.NOT. ncdf_read(varid))THEN
803
804	status = nf90_inquire_variable(ncdf_ncid, varid, xtype=TYPE, ndims=ndim)
805
806	IF (TYPE .NE. NF90_CHAR) THEN ! only read scalar variables
807
808	IF(ndim == 1)THEN
809	CALL ropp_io_read_ncdf_get_vlistD0d(varid, data%vlist%VlistD0d, irec)
810	ENDIF
811	IF(ndim == 2)THEN
812	CALL ropp_io_read_ncdf_get_vlistD1d(varid, data%vlist%VlistD1d, irec)
813	ENDIF
814	IF(ndim == 3)THEN
815	CALL ropp_io_read_ncdf_get_vlistD2d(varid, data%vlist%VlistD2d, irec)
816	ENDIF
817	ENDIF
818
819	ENDIF
820
821	ncdf_read(varid) = .FALSE. ! reset 'read variable' flag
822	ENDDO
823
824	! 1.17 Clean up
825	! -------------
826
827	CALL message_set_routine(routine)
828
829	END SUBROUTINE ropp_io_read_ncdf_get_rodata
830
831
832	!-------------------------------------------------------------------------------
833	! 2. Core RO data (two-dimensional meteorological data)
834	!-------------------------------------------------------------------------------
835
836	SUBROUTINE ropp_io_read_ncdf_get_rodata_2d(DATA, rec)
837
838	! 2.1 Declarations
839	! ----------------
840
841	USE ropp_utils
842	USE ncdf
843	USE ropp_io, not_this => ropp_io_read_ncdf_get_rodata_2d
844	USE ropp_io_types, ONLY: ROprof2d, &
845	ThisFmtVer
846
847	IMPLICIT NONE
848
849	TYPE(ROprof2d), INTENT(inout) :: DATA
850	INTEGER, OPTIONAL :: rec
851
852	INTEGER :: n, ierr
853	INTEGER :: irec
854	INTEGER, DIMENSION(2) :: n2d
855
856	INTEGER, DIMENSION(8) :: DT8 ! Date/time array
857
858	CHARACTER(LEN=15) :: cval
859
860	REAL(dp) :: time, start_time, dtocc_time, now_time !dp defined in DateTimeTypes
861	CHARACTER(len = 23) :: proc_date
862	CHARACTER(len = 256) :: routine
863
864	REAL :: fmtver
865
866	INTEGER :: status, varid, ndim, TYPE
867
868	! 2.2 Error handling
869	! ------------------
870
871	CALL message_get_routine(routine)
872	CALL message_set_routine('ropp_io_read_ncdf_get')
873
874	! 2.3 Default parameters
875	! ----------------------
876
877	IF (PRESENT(rec)) THEN
878	irec = rec
879	ELSE
880	irec = 1
881	ENDIF
882
883	! 2.4 (Global) Attributes
884	! ------------------------
885
886	data%FmtVersion = ' ' ; CALL ncdf_getatt('format_version', data%FmtVersion)
887	READ ( data%FmtVersion(11:), fmt=*, iostat=ierr ) fmtver
888	IF ( ierr /= 0 ) data%FmtVersion = ThisFmtVer
889	data%processing_centre = ' ' ; CALL ncdf_getatt('processing_centre', data%processing_centre)
890	IF (ncdf_isatt('processing_software')) THEN ! added at V8.0
891	data%processing_software = ' ' ; CALL ncdf_getatt('processing_software', data%processing_software)
892	ENDIF
893	proc_date = ' ' ; CALL ncdf_getatt('processing_date', proc_date)
894	data%pod_method = ' ' ; CALL ncdf_getatt('pod_method', data%pod_method)
895	data%phase_method = ' ' ; CALL ncdf_getatt('phase_method', data%phase_method)
896	data%bangle_method = ' ' ; CALL ncdf_getatt('bangle_method', data%bangle_method)
897	data%refrac_method = ' ' ; CALL ncdf_getatt('refrac_method', data%refrac_method)
898	data%meteo_method = ' ' ; CALL ncdf_getatt('meteo_method', data%meteo_method)
899	IF(ncdf_isatt('thin_method'))THEN ! added at V1.1
900	data%thin_method = ' ' ; CALL ncdf_getatt('thin_method', data%thin_method)
901	ENDIF
902	data%software_version = ' ' ; CALL ncdf_getatt('software_version', data%software_version)
903
904	IF (proc_date( 1: 4) /= ' ') READ(proc_date( 1: 4), *) data%DTpro%Year
905	IF (proc_date( 6: 7) /= ' ') READ(proc_date( 6: 7), *) data%DTpro%Month
906	IF (proc_date( 9:10) /= ' ') READ(proc_date( 9:10), *) data%DTpro%Day
907	IF (proc_date(12:13) /= ' ') READ(proc_date(12:13), *) data%DTpro%Hour
908	IF (proc_date(15:16) /= ' ') READ(proc_date(15:16), *) data%DTpro%Minute
909	IF (proc_date(18:19) /= ' ') READ(proc_date(18:19), *) data%DTpro%Second
910	IF (proc_date(21:23) /= ' ') READ(proc_date(21:23), *) data%DTpro%Msec
911
912	! 2.5 Header variables
913	! --------------------
914
915	CALL ncdf_getvar('occ_id', data%occ_id, rec = irec)
916	CALL ncdf_getvar('gns_id', data%gns_id, rec = irec)
917	CALL ncdf_getvar('leo_id', data%leo_id, rec = irec)
918	CALL ncdf_getvar('stn_id', data%stn_id, rec = irec)
919
920	! 2.6 Date and time
921	! -----------------
922
923	CALL ncdf_getvar('start_time', start_time, rec=irec)
924	CALL ncdf_getvar('year', data%DTocc%Year, &
925	units = data%DTocc%units%Year, &
926	range = data%DTocc%range%Year, &
927	rec = irec)
928	CALL ncdf_getvar('month', data%DTocc%Month, &
929	units = data%DTocc%units%Month, &
930	range = data%DTocc%range%Month, &
931	rec = irec)
932	CALL ncdf_getvar('day', data%DTocc%Day, &
933	units = data%DTocc%units%Day, &
934	range = data%DTocc%range%Day, &
935	rec = irec)
936	CALL ncdf_getvar('hour', data%DTocc%Hour, &
937	units = data%DTocc%units%Hour, &
938	range = data%DTocc%range%Hour, &
939	rec = irec)
940	CALL ncdf_getvar('minute', data%DTocc%Minute, &
941	units = data%DTocc%units%Minute, &
942	range = data%DTocc%range%Minute, &
943	rec = irec)
944	CALL ncdf_getvar('second', data%DTocc%Second, &
945	units = data%DTocc%units%Second, &
946	range = data%DTocc%range%Second, &
947	rec = irec)
948	CALL ncdf_getvar('msec', data%DTocc%Msec, &
949	units = data%DTocc%units%Msec, &
950	range = data%DTocc%range%Msec, &
951	rec = irec)
952
953	! 2.6.1 Check consistency: start_time and DTocc (if both are valid)
954	! should refer to the same epoch within 1ms. Issue a warning if not.
955
956	CALL Date_and_Time_UTC(Values=DT8)
957
958	CALL TimeSince(DT8, now_time, 1, Base="JS2000")
959
960	IF (isroppinrange(data%DTocc)) THEN
961	IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN
962	DT8 = (/ data%DTocc%Year, data%DTocc%Month, data%DTocc%Day, 0, &
963	data%DTocc%Hour, data%DTocc%Minute, data%DTocc%Second, &
964	data%DTocc%Msec /)
965	CALL TimeSince(DT8, dtocc_time, 1, Base="JS2000")
966	time = ABS(start_time - dtocc_time)
967	IF ( time > 0.0005_dp ) THEN
968	WRITE ( cval, FMT="(F15.3)" ) time
969	CALL message(msg_warn, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // &
970	"timestamps differ by " // &
971	TRIM(ADJUSTL(cval))//" seconds - using yr/../ms timestamp")
972	END IF
973	END IF
974
975	! If any DTocc element is invalid, substitute converted start_time
976	! (if that is valid) and issue a warning.
977	ELSE
978
979	IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN
980	CALL message(msg_warn, "One or more of yr/mo/dy/hr/mn/sc/ms times " // &
981	"are invalid - using 'start_time' instead")
982	CALL TimeSince(DT8, start_time, -1, Base="JS2000")
983	data%DTocc%Year = DT8(1)
984	data%DTocc%Month = DT8(2)
985	data%DTocc%Day = DT8(3)
986	data%DTocc%Hour = DT8(5)
987	data%DTocc%Minute = DT8(6)
988	data%DTocc%Second = DT8(7)
989	data%DTocc%Msec = DT8(8)
990	END IF
991
992	! NB if neither DTocc nor start_time are valid, do nothing here; missing
993	! coordinates should be dealt with as part of generic Q/C.
994	END IF
995
996	! 2.7 Overall quality
997	! -------------------
998
999	CALL ncdf_getvar('pcd', data%pcd, &
1000	units = data%units%pcd, &
1001	range = data%range%pcd, &
1002	rec = irec)
1003	CALL ncdf_getvar('overall_qual', data%overall_qual, &
1004	units = data%units%overall_qual, &
1005	range = data%range%overall_qual, &
1006	rec = irec)
1007
1008	! 2.8 Georeferencing
1009	! ------------------
1010
1011	CALL ncdf_getvar('time', time, rec=irec)
1012	CALL ncdf_getvar('lat', data%georef%lat, &
1013	units = data%georef%units%lat, &
1014	range = data%georef%range%lat, &
1015	rec = irec)
1016	CALL ncdf_getvar('lon', data%georef%lon, &
1017	units = data%georef%units%lon, &
1018	range = data%georef%range%lon, &
1019	rec = irec)
1020	CALL ncdf_getvar('time_offset', data%georef%time_offset, &
1021	units = data%georef%units%time_offset, &
1022	range = data%georef%range%time_offset, &
1023	rec = irec)
1024	CALL ncdf_getvar('undulation', data%georef%Undulation, &
1025	units = data%georef%units%Undulation, &
1026	range = data%georef%range%undulation, &
1027	rec = irec)
1028	CALL ncdf_getvar('roc', data%georef%roc, &
1029	units = data%georef%units%roc, &
1030	range = data%georef%range%roc, &
1031	rec = irec)
1032	CALL ncdf_getvar('r_coc', data%georef%r_coc, &
1033	units = data%georef%units%r_coc, &
1034	range = data%georef%range%r_coc, &
1035	rec = irec)
1036	CALL ncdf_getvar('azimuth', data%georef%azimuth, &
1037	units = data%georef%units%azimuth, &
1038	range = data%georef%range%azimuth, &
1039	rec = irec)
1040
1041	! 2.8.1 Other attributes
1042
1043	CALL ncdf_getatt('reference_frame', data%georef%reference_frame%r_coc, varname= 'r_coc')
1044
1045	! 2.9 Background characterisation (if any)
1046	! ----------------------------------------
1047
1048	IF (ncdf_isvar('bg_source')) THEN
1049	data%BG%Source = 'TBD'
1050	ELSE
1051	data%BG%Source = 'NONE'
1052	ENDIF
1053
1054	IF (data%BG%Source /= 'NONE') THEN
1055	CALL ncdf_getvar('bg_source', data%BG%Source, rec = irec)
1056	CALL ncdf_getvar('bg_year', data%BG%Year, &
1057	units = data%BG%units%Year, &
1058	range = data%BG%range%Year, &
1059	rec = irec)
1060	CALL ncdf_getvar('bg_month', data%BG%Month, &
1061	units = data%BG%units%Month, &
1062	range = data%BG%range%Month, &
1063	rec = irec)
1064	CALL ncdf_getvar('bg_day', data%BG%Day, &
1065	units = data%BG%units%Day, &
1066	range = data%BG%range%Day, &
1067	rec = irec)
1068	CALL ncdf_getvar('bg_hour', data%BG%Hour, &
1069	units = data%BG%units%Hour, &
1070	range = data%BG%range%Hour, &
1071	rec = irec)
1072	CALL ncdf_getvar('bg_minute', data%BG%Minute, &
1073	units = data%BG%units%Minute, &
1074	range = data%BG%range%Minute, &
1075	rec = irec)
1076	CALL ncdf_getvar('bg_fcperiod', data%BG%fcPeriod, &
1077	units = data%BG%units%fcPeriod, &
1078	range = data%BG%range%fcPeriod, &
1079	rec = irec)
1080	ENDIF
1081
1082	! 2.10 Level1a variables (if any)
1083	! ------------------------------
1084
1085	IF (ncdf_isvar('dtime')) THEN
1086	CALL ncdf_getsize('dtime', n, dim = 1)
1087	CALL ropp_io_init(data%Lev1a, n)
1088	ELSE
1089	data%Lev1a%Npoints = 0
1090	ENDIF
1091
1092	IF (data%Lev1a%Npoints > 0) THEN
1093
1094	CALL ncdf_getvar('dtime', data%Lev1a%dtime, &
1095	units = data%Lev1a%units%dtime, &
1096	range = data%Lev1a%range%dtime, &
1097	rec = irec)
1098	CALL ncdf_getvar('snr_L1ca', data%Lev1a%snr_L1ca, &
1099	units = data%Lev1a%units%snr, &
1100	range = data%Lev1a%range%snr, &
1101	rec = irec)
1102	CALL ncdf_getvar('snr_L1p', data%Lev1a%snr_L1p, &
1103	units = data%Lev1a%units%snr, &
1104	range = data%Lev1a%range%snr, &
1105	rec = irec)
1106	CALL ncdf_getvar('snr_L2p', data%Lev1a%snr_L2p, &
1107	units = data%Lev1a%units%snr, &
1108	range = data%Lev1a%range%snr, &
1109	rec = irec)
1110	CALL ncdf_getvar('phase_L1', data%Lev1a%phase_L1, &
1111	units = data%Lev1a%units%phase, &
1112	range = data%Lev1a%range%phase, &
1113	rec = irec)
1114	CALL ncdf_getvar('phase_L2', data%Lev1a%phase_L2, &
1115	units = data%Lev1a%units%phase, &
1116	range = data%Lev1a%range%phase, &
1117	rec = irec)
1118	CALL ncdf_getvar('r_gns', data%Lev1a%r_gns, &
1119	units = data%Lev1a%units%r_gns, &
1120	range = data%Lev1a%range%r_gns, &
1121	rec = irec)
1122	CALL ncdf_getvar('v_gns', data%Lev1a%v_gns, &
1123	units = data%Lev1a%units%v_gns, &
1124	range = data%Lev1a%range%v_gns, &
1125	rec = irec)
1126	CALL ncdf_getvar('r_leo', data%Lev1a%r_leo, &
1127	units = data%Lev1a%units%r_leo, &
1128	range = data%Lev1a%range%r_leo, &
1129	rec = irec)
1130	CALL ncdf_getvar('v_leo', data%Lev1a%v_leo, &
1131	units = data%Lev1a%units%v_leo, &
1132	range = data%Lev1a%range%v_leo, &
1133	rec = irec)
1134	CALL ncdf_getvar('phase_qual', data%Lev1a%phase_qual, &
1135	units = data%Lev1a%units%phase_qual, &
1136	range = data%Lev1a%range%phase_qual, &
1137	rec = irec)
1138
1139	! 2.10.1 Other attributes
1140
1141	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_gns, varname = 'r_gns')
1142	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_gns, varname = 'v_gns')
1143	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_leo, varname = 'r_leo')
1144	CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_leo, varname = 'v_leo')
1145
1146	ENDIF
1147
1148	! 2.11 Level1b variables (if any)
1149	! -------------------------------
1150
1151	IF (ncdf_isvar('lat_tp')) THEN
1152	CALL ncdf_getsize('lat_tp', n, dim = 1)
1153	CALL ropp_io_init(data%Lev1b, n)
1154	ELSE
1155	data%Lev1b%Npoints = 0
1156	ENDIF
1157
1158	IF (data%Lev1b%Npoints > 0) THEN
1159
1160	CALL ncdf_getvar('lat_tp', data%Lev1b%lat_tp, &
1161	units = data%Lev1b%units%lat_tp, &
1162	range = data%Lev1b%range%lat_tp, &
1163	rec = irec)
1164	CALL ncdf_getvar('lon_tp', data%Lev1b%lon_tp, &
1165	units = data%Lev1b%units%lon_tp, &
1166	range = data%Lev1b%range%lon_tp, &
1167	rec = irec)
1168	CALL ncdf_getvar('azimuth_tp', data%Lev1b%azimuth_tp, &
1169	units = data%Lev1b%units%azimuth_tp, &
1170	range = data%Lev1b%range%azimuth_tp, &
1171	rec = irec)
1172
1173	CALL ncdf_getvar('impact_L1', data%Lev1b%impact_L1, &
1174	units = data%Lev1b%units%impact, &
1175	range = data%Lev1b%range%impact, &
1176	rec = irec)
1177	CALL ncdf_getvar('impact_L2', data%Lev1b%impact_L2, &
1178	units = data%Lev1b%units%impact, &
1179	range = data%Lev1b%range%impact, &
1180	rec = irec)
1181	CALL ncdf_getvar('impact', data%Lev1b%impact, &
1182	units = data%Lev1b%units%impact, &
1183	range = data%Lev1b%range%impact, &
1184	rec = irec)
1185	IF (ncdf_isvar('impact_opt')) & ! added at v1.1
1186	CALL ncdf_getvar('impact_opt', data%Lev1b%impact_opt, &
1187	units = data%Lev1b%units%impact, &
1188	range = data%Lev1b%range%impact, &
1189	rec = irec)
1190
1191	CALL ncdf_getvar('bangle_L1', data%Lev1b%bangle_L1, &
1192	units = data%Lev1b%units%bangle, &
1193	range = data%Lev1b%range%bangle, &
1194	rec = irec)
1195	CALL ncdf_getvar('bangle_L2', data%Lev1b%bangle_L2, &
1196	units = data%Lev1b%units%bangle, &
1197	range = data%Lev1b%range%bangle, &
1198	rec = irec)
1199	CALL ncdf_getvar('bangle', data%Lev1b%bangle, &
1200	units = data%Lev1b%units%bangle, &
1201	range = data%Lev1b%range%bangle, &
1202	rec = irec)
1203	IF (ncdf_isvar('bangle_opt')) & ! added at v1.1
1204	CALL ncdf_getvar('bangle_opt', data%Lev1b%bangle_opt, &
1205	units = data%Lev1b%units%bangle, &
1206	range = data%Lev1b%range%bangle, &
1207	rec = irec)
1208
1209	CALL ncdf_getvar('bangle_L1_sigma', data%Lev1b%bangle_L1_sigma, &
1210	units = data%Lev1b%units%bangle_sigma, &
1211	range = data%Lev1b%range%bangle_sigma, &
1212	rec = irec)
1213	CALL ncdf_getvar('bangle_L2_sigma', data%Lev1b%bangle_L2_sigma, &
1214	units = data%Lev1b%units%bangle_sigma, &
1215	range = data%Lev1b%range%bangle_sigma, &
1216	rec = irec)
1217	CALL ncdf_getvar('bangle_sigma', data%Lev1b%bangle_sigma, &
1218	units = data%Lev1b%units%bangle_sigma, &
1219	range = data%Lev1b%range%bangle_sigma, &
1220	rec = irec)
1221	IF (ncdf_isvar('bangle_opt_sigma')) & ! added at v1.1
1222	CALL ncdf_getvar('bangle_opt_sigma', data%Lev1b%bangle_opt_sigma, &
1223	units = data%Lev1b%units%bangle_sigma, &
1224	range = data%Lev1b%range%bangle_sigma, &
1225	rec = irec)
1226
1227	CALL ncdf_getvar('bangle_L1_qual', data%Lev1b%bangle_L1_qual, &
1228	units = data%Lev1b%units%bangle_qual, &
1229	range = data%Lev1b%range%bangle_qual, &
1230	rec = irec)
1231	CALL ncdf_getvar('bangle_L2_qual', data%Lev1b%bangle_L2_qual, &
1232	units = data%Lev1b%units%bangle_qual, &
1233	range = data%Lev1b%range%bangle_qual, &
1234	rec = irec)
1235	CALL ncdf_getvar('bangle_qual', data%Lev1b%bangle_qual, &
1236	units = data%Lev1b%units%bangle_qual, &
1237	range = data%Lev1b%range%bangle_qual, &
1238	rec = irec)
1239	IF (ncdf_isvar('bangle_opt_qual')) & ! added at v1.1
1240	CALL ncdf_getvar('bangle_opt_qual', data%Lev1b%bangle_opt_qual, &
1241	units = data%Lev1b%units%bangle_qual, &
1242	range = data%Lev1b%range%bangle_qual, &
1243	rec = irec)
1244	ENDIF
1245
1246	! 2.12 Level2a variables (if any)
1247	! -------------------------------
1248
1249	IF (ncdf_isvar('alt_refrac')) THEN
1250	CALL ncdf_getsize('alt_refrac', n, dim = 1)
1251	CALL ropp_io_init(data%Lev2a, n)
1252	ELSE
1253	data%Lev2a%Npoints = 0
1254	ENDIF
1255
1256	IF (data%Lev2a%Npoints > 0) THEN
1257
1258	CALL ncdf_getvar('alt_refrac', data%Lev2a%alt_refrac, &
1259	units = data%Lev2a%units%alt_refrac, &
1260	range = data%Lev2a%range%alt_refrac, &
1261	rec = irec)
1262	CALL ncdf_getvar('geop_refrac', data%Lev2a%geop_refrac, &
1263	units = data%Lev2a%units%geop_refrac, &
1264	range = data%Lev2a%range%geop_refrac, &
1265	rec = irec)
1266	CALL ncdf_getvar('refrac', data%Lev2a%refrac, &
1267	units = data%Lev2a%units%refrac, &
1268	range = data%Lev2a%range%refrac, &
1269	rec = irec)
1270	CALL ncdf_getvar('refrac_sigma', data%Lev2a%refrac_sigma, &
1271	units = data%Lev2a%units%refrac_sigma, &
1272	range = data%Lev2a%range%refrac_sigma, &
1273	rec = irec)
1274	CALL ncdf_getvar('refrac_qual', data%Lev2a%refrac_qual, &
1275	units = data%Lev2a%units%refrac_qual, &
1276	range = data%Lev2a%range%refrac_qual, &
1277	rec = irec)
1278	IF (ncdf_isvar('dry_temp')) THEN !For backward compatibility
1279	CALL ncdf_getvar('dry_temp', data%Lev2a%dry_temp, &
1280	units = data%Lev2a%units%dry_temp, &
1281	range = data%Lev2a%range%dry_temp, &
1282	rec = irec)
1283	CALL ncdf_getvar('dry_temp_sigma', data%Lev2a%dry_temp_sigma, &
1284	units = data%Lev2a%units%dry_temp_sigma, &
1285	range = data%Lev2a%range%dry_temp_sigma, &
1286	rec = irec)
1287	CALL ncdf_getvar('dry_temp_qual', data%Lev2a%dry_temp_qual, &
1288	units = data%Lev2a%units%dry_temp_qual, &
1289	range = data%Lev2a%range%dry_temp_qual, &
1290	rec = irec)
1291	ENDIF
1292
1293	ENDIF
1294
1295	! 2.13 Level2b variables (if any)
1296	! -------------------------------
1297
1298	IF (ncdf_isvar('geop')) THEN
1299	n2d(1) = 0
1300	n2d(2) = 0
1301	CALL ncdf_getsize('geop', n2d)
1302	CALL ropp_io_init(data%Lev2b, n2d)
1303	ELSE
1304	data%Lev2b%Npoints = 0
1305	data%Lev2b%Nhoriz = 0
1306	ENDIF
1307
1308	IF (data%Lev2b%Npoints > 0) THEN
1309
1310	CALL ncdf_getvar('geop', data%Lev2b%geop, &
1311	units = data%Lev2b%units%geop, &
1312	range = data%Lev2b%range%geop, &
1313	rec = irec)
1314	CALL ncdf_getvar('geop_sigma', data%Lev2b%geop_sigma, &
1315	units = data%Lev2b%units%geop_sigma, &
1316	range = data%Lev2b%range%geop_sigma, &
1317	rec = irec)
1318	CALL ncdf_getvar('press', data%Lev2b%press, &
1319	units = data%Lev2b%units%press, &
1320	range = data%Lev2b%range%press, &
1321	rec = irec)
1322	CALL ncdf_getvar('press_sigma', data%Lev2b%press_sigma, &
1323	units = data%Lev2b%units%press_sigma, &
1324	range = data%Lev2b%range%press_sigma, &
1325	rec = irec)
1326	CALL ncdf_getvar('temp', data%Lev2b%temp, &
1327	units = data%Lev2b%units%temp, &
1328	range = data%Lev2b%range%temp, &
1329	rec = irec)
1330	CALL ncdf_getvar('temp_sigma', data%Lev2b%temp_sigma, &
1331	units = data%Lev2b%units%temp_sigma, &
1332	range = data%Lev2b%range%temp_sigma, &
1333	rec = irec)
1334	CALL ncdf_getvar('shum', data%Lev2b%shum, &
1335	units = data%Lev2b%units%shum, &
1336	range = data%Lev2b%range%shum, &
1337	rec = irec)
1338	CALL ncdf_getvar('shum_sigma', data%Lev2b%shum_sigma, &
1339	units = data%Lev2b%units%shum_sigma, &
1340	range = data%Lev2b%range%shum_sigma, &
1341	rec = irec)
1342	CALL ncdf_getvar('meteo_qual', data%Lev2b%meteo_qual, &
1343	units = data%Lev2b%units%meteo_qual, &
1344	range = data%Lev2b%range%meteo_qual, &
1345	rec = irec)
1346
1347	ENDIF
1348
1349	! 2.14 Level2c variables (if any)
1350	! -------------------------------
1351
1352	IF (ncdf_isvar('geop_sfc')) THEN
1353	n2d(1) = 1
1354	n2d(2) = 0
1355	CALL ncdf_getsize('geop_sfc', n2d(2), dim = 1)
1356	CALL ropp_io_init(data%Lev2c, n2d)
1357	ELSE
1358	data%Lev2c%Npoints = 0
1359	data%Lev2c%Nhoriz = 0
1360	ENDIF
1361
1362	IF (data%Lev2c%Npoints > 0) THEN
1363
1364	! new 2d code
1365
1366	CALL ncdf_getvar('dtheta', data%Lev2c%dtheta, &
1367	units = data%Lev2c%units%dtheta, &
1368	range = data%Lev2c%range%dtheta, &
1369	rec = irec)
1370	CALL ncdf_getvar('lat_2d', data%Lev2c%lat_2d, &
1371	units = data%Lev2c%units%lat_2d, &
1372	range = data%Lev2c%range%lat_2d, &
1373	rec = irec)
1374	CALL ncdf_getvar('lon_2d', data%Lev2c%lon_2d, &
1375	units = data%Lev2c%units%lon_2d, &
1376	range = data%Lev2c%range%lon_2d, &
1377	rec = irec)
1378	CALL ncdf_getvar('geop_sfc', data%Lev2c%geop_sfc, &
1379	units = data%Lev2c%units%geop_sfc, &
1380	range = data%Lev2c%range%geop_sfc, &
1381	rec = irec)
1382	CALL ncdf_getvar('press_sfc', data%Lev2c%press_sfc, &
1383	units = data%Lev2c%units%press_sfc, &
1384	range = data%Lev2c%range%press_sfc, &
1385	rec = irec)
1386	CALL ncdf_getvar('press_sfc_sigma', data%Lev2c%press_sfc_sigma, &
1387	units = data%Lev2c%units%press_sfc_sigma, &
1388	range = data%Lev2c%range%press_sfc_sigma, &
1389	rec = irec)
1390	CALL ncdf_getvar('press_sfc_qual', data%Lev2c%press_sfc_qual, &
1391	units = data%Lev2c%units%press_sfc_qual, &
1392	range = data%Lev2c%range%press_sfc_qual, &
1393	rec = irec)
1394	ENDIF
1395
1396	! 2.15 Level2d variables (if any)
1397	! -------------------------------
1398
1399	IF (ncdf_isvar('level_coeff_a')) THEN
1400	CALL ncdf_getsize('level_coeff_a', n, dim = 1)
1401	CALL ropp_io_init(data%Lev2d, n)
1402	ELSE
1403	data%Lev2d%Npoints = 0
1404	ENDIF
1405
1406	IF (data%Lev2d%Npoints > 0) THEN
1407
1408	CALL ncdf_getvar('level_type', data%Lev2d%level_type, &
1409	rec = irec)
1410	CALL ncdf_getvar('level_coeff_a', data%Lev2d%level_coeff_a, &
1411	units = data%Lev2d%units%level_coeff_a, &
1412	range = data%Lev2d%range%level_coeff_a, &
1413	rec = irec)
1414	CALL ncdf_getvar('level_coeff_b', data%Lev2d%level_coeff_b, &
1415	units = data%Lev2d%units%level_coeff_b, &
1416	range = data%Lev2d%range%level_coeff_b, &
1417	rec = irec)
1418	ENDIF
1419
1420	! 2.16 Additional variables (if any)
1421	! ----------------------------------
1422
1423	CALL ropp_io_init(data%vlist)
1424
1425	DO varid=1,ncdf_nvars
1426
1427	IF(.NOT. ncdf_read(varid))THEN
1428
1429	status = nf90_inquire_variable(ncdf_ncid, varid, xtype=TYPE, ndims=ndim)
1430
1431	IF (TYPE .NE. NF90_CHAR) THEN ! only read scalar variables
1432
1433	IF(ndim == 1)THEN
1434	CALL ropp_io_read_ncdf_get_vlistD0d(varid, data%vlist%VlistD0d, irec)
1435	ENDIF
1436	IF(ndim == 2)THEN
1437	CALL ropp_io_read_ncdf_get_vlistD1d(varid, data%vlist%VlistD1d, irec)
1438	ENDIF
1439	IF(ndim == 3)THEN
1440	CALL ropp_io_read_ncdf_get_vlistD2d(varid, data%vlist%VlistD2d, irec)
1441	ENDIF
1442	ENDIF
1443
1444	ENDIF
1445
1446	ncdf_read(varid) = .FALSE. ! reset 'read variable' flag
1447
1448	ENDDO
1449
1450	! 2.17 Clean up
1451	! -------------
1452
1453	CALL message_set_routine(routine)
1454
1455	END SUBROUTINE ropp_io_read_ncdf_get_rodata_2d
1456
1457
1458	!-------------------------------------------------------------------------------
1459	! 3. Error correlation and covariance matrices
1460	!-------------------------------------------------------------------------------
1461
1462	SUBROUTINE ropp_io_read_ncdf_get_rocorcov(DATA)
1463
1464	! 3.1 Declarations
1465	! ----------------
1466
1467	USE ropp_utils
1468	USE ncdf
1469	USE ropp_io, not_this => ropp_io_read_ncdf_get_rocorcov
1470	USE ropp_io_types, ONLY: ROcorcov
1471
1472	IMPLICIT NONE
1473
1474	TYPE(ROcorcov), DIMENSION(:), POINTER :: DATA
1475
1476	REAL(wp), DIMENSION(:), POINTER :: lat_min => null()
1477	REAL(wp), DIMENSION(:), POINTER :: lat_max => null()
1478	REAL(wp), DIMENSION(:,:), POINTER :: corr => null()
1479	REAL(wp), DIMENSION(:,:), POINTER :: sigma => null()
1480
1481	INTEGER :: i, m, n
1482	CHARACTER(len = 256) :: routine
1483
1484	! 3.2 Error handling
1485	! ------------------
1486
1487	CALL message_get_routine(routine)
1488	CALL message_set_routine('ropp_io_read_ncdf_get')
1489
1490	! 3.3 Latitude bins
1491	! -----------------
1492
1493	IF (ncdf_isvar('lat_min') .AND. ncdf_isvar('lat_max')) THEN
1494	CALL ncdf_getsize('lat_min', m)
1495	CALL ncdf_getsize('lat_max', n)
1496	ELSE
1497	CALL message(msg_fatal, &
1498	"NetCDF data file does not seem to contain an error correlation or covariance structure.")
1499	ENDIF
1500
1501	IF (m /= m) THEN
1502	CALL message(msg_fatal, &
1503	"Number of latitude bin boundaries in the netCDF data file is inconsistent.")
1504	ENDIF
1505
1506	CALL ncdf_getvar_alloc('lat_min', lat_min)
1507	CALL ncdf_getvar_alloc('lat_max', lat_max)
1508
1509	! 3.4 Error correlation matrices
1510	! ------------------------------
1511
1512	IF (ncdf_isvar('corr')) THEN
1513	CALL ncdf_getvar_alloc('corr', corr)
1514	ELSE
1515	CALL message(msg_fatal, &
1516	"NetCDF data file does not seem to contain an error correlation matrix.")
1517	ENDIF
1518
1519	! 3.5 Error standard deviations
1520	! -----------------------------
1521
1522	IF (ncdf_isvar('sigma')) THEN
1523	CALL ncdf_getvar_alloc('sigma', sigma)
1524	ENDIF
1525
1526	! 3.6 Allocate and fill ROPP structure
1527	! ------------------------------------
1528
1529	ALLOCATE(DATA(n))
1530
1531	DO i = 1, n
1532
1533	DATA(i)%lat_min = lat_min(i)
1534	DATA(i)%lat_max = lat_max(i)
1535
1536	ALLOCATE(DATA(i)%corr(SIZE(corr(:,i), 1)))
1537	DATA(i)%corr = corr(:,i)
1538
1539	IF (ASSOCIATED(sigma)) THEN
1540	ALLOCATE(DATA(i)%sigma(SIZE(sigma(:,i), 1)))
1541	DATA(i)%sigma = sigma(:, i)
1542	ENDIF
1543
1544	ENDDO
1545
1546	! 3.7 Clean up
1547	! ------------
1548
1549	DEALLOCATE(lat_min)
1550	DEALLOCATE(lat_max)
1551	DEALLOCATE(corr)
1552	IF (ASSOCIATED(sigma)) DEALLOCATE(sigma)
1553
1554	CALL message_set_routine(routine)
1555
1556	END SUBROUTINE ropp_io_read_ncdf_get_rocorcov
1557
1558	!-------------------------------------------------------------------------------
1559	! 4. wrapper for other centres' RO data
1560	!-------------------------------------------------------------------------------
1561
1562	SUBROUTINE ropp_io_read_ncdf_get_otherdata(DATA, file, centre, rec, resolution, getlevel1a, getbufr)
1563
1564	! 4.1 Declarations
1565	! ----------------
1566
1567	USE ropp_io, not_this => ropp_io_read_ncdf_get_otherdata
1568	USE ropp_io_types, ONLY: ROprof
1569
1570	IMPLICIT NONE
1571
1572	TYPE(ROprof), INTENT(inout) :: DATA
1573	CHARACTER (len = *), INTENT(in) :: file
1574	CHARACTER (len=20), INTENT(in) :: centre
1575	INTEGER, OPTIONAL :: rec
1576	CHARACTER (len=20), OPTIONAL :: resolution
1577	LOGICAL, OPTIONAL :: getlevel1a
1578	LOGICAL, OPTIONAL :: getbufr
1579
1580	CHARACTER (len=20) :: lresolution = 'thinned'
1581	LOGICAL :: lgetlevel1a = .FALSE.
1582	LOGICAL :: lgetbufr = .FALSE.
1583	LOGICAL :: ldummy = .FALSE.
1584
1585	! defaults
1586	IF (PRESENT(resolution)) lresolution=resolution
1587	IF (PRESENT(getlevel1a)) lgetlevel1a=getlevel1a
1588	IF (PRESENT(getbufr)) lgetbufr=getbufr
1589
1590	! call the appropriate data handling function, default is ROPP format
1591	SELECT CASE (centre)
1592	CASE('UCAR')
1593	CALL ropp_io_read_ncdf_get_ucardata(DATA, file)
1594	CASE('EUM')
1595	CALL ropp_io_read_ncdf_get_eumdata(DATA, file, lresolution, lgetlevel1a, lgetbufr, ldummy)
1596	CASE default
1597	CALL ropp_io_read_ncdf_get_rodata(DATA, rec)
1598	END SELECT
1599
1600	END SUBROUTINE ropp_io_read_ncdf_get_otherdata
1601
1602	!-------------------------------------------------------------------------------
1603	! 5. UCAR RO data
1604	!-------------------------------------------------------------------------------
1605
1606	SUBROUTINE ropp_io_read_ncdf_get_ucardata(DATA, file)
1607
1608	! 5.1 Declarations
1609	! ----------------
1610
1611	USE ncdf
1612	USE ropp_utils
1613	USE ropp_io_types, ONLY: ROprof
1614
1615	IMPLICIT NONE
1616
1617	TYPE(ROprof), INTENT(inout) :: DATA
1618	CHARACTER(len = *), INTENT(in) :: file
1619	CHARACTER(len = 256) :: routine
1620
1621	! 5.2 Error handling
1622	! ------------------
1623
1624	CALL message_get_routine(routine)
1625	CALL message_set_routine('ropp_io_read_ucardata')
1626
1627	! 5.3 Identify file type
1628	! ----------------------
1629
1630	IF (ncdf_isvar('Bend_ang') .AND. &
1631	(ncdf_isvar('Impact_parm') .OR. ncdf_isvar('Impact_height'))) THEN
1632	CALL ropp_io_read_ucardata_atmPrf(DATA, file)
1633	ELSE IF (ncdf_isvar('pL1Snr') .AND. ncdf_isvar('pL2Snr')) THEN
1634	CALL ropp_io_read_ucardata_atmPhs(DATA, file)
1635	ELSE IF (ncdf_isvar('MSL_alt') .AND. ncdf_isvar ('Pres')) THEN
1636	CALL ropp_io_read_ucardata_atmPrf(DATA, file)
1637	ELSE
1638	CALL message(msg_fatal, &
1639	"Routine ropp_io_read_ncdf_get_ucardata does not support this" // &
1640	"file type. Only atmPrf, ecmPrf, gfsPrf, ncpPrf, sonPrf and atmPhs "// &
1641	"files supported. Check input file type.")
1642	ENDIF
1643
1644	CALL message_set_routine(routine)
1645
1646	CONTAINS
1647
1648	!-------------------------------------------------------------------------------
1649	! 6. UCAR RO data - atmPrf files
1650	!-------------------------------------------------------------------------------
1651
1652	SUBROUTINE ropp_io_read_ucardata_atmPrf(DATA, file)
1653
1654	! NB this routine only supports UCAR 'Prf' format netCDF files
1655	! (i.e. atmPrf, ecmPrf, gfsPrf, ncpPrf and sonPrf)
1656	! See: http://cosmic-io.cosmic.ucar.edu/cdaac/fileFormats/atmPrf.html
1657
1658	! 6.1 Declarations
1659	! ----------------
1660
1661	USE DateTimeProgs, ONLY: Date_and_Time_UTC
1662	USE DateTimeTypes
1663	USE ropp_utils
1664	USE ncdf
1665	USE ropp_io
1666	USE ropp_io_types, ONLY: ROprof, &
1667	ThisFmtVer, &
1668	PCD_open_loop, &
1669	PCD_rising, &
1670	PCD_occultation
1671	USE geodesy, ONLY: geometric2geopotential
1672
1673	IMPLICIT NONE
1674
1675	TYPE(ROprof), INTENT(inout) :: DATA
1676	CHARACTER(len = *), INTENT(in) :: file
1677
1678	INTEGER :: n
1679	INTEGER :: readint
1680	REAL(wp) :: readreal
1681	CHARACTER (len = 256) :: readstr
1682
1683	REAL(wp), PARAMETER :: g_wmo = 9.80665_wp
1684	REAL(wp), PARAMETER :: epsilon_water = 0.621971_wp
1685	INTEGER, DIMENSION(8) :: DTnow
1686
1687	! holds where output
1688	! INTEGER, DIMENSION(:), POINTER :: idx => null()
1689	INTEGER :: nidx
1690
1691	! 6.3 Header variables
1692	! --------------------
1693
1694	readstr = ' '
1695	CALL ncdf_getatt('fileStamp', readstr)
1696	data%leo_id = readstr(1:4)
1697	CALL ncdf_getatt('occulting_sat_id', readint)
1698	WRITE(data%gns_id,'(A1,I3.3)') 'G', readint
1699	readstr = ' '
1700	IF (ncdf_isatt('fiducial_id')) THEN
1701	CALL ncdf_getatt('fiducial_id', readstr)
1702	IF(readstr /= " ") data%stn_id = readstr(1:4)
1703	ENDIF
1704
1705	! 6.4 Overall quality
1706	! -------------------
1707
1708	IF (ncdf_isvar('Bend_ang') .AND. &
1709	(ncdf_isvar('Impact_parm') .OR. ncdf_isvar('Impact_height'))) THEN
1710	CALL ncdf_getatt('bad', readstr)
1711	data%PCD = 0
1712	IF (TRIM(readstr) == "0") THEN
1713	data%PCD = 0
1714	data%overall_qual = 100.0
1715	END IF
1716	IF (TRIM(readstr) == "1") THEN
1717	data%PCD = 1 ! non nominal
1718	data%overall_qual = 0.0
1719	END IF
1720	data%units%overall_qual = "%"
1721	IF (ncdf_isatt('iol') .AND. ncdf_isatt('irs')) THEN
1722	CALL ncdf_getatt('iol', readint)
1723	IF (readint == 1) &
1724	data%PCD = IBSET(data%PCD, PCD_open_loop) ! open loop used
1725	CALL ncdf_getatt('irs', readint)
1726	IF (readint == -1) &
1727	data%PCD = IBSET(data%PCD, PCD_rising) ! rising occultation
1728	ENDIF
1729	ELSE
1730	data%PCD = 0
1731	data%PCD = IBSET(data%PCD, PCD_occultation) ! background data
1732	ENDIF
1733
1734	! 6.5 Date and time
1735	! -----------------
1736
1737	CALL gettime(data%DTocc)
1738
1739	! 6.6 Georeferencing
1740	! ------------------
1741
1742	CALL ncdf_getatt('lat', data%georef%lat)
1743	data%georef%units%lat = "degrees_north"
1744
1745	CALL ncdf_getatt('lon', data%georef%lon)
1746	data%georef%units%lon = "degrees_east"
1747
1748	IF (ncdf_isvar('Bend_ang') .AND. &
1749	(ncdf_isvar('Impact_parm') .OR. ncdf_isvar('Impact_height'))) THEN
1750
1751	CALL ncdf_getatt('occpt_offset', data%georef%time_offset)
1752	data%georef%units%time_offset = "seconds"
1753
1754	CALL ncdf_getatt('rgeoid', data%georef%Undulation)
1755	data%georef%Undulation = data%georef%Undulation * 1000.0
1756	data%georef%units%Undulation = "meters"
1757
1758	CALL ncdf_getatt('rfict', data%georef%roc)
1759	data%georef%roc = data%georef%roc * 1000.0
1760	data%georef%units%roc = "meters"
1761
1762	CALL ncdf_getatt('curv', data%georef%r_coc)
1763	data%georef%r_coc = data%georef%r_coc * 1000.0
1764	data%georef%units%r_coc = "meters"
1765
1766	CALL ncdf_getatt('azim', data%georef%azimuth)
1767	IF ((data%georef%azimuth < 0.0) .AND. &
1768	(data%georef%azimuth > -180.0)) &
1769	data%georef%azimuth = data%georef%azimuth + 360.0
1770	data%georef%units%azimuth = "degrees"
1771
1772	ENDIF
1773
1774	! 6.7 Level1b variables (if any)
1775	! -------------------------------
1776
1777	IF (ncdf_isvar('Bend_ang') .AND. &
1778	(ncdf_isvar('Impact_parm') .OR. ncdf_isvar('Impact_height'))) THEN
1779	CALL ncdf_getsize('Bend_ang', n, dim = 1)
1780	CALL ropp_io_init(data%Lev1b, n)
1781	ELSE
1782	data%Lev1b%Npoints = 0
1783	ENDIF
1784
1785	IF (data%Lev1b%Npoints > 0) THEN
1786
1787	CALL ncdf_getvar('Lat', data%Lev1b%lat_tp)
1788	data%Lev1b%units%lat_tp = "degrees_north"
1789
1790	CALL ncdf_getvar('Lon', data%Lev1b%lon_tp)
1791	data%Lev1b%units%lon_tp = "degrees_east"
1792
1793	CALL ncdf_getvar('Azim', data%Lev1b%azimuth_tp)
1794	! idx => WHERE( (data%Lev1b%azimuth_tp < 0.0) .AND. &
1795	! (data%Lev1b%azimuth_tp > -180.0), nidx)
1796	! IF (nidx > 0) data%Lev1b%azimuth_tp(idx) = data%Lev1b%azimuth_tp(idx) + 360.0
1797	WHERE ( (data%Lev1b%azimuth_tp < 0.0) .AND. &
1798	(data%Lev1b%azimuth_tp > -180.0) ) &
1799	data%Lev1b%azimuth_tp = data%Lev1b%azimuth_tp + 360.0
1800	data%Lev1b%units%azimuth_tp = "degrees"
1801
1802	IF (ncdf_isvar('Impact_parm')) THEN
1803
1804	CALL ncdf_getvar('Impact_parm', data%Lev1b%impact, &
1805	units = data%Lev1b%units%impact)
1806
1807	ELSE ! If both are present, use Impact_parm
1808
1809	IF (ncdf_isvar('Impact_height')) THEN
1810
1811	CALL ncdf_getvar('Impact_height', data%Lev1b%impact, &
1812	units = data%Lev1b%units%impact)
1813	data%Lev1b%impact = data%Lev1b%impact + data%georef%roc ! Both are in m at this point
1814
1815	END IF
1816
1817	END IF
1818
1819	data%lev1b%impact_opt = data%Lev1b%impact
1820
1821	CALL ncdf_getvar('Bend_ang', data%Lev1b%bangle)
1822	data%Lev1b%units%bangle = "radians"
1823
1824	CALL ncdf_getvar('Bend_ang_stdv', data%Lev1b%bangle_sigma)
1825	data%Lev1b%units%bangle_sigma = "radians"
1826
1827	data%Lev1b%bangle_opt_sigma = data%Lev1b%bangle_sigma
1828
1829	! set the quality for bangle
1830	! CALL ncdf_getatt('_FillValue', readreal, 'Opt_bend_ang')
1831	! idx => WHERE( data%Lev1b%bangle > readreal, nidx)
1832	! IF (nidx > 0) data%Lev1b%bangle_qual(idx) = 100.0
1833	WHERE (data%Lev1b%bangle > readreal) &
1834	data%Lev1b%bangle_qual = 100.0
1835
1836	IF (ncdf_isvar('Opt_bend_ang')) THEN
1837	CALL ncdf_getvar('Opt_bend_ang', data%Lev1b%bangle_opt)
1838	CALL ncdf_getatt('_FillValue', readreal, 'Opt_bend_ang')
1839	data%Lev1b%bangle_opt_qual = data%Lev1b%bangle_qual
1840	ENDIF
1841
1842	ENDIF
1843
1844	! 6.8 Level2a variables (if any)
1845	! -------------------------------
1846
1847	IF (ncdf_isvar('MSL_alt') .AND. ncdf_isvar('Ref')) THEN
1848	CALL ncdf_getsize('MSL_alt', n, dim = 1)
1849	CALL ropp_io_init(data%Lev2a, n)
1850	ELSE
1851	data%Lev2a%Npoints = 0
1852	ENDIF
1853
1854	IF (data%Lev2a%Npoints > 0) THEN
1855
1856	CALL ncdf_getvar('MSL_alt', data%Lev2a%alt_refrac, &
1857	units = data%Lev2a%units%alt_refrac)
1858	! geopotential not in UCAR files - generate from altitude
1859	! idx => WHERE((data%Lev2a%alt_refrac > -999.0), nidx)
1860	! IF (nidx > 0 .AND. &
1861	! data%georef%lat > -999.0) &
1862	! data%Lev2a%geop_refrac(idx) = geometric2geopotential(data%georef%lat, &
1863	! data%Lev2a%alt_refrac(idx))
1864	IF (data%georef%lat > -999.0) THEN
1865	WHERE(data%Lev2a%alt_refrac > -999.0) &
1866	data%Lev2a%geop_refrac = geometric2geopotential(data%georef%lat, &
1867	data%Lev2a%alt_refrac)
1868	ENDIF
1869
1870	CALL ncdf_getvar('Ref', data%Lev2a%refrac)
1871	data%Lev2a%units%refrac = "1"
1872
1873	IF (ncdf_isvar('Ref_stdv')) THEN
1874	CALL ncdf_getvar('Ref_stdv', data%Lev2a%refrac_sigma)
1875	data%Lev2a%units%refrac_sigma = "1"
1876
1877	! set the quality for ref
1878	CALL ncdf_getatt('_FillValue', readreal, 'Ref' )
1879	! idx => WHERE( data%Lev2a%refrac > readreal, nidx)
1880	! IF (nidx > 0) data%Lev2a%refrac_qual(idx) = 100.0
1881	WHERE( data%Lev2a%refrac > readreal) &
1882	data%Lev2a%refrac_qual = 100.0
1883	ENDIF
1884
1885	!! include dry temperature in Level 2a - if variable 'Bend_ang' exists, the
1886	!! UCAR file is an atmPrf file, and the 'Temp' variable is actually dry temperature.
1887	IF (ncdf_isvar('Temp') .AND. ncdf_isvar('Bend_ang')) THEN
1888
1889	CALL ncdf_getvar('Temp', data%Lev2a%dry_temp)
1890
1891	! idx => WHERE((data%Lev2a%dry_temp > -999.0), nidx)
1892	! IF (nidx > 0) THEN
1893	! data%Lev2a%dry_temp(idx) = data%Lev2a%dry_temp(idx) + 273.15_wp
1894	! ENDIF
1895	WHERE(data%Lev2a%dry_temp > -999.0) &
1896	data%Lev2a%dry_temp = data%Lev2a%dry_temp + 273.15_wp
1897	data%Lev2a%units%dry_temp = "kelvin"
1898
1899	END IF
1900
1901	ENDIF
1902
1903	! 6.9 Level2b variables (if any)
1904	! -------------------------------
1905
1906	IF (ncdf_isvar('MSL_alt') .AND. ncdf_isvar('Pres')) THEN
1907	CALL ncdf_getsize('MSL_alt', n, dim = 1)
1908	CALL ropp_io_init(data%Lev2b, n)
1909	ELSE
1910	data%Lev2b%Npoints = 0
1911	ENDIF
1912
1913	IF (data%Lev2b%Npoints > 0) THEN
1914
1915	CALL ncdf_getvar('Pres', data%Lev2b%press)
1916	data%Lev2b%units%press = "hPa"
1917	CALL ncdf_getvar('Temp', data%Lev2b%temp)
1918	! idx => WHERE(data%Lev2b%temp > -999.0, nidx)
1919	! IF (nidx > 0) data%Lev2b%temp(idx) = data%Lev2b%temp(idx) + 273.15_wp
1920	WHERE (data%Lev2b%temp > -999.0) &
1921	data%Lev2b%temp = data%Lev2b%temp + 273.15_wp
1922	data%Lev2b%units%temp = "kelvin"
1923
1924	IF (ncdf_isvar('Vp')) THEN
1925	CALL ncdf_getvar('Vp', data%Lev2b%shum)
1926	! idx => WHERE(data%Lev2b%shum > -999.0, nidx)
1927	! IF (nidx > 0) data%Lev2b%shum(idx) = &
1928	! (data%Lev2b%shum(idx)*epsilon_water) / &
1929	! (data%Lev2b%press(idx) - &
1930	! (data%Lev2b%shum(idx)*(1.0_wp - epsilon_water)))
1931	WHERE (data%Lev2b%shum > -999.0) &
1932	data%Lev2b%shum = &
1933	(data%Lev2b%shum*epsilon_water) / &
1934	(data%Lev2b%press - &
1935	(data%Lev2b%shum*(1.0_wp - epsilon_water)))
1936	data%Lev2b%units%shum = "kilogram/kilogram"
1937	ENDIF
1938
1939	! get the geopotential from L2a
1940	CALL ncdf_getvar('MSL_alt', data%Lev2b%geop, units = 'metres')
1941	data%Lev2b%geop = geometric2geopotential(data%georef%lat,data%Lev2b%geop)
1942
1943	! background time
1944	data%bg%Year = data%DTocc%Year
1945	data%bg%Month = data%DTocc%Month
1946	data%bg%Day = data%DTocc%Day
1947	data%bg%Hour = data%DTocc%Hour
1948	data%bg%Minute = data%DTocc%Minute
1949	nidx = INDEX(file, 'Prf_', .TRUE.)
1950	IF ( nidx > 3 ) data%bg%source = file(nidx-3:nidx-1) ! read from filename
1951
1952	ENDIF
1953
1954	! 6.10 (Global) Attributes
1955	! ------------------------
1956
1957	data%FmtVersion = ' ' ; data%FmtVersion = ThisFmtVer
1958	data%processing_centre = ' ' ; CALL ncdf_getatt('center', data%processing_centre)
1959	data%pod_method = ' ' ; data%pod_method = "UNKNOWN"
1960	data%phase_method = ' ' ; data%phase_method = "UNKNOWN"
1961	data%bangle_method = ' ' ; data%bangle_method = "UNKNOWN"
1962	data%refrac_method = ' ' ; data%refrac_method = "UNKNOWN"
1963	data%meteo_method = ' ' ; data%meteo_method = "UNKNOWN"
1964	data%thin_method = ' ' ; data%thin_method = "UNKNOWN"
1965	data%software_version = ' ' ; data%software_version = "UNKNOWN"
1966	nidx = INDEX(file, 'Prf_', .TRUE.)
1967	IF ( nidx > 0 ) & ! VNN.nnn from file name
1968	data%software_version = "V" // file(nidx+38:nidx+39) // &
1969	"." // file(nidx+32:nidx+34)
1970
1971	! COSMIC has no processing date, set to current utc date/time
1972	CALL Date_and_Time_UTC ( Values=DTnow )
1973	data%DTpro%Year = DTnow(1)
1974	data%DTpro%Month = DTnow(2)
1975	data%DTpro%Day = DTnow(3)
1976	data%DTpro%Hour = DTnow(5)
1977	data%DTpro%Minute = DTnow(6)
1978	data%DTpro%Second = DTnow(7)
1979	data%DTpro%Msec = DTnow(8)
1980
1981	! 6.11 Occultation ID
1982	! -------------------
1983
1984	CALL ropp_io_occid(DATA)
1985
1986	! 6.12 Clean up
1987	! -------------
1988
1989	! IF (ASSOCIATED(idx)) DEALLOCATE(idx)
1990
1991	END SUBROUTINE ropp_io_read_ucardata_atmPrf
1992
1993	!-------------------------------------------------------------------------------
1994	! 7. UCAR RO data - atmPhs files
1995	!-------------------------------------------------------------------------------
1996
1997	SUBROUTINE ropp_io_read_ucardata_atmPhs(DATA, file)
1998
1999	! NB this routine only supports UCAR atmPhs netCDF files,
2000	! See: http://cosmic-io.cosmic.ucar.edu/cdaac/fileFormats/atmPhs.html
2001
2002	! 7.1 Declarations
2003	! ----------------
2004
2005	USE DateTimeProgs, ONLY: Date_and_Time_UTC
2006	USE DateTimeTypes
2007	USE ropp_utils
2008	USE ncdf
2009	USE ropp_io
2010	USE ropp_io_types, ONLY: ROprof, &
2011	ThisFmtVer
2012
2013	IMPLICIT NONE
2014
2015	TYPE(ROprof), INTENT(inout) :: DATA
2016	CHARACTER(len = *), INTENT(in) :: file
2017
2018	INTEGER :: n
2019	INTEGER :: readint
2020	REAL(wp) :: readreal
2021	CHARACTER (len = 256) :: readstr
2022
2023	REAL(wp), PARAMETER :: g_wmo = 9.80665
2024	INTEGER, DIMENSION(8) :: DTnow
2025
2026	REAL(wp), DIMENSION(:), ALLOCATABLE :: workdata
2027
2028	! holds where output
2029	! INTEGER, DIMENSION(:), POINTER :: idx => null()
2030	INTEGER :: nidx
2031
2032	! 7.3 Header variables
2033	! --------------------
2034
2035	readstr = ' '
2036	CALL ncdf_getatt('fileStamp', readstr)
2037	data%leo_id = readstr(1:4)
2038	CALL ncdf_getatt('occsatId', readint)
2039	WRITE(data%gns_id,'(A1,I3.3)') 'G', readint
2040	readstr = ' '
2041	data%stn_id = ' '
2042
2043	! 7.4 Date and time
2044	! -----------------
2045
2046	CALL gettime(data%DTocc)
2047
2048	! 7.5 Overall quality
2049	! -------------------
2050
2051	CALL ncdf_getatt('bad', readstr)
2052	data%PCD = 0
2053	IF (TRIM(readstr) == "0") THEN
2054	data%PCD = 0
2055	data%overall_qual = 100.0
2056	END IF
2057	IF (TRIM(readstr) == "1") THEN
2058	data%PCD = 1 ! non nominal
2059	data%overall_qual = 0.0
2060	END IF
2061	data%units%overall_qual = "%"
2062
2063	! 7.6 Level1a variables (if any)
2064	! -------------------------------
2065
2066	IF (ncdf_isvar('time')) THEN
2067	CALL ncdf_getsize('time', n, dim = 1)
2068	CALL ropp_io_init(data%Lev1a, n)
2069	ELSE
2070	data%Lev1a%Npoints = 0
2071	ENDIF
2072
2073	IF (data%Lev1a%Npoints > 0) THEN
2074
2075	CALL ncdf_getvar('time', data%Lev1a%dtime)
2076	CALL ncdf_getvar('caL1Snr', data%Lev1a%snr_L1ca)
2077	CALL ncdf_getvar('pL1Snr', data%Lev1a%snr_L1p)
2078	CALL ncdf_getvar('pL2Snr', data%Lev1a%snr_L2p)
2079	data%Lev1a%units%snr = "volt/volt"
2080
2081	CALL ncdf_getvar('xLeo', data%Lev1a%r_leo(:,1))
2082	CALL ncdf_getvar('yLeo', data%Lev1a%r_leo(:,2))
2083	CALL ncdf_getvar('zLeo', data%Lev1a%r_leo(:,3))
2084	data%Lev1a%r_leo(:,:) = data%Lev1a%r_leo(:,:) * 1000.0_wp
2085	data%Lev1a%units%r_leo = "metres"
2086	data%Lev1a%reference_frame%r_leo = "ECI"
2087
2088	CALL ncdf_getvar('xdLeo', data%Lev1a%v_leo(:,1))
2089	CALL ncdf_getvar('ydLeo', data%Lev1a%v_leo(:,2))
2090	CALL ncdf_getvar('zdLeo', data%Lev1a%v_leo(:,3))
2091	data%Lev1a%v_leo(:,:) = data%Lev1a%v_leo(:,:) * 1000.0_wp
2092	data%Lev1a%units%v_leo = "metres / seconds"
2093
2094	CALL ncdf_getvar('xGps', data%Lev1a%r_gns(:,1))
2095	CALL ncdf_getvar('yGps', data%Lev1a%r_gns(:,2))
2096	CALL ncdf_getvar('zGps', data%Lev1a%r_gns(:,3))
2097	data%Lev1a%r_gns(:,:) = data%Lev1a%r_gns(:,:) * 1000.0_wp
2098	data%Lev1a%units%r_gns = "metres"
2099	data%Lev1a%reference_frame%r_gns = "ECI"
2100
2101	CALL ncdf_getvar('xdGps', data%Lev1a%v_gns(:,1))
2102	CALL ncdf_getvar('ydGps', data%Lev1a%v_gns(:,2))
2103	CALL ncdf_getvar('zdGps', data%Lev1a%v_gns(:,3))
2104	data%Lev1a%v_gns(:,:) = data%Lev1a%v_gns(:,:) * 1000.0_wp
2105	data%Lev1a%units%v_gns = "metres / seconds"
2106
2107	CALL ncdf_getvar('exL1', data%Lev1a%phase_L1)
2108	CALL ncdf_getvar('exL2', data%Lev1a%phase_L2)
2109	data%Lev1a%units%phase = "metres"
2110
2111	! open loop phase model data
2112	ALLOCATE(workdata(data%Lev1a%Npoints))
2113	CALL ncdf_getvar('xmdl', workdata)
2114	CALL ropp_io_addvar_rodataD1d(DATA, &
2115	name = "open_loop_lcf", &
2116	long_name= "Lost Carrier Flag", &
2117	units = "metres", &
2118	range = (/-1000000.0_wp, 1000000.0_wp/), &
2119	DATA = workdata)
2120	DEALLOCATE(workdata)
2121
2122	! set the quality for phase
2123	CALL ncdf_getatt('_FillValue', readreal, 'exL1')
2124	! idx => WHERE( data%Lev1a%phase_L1 > readreal, nidx)
2125	! IF (nidx > 0) data%Lev1a%phase_qual(idx) = 100.0
2126	WHERE (data%Lev1a%phase_L1 > readreal) &
2127	data%Lev1a%phase_qual = 100.0
2128	ENDIF
2129
2130	! 7.7 (Global) Attributes
2131	! ------------------------
2132
2133	data%FmtVersion = ' ' ; data%FmtVersion = ThisFmtVer
2134	data%processing_centre = ' ' ; CALL ncdf_getatt('center', data%processing_centre)
2135	data%pod_method = ' ' ; data%pod_method = "UNKNOWN"
2136	data%phase_method = ' ' ; data%phase_method = "UNKNOWN"
2137	data%bangle_method = ' ' ; data%bangle_method = "UNKNOWN"
2138	data%refrac_method = ' ' ; data%refrac_method = "UNKNOWN"
2139	data%meteo_method = ' ' ; data%meteo_method = "UNKNOWN"
2140	data%thin_method = ' ' ; data%thin_method = "UNKNOWN"
2141	data%software_version = ' ' ; data%software_version = "UNKNOWN"
2142	nidx = INDEX(file, 'atmPhs_', .TRUE.)
2143	IF ( nidx > 0 ) & ! VNN.nnn from file name
2144	data%software_version = "V" // file(nidx+38:nidx+39) // &
2145	"." // file(nidx+32:nidx+34)
2146
2147	! COSMIC has no processing date, set to current utc date/time
2148	CALL Date_and_Time_UTC ( Values=DTnow )
2149	data%DTpro%Year = DTnow(1)
2150	data%DTpro%Month = DTnow(2)
2151	data%DTpro%Day = DTnow(3)
2152	data%DTpro%Hour = DTnow(5)
2153	data%DTpro%Minute = DTnow(6)
2154	data%DTpro%Second = DTnow(7)
2155	data%DTpro%Msec = DTnow(8)
2156
2157	! 7.8 Occultation ID
2158	! -------------------
2159
2160	CALL ropp_io_occid(DATA)
2161
2162	! 7.9 Clean up
2163	! -------------
2164
2165	! IF (ASSOCIATED(idx)) DEALLOCATE(idx)
2166	CALL message_set_routine(routine)
2167
2168	END SUBROUTINE ropp_io_read_ucardata_atmPhs
2169
2170	SUBROUTINE gettime(DTocc)
2171
2172	USE ropp_io_types, ONLY: DT7type
2173	USE ncdf
2174
2175	IMPLICIT NONE
2176
2177	TYPE(DT7type), INTENT(inout) :: DTocc
2178	REAL :: sec
2179
2180	CALL ncdf_getatt('year', DTocc%Year)
2181	CALL ncdf_getatt('month', DTocc%Month)
2182	CALL ncdf_getatt('day', DTocc%Day)
2183	CALL ncdf_getatt('hour', DTocc%Hour)
2184	CALL ncdf_getatt('minute', DTocc%Minute)
2185	CALL ncdf_getatt('second', sec)
2186	DTocc%Second = NINT(sec)
2187	DTocc%Msec = 0
2188
2189	END SUBROUTINE gettime
2190
2191	END SUBROUTINE ropp_io_read_ncdf_get_ucardata
2192
2193	!-------------------------------------------------------------------------------
2194	! 8. EUM RO data
2195	!-------------------------------------------------------------------------------
2196
2197	SUBROUTINE ropp_io_read_ncdf_get_eumdata(DATA, file, resolution, getlevel1a, getbufr, ldummy)
2198
2199	! 8.1 Declarations
2200	! ----------------
2201
2202	USE ncdf
2203	USE ropp_utils
2204	USE ropp_io_types, ONLY: ROprof
2205
2206	IMPLICIT NONE
2207
2208	TYPE(ROprof), INTENT(inout) :: DATA
2209	CHARACTER(len = *), INTENT(in) :: file
2210	CHARACTER(len = 20), INTENT(in) :: resolution
2211	LOGICAL, INTENT(IN) :: getlevel1a
2212	LOGICAL, INTENT(IN) :: getbufr
2213	CHARACTER(len = 256) :: routine
2214	LOGICAL, INTENT(IN) :: ldummy ! needed to differentiate between
2215	! ropp_io_read_ncdf_get_eumdata and
2216	! ropp_io_read_ncdf_get_otherdata
2217
2218	! 8.2 Error handling
2219	! ------------------
2220
2221	CALL message_get_routine(routine)
2222	CALL message_set_routine('ropp_io_read_eumdata')
2223
2224	CALL ropp_io_read_eumdata(DATA, file, resolution, getlevel1a, getbufr)
2225
2226	CALL message_set_routine(routine)
2227
2228
2229	CONTAINS
2230
2231
2232	!-------------------------------------------------------------------------------
2233	! 9. EUMETSAT RO data - level 1B
2234	!-------------------------------------------------------------------------------
2235
2236	SUBROUTINE ropp_io_read_eumdata(DATA, file, resolution, getlevel1a, getbufr)
2237
2238	! This routine reads the netCDF4 EUMETSAT format into an internal ROPP structure.
2239
2240	! 9.1 Declarations
2241	! ----------------
2242
2243	USE DateTimeTypes
2244	USE ropp_utils
2245	USE coordinates
2246	USE ncdf
2247	USE ropp_io
2248	USE ropp_io_types, ONLY: ROprof, &
2249	ThisFmtVer, &
2250	PCD_open_loop, &
2251	PCD_rising
2252
2253	IMPLICIT NONE
2254
2255	TYPE(ROprof), INTENT(inout) :: DATA
2256	TYPE(ROprof) :: DATA_CL, DATA_RS
2257	CHARACTER(len = 20), INTENT(in) :: resolution
2258	CHARACTER(len = *), INTENT(in) :: file
2259	LOGICAL, INTENT(IN) :: getlevel1a
2260	LOGICAL, INTENT(IN) :: getbufr
2261
2262	REAL(wp), DIMENSION(:), ALLOCATABLE :: rs_navbit_ext ! RS navbits (external)
2263	REAL(wp), DIMENSION(:), ALLOCATABLE :: rs_navbit_int ! RS navbits (internal)
2264	REAL(wp), DIMENSION(:), ALLOCATABLE :: temparray ! Temporary array
2265	INTEGER, DIMENSION(:), ALLOCATABLE :: rs_lcf, cl_lcf, lcf ! Lost carrier flag
2266
2267	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_gns, v_gns ! Temporary pos/velocity arrays
2268	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_leo, v_leo ! Temporary pos/velocity arrays
2269
2270	REAL(wp), DIMENSION(:), ALLOCATABLE :: rs_i_ca_uncorr ! RS I component [V]
2271	REAL(wp), DIMENSION(:), ALLOCATABLE :: rs_q_ca_uncorr ! RS Q component [V]
2272	REAL(wp), DIMENSION(:), ALLOCATABLE :: rs_exphase_l1_nco ! RS NCO excess phase [m]
2273	REAL(wp), DIMENSION(:), ALLOCATABLE :: rs_phase_l1_iq ! I/Q contribution
2274	Integer, DIMENSION(:), ALLOCATABLE :: tracking_state ! Tracking state
2275
2276	INTEGER :: n, n_cl, n_rs, first_valid, i, j
2277	INTEGER :: readint, readint2, readint3
2278	REAL(EightByteReal) :: readreal, readreal2, ts, ts1
2279	CHARACTER (len = 256) :: readstr, readstr2
2280
2281	INTEGER(OneByteInt) :: readbyte1, readbyte2, readbyte3
2282
2283	CHARACTER (len = 256) :: sdir ! directory to science info in netCDF4
2284	CHARACTER (len = 256) :: ddir ! directory to science/data in netCDF4
2285	CHARACTER (len = 256) :: tdir ! tmp directory for use
2286
2287	REAL(wp), PARAMETER :: g_wmo = 9.80665_wp
2288	REAL(wp), PARAMETER :: epsilon_water = 0.621971_wp
2289	REAL(wp), PARAMETER :: f_L1 = 1.57542e9_wp
2290	REAL(wp), PARAMETER :: c_light = 299792458.0_wp
2291	! Pi already defined as parameter in coordinates module. Confuses pgf95.
2292	REAL(wp), PARAMETER :: pi1 = 3.141592653589793238_wp
2293
2294	INTEGER :: have_nb = 0
2295
2296	! INTEGER, DIMENSION(:), POINTER :: idx => null() ! Holds 'where' output
2297	! INTEGER :: nidx
2298
2299	! 9.2 Group path for science, level 1B data
2300	! -----------------------------------------
2301
2302	sdir = '/data/'
2303
2304	ddir = TRIM(sdir) // 'level_1b/' // TRIM(resolution) // '/'
2305
2306	! 9.3 Header variables
2307	! --------------------
2308
2309	! receiver satellite ID in ROPP format
2310	readstr = ' '
2311	CALL ncdf_getatt('spacecraft', readstr)
2312	IF ( readstr == 'M02' ) data%leo_id = 'META'
2313	IF ( readstr == 'M01' ) data%leo_id = 'METB'
2314	IF ( readstr == 'M03' ) data%leo_id = 'METC'
2315
2316	! occulting GNSS satellite
2317	CALL ncdf_getvar(TRIM(sdir)//'occultation/prn', readint)
2318	WRITE(data%gns_id,'(A1,I3.3)') 'G', readint
2319
2320	! station ID unused for EUMETSAT data
2321	! readstr = ' '
2322	! CALL ncdf_getatt('fiducial_id', readstr)
2323	! IF(readstr /= ' ') data%stn_id = readstr(1:4)
2324
2325	! 9.4 Overall quality
2326	! -------------------
2327
2328	! Overall nominal / non-nominal
2329	tdir = '/quality/'
2330	CALL ncdf_getvar(TRIM(tdir)//'overall_quality_ok', readbyte1)
2331
2332	data%PCD = 0
2333	IF (readbyte1 == 0) data%PCD = 1
2334
2335	! BA nominal (set to same value as overall)
2336	IF (readbyte1 == 0) data%PCD = IBSET(data%PCD, PCD_bangle) ! set to same value as overall
2337
2338	! quality indicator, currently just set to 100% or 0%
2339	data%overall_qual = 100.0_wp * readbyte1
2340	data%units%overall_qual = '%'
2341
2342	! open loop data used?
2343	CALL ncdf_getvar(TRIM(tdir)//'ol_data_used', readbyte1)
2344	IF (readbyte1 /= 0) data%PCD = IBSET(data%PCD, PCD_open_loop)
2345
2346	! setting or rising
2347	CALL ncdf_getatt(TRIM(sdir)//'occultation/occultation_type', readstr)
2348	IF (TRIM(readstr) == 'rising') data%PCD = IBSET(data%PCD, PCD_rising)
2349
2350	! closed loop phase measurements okay?
2351	CALL ncdf_getvar(TRIM(tdir)//'cl_snr_ca_ok', readbyte1)
2352	CALL ncdf_getvar(TRIM(tdir)//'cl_snr_p1_ok', readbyte2)
2353	CALL ncdf_getvar(TRIM(tdir)//'cl_snr_p2_ok', readbyte3)
2354	IF ((readbyte1 == 0) .AND. (readbyte2 == 0) .AND. (readbyte3 == 0)) &
2355	data%PCD = IBSET(data%PCD, PCD_phase)
2356
2357	! NRT or offline processing?
2358	CALL ncdf_getatt('environment', readstr)
2359	IF (TRIM(readstr) == 'Operational' .OR. &
2360	TRIM(readstr) == 'Validation') THEN
2361	data%PCD = IBCLR(data%PCD, PCD_offline)
2362	ELSE
2363	data%PCD = IBSET(data%PCD, PCD_offline)
2364	ENDIF
2365
2366	! 9.5 Date and time
2367	! -----------------
2368
2369	! start time of occultation
2370	tdir = TRIM(sdir) // 'level_1b/'
2371
2372	CALL ncdf_getvar(TRIM(tdir)//'utc_start_absdate', readint)
2373	CALL ncdf_getvar(TRIM(tdir)//'utc_start_abstime', readreal)
2374
2375	! check that we have the right time
2376	CALL ncdf_getatt(TRIM(tdir)//'units', readstr, 'utc_start_absdate')
2377	IF ( readstr == 'days since 2000-01-01 00:00:00.00' ) THEN
2378	CALL abstimetoDT(readint, readreal, data%DTocc)
2379	ELSE
2380	CALL message(msg_fatal, &
2381	'Time units found for utc_start_absdate are inconsistent.')
2382	ENDIF
2383
2384	! difference between start time and georef time
2385	tdir = TRIM(sdir) // 'occultation/'
2386	CALL ncdf_getvar(TRIM(tdir)//'utc_georef_absdate', readint2)
2387	CALL ncdf_getvar(TRIM(tdir)//'utc_georef_abstime', readreal2)
2388
2389	! check that we have the right time
2390	CALL ncdf_getatt(TRIM(tdir)//'units', readstr, 'utc_georef_absdate')
2391	IF ( .NOT. (readstr == 'days since 2000-01-01 00:00:00.00') ) &
2392	CALL message(msg_fatal, &
2393	'Time units found for utc_georef_absdate are inconsistent.')
2394
2395	data%georef%time_offset = 86400.d0*(readint2 - readint) + (readreal2 - readreal)
2396
2397	! 9.6 Georeferencing
2398	! ------------------
2399
2400	CALL ncdf_getvar(TRIM(tdir)//'latitude', data%georef%lat, &
2401	units=data%georef%units%lat)
2402	CALL ncdf_getvar(TRIM(tdir)//'longitude', data%georef%lon, &
2403	units=data%georef%units%lon)
2404	CALL ncdf_getvar(TRIM(tdir)//'undulation', data%georef%Undulation, &
2405	units=data%georef%units%Undulation)
2406	CALL ncdf_getvar(TRIM(tdir)//'r_curve', data%georef%roc, &
2407	units=data%georef%units%roc)
2408	CALL ncdf_getvar(TRIM(tdir)//'r_curve_centre_fixed', data%georef%r_coc, &
2409	units=data%georef%units%r_coc)
2410	CALL ncdf_getvar(TRIM(tdir)//'azimuth_north', data%georef%azimuth, &
2411	units=data%georef%units%azimuth)
2412
2413	! FIXME: Do we need this?
2414	data%georef%reference_frame%r_coc = 'ECF'
2415
2416	! get the location, velocity of leo and gnss, assure that the getbufr flag
2417	! is set, not the getlevel1a. The getbufr only allocates one level 1a for
2418	! this info and reads the value from the correct group occultation
2419	! FIXME: Is this always done like this? What about when the level 1a has
2420	! a dimension > 1, what is actually chosen for the one bufr value?
2421	! It would be better to store the reference values for e.g. the bufr
2422	! in a dedicated dimension, not in the level 1a one! This dimension
2423	! would always have the value 1 when used, but allows to separate data
2424	! at the reference point from the general level 1a data, which covers the
2425	! whole occultation
2426
2427	IF (getbufr .AND. getlevel1a) THEN
2428
2429	CALL message(msg_fatal, &
2430	'Extraction of both BUFR and level1a data is not allowed.')
2431
2432	ELSE
2433
2434	IF (getbufr) THEN
2435
2436	CALL ropp_io_init(data%Lev1a, 1)
2437	CALL ncdf_getvar(TRIM(tdir)//'position_rec_fixed', data%Lev1a%r_leo(1,:), &
2438	units=data%Lev1a%units%r_leo, range = data%Lev1a%range%r_leo)
2439	CALL ncdf_getvar(TRIM(tdir)//'velocity_rec', data%Lev1a%v_leo(1,:), &
2440	units=data%Lev1a%units%v_leo, range = data%Lev1a%range%v_leo)
2441	CALL ncdf_getvar(TRIM(tdir)//'position_gns_fixed', data%Lev1a%r_gns(1,:), &
2442	units=data%Lev1a%units%r_gns, range = data%Lev1a%range%r_gns)
2443	CALL ncdf_getvar(TRIM(tdir)//'velocity_gns', data%Lev1a%v_gns(1,:), &
2444	units=data%Lev1a%units%v_gns, range = data%Lev1a%range%v_gns)
2445
2446	! FIXME: Do we need this?
2447	data%Lev1a%reference_frame%r_leo = 'ECF'
2448	data%Lev1a%reference_frame%r_gns = 'ECF'
2449	data%Lev1a%reference_frame%v_leo = 'ECI'
2450	data%Lev1a%reference_frame%v_leo = 'ECI'
2451
2452	! dtime set to zero for only 1 level 1A field
2453	data%Lev1a%dtime = 0.0_wp
2454
2455	! mean SNR values for CA, P1, P2
2456	! FIXME: - this is the EUMETSAT mean value for SLTA > 60km,
2457	! should be mentioned in ROPP netCDF file if required!
2458	! - data appears when doing a bufr2ropp, however seems not to be
2459	! part of the bufr content.
2460	!tdir = TRIM(sdir)//'level_1a/closed_loop/'
2461	!CALL ncdf_getvar(TRIM(tdir)//'snr_ca_mean', data%Lev1a%snr_L1ca(1), &
2462	! units=data%Lev1a%units%snr, range = data%Lev1a%range%snr)
2463	!CALL ncdf_getvar(TRIM(tdir)//'snr_p1_mean', data%Lev1a%snr_L1p(1), &
2464	! units=data%Lev1a%units%snr, range = data%Lev1a%range%snr)
2465	!CALL ncdf_getvar(TRIM(tdir)//'snr_p2_mean', data%Lev1a%snr_L2p(1), &
2466	! units=data%Lev1a%units%snr, range = data%Lev1a%range%snr)
2467
2468	! data are available
2469	data%Lev1a%Missing = .FALSE.
2470
2471	ENDIF ! getbufr
2472
2473	ENDIF ! getbufr .AND. getlevel1a not both true
2474
2475
2476	! 9.7 Level1a variables (if any)
2477	! ------------------------------
2478
2479	! 9.7.1 Closed Loop Level1a variables (if requested)
2480	! --------------------------------------------------
2481
2482	IF (getlevel1a) THEN
2483
2484	! FIXME: there should be a switch, similar to the level 1b resolution, to
2485	! determine what lev1a data to read (cl only, cl+rs, cl+ol). For now the
2486	! approach is to read and combine closed loop and raw sampling records.
2487	! This and related stuff is Kjartan's update.
2488
2489	tdir = TRIM(sdir) // 'level_1a/closed_loop/'
2490
2491	IF (ncdf_isvar(TRIM(tdir)//'dtime')) THEN ! read CL data
2492
2493	CALL ncdf_getsize(TRIM(tdir)//'dtime', n_cl, dim = 1)
2494	CALL ropp_io_init(data_cl%Lev1a, n_cl)
2495
2496	! Time
2497
2498	CALL ncdf_getvar(TRIM(tdir)//'dtime', data_cl%Lev1a%dtime, units=data_cl%Lev1a%units%dtime)
2499
2500	! Position and velocity
2501
2502	ALLOCATE(r_gns(3,n_cl), v_gns(3,n_cl), r_leo(3,n_cl), v_leo(3,n_cl))
2503
2504	CALL ncdf_getvar(TRIM(tdir)//'r_transmitter', r_gns, &
2505	units=data_cl%Lev1a%units%r_gns)
2506	CALL ncdf_getvar(TRIM(tdir)//'v_transmitter', v_gns, &
2507	units=data_cl%Lev1a%units%v_gns)
2508	CALL ncdf_getvar(TRIM(tdir)//'r_receiver', r_leo, &
2509	units=data_cl%Lev1a%units%r_leo)
2510	CALL ncdf_getvar(TRIM(tdir)//'v_receiver', v_leo, &
2511	units=data_cl%Lev1a%units%v_leo)
2512
2513	data_cl%Lev1a%r_gns = TRANSPOSE(r_gns)
2514	data_cl%Lev1a%v_gns = TRANSPOSE(v_gns)
2515	data_cl%Lev1a%r_leo = TRANSPOSE(r_leo)
2516	data_cl%Lev1a%v_leo = TRANSPOSE(v_leo)
2517
2518	DEALLOCATE(r_gns, v_gns, r_leo, v_leo)
2519
2520	! Phase and amplitude
2521
2522	CALL ncdf_getvar(TRIM(tdir)//'exphase_ca', data_cl%Lev1a%phase_L1, &
2523	units=data_cl%Lev1a%units%phase)
2524	CALL ncdf_getvar(TRIM(tdir)//'exphase_p2', data_cl%Lev1a%phase_L2, &
2525	units=data_cl%Lev1a%units%phase)
2526	CALL ncdf_getvar(TRIM(tdir)//'snr_ca', data_cl%Lev1a%snr_L1ca, &
2527	units=data_cl%Lev1a%units%snr)
2528	CALL ncdf_getvar(TRIM(tdir)//'snr_p1', data_cl%Lev1a%snr_L1p, &
2529	units=data_cl%Lev1a%units%snr)
2530	CALL ncdf_getvar(TRIM(tdir)//'snr_p2', data_cl%Lev1a%snr_L2p, &
2531	units=data_cl%Lev1a%units%snr)
2532
2533	! Closed loop lost carrier flag (looking for data gaps)
2534
2535	ALLOCATE(cl_LCF(n_cl))
2536	cl_LCF(:) = 0
2537
2538	ts = MINVAL(ABS(data_cl%Lev1a%dtime(2:n_cl) - data_cl%Lev1a%dtime(1:n_cl-1)))
2539
2540	IF (BTEST(data_cl%PCD, PCD_rising)) THEN ! Rising occultation
2541	DO j=n_cl-1,1,-1
2542	ts1 = ABS(data_cl%Lev1a%dtime(j) - data_cl%Lev1a%dtime(j+1))
2543	IF (ts1 > 1.05*ts) THEN
2544	cl_LCF(j) = IBSET(cl_LCF(j), 3)
2545	cl_LCF(j+1) = IBSET(cl_LCF(j+1),3)
2546	ENDIF
2547	ENDDO
2548	ELSE ! Setting occultation
2549	DO j=2,n_cl
2550	ts1 = ABS(data_cl%Lev1a%dtime(j) - data_cl%Lev1a%dtime(j-1))
2551	IF (ts1 > 1.05*ts) THEN
2552	cl_LCF(j) = IBSET(cl_LCF(j), 3)
2553	cl_LCF(j-1) = IBSET(cl_LCF(j-1),3)
2554	ENDIF
2555	ENDDO
2556
2557	ENDIF ! Rising occultation
2558
2559	ENDIF ! /level_1a/closed_loop/dtime = true
2560
2561	! 9.7.2 Raw Sampling Level1a variables (if requested)
2562	! ---------------------------------------------------
2563
2564	n_rs = 0 ! will be set differently if RS data is available
2565
2566	CALL ncdf_getvar('/quality/rs_data_available', readbyte1)
2567
2568	IF (readbyte1 == 1) THEN
2569
2570	tdir = TRIM(sdir) // 'level_1a/raw_sampling/'
2571
2572	IF (ncdf_isvar(TRIM(tdir)//'dtime')) THEN ! read RS data
2573
2574	data%PCD = IBSET(data%PCD, PCD_open_loop)
2575	CALL ncdf_getsize(TRIM(tdir)//'dtime', n_rs, dim = 1)
2576	CALL ropp_io_init(data_rs%Lev1a, n_rs)
2577
2578	! Time
2579
2580	CALL ncdf_getvar(TRIM(tdir)//'dtime', data_rs%Lev1a%dtime, units=data_rs%Lev1a%units%dtime)
2581
2582	! Position and velocity
2583
2584	ALLOCATE(r_leo(3,n_rs), v_leo(3,n_rs), r_gns(3,n_rs), v_gns(3,n_rs))
2585
2586	CALL ncdf_getvar(TRIM(tdir)//'r_transmitter', r_gns, &
2587	units=data_rs%Lev1a%units%r_gns)
2588	CALL ncdf_getvar(TRIM(tdir)//'v_transmitter', v_gns, &
2589	units=data_rs%Lev1a%units%v_gns)
2590	CALL ncdf_getvar(TRIM(tdir)//'r_receiver', r_leo, &
2591	units=data_rs%Lev1a%units%r_leo)
2592	CALL ncdf_getvar(TRIM(tdir)//'v_receiver', v_leo, &
2593	units=data_rs%Lev1a%units%v_leo)
2594
2595	data_rs%Lev1a%r_gns = TRANSPOSE(r_gns)
2596	data_rs%Lev1a%v_gns = TRANSPOSE(v_gns)
2597	data_rs%Lev1a%r_leo = TRANSPOSE(r_leo)
2598	data_rs%Lev1a%v_leo = TRANSPOSE(v_leo)
2599
2600	DEALLOCATE(r_gns, v_gns, r_leo, v_leo)
2601
2602	! Phase and amplitude
2603
2604	ALLOCATE(rs_navbit_int(n_rs))
2605	ALLOCATE(rs_navbit_ext(n_rs))
2606	ALLOCATE(tracking_state(n_rs))
2607
2608	CALL ncdf_getvar(TRIM(tdir)//'navbits_internal', rs_navbit_int)
2609	CALL ncdf_getvar('/quality/rs_external_navbits_applied', have_nb)
2610	CALL ncdf_getvar(TRIM(tdir)//'tracking_state', tracking_state)
2611
2612	IF (have_nb > 0) THEN
2613	CALL ncdf_getvar(TRIM(tdir)//'navbits_external', rs_navbit_ext)
2614	ELSE
2615	rs_navbit_ext(:) = rs_navbit_int(:)
2616	ENDIF
2617
2618	ALLOCATE(rs_phase_l1_iq(n_rs))
2619	ALLOCATE(rs_i_ca_uncorr(n_rs))
2620	ALLOCATE(rs_q_ca_uncorr(n_rs))
2621	ALLOCATE(rs_exphase_l1_nco(n_rs))
2622
2623	CALL ncdf_getvar(TRIM(tdir)//'i_ca_uncorr', rs_i_ca_uncorr)
2624	CALL ncdf_getvar(TRIM(tdir)//'q_ca_uncorr', rs_q_ca_uncorr)
2625	CALL ncdf_getvar(TRIM(tdir)//'exphase_l1_nco', rs_exphase_l1_nco, &
2626	units=data_rs%Lev1a%units%phase)
2627	CALL ncdf_getvar(TRIM(tdir)//'snr_ca', data_rs%Lev1a%snr_L1ca, &
2628	units=data_rs%Lev1a%units%snr)
2629
2630	DO j=1,n_rs
2631	rs_phase_l1_iq(j) = ATAN2(rs_q_ca_uncorr(j), rs_i_ca_uncorr(j)+TINY(1.0_wp))
2632	ENDDO
2633
2634	CALL Accumulate_Phase(rs_phase_l1_iq)
2635
2636	data_rs%lev1a%phase_l1(:) = rs_exphase_l1_nco(:) + &
2637	(c_light/(2.0_wppi1f_L1))*rs_phase_l1_iq(:)
2638
2639	! Comment in the line below if you wish to read the excess phase directly
2640	! from the EUMETSAT file instead of deriving from I and Q
2641	! CALL ncdf_getvar(TRIM(tdir)//'exphase_ca', &
2642	! data_rs%Lev1a%phase_L1, units=data_rs%Lev1a%units%phase)
2643
2644	! For rising profiles remove all data with tracking state = 2 in the
2645	! beginning of the record as these are not valid.
2646
2647	IF (BTEST(data_rs%PCD, PCD_rising)) THEN ! Rising occultation
2648
2649	first_valid = -1 ! value to mark that this hasn't been set yet
2650
2651	DO j=1,n_rs
2652	IF ((tracking_state(j) /= 2) .AND. (first_valid == -1)) THEN
2653	first_valid = j
2654	ENDIF
2655	ENDDO
2656
2657	CALL ropp_io_shrink(data_rs%Lev1a, first_valid, n_rs, 1)
2658
2659	ALLOCATE(temparray(data_rs%Lev1a%Npoints))
2660
2661	temparray = rs_navbit_int(first_valid:n_rs)
2662
2663	DEALLOCATE(rs_navbit_int)
2664	ALLOCATE(rs_navbit_int(data_rs%Lev1a%Npoints))
2665	rs_navbit_int = temparray
2666
2667	temparray = rs_navbit_ext(first_valid:n_rs)
2668	DEALLOCATE(rs_navbit_ext)
2669	ALLOCATE(rs_navbit_ext(data_rs%Lev1a%Npoints))
2670	rs_navbit_ext = temparray
2671
2672	DEALLOCATE(temparray)
2673
2674	n_rs = data_rs%Lev1a%Npoints
2675
2676	ENDIF ! Rising occultation
2677
2678	! data flag
2679
2680	ALLOCATE(rs_LCF(n_rs))
2681	rs_LCF(:) = 0
2682
2683	ts = MINVAL(ABS(data_rs%Lev1a%dtime(2:n_rs) - data_rs%Lev1a%dtime(1:n_rs-1)))
2684
2685	IF (BTEST(data_rs%PCD, PCD_rising)) THEN ! Rising occultation
2686	DO j=n_rs-1,1,-1
2687	ts1 = ABS(data_rs%Lev1a%dtime(j) - data_rs%Lev1a%dtime(j+1))
2688	IF (ts1 > 1.05*ts) THEN
2689	rs_LCF(j) = IBSET(rs_LCF(j), 3)
2690	rs_LCF(j+1) = IBSET(rs_LCF(j+1),3)
2691	ENDIF
2692	ENDDO
2693	ELSE ! Setting occultation
2694	DO j=2,n_rs
2695	ts1 = ABS(data_rs%Lev1a%dtime(j) - data_rs%Lev1a%dtime(j-1))
2696	IF (ts1 > 1.05*ts) THEN
2697	rs_LCF(j) = IBSET(rs_LCF(j), 3)
2698	rs_LCF(j-1) = IBSET(rs_LCF(j-1),3)
2699	ENDIF
2700	ENDDO
2701	ENDIF
2702
2703	rs_LCF(:) = IBSET(rs_LCF(:), 0) ! Open loop mode
2704
2705	WHERE (NINT(rs_navbit_ext(:)) == 1) ! External navbit
2706	rs_LCF(:) = IBSET(rs_LCF(:), 1)
2707	ENDWHERE
2708
2709	rs_LCF(:) = IBSET(rs_LCF(:), 2) ! Navbit quality OK
2710
2711	WHERE (NINT(rs_navbit_int(:)) == 1) ! Alternative navbit
2712	rs_LCF(:) = IBSET(rs_LCF(:), 4)
2713	ENDWHERE
2714
2715	rs_LCF(:) = IBSET(rs_LCF(:), 5) ! Alternative navbit quality
2716
2717	! Setting flag for duplicated CL record
2718	WHERE((data_rs%Lev1a%dtime(1) < data_cl%Lev1a%dtime(:)) .AND. &
2719	(data_cl%Lev1a%dtime(:) < data_rs%Lev1a%dtime(n_rs)))
2720	cl_LCF(:) = IBSET(cl_LCF(:), 6)
2721	ENDWHERE
2722
2723	ENDIF ! /level_1a/raw_sampling/dtime exists
2724
2725	ENDIF ! /quality_control/rs_data_available is true
2726
2727	! 9.7.3 Combine CL and RS Level1a variables (if requested)
2728	! ---------------------------------------------------------
2729
2730	n = n_cl + n_rs
2731
2732	! Initialise ro_data structure variables
2733
2734	CALL ropp_io_init(data%Lev1a, n)
2735	ALLOCATE(lcf(n))
2736
2737	! Set output variables
2738
2739	IF (BTEST(data%PCD, PCD_rising)) THEN ! Rising occultation
2740
2741	DO j=1,n_cl
2742	data%lev1a%dtime(n_rs+j) = data_cl%lev1a%dtime(j)
2743	data%lev1a%r_gns(n_rs+j,:) = data_cl%lev1a%r_gns(j,:)
2744	data%lev1a%v_gns(n_rs+j,:) = data_cl%lev1a%v_gns(j,:)
2745	data%lev1a%r_leo(n_rs+j,:) = data_cl%lev1a%r_leo(j,:)
2746	data%lev1a%v_leo(n_rs+j,:) = data_cl%lev1a%v_leo(j,:)
2747	data%lev1a%snr_L1ca(n_rs+j) = data_cl%lev1a%snr_L1ca(j)
2748	data%lev1a%snr_L1p(n_rs+j) = data_cl%lev1a%snr_L1p(j)
2749	data%lev1a%snr_L2p(n_rs+j) = data_cl%lev1a%snr_L2p(j)
2750	data%lev1a%phase_L1(n_rs+j) = data_cl%lev1a%phase_L1(j)
2751	data%lev1a%phase_L2(n_rs+j) = data_cl%lev1a%phase_L2(j)
2752	lcf(n_rs+j) = cl_LCF(j)
2753	ENDDO
2754
2755	IF (n_rs /= 0) THEN
2756
2757	DO j=1,n_rs
2758	data%lev1a%dtime(j) = data_rs%lev1a%dtime(j)
2759	data%lev1a%r_gns(j,:) = data_rs%lev1a%r_gns(j,:)
2760	data%lev1a%v_gns(j,:) = data_rs%lev1a%v_gns(j,:)
2761	data%lev1a%r_leo(j,:) = data_rs%lev1a%r_leo(j,:)
2762	data%lev1a%v_leo(j,:) = data_rs%lev1a%v_leo(j,:)
2763	data%lev1a%snr_L1ca(j) = data_rs%lev1a%snr_L1ca(j)
2764	data%lev1a%snr_L1p(j) = ropp_MDFV
2765	data%lev1a%snr_L2p(j) = ropp_MDFV
2766	data%lev1a%phase_L1(j) = data_rs%lev1a%phase_L1(j)
2767	data%lev1a%phase_L2(j) = ropp_MDFV
2768	lcf(j) = rs_LCF(j)
2769	ENDDO
2770
2771	! test for gap between cl and rs records
2772	ts = MINVAL(ABS(data_cl%lev1a%dtime(2:n_cl) - data_cl%lev1a%dtime(1:n_cl-1)))
2773
2774	IF (data_rs%lev1a%dtime(n_rs) < data_cl%lev1a%dtime(1)-1.5*ts) THEN
2775	LCF(n_rs:n_rs+1) = IBSET(LCF(n_rs:n_rs+1),3)
2776	ENDIF
2777
2778	ENDIF ! n_rs /= 0
2779
2780	ELSE ! Setting occultation
2781
2782	DO j=1,n_cl
2783	data%lev1a%dtime(j) = data_cl%lev1a%dtime(j)
2784	data%lev1a%r_gns(j,:) = data_cl%lev1a%r_gns(j,:)
2785	data%lev1a%v_gns(j,:) = data_cl%lev1a%v_gns(j,:)
2786	data%lev1a%r_leo(j,:) = data_cl%lev1a%r_leo(j,:)
2787	data%lev1a%v_leo(j,:) = data_cl%lev1a%v_leo(j,:)
2788	data%lev1a%snr_L1ca(j) = data_cl%lev1a%snr_L1ca(j)
2789	data%lev1a%snr_L1p(j) = data_cl%lev1a%snr_L1p(j)
2790	data%lev1a%snr_L2p(j) = data_cl%lev1a%snr_L2p(j)
2791	data%lev1a%phase_L1(j) = data_cl%lev1a%phase_L1(j)
2792	data%lev1a%phase_L2(j) = data_cl%lev1a%phase_L2(j)
2793	lcf(j) = cl_LCF(j)
2794	ENDDO
2795
2796	IF (n_rs /= 0) THEN
2797
2798	DO j=1,n_rs
2799	data%lev1a%dtime(n_cl+j) = data_rs%lev1a%dtime(j)
2800	data%lev1a%r_gns(n_cl+j,:) = data_rs%lev1a%r_gns(j,:)
2801	data%lev1a%v_gns(n_cl+j,:) = data_rs%lev1a%v_gns(j,:)
2802	data%lev1a%r_leo(n_cl+j,:) = data_rs%lev1a%r_leo(j,:)
2803	data%lev1a%v_leo(n_cl+j,:) = data_rs%lev1a%v_leo(j,:)
2804	data%lev1a%snr_L1ca(n_cl+j) = data_rs%lev1a%snr_L1ca(j)
2805	data%lev1a%snr_L1p(n_cl+j) = ropp_MDFV
2806	data%lev1a%snr_L2p(n_cl+j) = ropp_MDFV
2807	data%lev1a%phase_L1(n_cl+j) = data_rs%lev1a%phase_L1(j)
2808	data%lev1a%phase_L2(n_cl+j) = ropp_MDFV
2809	lcf(n_cl+j) = rs_LCF(j)
2810	ENDDO
2811
2812	ts = MINVAL(ABS(data_cl%lev1a%dtime(2:n_cl) - data_cl%lev1a%dtime(1:n_cl-1)))
2813
2814	IF (data_cl%lev1a%dtime(n_cl) < data_rs%lev1a%dtime(1)-1.5*ts) THEN
2815	LCF(n_cl:n_cl+1) = IBSET(LCF(n_cl:n_cl+1),3)
2816	ENDIF
2817
2818	ENDIF ! n_rs /= 0
2819
2820	ENDIF ! Rising occultation
2821
2822	! 9.7.4 Missing/invalid data checks
2823	! ---------------------------------
2824
2825	WHERE (ropp_io_isnan(data%Lev1a%phase_L1))
2826	data%Lev1a%phase_L1 = ropp_MDFV
2827	LCF = IBSET(LCF, 3)
2828	ENDWHERE
2829
2830	WHERE (ropp_io_isnan(data%Lev1a%phase_L2))
2831	data%Lev1a%phase_L2 = ropp_MDFV
2832	ENDWHERE
2833
2834	WHERE (ropp_io_isnan(data%Lev1a%snr_L1ca))
2835	data%Lev1a%snr_L1ca = ropp_MDFV
2836	LCF = IBSET(LCF, 3)
2837	ENDWHERE
2838
2839	WHERE (ropp_io_isnan(data%Lev1a%snr_L1p))
2840	data%Lev1a%snr_L1p = ropp_MDFV
2841	ENDWHERE
2842
2843	WHERE (ropp_io_isnan(data%Lev1a%snr_L2p))
2844	data%Lev1a%snr_L2p = ropp_MDFV
2845	ENDWHERE
2846
2847	data%Lev1a%reference_frame%r_gns = 'ECI'
2848	data%Lev1a%reference_frame%r_leo = 'ECI'
2849
2850	data%Lev1a%range%phase = (/ &
2851	MIN(MINVAL(data%Lev1a%phase_L1), data%Lev1a%range%phase(1)), &
2852	MAX(MAXVAL(data%Lev1a%phase_L1), data%Lev1a%range%phase(2)) /)
2853
2854	! Add lost carrier information to file
2855
2856	CALL ropp_io_addvar_rodataD1d(data, &
2857	name = 'open_loop_lcf', &
2858	long_name= 'Lost Carrier Flag', &
2859	units = '', &
2860	range = (/-1000000.0_wp, 1000000.0_wp/), &
2861	DATA = REAL(LCF,wp) )
2862
2863	DEALLOCATE(LCF, CL_LCF)
2864
2865	CALL ropp_io_free(data_cl)
2866
2867	IF (n_rs /= 0) THEN
2868	DEALLOCATE(RS_LCF)
2869	CALL ropp_io_free(data_rs)
2870	ENDIF
2871
2872	ENDIF ! getlevel1a is true
2873
2874
2875	! 9.8 Level1b variables (if any)
2876	! ------------------------------
2877
2878	IF (ncdf_isvar(TRIM(ddir)//'impact')) THEN
2879	CALL ncdf_getsize(TRIM(ddir)//'impact', n, dim = 1)
2880	CALL ropp_io_init(data%Lev1b, n)
2881	ELSE
2882	data%Lev1b%Npoints = 0
2883	ENDIF
2884
2885	IF (data%Lev1b%Npoints > 0) THEN
2886
2887	CALL ncdf_getvar(TRIM(ddir)//'lat_tp', data%Lev1b%lat_tp, &
2888	units=data%Lev1b%units%lat_tp)
2889	CALL ncdf_getvar(TRIM(ddir)//'lon_tp', data%Lev1b%lon_tp, &
2890	units=data%Lev1b%units%lon_tp)
2891	CALL ncdf_getvar(TRIM(ddir)//'azimuth_tp', data%Lev1b%azimuth_tp, &
2892	units=data%Lev1b%units%azimuth_tp)
2893	CALL ncdf_getvar(TRIM(ddir)//'impact', data%Lev1b%impact, &
2894	units=data%Lev1b%units%impact)
2895	CALL ncdf_getvar(TRIM(ddir)//'impact', data%Lev1b%impact_L1, &
2896	units=data%Lev1b%units%impact)
2897	CALL ncdf_getvar(TRIM(ddir)//'impact', data%Lev1b%impact_L2, &
2898	units=data%Lev1b%units%impact)
2899	CALL ncdf_getvar(TRIM(ddir)//'bangle', data%Lev1b%bangle, &
2900	units=data%Lev1b%units%bangle)
2901	CALL ncdf_getvar(TRIM(ddir)//'bangle_ca', data%Lev1b%bangle_L1, &
2902	units=data%Lev1b%units%bangle)
2903	CALL ncdf_getvar(TRIM(ddir)//'bangle_p2', data%Lev1b%bangle_L2, &
2904	units=data%Lev1b%units%bangle)
2905
2906	! remove all NaN from EUM fields
2907	! FIXME: maybe do this in ncdf_getvar?
2908
2909	WHERE( .NOT. ropp_io_isnan(data%Lev1b%bangle) ) data%Lev1b%bangle_qual = 100.0_wp
2910
2911	WHERE( ropp_io_isnan(data%Lev1b%lat_tp) ) data%Lev1b%lat_tp = ropp_MDFV
2912	WHERE( ropp_io_isnan(data%Lev1b%lon_tp) ) data%Lev1b%lon_tp = ropp_MDFV
2913	WHERE( ropp_io_isnan(data%Lev1b%azimuth_tp) ) data%Lev1b%azimuth_tp = ropp_MDFV
2914	WHERE( ropp_io_isnan(data%Lev1b%bangle) ) data%Lev1b%bangle = ropp_MDFV
2915	WHERE( ropp_io_isnan(data%Lev1b%bangle_L1) ) data%Lev1b%bangle_L1 = ropp_MDFV
2916	WHERE( ropp_io_isnan(data%Lev1b%bangle_L2) ) data%Lev1b%bangle_L2 = ropp_MDFV
2917
2918	! set the quality for bangle, EUMETSAT data uses NaN for missing
2919
2920	data%Lev1b%Missing = .FALSE.
2921
2922	ENDIF ! data%Lev1b%Npoints > 0
2923
2924
2925	! 9.9 (Global) Attributes
2926	! -----------------------
2927
2928	data%FmtVersion = ThisFmtVer
2929
2930	data%processing_centre = ' '
2931	CALL ncdf_getatt('institution', data%processing_centre)
2932
2933	data%pod_method = ' '
2934	CALL ncdf_getatt(TRIM(sdir)//'occultation/pod_method', data%pod_method)
2935
2936	data%phase_method = ' '
2937	CALL ncdf_getatt(TRIM(sdir)//'occultation/phase_method', data%phase_method)
2938
2939	data%bangle_method = ' '
2940	CALL ncdf_getatt(TRIM(sdir)//'occultation/retrieval_method', data%bangle_method)
2941
2942	data%refrac_method = 'UNKNOWN'
2943
2944	data%meteo_method = 'UNKNOWN'
2945
2946	IF (TRIM(resolution) == 'high_resolution') THEN
2947	data%thin_method = 'Unthinned data'
2948	ELSE
2949	data%thin_method = ' '
2950	CALL ncdf_getatt(TRIM(ddir)//'thinner_method', data%thin_method)
2951	ENDIF
2952
2953	! Software version is the ROPP software version, in the format 'vnn.mmm'.
2954	! The last digit of the fractional part (mmm) is reserved for indicating
2955	! the bending angle processing method in case of bufr encoding.
2956	! At ROPP8.0, a new variable, 'processing_software', has been introduced,
2957	! which can hold information about other software - in this case, the EUM
2958	! processing code that generated the data in the first place.
2959
2960	data%software_version = ' ' ; readstr = ' ' ; readstr2 = ' '
2961
2962	CALL ncdf_getatt('/status/processing/processor_name' , readstr)
2963	CALL ncdf_getatt('/status/processing/processor_version' , readstr2)
2964
2965	i = SCAN(readstr2, ".", .FALSE.)
2966	j = SCAN(readstr2, ".", .TRUE.)
2967
2968	IF (i < 2) THEN
2969	CALL message(msg_warn, "Processor version cannot be decoded")
2970	ELSE
2971	IF (i == j) THEN
2972	j = LEN_TRIM(readstr2) + 1
2973	READ(readstr2(1:i-1), FMT=*) readint
2974	READ(readstr2(i+1:j-1), FMT=*) readint2
2975	ELSE
2976	READ(readstr2(1:i-1), FMT=*) readint
2977	READ(readstr2(i+1:j-1), FMT=*) readint2
2978	ENDIF
2979	IF (getbufr) THEN
2980	SELECT CASE (TRIM(data%bangle_method))
2981	CASE ("FSI")
2982	readint3 = 1
2983	CASE ("GO")
2984	readint3 = 0
2985	CASE DEFAULT
2986	readint3 = 9
2987	END SELECT
2988	WRITE(data%software_version, '("v",i2.2,".",i2.2,i1)') readint, readint2, readint3
2989	ELSE
2990	WRITE(data%software_version, '("v",i2.2,".",i2.2)') readint, readint2
2991	ENDIF
2992	ENDIF
2993
2994	data%processing_software = TRIM(readstr) // ' ' // TRIM(readstr2)
2995
2996	! Processing time, split up string
2997	! FIXME: Is there a library function for this that does some more format checks?
2998	CALL ncdf_getatt('/status/processing/creation_time', readstr)
2999
3000	SELECT CASE (LEN(TRIM(readstr)))
3001	CASE (14)
3002	READ(readstr, '(i4,i2,i2,i2,i2,i2)') &
3003	data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, &
3004	data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second
3005	CASE (17)
3006	READ(readstr, '(i4,i2,i2,i2,i2,i2,i3)') &
3007	data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, &
3008	data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second, data%DTpro%Msec
3009	CASE (19)
3010	READ(readstr, '(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
3011	data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, &
3012	data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second
3013	CASE (23)
3014	READ(readstr, '(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2,1x,i3)') &
3015	data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, &
3016	data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second, data%DTpro%Msec
3017	CASE DEFAULT
3018	CALL message( msg_warn, &
3019	'Invalid /status/processing/creation_time of ' // &
3020	TRIM(readstr) // ' ... defaulting DTpro%year etc.' )
3021	END SELECT
3022
3023
3024	! 9.10 Occultation ID
3025	! -------------------
3026
3027	CALL ropp_io_occid(DATA)
3028
3029
3030	! 9.11 Clean up
3031	! -------------
3032
3033	! IF (ASSOCIATED(idx)) DEALLOCATE(idx)
3034
3035	END SUBROUTINE ropp_io_read_eumdata
3036
3037
3038	!-------------------------------------------------------------------------------
3039	! 10. Accumulate phase
3040	!-------------------------------------------------------------------------------
3041
3042	SUBROUTINE Accumulate_Phase(Ph, Sign) ! (Array of (accumulated) phase, dir)
3043
3044	! Method:
3045	! Sign = 0 or no Sign:
3046	! Adding +-2*Pi where phase jumps from
3047	! +-Pi to -+Pi,
3048	! Sign > 0:
3049	! Adding +2*Pi where phase jumps from
3050	! - to +
3051	! Sign < 0
3052	! Adding -2*Pi where phase jumps from
3053	! + to -
3054
3055	! 10.1 Declarations
3056
3057	USE typesizes, ONLY: wp => EightByteReal
3058
3059	IMPLICIT NONE
3060
3061	REAL(wp), DIMENSION(:), INTENT(inout) :: Ph ! Phase --> accumulated phase
3062	INTEGER, OPTIONAL, INTENT(in) :: Sign ! Phase change sign
3063	! Pi already defined as parameter in coordinates module. Confuses pgf95.
3064	REAL(wp), PARAMETER :: pi1 = 3.141592653589793238_wp
3065	INTEGER :: i ! Array index
3066	INTEGER :: PSign ! Phase change sign
3067
3068	! 10.2 Determine phase change sign
3069
3070	IF (.NOT. PRESENT(Sign)) THEN
3071	PSign = 0
3072	ELSE
3073	PSign = Sign
3074	ENDIF
3075
3076	! 10.3 Accumulate phase
3077
3078	IF (PSign == 0) THEN
3079	DO i=2,SIZE(Ph)
3080	Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1)+pi1, 2.0_wp*pi1) - pi1
3081	ENDDO
3082	ELSEIF (PSign > 0) THEN
3083	DO i=2,SIZE(Ph)
3084	Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1), 2.0_wp*pi1)
3085	ENDDO
3086	ELSEIF (PSign < 0) THEN
3087	DO i=2,SIZE(Ph)
3088	Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1)+2.0_wppi1, 2.0_wppi1) - 2.0_wp*pi1
3089	ENDDO
3090	ENDIF
3091
3092	END SUBROUTINE Accumulate_Phase
3093
3094	!-------------------------------------------------------------------------------
3095	! 11. Handle time format of EUM files
3096	!-------------------------------------------------------------------------------
3097
3098	SUBROUTINE abstimetoDT(i, r, DTocc)
3099
3100	! Based on 'Practical Ephemeris Calculations' by Oliver Montenbruck
3101	! (Springer-Verlag, 1989). Added handling of hour, seconds as well.
3102	!
3103	USE ropp_io_types, ONLY: DT7type
3104
3105	IMPLICIT NONE
3106
3107	INTEGER, INTENT(in) :: i ! number of days since 2000-01-01 00:00:00
3108	REAL(KIND(1.0D0)), INTENT(in) :: r ! seconds since start of day
3109	TYPE(DT7type), INTENT(inout) :: DTocc
3110	REAL(KIND(1.0D0)) :: h, mi, s, ms
3111	INTEGER :: a, b, c, d, e, f, y, m, dd, jd
3112
3113	! calculate JD from input, 2451545 is JD of 2000-01-01
3114	jd = i + 2451545
3115	a = INT(jd+0.5D0)
3116
3117	IF (a < 2299161) THEN
3118	c = a + 1524
3119	ELSE
3120	b = INT( (a - 1867216.25D0) / 36524.25D0 )
3121	c = a + b - INT(b/4D0) + 1525
3122	ENDIF
3123
3124	d = INT( ( c-122.1D0 ) / 365.25D0 )
3125	e = INT( 365.25D0*d)
3126	f = INT( (c-e) / 30.6001D0 )
3127	dd = c - e - INT( 30.6001*f ) + MOD( ( jd+0.5D0 ), DBLE(a) )
3128	m = f - 1 - 12*INT( f/14D0 )
3129	y = d - 4715 - INT( ( 7+m )/10.D0 )
3130
3131
3132	DTocc%Year = y
3133	DTocc%Month = m
3134	DTocc%Day = dd
3135	h = FLOOR(r/3600.D0)
3136	DTocc%Hour = INT(h)
3137	mi = FLOOR( (r - h*3600.D0) / 60.D0 )
3138	s = FLOOR( (r - h3600.D0 - mi60.D0 ) )
3139	ms = (r - h3600.D0 - mi60.D0 - s) * 1000.D0
3140	DTocc%Hour = INT(h)
3141	DTocc%Minute = INT(mi)
3142	DTocc%Second = INT(s)
3143	DTocc%Msec = INT(ms)
3144
3145	END SUBROUTINE abstimetoDT
3146
3147	!-------------------------------------------------------------------------------
3148	! 12. Find where reals are NaN (used by EUMETSAT to indicate missing data)
3149	!-------------------------------------------------------------------------------
3150
3151	FUNCTION ropp_io_isnan(x) RESULT(lnan) ! Says where reals are NaNs.
3152
3153	USE typesizes, ONLY: wp => EightByteReal
3154
3155	IMPLICIT NONE
3156
3157	REAL(wp), DIMENSION(:), INTENT(IN) :: x
3158	LOGICAL, DIMENSION(SIZE(x)) :: lnan
3159	INTEGER :: k
3160
3161	lnan(:) = .FALSE.
3162
3163	! g95 doesn't like this.
3164	! WHERE ( &
3165	! (x /= x) .OR. &
3166	! (x + 1.0_wp == x) .OR. &
3167	! ((x > 0) .EQV. (x <= 0)) &
3168	! ) lnan = .TRUE.
3169
3170	DO k=1,SIZE(x)
3171	IF ( &
3172	(x(k) /= x(k)) .OR. &
3173	(x(k) + 1.0_wp == x(k)) .OR. &
3174	((x(k) > 0) .EQV. (x(k) <= 0)) &
3175	) lnan(k) = .TRUE.
3176	END DO
3177
3178	END FUNCTION ropp_io_isnan
3179
3180
3181	END SUBROUTINE ropp_io_read_ncdf_get_eumdata

Download in other formats:

Original format