Ticket #262: ropp_pp_spectra_tool.f90

File ropp_pp_spectra_tool.f90, 16.9 KB (added by kmk, 14 years ago)

Suggested modified version of ropp_pp_spectra_tool

Line 
1! $Id: ropp_pp_spectra_tool.f90 2021 2009-01-16 10:49:04Z frhl $
2
3PROGRAM ropp_pp_spectra_tool
4
5!****p* Programs/ropp_pp_spectra_tool *
6!
7! NAME
8! ropp_pp_spectra_tool
9!
10! SYNOPSIS
11! Pre-processing tool to compute local spatial spectra of the
12! complex wave field.
13!
14! > ropp_pp_spectra_tool <infile(s)> [-c <cnfgfile>]
15! [-o <outfile>] [-h] [-v]
16!
17! ARGUMENTS
18! <infile(s)> One (or more) input file names
19!
20! OPTIONS
21! -c <cnfgfile> name of configuration file
22! -o <outfile> name of spectra output file
23! (default: ropp_pp_spectra_<type>.dat)
24! -h help
25! -v version information
26!
27! DESCRIPTION
28! This program reads RO L1 and L2 excess phase data as a function of time
29! from the input data files and calculates local spatial spectra as a
30! function of doppler frequency shift and time or as a function of bending
31! angle and impact parameter. The results are written to ASCII files for
32! plotting (see tool ropp_pp_plot_spectra.pro)
33!
34! NOTES
35! If the input file is a multifile, or more than one input files are
36! specified, output for all profiles are written out in separate files.
37!
38! Already existing output files will be overwritten.
39!
40! EXAMPLE
41! To calculate spectra for one of the example files in the data directory:
42!
43! > ropp_pp_spectra_tool ../data/input.nc
44!
45! SEE ALSO
46! For an example plotting tool to view the resulting spectra see
47! ropp_pp_plot_spectra.pro
48!
49! REFERENCES
50! Gorbunov M.E., Lauritsen K.B., Rhodin A., Tomassini M. and
51! Kornblueh L. 2006
52! Radio holographic filtering, error estimation, and quality control of
53! radio occultation data
54! Journal of Geophysical Research (111) D10105
55!
56! AUTHOR
57! Met Office, Exeter, UK.
58! Any comments on this software should be given via the GRAS SAF
59! Helpdesk at http://www.grassaf.org
60!
61! COPYRIGHT
62! (c) EUMETSAT. All rights reserved.
63! For further details please refer to the file COPYRIGHT
64! which you should have received as part of this distribution.
65!
66!****
67
68!-------------------------------------------------------------------------------
69! 1. Declarations
70!-------------------------------------------------------------------------------
71
72 USE typesizes, ONLY: wp => EightByteReal
73 USE ropp_utils
74 USE ropp_io
75 USE ropp_io_types, ONLY: ROprof, L1atype, PCD_rising, PCD_open_loop
76 USE ropp_pp
77 USE ropp_pp_preproc
78 USE ropp_pp_types, ONLY: PPConfig
79
80
81
82 IMPLICIT NONE
83
84 TYPE(ROprof) :: ro_data ! Input RO data
85 TYPE(ppConfig) :: config ! Configuration options
86 TYPE(L1atype) :: Lev1a ! Temporary Level1a structure for storage
87
88 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: A ! L1 and L2 amplitude
89 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: S ! L1 and L2 excess phase
90 REAL(wp), DIMENSION(:), POINTER :: phase => null() ! Work array
91 REAL(wp), DIMENSION(:), POINTER :: phase_LM => null() ! Model excess phase (m)
92 REAL(wp), DIMENSION(:), POINTER :: impact_LM => null() ! Model impact (m)
93
94 INTEGER, DIMENSION(:), POINTER :: LCF => null() ! Lost carrier flag
95! INTEGER, DIMENSION(:), ALLOCATABLE :: LCF ! Lost carrier flag
96 REAL(wp) :: secs_past_hour ! No. of seconds past hour
97 REAL(wp) :: sb ! Phase base
98 INTEGER :: iol ! Open loop discont index
99 INTEGER :: ic ! Index counter
100 INTEGER :: n ! Number of data points
101
102 INTEGER :: idummy
103 INTEGER :: i, iargc, argc, k
104 INTEGER :: n_files, n_profiles
105 LOGICAL :: give_help
106
107 CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
108 CHARACTER(len = 4096) :: ofile
109 CHARACTER(len = 4096) :: cfg_file = " "
110 CHARACTER(len = 256) :: buffer
111 CHARACTER(len = 80) :: mfile = " "
112 CHARACTER(len = 80) :: navfile = " "
113 CHARACTER(len = 10) :: method = " "
114 CHARACTER(len = 4) :: istr
115 CHARACTER(len = 6) :: nstr
116
117!-------------------------------------------------------------------------------
118! 2. Default settings
119!-------------------------------------------------------------------------------
120
121 give_help = .FALSE.
122 ofile = "ropp_pp_spectra"
123
124 CALL message(msg_noin, '')
125 CALL message(msg_noin, &
126 '----------------------------------------------------------------------')
127 CALL message(msg_noin, &
128 ' ROPP Pre-processor Spectra Tool' )
129 CALL message(msg_noin, &
130 '----------------------------------------------------------------------')
131 CALL message(msg_noin, '')
132
133!-------------------------------------------------------------------------------
134! 3. Command line arguments
135!-------------------------------------------------------------------------------
136
137 argc = iargc()
138 i = 1
139 n_files = 0
140 ALLOCATE(ifiles(argc))
141
142 DO WHILE(i <= argc)
143 CALL getarg(i, buffer)
144 SELECT CASE (buffer)
145 CASE('-o') ! Output file name (netCDF output)
146 CALL getarg(i+1, buffer)
147 ofile = buffer
148 i = i + 1
149 CASE('-c') ! Configuration file name
150 CALL getarg(i+1, buffer)
151 cfg_file = buffer
152 i = i + 1
153 CASE('-m') ! Ionospheric correction method
154 CALL getarg(i+1, buffer)
155 method = buffer
156 i = i + 1
157 CASE('-mfile') ! Model refractivity file
158 CALL getarg(i+1, buffer)
159 mfile = buffer
160 i = i + 1
161 CASE('-navfile') ! External navigation bit file
162 CALL getarg(i+1, buffer)
163 navfile = buffer
164 i = i + 1
165 CASE('-h', '--help', '?') ! help
166 give_help = .TRUE.
167 CASE('-v', '-V', '--version') ! version information
168 CALL version_info()
169 CALL EXIT(0)
170 CASE default ! Input file name
171 IF ( buffer(1:1) /= "-" ) THEN
172 n_files = n_files + 1
173 ifiles(n_files) = buffer
174 END IF
175 END SELECT
176 i = i + 1
177 END DO
178
179 IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN
180 CALL usage()
181 CALL EXIT(0)
182 ENDIF
183
184 ! 3.1 Read configuration file (if exists)
185
186 IF (method /= " ") config%method = method
187 IF (mfile /= " ") config%mfile = mfile
188 IF (navfile /= " ") config%navbit_file = navfile
189
190 IF (cfg_file /= " ") THEN
191 CALL message(msg_info, &
192 "Reading configuration file " // TRIM(cfg_file) // ".\n")
193 CALL ropp_pp_read_config(cfg_file, config)
194 IF (method /= " ") config%method = method
195 IF (mfile /= " ") config%mfile = mfile
196 IF (navfile /= " ") config%navbit_file = navfile
197 ENDIF
198
199!-------------------------------------------------------------------------------
200! 4. Remove pre-existing output file
201!-------------------------------------------------------------------------------
202
203 CALL file_delete(ofile, idummy)
204
205!-------------------------------------------------------------------------------
206! 5. Loop over all input files and profiles
207!-------------------------------------------------------------------------------
208
209 DO k = 1, n_files
210
211 n_profiles = ropp_io_nrec(ifiles(k))
212
213 DO i = 1, n_profiles
214
215 WRITE(istr, '(i4)') i
216 WRITE(nstr, '(i6)') n_profiles
217 CALL message(msg_info, "Processing profile " // istr // " of " // nstr )
218
219!-------------------------------------------------------------------------------
220! 6. Read data
221!-------------------------------------------------------------------------------
222
223! 6.1 Read input file
224
225 CALL ropp_io_read(ro_data, ifiles(k), rec=i, ranchk=.TRUE.)
226 CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")
227
228 IF (ro_data%Lev1a%Npoints == 0) THEN
229 config%obs_ok = .FALSE.
230 CALL message(msg_fatal, "FAILURE: No Level1a data in file " // &
231 ifiles(k) // "\n")
232 ENDIF
233
234 ! 6.2 Shrink ro_data to correct size (for multiple profiles)
235
236 CALL ropp_io_init(Lev1a, ro_data%lev1a%npoints)
237 Lev1a = ro_data%lev1a
238 ro_data%lev1a = Lev1a
239
240!-------------------------------------------------------------------------------
241! 7. Check coordinate frames
242!-------------------------------------------------------------------------------
243
244 CALL ropp_pp_set_coordinates(ro_data)
245
246!-------------------------------------------------------------------------------
247! 8. Determine occultation point georeference information
248!-------------------------------------------------------------------------------
249
250 CALL occ_point( ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
251 ro_data%georef%lat, ro_data%georef%lon, &
252 ro_data%georef%r_coc, ro_data%georef%roc, &
253 ro_data%georef%azimuth, ro_data%georef%undulation , &
254 config%egm96, config%corr_egm96)
255
256!-------------------------------------------------------------------------------
257! 9. Mission-specific pre-processing and input data cut-off, compute spectra
258!-------------------------------------------------------------------------------
259
260 ! 9.1 Set default configuration options
261
262 IF (ANY(ro_data%Lev1a%r_leo == ropp_MDFV) .OR. &
263 ANY(ro_data%Lev1a%r_gns == ropp_MDFV)) THEN
264 CALL message(msg_warn, 'Invalid coordinate values. Exiting processing.')
265 config%obs_ok = .false.
266 EXIT
267 ENDIF
268
269 config%opt_spectra = .true.
270 config%Acut = 0.0_wp
271 config%Pcut = -9999.9_wp
272 config%Bcut = 0.0_wp
273 config%obs_ok = .TRUE.
274
275 n = ro_data%Lev1a%Npoints
276
277 ALLOCATE(LCF(n))
278 LCF(:) = 0
279
280 CALL ropp_pp_cutoff_amplitude(ro_data, LCF, config)
281
282 SELECT CASE (ro_data%leo_id(1:2))
283
284 ! 9.1 COSMIC
285
286 CASE ('C0','CO')
287 CALL message(msg_info, 'COSMIC data preprocessing')
288 CALL ropp_pp_preprocess_COSMIC(ro_data, config, LCF)
289
290 ! 9.2 CHAMP and GRACE
291
292 CASE ('CH','GR')
293 CALL message(msg_info, 'CHAMP data preprocessing')
294
295 ! 9.2.1 Correct L2 amplitude if not defined
296
297 IF ( ANY(ro_data%Lev1a%snr_L2p == 0.0_wp) ) THEN
298 ro_data%Lev1a%snr_L2p(:) = ro_data%Lev1a%snr_L1ca(:)
299 ENDIF
300
301 ! 9.3 GRAS
302
303 CASE ('ME')
304
305 CALL message(msg_info, 'GRAS data preprocessing')
306
307 CALL ropp_pp_preprocess_GRASRS(ro_data, config, LCF)
308
309 CALL ropp_pp_cutoff_amplitude(ro_data, LCF, config)
310
311
312 WRITE(nstr, '(i6)') ro_data%lev1a%npoints
313 CALL message(msg_diag, 'Merged RS+CL data size: ' // nstr)
314
315 CASE default
316
317 CALL message(msg_warn, 'Occultation LEO id '// TRIM(ro_data%leo_id) // &
318 ' not recognised \n. ' // &
319 ' No mission-specific data pre-processing conducted.')
320
321 END SELECT
322
323!-------------------------------------------------------------------------------
324! 10. Compute model excess phase
325!-------------------------------------------------------------------------------
326 IF (ro_data%Lev1a%Npoints .lt. 100) THEN
327 CALL message(msg_warn, "Error: Too few data")
328 config%obs_ok = .FALSE.
329 RETURN
330 ENDIF
331
332 n = ro_data%Lev1a%Npoints
333 ALLOCATE(phase_LM(n))
334 ALLOCATE(impact_LM(n))
335
336
337 CALL ropp_pp_modelphase(ro_data%dtocc%month, ro_data%georef%lat, &
338 ro_data%georef%lon, &
339 ro_data%Lev1a%dtime, ro_data%lev1a%r_leo, &
340 ro_data%lev1a%r_gns, ro_data%georef%r_coc, &
341 ro_data%georef%roc, phase_LM, impact_LM, config)
342
343write(*,*) 'here'
344!-------------------------------------------------------------------------------
345! 11. Open loop processing (GRAS and COSMIC missions only)
346!-------------------------------------------------------------------------------
347
348 SELECT CASE (ro_data%leo_id(1:2))
349
350 CASE ('ME')
351
352 CALL message(msg_info,'GRAS data: openloop preprocessing')
353
354 WHERE(MOD(LCF(:),2) /= 0)
355 ro_data%Lev1a%phase_L2(:) = phase_LM(:)
356 ENDWHERE
357
358 WHERE(ro_data%Lev1a%phase_L2(:) < ropp_MDTV)
359 ro_data%Lev1a%phase_L2(:) = -1.0_wp
360 ro_data%Lev1a%phase_L2(:) = phase_LM(:)
361 ENDWHERE
362
363 secs_past_hour = 0.0_wp
364 IF(ro_data%DTocc%Minute .LE. 59) &
365 secs_past_hour = secs_past_hour + ro_data%DTocc%Minute*60.0_wp
366 IF(ro_data%DTocc%Second .LE. 59) &
367 secs_past_hour = secs_past_hour + ro_data%DTocc%Second*1.0_wp
368 IF(ro_data%DTocc%Msec .LT. 999) & ! Sometimes msec not set, so defaults to 999
369 secs_past_hour = secs_past_hour + ro_data%DTocc%Msec*1.0e-3_wp
370
371 CALL ropp_pp_openloop(ro_data%Lev1a%dtime+secs_past_hour, &
372 ro_data%Lev1a%phase_L1, ro_data%Lev1a%phase_L2, &
373 phase_LM, ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
374 ro_data%georef%r_coc, LCF)
375
376 ! Set PCD flag
377 ro_data%PCD = IBSET(ro_data%PCD, PCD_open_loop)
378
379 CASE ('C0','CO')
380 CALL message(msg_info,'COSMIC data: openloop preprocessing')
381
382 CALL ropp_pp_openloop(ro_data%Lev1a%dtime, ro_data%Lev1a%phase_L1, &
383 ro_data%Lev1a%phase_L2, phase_LM, &
384 ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
385 ro_data%georef%r_coc, LCF)
386
387 ! Set PCD flag
388 ro_data%PCD = IBSET(ro_data%PCD, PCD_open_loop)
389 END SELECT
390
391!-------------------------------------------------------------------------------
392! 12. Data cutoff
393!-------------------------------------------------------------------------------
394
395! CALL ropp_pp_cutoff(ro_data, config, phase_LM, impact_LM)
396
397!-------------------------------------------------------------------------------
398! 13. Calculate spectra
399!-------------------------------------------------------------------------------
400
401 CALL message(msg_info, 'Compute frequency-time spectra')
402 CALL ropp_pp_spectra(ro_data%Lev1a%dtime, ro_data%Lev1a%phase_L1, &
403 ro_data%Lev1a%phase_L2, phase_LM, impact_LM, config, &
404 OutRO=.true., filnam=ofile)
405
406 ALLOCATE(A(2,ro_data%Lev1a%npoints))
407 ALLOCATE(S(2,ro_data%Lev1a%npoints))
408 A(1,:) = ro_data%Lev1a%snr_L1ca(:)
409 A(2,:) = ro_data%Lev1a%snr_L2p(:)
410 S(1,:) = ro_data%Lev1a%phase_L1(:)
411 S(2,:) = ro_data%Lev1a%phase_L2(:)
412
413 CALL message(msg_info, 'Compute bending-impact spectra')
414 CALL ropp_pp_radiooptic_analysis(ro_data%Lev1a%dtime, ro_data%lev1a%r_leo, &
415 ro_data%lev1a%r_gns, ro_data%georef%r_coc, ro_data%georef%roc, &
416 phase_LM, S, A, OutRO=.true., filnam=ofile)
417
418 DEALLOCATE(A)
419 DEALLOCATE(S)
420
421!-------------------------------------------------------------------------------
422! 14. Clean up
423!-------------------------------------------------------------------------------
424
425 IF (config%obs_ok) THEN
426 CALL message(msg_info, &
427 "Computed spectra successfully \n")
428 CALL message(msg_info, &
429 "Output spectra to files " // TRIM(ro_data%occ_id) // "* \n")
430 CALL message(msg_info, &
431 "See tests/it_pp_plot_spectra_dt.pro and tests/it_pp_plot_spectra_ep.pro for plotting \n")
432
433 ELSE
434 CALL message(msg_warn, &
435 "Error with input data - cannot compute spectra. " // &
436 "Program terminated. \n")
437
438 ENDIF
439
440 CALL ropp_io_free(ro_data)
441 DEALLOCATE(LCF)
442 DEALLOCATE(impact_LM)
443 DEALLOCATE(phase_LM)
444
445 END DO
446 END DO
447
448 CONTAINS
449
450!-------------------------------------------------------------------------------
451! 23. Usage information
452!-------------------------------------------------------------------------------
453
454 SUBROUTINE usage()
455 PRINT *, 'Purpose:'
456 PRINT *, ' Pre-processor spectra tool - calculate local spatial spectra'
457 PRINT *, ' from L1 and L2 phase and amplitude data'
458 PRINT *, 'Usage:'
459 PRINT *, ' > ropp_pp_spectra_tool [<options>] <input_file(s)>'
460 PRINT *, 'Options:'
461 PRINT *, ' -c <config_file> name of configuration file'
462 PRINT *, ' -o <output_file> name of spectra output file'
463 PRINT *, ' (default: ropp_pp_spectra_<type>.dat)'
464 PRINT *, ' -h this help'
465 PRINT *, ' -v version information'
466 PRINT *, ''
467 END SUBROUTINE usage
468
469!-------------------------------------------------------------------------------
470! 24. Version information
471!-------------------------------------------------------------------------------
472
473 SUBROUTINE version_info()
474 CHARACTER (LEN=40) :: version
475 version = ropp_pp_version()
476 PRINT *, 'ropp_pp_spectra_tool - Pre-processor spectra tool: '
477 PRINT *, ' Calculate local spatial spectra from'
478 PRINT *, ' L1 and L2 phase and amplitude data'
479 PRINT *, ''
480 PRINT *, 'This program is part of ROPP (PP) Release ' // TRIM(version)
481 PRINT *, ''
482 END SUBROUTINE version_info
483
484END PROGRAM ropp_pp_spectra_tool