Ticket #349: ropp_pp_preprocess_grasrs.f90_SSY

File ropp_pp_preprocess_grasrs.f90_SSY, 28.7 KB (added by Ian Culverwell, 12 years ago)

ropp_pp_preprocess_grasrs.f90_SSY

Line 
1! $Id: ropp_pp_preprocess_GRASRS.f90 2021 2009-01-16 10:49:04Z frhl $
2
3SUBROUTINE ropp_pp_preprocess_GRASRS(ro_data, config, LCF)
4
5!****s* Preprocessing/ropp_pp_preprocess_GRASRS *
6!
7! NAME
8! ropp_pp_preprocess_GRASRS - Mission-specific Level1a data preprocessing
9! for GRAS Raw Sampling data
10!
11! SYNOPSIS
12! call ropp_pp_preprocess_GRASRS(ro_data, config, LCF)
13!
14! DESCRIPTION
15! Merge and upsample CL and RS data
16! 1. Select CL and RS records by LCF flag
17! 2. Generate merged time grid anchored at highest point of RS record
18! 3. Interpolate CL and RS data to merged time grid
19! 4. Merge data
20! 5. Generate phase model and connecting phase
21! 6. Interpolate data on merged time grid within small gaps
22! 7. Restore phase variation
23!
24! INPUTS
25! type(ROprof) :: ro_data ! Radio occultation data strucuture
26! type(PPConfig) :: config ! Configuration options
27! integer :: LCF ! Lost carrier flag
28!
29! OUTPUT
30! type(ROprof) :: ro_data ! Corrected radio occultation data
31! type(PPConfig) :: config ! Configuration options
32! integer :: LCF ! Lost carrier flag
33!
34! NOTES
35! Requires ROprof data structure type, defined in ropp_io module. This
36! routine therefore requires that the ropp_io module is pre-installed before
37! compilation.
38!
39! REFERENCES
40!
41! AUTHOR
42! M Gorbunov, Russian Academy of Sciences, Russia.
43! Any comments on this software should be given via the ROM SAF
44! Helpdesk at http://www.romsaf.org
45!
46! COPYRIGHT
47! Copyright (c) 1998-2010 Michael Gorbunov <michael.gorbunov@zmaw.de>
48! For further details please refer to the file COPYRIGHT
49! which you should have received as part of this distribution.
50!
51!****
52
53!-------------------------------------------------------------------------------
54! 1. Declarations
55!-------------------------------------------------------------------------------
56
57 USE typesizes, ONLY: wp => EightByteReal
58 USE arrays, ONLY: WHERE
59 USE ropp_utils, ONLY: impact_parameter, ropp_MDFV, occ_point
60 USE ropp_io, ONLY: ropp_io_init, ropp_io_free, ropp_io_addvar
61 USE ropp_io_types, ONLY: ROprof, L1atype
62 USE ropp_pp_preproc, not_this => ropp_pp_preprocess_GRASRS
63 USE ropp_pp
64 USE ropp_pp_spline
65 USE ropp_pp_constants
66 USE ropp_pp_types, ONLY: PPConfig
67
68 IMPLICIT NONE
69
70 TYPE(ROprof), INTENT(inout) :: ro_data ! Occultation data struct
71 TYPE(PPconfig), INTENT(inout) :: config ! Configuration options
72 INTEGER, DIMENSION(:), POINTER :: LCF ! Lost carrier flag
73
74 INTEGER, PARAMETER :: np = 500 ! Dimension phase model
75 INTEGER, PARAMETER :: nt = 300 ! RS grid interval
76 INTEGER, PARAMETER :: nv = 5 ! Regression order
77! REAL(wp), PARAMETER :: gap = 0.8 ! Maximum data gap size
78 REAL(wp), PARAMETER :: gap = 0.04 ! Maximum data gap size
79
80 TYPE(L1atype) :: mg_data ! Merged data
81 TYPE(L1atype) :: rs_data ! Raw sampling data
82 TYPE(L1atype) :: cl_data ! Closed loop data
83 INTEGER, DIMENSION(:), ALLOCATABLE :: mg_lcf ! Merged LCF
84 INTEGER, DIMENSION(:), ALLOCATABLE :: rs_lcf ! Raw sampling LCF
85 INTEGER, DIMENSION(:), ALLOCATABLE :: cl_lcf ! Closed loop LCF
86
87 REAL(wp), DIMENSION(:), ALLOCATABLE :: mg_msis ! MSIS phase model (merged)
88 REAL(wp), DIMENSION(:), ALLOCATABLE :: mg_impact ! MSIS phase model (merged)
89 REAL(wp), DIMENSION(:), ALLOCATABLE :: mg_ds ! phase deviation (merged)
90 REAL(wp), DIMENSION(:), ALLOCATABLE :: t_norm ! Normalised time
91 REAL(wp), DIMENSION(:,:),ALLOCATABLE :: KV ! Regression matrix
92 REAL(wp), DIMENSION(0:nv,3) :: coeff_vleo ! Regression coeffs
93 REAL(wp), DIMENSION(0:nv,3) :: coeff_vgns ! Regression coeffs
94
95 INTEGER :: ocd
96 INTEGER :: i, j, n, n1, n_cl, n_rs
97 INTEGER :: icl_min, icl_max
98 INTEGER :: irs_min, irs_max
99 INTEGER :: ig1, ig2, igl
100 INTEGER :: i_int
101 LOGICAL :: mrg_Mode ! Merging mode
102
103 REAL(wp) :: p1, pN, sb
104 REAL(wp) :: ts, t1, tgl, tmin, tmax
105
106 INTEGER, DIMENSION(:), POINTER :: idx
107 INTEGER :: nidx
108 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dtime_idc
109
110!-------------------------------------------------------------------------------
111! 2. Retrieve lost carrier flag information from input data file
112!-------------------------------------------------------------------------------
113
114 IF (ASSOCIATED(ro_data%vlist%VlistD1d)) THEN
115 LCF(:) = NINT(ro_data%vlist%VlistD1d%data(1:size(LCF)))
116 DEALLOCATE(ro_data%vlist%VlistD1d%data)
117 ENDIF
118
119!-------------------------------------------------------------------------------
120! 3. Data quality checks
121!-------------------------------------------------------------------------------
122
123 WHERE(ro_data%Lev1a%phase_L1(:) == ropp_MDFV)
124 ro_data%Lev1a%phase_L1(:) = -1.0_wp
125 ENDWHERE
126 WHERE(ro_data%Lev1a%phase_L2(:) == ropp_MDFV)
127 ro_data%Lev1a%phase_L2(:) = -1.0_wp
128 ENDWHERE
129
130!-------------------------------------------------------------------------------
131! 4. Initialisation
132!-------------------------------------------------------------------------------
133
134 mrg_Mode = .TRUE.
135
136 ! 4.1 Select CL data
137
138 icl_min = SUM(MINLOC(ro_data%lev1a%dtime(:), MASK = .NOT. BTEST(LCF(:),0)))
139 icl_max = SUM(MAXLOC(ro_data%lev1a%dtime(:), MASK = .NOT. BTEST(LCF(:),0)))
140
141 n_cl = icl_max-icl_min+1
142 CALL ropp_io_init(cl_data, n_cl)
143 ALLOCATE(cl_lcf(n_cl))
144
145 cl_data%dtime = ro_data%lev1a%dtime(icl_min:icl_max)
146 DO j=1,3
147 cl_data%r_gns(:,j) = ro_data%lev1a%r_gns(icl_min:icl_max,j)
148 cl_data%v_gns(:,j) = ro_data%lev1a%v_gns(icl_min:icl_max,j)
149 cl_data%r_leo(:,j) = ro_data%lev1a%r_leo(icl_min:icl_max,j)
150 cl_data%v_leo(:,j) = ro_data%lev1a%v_leo(icl_min:icl_max,j)
151 ENDDO
152 cl_data%snr_L1ca = ro_data%lev1a%snr_L1ca(icl_min:icl_max)
153 cl_data%snr_L1p = ro_data%lev1a%snr_L1p(icl_min:icl_max)
154 cl_data%snr_L2p = ro_data%lev1a%snr_L2p(icl_min:icl_max)
155 cl_data%phase_L1 = ro_data%lev1a%phase_L1(icl_min:icl_max)
156 cl_data%phase_L2 = ro_data%lev1a%phase_L2(icl_min:icl_max)
157 cl_lcf = lcf(icl_min:icl_max)
158
159 ! 4.2 Select RS data
160
161 irs_min = SUM(MINLOC(ro_data%lev1a%dtime(:), MASK = BTEST(LCF(:),0)))
162! irs_min = SUM(MINLOC(ro_data%lev1a%dtime(:), MASK = (BTEST(LCF(:),0).and..not.BTEST(LCF(:),3))))
163 irs_max = SUM(MAXLOC(ro_data%lev1a%dtime(:), MASK = BTEST(LCF(:),0)))
164! irs_max = SUM(MAXLOC(ro_data%lev1a%dtime(:), MASK = (BTEST(LCF(:),0).and..not.BTEST(LCF(:),3))))
165
166 IF (irs_min >= 1 .AND. irs_min <= ro_data%lev1a%npoints) THEN
167
168 n_rs = irs_max-irs_min+1
169 CALL ropp_io_init(rs_data, n_rs)
170 ALLOCATE(rs_lcf(n_rs))
171 rs_data%dtime = ro_data%lev1a%dtime(irs_min:irs_max)
172 DO j=1,3
173 rs_data%r_gns(:,j) = ro_data%lev1a%r_gns(irs_min:irs_max,j)
174 rs_data%v_gns(:,j) = ro_data%lev1a%v_gns(irs_min:irs_max,j)
175 rs_data%r_leo(:,j) = ro_data%lev1a%r_leo(irs_min:irs_max,j)
176 rs_data%v_leo(:,j) = ro_data%lev1a%v_leo(irs_min:irs_max,j)
177 ENDDO
178 rs_data%snr_L1ca = ro_data%lev1a%snr_L1ca(irs_min:irs_max)
179 rs_data%snr_L1p = ro_data%lev1a%snr_L1p(irs_min:irs_max)
180 rs_data%snr_L2p = ro_data%lev1a%snr_L2p(irs_min:irs_max)
181 rs_data%phase_L1 = ro_data%lev1a%phase_L1(irs_min:irs_max)
182 rs_data%phase_L2 = ro_data%lev1a%phase_L2(irs_min:irs_max)
183 rs_lcf = lcf(irs_min:irs_max)
184
185
186! idx => WHERE((BTEST(LCF(:),0)) .AND. (.NOT.BTEST(LCF(:),3)), nidx)
187
188
189! IF (nidx > 0) THEN
190
191! n_rs = nidx
192! CALL ropp_io_init(rs_data, n_rs)
193! ALLOCATE(rs_lcf(n_rs))
194! rs_data%dtime = ro_data%lev1a%dtime(idx)
195! DO j=1,3
196! rs_data%r_gns(:,j) = ro_data%lev1a%r_gns(idx,j)
197! rs_data%v_gns(:,j) = ro_data%lev1a%v_gns(idx,j)
198! rs_data%r_leo(:,j) = ro_data%lev1a%r_leo(idx,j)
199! rs_data%v_leo(:,j) = ro_data%lev1a%v_leo(idx,j)
200! END DO
201! rs_data%snr_L1ca = ro_data%lev1a%snr_L1ca(idx)
202! rs_data%snr_L1p = ro_data%lev1a%snr_L1p(idx)
203! rs_data%snr_L2p = ro_data%lev1a%snr_L2p(idx)
204! rs_data%phase_L1 = ro_data%lev1a%phase_L1(idx)
205! rs_data%phase_L2 = ro_data%lev1a%phase_L2(idx)
206! rs_lcf = lcf(idx)
207
208 ELSE
209
210 n_rs = 0
211 mrg_Mode = .FALSE. ! Closed loop data only
212
213 ENDIF
214
215!-------------------------------------------------------------------------------
216! 5. Generate merged time grid
217!-------------------------------------------------------------------------------
218
219 ! 5.1 Determine occultation direction
220
221 p1 = impact_parameter(cl_data%r_leo(1,:), cl_data%r_gns(1,:))
222 pN = impact_parameter(cl_data%r_leo(n_cl,:), cl_data%r_gns(n_cl,:))
223 ocd = NINT(SIGN(1.0_wp, pN - p1))
224
225 ! 5.2 Determine reference points and grid time step
226
227 IF (mrg_Mode) THEN ! Merge RS+CL
228
229 ts = MINVAL(ABS(rs_data%dtime(1+nt:n_rs)-rs_data%dtime(1:n_rs-nt)))/nt
230 IF (ocd < 0) THEN
231 tmin = MINVAL(cl_data%dtime(:))
232 tmax = MAXVAL(rs_data%dtime(:))
233 t1 = MINVAL(rs_data%dtime(:))
234 ELSE
235 tmin = MINVAL(rs_data%dtime(:))
236 tmax = MAXVAL(cl_data%dtime(:))
237 t1 = MAXVAL(rs_data%dtime(:))
238 ENDIF
239
240 ELSE
241
242 ts = MINVAL(ABS(cl_data%dtime(2:n_cl)-cl_data%dtime(1:n_cl-1)))
243 tmin = MINVAL(cl_data%dtime(:))
244 tmax = MAXVAL(cl_data%dtime(:))
245 IF (ocd < 0) THEN
246 t1 = MINVAL(cl_data%dtime(:))
247 ELSE
248 t1 = MAXVAL(cl_data%dtime(:))
249 ENDIF
250
251 ENDIF
252
253 ! 5.3 Determine anchor point index and grid dimension
254
255 n1 = FLOOR((t1 - tmin)/ts) + 1
256 n = FLOOR((tmax - t1)/ts) + n1
257
258 ! 5.4 Generate merged time grid
259 CALL ropp_io_init(mg_data, n)
260 ALLOCATE(mg_LCF(n))
261 mg_data%reference_frame%r_leo = ro_data%lev1a%reference_frame%r_leo
262 mg_data%reference_frame%r_gns = ro_data%lev1a%reference_frame%r_gns
263
264 IF (ocd < 0) THEN
265 DO i=1,n
266 mg_data%dtime(i) = ((n-i)*t1 + (i-n1)*tmax)/(n-n1)
267 ENDDO
268 ELSE
269 DO i=1,n
270 mg_data%dtime(i) = ((n1-i)*tmin + (i-1)*t1)/(n1-1)
271 ENDDO
272 ENDIF
273print*,'n, minval(mg_data%dtime), maxval(mg_data%dtime), ro_data%Lev1a%npoints = ', n, minval(mg_data%dtime), maxval(mg_data%dtime), ro_data%Lev1a%npoints
274
275!-------------------------------------------------------------------------------
276! 6. Interpolate CL and RS data to merged time grid
277!-------------------------------------------------------------------------------
278
279 ! 6.1 Interpolate trajectories
280
281 ALLOCATE(t_norm(ro_data%Lev1a%npoints))
282! ALLOCATE(t_norm(n_cl+n_rs))
283 ALLOCATE(KV(ro_data%Lev1a%npoints, 0:nv))
284! ALLOCATE(KV(n_cl+n_rs, 0:nv))
285 t_norm(:) = (ro_data%Lev1a%dtime(:) - ro_data%Lev1a%dtime(1))/ &
286 (ro_data%Lev1a%dtime(ro_data%Lev1a%npoints) - ro_data%Lev1a%dtime(1))
287! t_norm(:) = (ro_data%Lev1a%dtime(:n_cl+n_rs) - ro_data%Lev1a%dtime(1))/ &
288! (ro_data%Lev1a%dtime(n_cl+n_rs) - ro_data%Lev1a%dtime(1))
289 CALL ropp_pp_init_polynomial(t_norm, KV)
290 DO j=1,3
291 CALL ropp_pp_regression(KV,ro_data%Lev1a%r_leo(:,j),coeff_vleo(:,j))
292! CALL ropp_pp_regression(KV,ro_data%Lev1a%r_leo(:n_cl+n_rs,j),coeff_vleo(:,j))
293 CALL ropp_pp_regression(KV,ro_data%Lev1a%r_gns(:,j),coeff_vgns(:,j))
294! CALL ropp_pp_regression(KV,ro_data%Lev1a%r_gns(:n_cl+n_rs,j),coeff_vgns(:,j))
295 ENDDO
296
297! For diagnostics, plot out residuals evaluated at original times.
298
299 DO i=1,ro_data%Lev1a%npoints
300 CALL ropp_pp_interpolate_trajectory(ro_data%Lev1a%dtime, &
301 coeff_vleo, coeff_vgns, &
302 ro_data%georef%r_coc, &
303 ro_data%Lev1a%dtime(i), &
304 mg_data%r_leo(i,:), &
305 mg_data%v_leo(i,:), &
306 mg_data%r_gns(i,:), &
307 mg_data%v_gns(i,:))
308 ENDDO
309
310print*, 'SIZE(LCF) = ', SIZE(LCF)
311print*, 'LCF(1:5) = ', LCF(1:5)
312
313 ALLOCATE(dtime_idc(ro_data%Lev1a%npoints, 3))
314 dtime_idc(:, 1) = mg_data%dtime(1:ro_data%Lev1a%npoints)
315 dtime_idc(:, 2) = ro_data%Lev1a%dtime(1:ro_data%Lev1a%npoints)
316 dtime_idc(:, 3) = REAL(LCF, wp)
317 CALL ropp_io_addvar( &
318 ro_data, &
319 name = "dtime_idc", &
320 long_name = "mg_data dtime, ro_data dtime and lcf", &
321 units = "s", &
322 range = (/-1.0E6_wp, 1.0E6_wp/), &
323 DATA = dtime_idc)
324 DEALLOCATE(dtime_idc)
325
326print*,'SHAPE(ro_data%Lev1a%r_gns) = ', SHAPE(ro_data%Lev1a%r_gns)
327print*,'SHAPE(mg_data%r_gns) = ', SHAPE(mg_data%r_gns)
328print*,'SHAPE(mg_data%r_gns-ro_data%Lev1a%r_gns) = ', SHAPE(mg_data%r_gns-ro_data%Lev1a%r_gns)
329
330print*,'SHAPE(ropp_io_isnan(ro_data%Lev1a%r_gns)) = ', SHAPE(ropp_io_isnan(ro_data%Lev1a%r_gns))
331print*,'COUNT(ropp_io_isnan(ro_data%Lev1a%r_gns)) = ', COUNT(ropp_io_isnan(ro_data%Lev1a%r_gns))
332
333print*,'SHAPE(ropp_io_isnan(mg_data%r_gns)) = ', SHAPE(ropp_io_isnan(mg_data%r_gns))
334print*,'COUNT(ropp_io_isnan(mg_data%r_gns)) = ', COUNT(ropp_io_isnan(mg_data%r_gns))
335
336print*,'SHAPE(ropp_io_isnan(ro_data%Lev1a%r_gns) .OR. ropp_io_isnan(mg_data%r_gns)) = ', SHAPE(ropp_io_isnan(ro_data%Lev1a%r_gns) .OR. ropp_io_isnan(mg_data%r_gns))
337print*,'COUNT(ropp_io_isnan(ro_data%Lev1a%r_gns) .OR. ropp_io_isnan(mg_data%r_gns)) = ', COUNT(ropp_io_isnan(ro_data%Lev1a%r_gns) .OR. ropp_io_isnan(mg_data%r_gns))
338
339! WHERE (ropp_io_isnan(ro_data%Lev1a%r_gns(1:ro_data%Lev1a%npoints)) .OR. ropp_io_isnan(mg_data%r_gns(1:ro_data%Lev1a%npoints)))
340! mg_data%r_gns = ropp_MDFV
341! ELSEWHERE
342! mg_data%r_gns = ro_data%Lev1a%r_gns - mg_data%r_gns
343! ENDWHERE
344
345 CALL ropp_io_addvar( &
346 ro_data, &
347 name = "gns_residual_1", &
348 long_name = "1st residual (data-fit) in GNS location", &
349 units = "m", &
350 range = (/-1.0_wp, 1.0_wp/), &
351 DATA = ro_data%Lev1a%r_gns(1:ro_data%Lev1a%npoints, :) - mg_data%r_gns(1:ro_data%Lev1a%npoints, :))
352
353! WHERE (ropp_io_isnan(ro_data%Lev1a%r_leo) .OR. ropp_io_isnan(mg_data%r_leo))
354! mg_data%r_leo = ropp_MDFV
355! ELSEWHERE
356! mg_data%r_leo = ro_data%Lev1a%r_leo - mg_data%r_leo
357! ENDWHERE
358
359 CALL ropp_io_addvar( &
360 ro_data, &
361 name = "leo_residual_1", &
362 long_name = "1st residual (data-fit) in LEO location", &
363 units = "m", &
364 range = (/-1.0_wp, 1.0_wp/), &
365 DATA = ro_data%Lev1a%r_leo(1:ro_data%Lev1a%npoints, :) - mg_data%r_leo(1:ro_data%Lev1a%npoints, :))
366
367
368! Perform regression on residual positions to gain higher accuracy
369 DO j=1,3
370 CALL ropp_pp_residual_regression(KV,t_norm, &
371 ro_data%Lev1a%r_leo(:,j),coeff_vleo(:,j))
372! ro_data%Lev1a%r_leo(:n_cl+n_rs,j),coeff_vleo(:,j))
373 CALL ropp_pp_residual_regression(KV,t_norm, &
374 ro_data%Lev1a%r_gns(:,j),coeff_vgns(:,j))
375! ro_data%Lev1a%r_gns(:n_cl+n_rs,j),coeff_vgns(:,j))
376 ENDDO
377 DEALLOCATE(t_norm)
378 DEALLOCATE(KV)
379
380 DO i=1,ro_data%Lev1a%npoints
381 CALL ropp_pp_interpolate_trajectory(ro_data%Lev1a%dtime, &
382 coeff_vleo, coeff_vgns, &
383 ro_data%georef%r_coc, &
384 ro_data%Lev1a%dtime(i), &
385 mg_data%r_leo(i,:), &
386 mg_data%v_leo(i,:), &
387 mg_data%r_gns(i,:), &
388 mg_data%v_gns(i,:))
389 ENDDO
390
391! WHERE (ropp_io_isnan(ro_data%Lev1a%r_gns(1:ro_data%Lev1a%npoints, :)) .OR. ropp_io_isnan(mg_data%r_gns(1:ro_data%Lev1a%npoints, :)))
392! mg_data%r_gns = ropp_MDFV
393! ELSEWHERE
394! mg_data%r_gns = ro_data%Lev1a%r_gns - mg_data%r_gns
395! ENDWHERE
396
397 CALL ropp_io_addvar( &
398 ro_data, &
399 name = "gns_residual_2", &
400 long_name = "2nd residual (data-fit) in GNS location", &
401 units = "m", &
402 range = (/-1.0_wp, 1.0_wp/), &
403 DATA = ro_data%Lev1a%r_gns(1:ro_data%Lev1a%npoints, :) - mg_data%r_gns(1:ro_data%Lev1a%npoints, :))
404
405! WHERE (ropp_io_isnan(ro_data%Lev1a%r_leo) .OR. ropp_io_isnan(mg_data%r_leo))
406! mg_data%r_leo = ropp_MDFV
407! ELSEWHERE
408! mg_data%r_leo = ro_data%Lev1a%r_leo - mg_data%r_leo
409! ENDWHERE
410
411 CALL ropp_io_addvar( &
412 ro_data, &
413 name = "leo_residual_2", &
414 long_name = "2nd residual (data-fit) in LEO location", &
415 units = "m", &
416 range = (/-1.0_wp, 1.0_wp/), &
417 DATA = ro_data%Lev1a%r_leo(1:ro_data%Lev1a%npoints, :) - mg_data%r_leo(1:ro_data%Lev1a%npoints, :))
418
419print*,'ro_data%vlist%vlistd2d%name = ', ro_data%vlist%vlistd2d%name
420print*,'ro_data%vlist%vlistd2d%next%name = ', ro_data%vlist%vlistd2d%next%name
421print*,'ro_data%vlist%vlistd2d%next%next%name = ', ro_data%vlist%vlistd2d%next%next%name
422print*,'ro_data%vlist%vlistd2d%next%next%next%name = ', ro_data%vlist%vlistd2d%next%next%next%name
423print*,'ro_data%vlist%vlistd2d%next%next%next%next%name = ', ro_data%vlist%vlistd2d%next%next%next%next%name
424
425 DO i=1,n
426 CALL ropp_pp_interpolate_trajectory(ro_data%Lev1a%dtime, &
427! CALL ropp_pp_interpolate_trajectory(ro_data%Lev1a%dtime(:n_cl+n_rs), &
428 coeff_vleo, coeff_vgns, &
429 ro_data%georef%r_coc, &
430 mg_data%dtime(i), &
431 mg_data%r_leo(i,:), &
432 mg_data%v_leo(i,:), &
433 mg_data%r_gns(i,:), &
434 mg_data%v_gns(i,:))
435 ENDDO
436
437 ! 6.2 Interpolate LCF, phase and amplitude
438
439 DO i=1,n
440
441 IF (mrg_Mode) THEN ! Merge LC + RS data
442
443 i_int = ropp_pp_seek_index(rs_data%dtime, mg_data%dtime(i))
444
445 IF (i_int == 0 .OR. i_int == rs_data%npoints) THEN
446
447 ! RS data range outside current time period
448 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
449 cl_LCF, mg_LCF(i))
450 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
451 cl_data%snr_L1ca, mg_data%snr_L1ca(i))
452 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
453 cl_data%snr_L1p, mg_data%snr_L1p(i))
454 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
455 cl_data%snr_L2p, mg_data%snr_L2p(i))
456 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
457 cl_data%phase_L1, mg_data%phase_L1(i))
458 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
459 cl_data%phase_L2, mg_data%phase_L2(i))
460
461 ELSE
462
463 ! RS data range within current time period
464 CALL ropp_pp_interpol(rs_data%dtime, mg_data%dtime(i), &
465 rs_LCF, mg_LCF(i))
466 CALL ropp_pp_interpol(rs_data%dtime, mg_data%dtime(i), &
467 rs_data%snr_L1ca, mg_data%snr_L1ca(i))
468 CALL ropp_pp_interpol(rs_data%dtime, mg_data%dtime(i), &
469 rs_data%snr_L1p, mg_data%snr_L1p(i))
470 CALL ropp_pp_interpol(rs_data%dtime, mg_data%dtime(i), &
471 rs_data%snr_L2p, mg_data%snr_L2p(i))
472 CALL ropp_pp_interpol(rs_data%dtime, mg_data%dtime(i), &
473 rs_data%phase_L1, mg_data%phase_L1(i))
474 CALL ropp_pp_interpol(rs_data%dtime, mg_data%dtime(i), &
475 rs_data%phase_L2, mg_data%phase_L2(i))
476
477 ENDIF
478
479 ELSE ! Use only CL data
480
481 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
482 cl_LCF, mg_LCF(i))
483 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
484 cl_data%snr_L1ca, mg_data%snr_L1ca(i))
485 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
486 cl_data%snr_L1p, mg_data%snr_L1p(i))
487 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
488 cl_data%snr_L2p, mg_data%snr_L2p(i))
489 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
490 cl_data%phase_L1, mg_data%phase_L1(i))
491 CALL ropp_pp_interpol(cl_data%dtime, mg_data%dtime(i), &
492 cl_data%phase_L2, mg_data%phase_L2(i))
493
494 ENDIF
495
496
497 ENDDO
498
499!-------------------------------------------------------------------------------
500! 7. Connect phase
501!-------------------------------------------------------------------------------
502
503 ! 7.1 Calculate MSIS excess phase on merged time grid
504 CALL occ_point( mg_data%r_leo, mg_data%r_gns, &
505 ro_data%georef%lat, ro_data%georef%lon, &
506 ro_data%georef%r_coc, ro_data%georef%roc, &
507 ro_data%georef%azimuth, &
508 ro_data%georef%undulation, &
509 config%egm96, config%corr_egm96)
510
511 ALLOCATE(mg_msis(n))
512 ALLOCATE(mg_impact(n))
513 CALL ropp_pp_modelphase(ro_data%dtocc%month, ro_data%georef%lat, &
514 ro_data%georef%lon, mg_data%dtime, &
515 mg_data%r_leo, mg_data%r_gns, &
516 ro_data%georef%r_coc, ro_data%georef%roc, &
517 mg_msis, mg_impact, config)
518 DEALLOCATE(mg_impact)
519
520 ! 7.2 Phase re-accumulation
521
522 ALLOCATE(mg_ds(n))
523 WHERE (.NOT. BTEST(mg_LCF(:),3))
524 mg_ds(:) = (2.0_wp*Pi*f_L1/c_light)*(mg_data%phase_L1(:)-mg_msis(:))
525 ELSEWHERE
526 mg_ds(:) = 0.0_wp
527 ENDWHERE
528
529 CALL Accumulate_Phase(mg_DS)
530
531!-------------------------------------------------------------------------------
532! 8. Fill in data gaps
533!-------------------------------------------------------------------------------
534
535 ts = MINVAL(ABS(cl_data%dtime(2:n_cl) - cl_data%dtime(1:n_cl-1)))
536
537 ig1 = 1
538
539 DO i=2,n
540
541 IF (BTEST(mg_LCF(i),3) .AND. .NOT. BTEST(mg_LCF(i-1),3)) THEN
542 ig1 = i-1
543 ENDIF
544
545 IF (.NOT. BTEST(mg_LCF(i),3) .AND. BTEST(mg_LCF(i-1),3)) THEN
546
547 ig2 = i
548 igl = ig2-ig1
549 tgl = mg_data%dtime(ig2) - mg_data%dtime(ig1)
550
551 IF (tgl < gap) THEN ! Fill in data gap
552
553 DO j=ig1+1,ig2-1
554 CALL ropp_pp_interpol(mg_data%dtime(ig1:ig2:igl), mg_data%dtime(j), &
555 mg_ds(ig1:ig2:igl), mg_ds(j))
556 mg_LCF(j) = IBCLR(mg_LCF(j),3)
557 ENDDO
558 ENDIF
559
560 ENDIF
561
562 ENDDO
563
564!-------------------------------------------------------------------------------
565! 9. Restore and normalise phase
566!-------------------------------------------------------------------------------
567
568 ! 9.1 Restore phase variation
569
570 mg_data%phase_L1(:) = mg_msis(:) + (c_light/(2.0_wp*pi*f_L1))*mg_ds(:)
571
572 ! 9.2 Normalise to zero minimum value
573
574 sb = MINVAL(mg_data%phase_L1(:))
575 mg_data%phase_L1(:) = mg_data%phase_L1(:) - sb
576
577!-------------------------------------------------------------------------------
578! 7. Setup merged output data
579!-------------------------------------------------------------------------------
580
581 CALL ropp_io_free(ro_data%Lev1a)
582 DEALLOCATE(LCF)
583
584 CALL ropp_io_init(ro_data%Lev1a, n)
585 ALLOCATE(LCF(n))
586
587 ro_data%Lev1a = mg_data
588 LCF(:) = mg_LCF(:)
589
590!-------------------------------------------------------------------------------
591! 8. Re-compute occultation point
592!-------------------------------------------------------------------------------
593
594 CALL occ_point( ro_data%Lev1a%r_leo, ro_data%Lev1a%r_gns, &
595 ro_data%georef%lat, ro_data%georef%lon, &
596 ro_data%georef%r_coc, ro_data%georef%roc, &
597 ro_data%georef%azimuth, &
598 ro_data%georef%undulation, &
599 config%egm96, config%corr_egm96)
600
601!-------------------------------------------------------------------------------
602! 9. Quality control
603!-------------------------------------------------------------------------------
604
605 WHERE(ro_data%Lev1a%phase_L1(:) == ropp_MDFV)
606 ro_data%Lev1a%phase_L1(:) = -1.0_wp
607 ENDWHERE
608 WHERE(ro_data%Lev1a%phase_L2(:) == ropp_MDFV)
609 ro_data%Lev1a%phase_L2(:) = -1.0_wp
610 ENDWHERE
611
612!-------------------------------------------------------------------------------
613! 10. Clean up
614!-------------------------------------------------------------------------------
615
616 DEALLOCATE(mg_msis)
617 CALL ropp_io_free(cl_data)
618 CALL ropp_io_free(rs_data)
619
620print*,'ro_data%vlist%vlistd2d%name = ', ro_data%vlist%vlistd2d%name
621print*,'ro_data%vlist%vlistd2d%data(1:5, 1:2) = ', ro_data%vlist%vlistd2d%data(1:5, 1:2)
622
623print*,'ro_data%vlist%vlistd2d%next%name = ', ro_data%vlist%vlistd2d%next%name
624print*,'ro_data%vlist%vlistd2d%next%data(1:5, 1:3) = ', ro_data%vlist%vlistd2d%next%data(1:5, 1:3)
625
626print*,'ro_data%vlist%vlistd2d%next%next%name = ', ro_data%vlist%vlistd2d%next%next%name
627print*,'ro_data%vlist%vlistd2d%next%next%data(1:5, 1:3) = ', ro_data%vlist%vlistd2d%next%next%data(1:5, 1:3)
628
629print*,'ro_data%vlist%vlistd2d%next%next%next%name = ', ro_data%vlist%vlistd2d%next%next%next%name
630print*,'ro_data%vlist%vlistd2d%next%next%next%data(1:5, 1:3) = ', ro_data%vlist%vlistd2d%next%next%next%data(1:5, 1:3)
631
632print*,'ro_data%vlist%vlistd2d%next%next%next%next%name = ', ro_data%vlist%vlistd2d%next%next%next%next%name
633print*,'ro_data%vlist%vlistd2d%next%next%next%next%data(1:5, 1:3) = ', ro_data%vlist%vlistd2d%next%next%next%next%data(1:5, 1:3)
634
635IF(ASSOCIATED(ro_data%vlist%VlistD0d)) print*,'SHAPE(ro_data%vlist%VlistD0d) = ', SHAPE(ro_data%vlist%VlistD0d)
636IF(ASSOCIATED(ro_data%vlist%VlistD1d)) print*,'SHAPE(ro_data%vlist%VlistD1d) = ', SHAPE(ro_data%vlist%VlistD1d)
637IF(ASSOCIATED(ro_data%vlist%VlistD2d)) print*,'SHAPE(ro_data%vlist%VlistD2d) = ', SHAPE(ro_data%vlist%VlistD2d)
638
639CONTAINS
640
641!-------------------------------------------------------------------------------
642! 11. Transform phase to accumulated phase
643!-------------------------------------------------------------------------------
644
645 SUBROUTINE Accumulate_Phase(Ph, Sign) ! (Array of (accumulated) phase, dir)
646
647! Method:
648! Sign = 0 or no Sign:
649! Adding +-2*Pi where phase jumps from
650! +-Pi to -+Pi,
651! Sign > 0:
652! Adding +2*Pi where phase jumps from
653! - to +
654! Sign < 0
655! Adding -2*Pi where phase jumps from
656! + to -
657
658 ! 11.1 Declarations
659
660 USE typesizes, ONLY: wp => EightByteReal
661 USE ropp_pp_constants, ONLY: pi
662 IMPLICIT NONE
663
664 REAL(wp), DIMENSION(:), INTENT(inout) :: Ph ! Phase --> accumulated phase
665 INTEGER, OPTIONAL, INTENT(in) :: Sign ! Phase change sign
666
667 INTEGER :: i ! Array index
668 INTEGER :: PSign ! Phase change sign
669
670 ! 11.2 Determine phase change sign
671
672 IF (.NOT. PRESENT(Sign)) THEN
673 PSign = 0
674 ELSE
675 PSign = Sign
676 ENDIF
677
678 ! 11.3 Accumulate phase
679
680 IF (PSign == 0) THEN
681 DO i=2,SIZE(Ph)
682 Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1)+pi, 2*pi) - pi
683 ENDDO
684 ELSEIF (PSign > 0) THEN
685 DO i=2,SIZE(Ph)
686 Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1), 2*pi)
687 ENDDO
688 ELSEIF (PSign < 0) THEN
689 DO i=2,SIZE(Ph)
690 Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1)+2*pi, 2*pi) - 2*pi
691 ENDDO
692 ENDIF
693
694 END SUBROUTINE Accumulate_Phase
695
696
697
698
699FUNCTION ropp_io_isnan(x) RESULT(lnan) ! Says where reals are NaNs.
700
701 USE typesizes, ONLY: wp => EightByteReal
702
703 IMPLICIT NONE
704
705 INTEGER :: IPInf, IMinf
706 REAL(wp) :: PInf, MInf
707 data IPInf/B'01111111100000000000000000000000'/ ! +Infinity
708 data IMInf/B'11111111100000000000000000000000'/ ! -Infinity
709
710 REAL(wp), DIMENSION(:,:), INTENT(IN) :: x
711 LOGICAL, DIMENSION(SIZE(x,1), SIZE(x,2)) :: lnan
712
713 PInf = TRANSFER(IPinf, Pinf)
714 Minf = TRANSFER(IMinf, Minf)
715
716 lnan = .FALSE.
717
718! WHERE ((x /= x) .OR. (x + 1.0_wp == x) .OR. ((x > 0) .EQV. (x <= 0)) .OR. (x == Pinf) .OR. (x == Minf)) &
719 WHERE ((x /= x) .OR. (x + 1.0_wp == x) .OR. ((x > 0) .EQV. (x <= 0))) &
720 lnan = .TRUE.
721
722END FUNCTION ropp_io_isnan
723
724END SUBROUTINE ropp_pp_preprocess_GRASRS