Ticket #255: ropp_pp_occ_tool.f90

File ropp_pp_occ_tool.f90, 41.7 KB (added by Ian Culverwell, 14 years ago)
Line 
1! $Id: ropp_pp_occ_tool.f90 2021 2009-01-16 10:49:04Z frhl $
2
3PROGRAM ropp_pp_occ_tool
4
5!****p* Programs/ropp_pp_occ_tool *
6!
7! NAME
8! ropp_pp_occ_tool - Pre-processing tool to calculate refractivity profile
9! from L1 and L2 excess phase radio occultation data using
10! geometric optics or wave optics, ionospheric correction,
11! statistical optimization and Abel transform.
12!
13! SYNOPSIS
14! ropp_pp_occ_tool <infile(s)> -o <outfile> [-d]
15!
16! DESCRIPTION
17! This program reads RO L1 and L2 excess phase data as a function of time
18! from the input data files and calculates vertical profiles L1 and L2
19! bending angles, ionospheric corrected bending angle and refractivity. The
20! result is written to an ROPP formatted output file.
21!
22! ARGUMENTS
23! <infile(s)> One (or more) input file names.
24!
25! -o <outfile> Name of output file (default: ropp_pp_occ.nc).
26!
27! NOTES
28! If the input file is a multifile, or more than one input files are
29! specified, the output file is a multifile.
30!
31! Already existing output files will be overwritten.
32!
33! EXAMPLE
34! To calculate bending angle and refractivity from one of the example
35! (single-) files in the data directory:
36!
37! ropp_pp_occ_tool ../data/input.nc -o example_01.nc
38!
39! To calculate bending angle and refractivity profiles from all singlefiles
40! in the data directory:
41!
42! ropp_pp_occ_tool ../data/*.nc -o example_02.nc
43!
44! Note that the resulting example_02.nc file contains processed data from
45! all example profiles.
46!
47! To calculate bending angle and refractivity profiles from
48! all profiles contained in the multifile multi.nc:
49!
50! ropp_pp_occ_tool ../data/multi.nc -o example_03.nc
51!
52! Since the multi_* file was generated by concatenating the other files
53! in the data directory, example_02.nc and example_03.nc should be identical
54! apart from the file names.
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, L1btype, L2atype
76 USE ropp_pp
77 USE ropp_pp_preproc
78 USE ropp_pp_types, ONLY: PPConfig, PPDiag
79
80 IMPLICIT NONE
81
82 TYPE(ROprof) :: ro_data ! Input RO data
83 TYPE(L1atype) :: Lev1a ! Temporary Level1a structure for storage
84 TYPE(L1btype) :: bangle ! Output bending angle profiles
85 TYPE(L1btype) :: out_ba ! Corrected bending angles
86 TYPE(L1btype) :: smt_ba ! Smoothed bending angle
87 TYPE(L1btype) :: mod_ba ! Model bending angle
88 TYPE(L2atype) :: out_refrac ! Retrieved refractivity profile
89
90 TYPE(ppConfig) :: config ! Configuration options
91 TYPE(ppDiag) :: diag ! Diagnostic variables
92
93 INTEGER :: idummy, imin, imax
94 INTEGER :: i, iargc, argc, k, j
95 INTEGER :: n_files, n_profiles, nbi
96 INTEGER :: ws_go_smooth, ws_go_full, ws_wo, ws_low
97 REAL(wp) :: Pmin, Pmax
98 LOGICAL :: give_help
99 LOGICAL :: earth_ellipsoid = .FALSE.
100 LOGICAL :: output_full = .FALSE.
101 LOGICAL :: twofit = .FALSE.
102 LOGICAL :: ranchk = .TRUE.
103
104 REAL(wp), allocatable :: WLC(:) ! SO-weight of linear combination !HGL!
105
106 CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
107 CHARACTER(len = 4096) :: ofile
108 CHARACTER(len = 4096) :: cfg_file = " "
109 CHARACTER(len = 256) :: buffer
110 CHARACTER(len = 64) :: version
111 CHARACTER(len = 64) :: outstr
112 CHARACTER(len = 80) :: mfile = " "
113 CHARACTER(len = 80) :: bfile = " "
114 CHARACTER(len = 80) :: navfile = " "
115 CHARACTER(len = 10) :: method = " "
116 CHARACTER(len = 10) :: occmethod = " "
117 CHARACTER(len = 10) :: filtmethod = " "
118 CHARACTER(len = 4) :: istr
119 CHARACTER(len = 10) :: nstr
120 CHARACTER(len = 10) :: pstr1
121 CHARACTER(len = 10) :: pstr2
122
123!-------------------------------------------------------------------------------
124! 2. Default settings
125!-------------------------------------------------------------------------------
126
127 version = ADJUSTL("5.0")
128 give_help = .FALSE.
129 ofile = "ropp_pp_occ.nc"
130
131 CALL message(msg_noin, '')
132 CALL message(msg_noin, &
133 '----------------------------------------------------------------------')
134 CALL message(msg_noin, &
135 ' ROPP Occultation Pre-processor Tool - v' // TRIM(version))
136 CALL message(msg_noin, &
137 '----------------------------------------------------------------------')
138 CALL message(msg_noin, '')
139
140!-------------------------------------------------------------------------------
141! 3. Command line arguments
142!-------------------------------------------------------------------------------
143
144 argc = iargc()
145 i = 1
146 n_files = 0
147 ALLOCATE(ifiles(argc))
148
149 DO WHILE(i <= argc)
150 CALL getarg(i, buffer)
151 SELECT CASE (buffer)
152 CASE('-o') ! Output file name (netCDF output)
153 CALL getarg(i+1, buffer)
154 ofile = buffer
155 i = i + 1
156 CASE('-c') ! Configuration file name
157 CALL getarg(i+1, buffer)
158 cfg_file = buffer
159 i = i + 1
160 CASE('-m') ! Ionospheric correction method
161 CALL getarg(i+1, buffer)
162 method = buffer
163 i = i + 1
164 CASE('-mfile') ! Model refractivity file
165 CALL getarg(i+1, buffer)
166 mfile = buffer
167 i = i + 1
168 CASE('-bfile') ! Background atmosphere profile file
169 CALL getarg(i+1, buffer)
170 bfile = buffer
171 i = i + 1
172 CASE('-navfile') ! External navigation bit file
173 CALL getarg(i+1, buffer)
174 navfile = buffer
175 i = i + 1
176 CASE('-occ', '--occ') ! Occultation processing method
177 CALL getarg(i+1, buffer)
178 occmethod = buffer
179 i = i + 1
180 CASE('-filter', '--filter') ! Filtering method
181 CALL getarg(i+1, buffer)
182 filtmethod = buffer
183 i = i + 1
184 CASE('-fit', '--fit') ! 2 parameter fitting method
185 twofit = .true.
186 CASE('-ellipsoid') ! Output height wrt reference ellipsoid
187 earth_ellipsoid = .TRUE.
188 CASE('-full') ! Output full height grid (IC plus MSIS)
189 output_full = .TRUE.
190 CASE('-d') ! 'Diagnostic' output mode
191 msg_MODE = VerboseMode
192 CASE ('-no_ranchk') ! Use no rangecheck on output
193 ranchk = .FALSE.
194 CASE('-h', '-help', '--help') ! Give some help
195 give_help = .TRUE.
196 CASE('-V', '-version', '--version') ! Give some version information
197 CALL version_info(version)
198 CALL EXIT(0)
199 CASE default ! Input file name
200 n_files = n_files + 1
201 ifiles(n_files) = buffer
202 END SELECT
203 i = i + 1
204 END DO
205
206 IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN
207 CALL usage()
208 CALL EXIT(0)
209 ENDIF
210
211 ! 3.1 Read configuration file (if exists), preserving command-line options
212
213 IF (method /= " ") config%method = method
214 IF (mfile /= " ") config%mfile = mfile
215 IF (bfile /= " ") config%bfile = bfile
216 IF (navfile /= " ") config%navbit_file = navfile
217 IF (occmethod /= " ") config%occ_method = occmethod
218 IF (filtmethod /= " ") config%filter_method = filtmethod
219
220 IF (cfg_file /= " ") THEN
221 CALL message(msg_info, &
222 "Reading configuration file " // TRIM(cfg_file) // ".\n")
223
224 CALL ropp_pp_read_config(cfg_file, config)
225 IF (method /= " ") config%method = method
226 IF (mfile /= " ") config%mfile = mfile
227 IF (bfile /= " ") config%bfile = bfile
228 IF (navfile /= " ") config%navbit_file = navfile
229 IF (occmethod /= " ") config%occ_method = occmethod
230 IF (filtmethod /= " ") config%filter_method = filtmethod
231 IF (twofit) config%nparm_fit = 2
232 ENDIF
233
234!-------------------------------------------------------------------------------
235! 4. Remove pre-existing output file
236!-------------------------------------------------------------------------------
237
238 CALL file_delete(ofile, idummy)
239
240!-------------------------------------------------------------------------------
241! 5. Loop over all input files and profiles
242!-------------------------------------------------------------------------------
243
244 DO k = 1, n_files
245
246 n_profiles = ropp_io_nrec(ifiles(k))
247
248 DO i = 1, n_profiles
249
250 WRITE(istr, '(i4)') i
251 WRITE(nstr, '(i6)') n_profiles
252 CALL message(msg_info, "Processing profile " // istr // " of " // nstr )
253
254!-------------------------------------------------------------------------------
255! 6. Read data
256!-------------------------------------------------------------------------------
257
258! 6.1 Read input file
259
260 CALL message(msg_info, &
261 "Reading input data file " // TRIM(ifiles(k)) // ".\n")
262 CALL ropp_io_read(ro_data, ifiles(k), rec = i, ranchk = .TRUE.)
263 CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")
264
265 IF (ro_data%Lev1a%Npoints == 0) THEN
266 config%obs_ok = .FALSE.
267 CALL ropp_io_free(ro_data%vlist)
268 CALL ropp_io_write(ro_data, ofile, append = .TRUE., ranchk = ranchk )
269 CALL message(msg_fatal, "FAILURE: No Level1a data in file " // &
270 ifiles(k) // "\n")
271 ENDIF
272
273! 6.2 Shrink ro_data to correct size (for multiple profiles)
274
275 CALL ropp_io_init(Lev1a, ro_data%lev1a%npoints)
276 Lev1a = ro_data%lev1a
277 ro_data%lev1a = Lev1a
278 CALL ropp_io_free(Lev1a)
279 CALL ropp_io_free(ro_data%Lev1b)
280
281!-------------------------------------------------------------------------------
282! 7. Check coordinate frames
283!-------------------------------------------------------------------------------
284
285 CALL ropp_pp_set_coordinates(ro_data)
286
287!-------------------------------------------------------------------------------
288! 8. Determine occultation point georeference information
289!-------------------------------------------------------------------------------
290
291 CALL occ_point( ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
292 ro_data%georef%lat, ro_data%georef%lon, &
293 ro_data%georef%r_coc, ro_data%georef%roc, &
294 ro_data%georef%azimuth, &
295 ro_data%georef%undulation, &
296 config%egm96, config%corr_egm96)
297 WRITE(outstr,'(2(A,1x,f7.2))') &
298 "Occultation point: Lat= ", ro_data%georef%lat, &
299 ", Lon= ", ro_data%georef%lon
300 CALL message(msg_diag, outstr)
301
302!-------------------------------------------------------------------------------
303! 9. Mission-specific pre-processing and input data cut-off
304!-------------------------------------------------------------------------------
305
306 CALL ropp_pp_preprocess(ro_data, config, diag)
307
308!-------------------------------------------------------------------------------
309! 10. Initialise bending angle structures
310!-------------------------------------------------------------------------------
311
312 IF (config%obs_ok) THEN
313
314 CALL ropp_io_init(smt_ba, ro_data%Lev1a%Npoints)
315 CALL ropp_io_init(bangle, ro_data%Lev1a%Npoints)
316
317 CALL tangent_point(ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
318 bangle%lat_tp, bangle%lon_tp, bangle%azimuth_tp)
319
320!-------------------------------------------------------------------------------
321! 11. Geometric optics processing
322!-------------------------------------------------------------------------------
323
324 CALL message(msg_info, &
325 "Retrieving bending angle profiles by GEOMETRIC OPTICS \n")
326
327
328! 11.1 Calculate smoothed bending angle
329
330 DO j=1,ro_data%Lev1a%Npoints
331 smt_ba%impact(j) = &
332 impact_parameter(ro_data%Lev1a%r_leo(j,:)-ro_data%georef%r_coc, &
333 ro_data%Lev1a%r_gns(j,:)-ro_data%georef%r_coc )
334 ENDDO
335 Pmax = MAXVAL(smt_ba%impact)
336 Pmin = MAX(MINVAL(smt_ba%impact), ro_data%georef%roc+2000.0_wp)
337 ws_go_smooth = CEILING( config%fw_go_smooth * &
338 (ro_data%Lev1a%Npoints-1)/ABS(Pmax-Pmin))
339
340 WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc
341 WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc
342 CALL message(msg_diag, "Smoothed bending angle. ")
343 CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
344 WRITE(nstr, '(i8)') ws_go_smooth
345 CALL message(msg_diag, "ws_go_smooth " // nstr)
346
347 IF (Pmax > Pmin) THEN
348
349 CALL ropp_pp_bending_angle_go( ro_data%Lev1a%dtime, &
350 ro_data%Lev1a%r_leo, &
351 ro_data%Lev1a%r_gns, &
352 ro_data%georef%r_coc, &
353 ro_data%Lev1a%phase_L1, &
354 ro_data%Lev1a%phase_L2, &
355 ws_go_smooth, &
356 config%filter_method, &
357 smt_ba%impact_L1, &
358 smt_ba%bangle_L1, &
359 smt_ba%impact_L2, &
360 smt_ba%bangle_L2 )
361
362
363 CALL ropp_pp_linear_combination( smt_ba%impact_L1, &
364 smt_ba%bangle_L1, &
365 smt_ba%impact_L2, &
366 smt_ba%bangle_L2, &
367 smt_ba%impact, &
368 smt_ba%bangle )
369
370 ELSE
371
372 config%obs_ok = .FALSE.
373 CALL message(msg_warn, "Cannot process profile (smt): Pmin > Pmax")
374
375 ENDIF
376
377! 11.2 Calculate bending angle
378
379 bangle%impact = smt_ba%impact_L1
380 Pmax = MAXVAL(smt_ba%impact_L1)
381 Pmin = MAX(MINVAL(smt_ba%impact_L1), ro_data%georef%roc+2000.0_wp)
382 ws_go_full = CEILING(config%fw_go_full*(ro_data%Lev1a%Npoints - 1)/ &
383 ABS(Pmax - Pmin))
384
385 WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc
386 WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc
387 CALL message(msg_diag, "Full resolution bending angle. ")
388 CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
389 WRITE(nstr, '(i8)') ws_go_full
390 CALL message(msg_diag, "ws_go_full " // nstr)
391
392 IF (Pmax > Pmin) THEN
393
394 CALL ropp_pp_bending_angle_go(ro_data%Lev1a%dtime, &
395 ro_data%Lev1a%r_leo, &
396 ro_data%Lev1a%r_gns, &
397 ro_data%georef%r_coc, &
398 ro_data%Lev1a%phase_L1, &
399 ro_data%Lev1a%phase_L2, &
400 ws_go_full, &
401 config%filter_method, &
402 bangle%impact_L1, &
403 bangle%bangle_L1, &
404 bangle%impact_L2, &
405 bangle%bangle_L2 )
406
407 WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints
408 WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp
409 WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp
410 CALL message(msg_info, &
411 TRIM(nstr) // " data points in output between " // &
412 TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n")
413
414 ELSE
415
416 config%obs_ok = .FALSE.
417 CALL message(msg_warn, "Cannot process profile (full): Pmin > Pmax")
418
419 ENDIF
420
421!-------------------------------------------------------------------------------
422! 12. Wave optics processing
423!-------------------------------------------------------------------------------
424
425 IF ( INDEX(config%occ_method, "WO" ) == 1 ) THEN
426
427 CALL message(msg_info, &
428 "Retrieving bending angle profiles by WAVE OPTICS \n")
429
430 ws_wo = CEILING(config%fw_wo * (bangle%npoints-1.0_wp) / &
431 ABS(Pmax - Pmin))
432 ws_low = CEILING(config%fw_low * (bangle%npoints-1.0_wp) / &
433 ABS(Pmax - Pmin))
434
435 WRITE(nstr,'(I5)') ws_go_smooth
436 CALL message(msg_diag, 'WM = ' // nstr)
437 WRITE(nstr,'(I5)') ws_wo
438 CALL message(msg_diag, 'W = ' // nstr)
439 WRITE(nstr,'(I5)') ws_low
440 CALL message(msg_diag, 'WL = ' // nstr)
441
442 IF (Pmax > Pmin) THEN
443
444 CALL ropp_pp_bending_angle_wo( ro_data%Lev1a%dtime, &
445 ro_data%Lev1a%r_leo, &
446 ro_data%Lev1a%r_gns, &
447 ro_data%georef%r_coc, &
448 ro_data%georef%roc, &
449 ro_data%Lev1a%phase_L1, &
450 ro_data%Lev1a%phase_L2, &
451 ro_data%Lev1a%snr_L1ca, &
452 ro_data%Lev1a%snr_L2p, &
453 ws_go_smooth, &
454 ws_wo, &
455 ws_low, &
456 config%hmax_wo, &
457 config%filter_method, &
458 config%opt_DL2, &
459 config%cff, &
460 config%dsh, &
461 bangle%impact_L1, &
462 bangle%bangle_L1, &
463 bangle%bangle_L1_sigma, &
464 bangle%impact_L2, &
465 bangle%bangle_L2, &
466 bangle%bangle_L2_sigma, &
467 diag )
468
469 WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints
470 WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp
471 WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp
472 CALL message(msg_info, &
473 TRIM(nstr) // " data points in output between " // &
474 TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n")
475
476 ELSE
477
478 config%obs_ok = .FALSE.
479 CALL message(msg_warn, "Cannot process profile (wo): Pmin > Pmax")
480
481 ENDIF
482
483
484 ENDIF
485
486!-------------------------------------------------------------------------------
487! 13. Check that profiles are monotonously increasing in height
488! (index 1 towards surface) and valid bending angles
489!-------------------------------------------------------------------------------
490
491 CALL ropp_pp_monotonous(bangle%impact_L1, -1)
492 CALL ropp_pp_monotonous(bangle%impact_L2, -1)
493 CALL ropp_pp_monotonous(smt_ba%impact, -1)
494
495 IF ( ANY(bangle%bangle_L1 <= ropp_MDTV) .OR. &
496 ANY(bangle%bangle_L2 <= ropp_MDTV) ) THEN
497 config%obs_ok = .FALSE.
498 CALL message(msg_warn, "Cannot process profile: " // &
499 "Invalid bending angles computed \n")
500 ENDIF
501
502 ENDIF
503
504!-------------------------------------------------------------------------------
505! 14. Compute un-optimised ionospheric corrected bending angle profile (LC)
506!-------------------------------------------------------------------------------
507
508 IF (config%obs_ok) THEN
509
510 CALL ropp_pp_linear_combination( bangle%impact_L1, bangle%bangle_L1, &
511 bangle%impact_L2, bangle%bangle_L2, &
512 bangle%impact, bangle%bangle )
513
514!-------------------------------------------------------------------------------
515! 15. Define configuration parameters
516!-------------------------------------------------------------------------------
517
518 config%r_curve = ro_data%georef%roc
519 config%npoints = bangle%npoints
520
521 config%Pmax = MAXVAL(bangle%impact_L1(:))
522 config%Pmin = MINVAL(bangle%impact_L1(:))
523
524!-------------------------------------------------------------------------------
525! 16. Interpolate input data to standard grid
526!-------------------------------------------------------------------------------
527
528 ! 16.1 Calculate size of standard grid
529
530 IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
531 Pmax = MAXVAL(bangle%impact_L1(:)) - 5000.0_wp
532 ELSE
533 Pmax = config%ztop_invert + config%r_curve
534 ENDIF
535 Pmin = MINVAL(bangle%impact_L1(:))
536 nbi = 1 + CEILING((Pmax - Pmin)/config%dpi)
537
538 WRITE(pstr2, '(f10.3)') Pmax-config%r_curve
539 WRITE(pstr1, '(f10.3)') Pmin-config%r_curve
540 CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
541 WRITE(nstr, '(I6)') nbi
542 CALL message(msg_diag, "Standard grid size nbi = " // nstr)
543
544! 16.2 Initialise standard grid data structures for output
545
546 CALL ropp_io_init(out_ba, nbi)
547 CALL ropp_io_init(out_refrac, nbi)
548
549! 16.3 Interpolate input data onto standard grid
550
551 CALL ropp_pp_merge_profile(bangle%impact_L1, bangle%bangle_L1, &
552 bangle%impact_L2, bangle%bangle_L2, &
553 out_ba%impact_L1, out_ba%bangle_L1, &
554 out_ba%impact_L2, out_ba%bangle_L2, &
555 Pmin, Pmax)
556
557 CALL ropp_pp_merge_profile(bangle%impact_L1,bangle%bangle_L1_sigma, &
558 bangle%impact_L2,bangle%bangle_L2_sigma, &
559 out_ba%impact_L1,out_ba%bangle_L1_sigma, &
560 out_ba%impact_L2,out_ba%bangle_L2_sigma, &
561 Pmin, Pmax)
562
563 ALLOCATE(WLC(nbi)) !HGL!
564 WLC(:) = 1.0 !HGL!
565
566
567!-------------------------------------------------------------------------------
568! 17. Ionospheric correction of bending angle profile by linear combination
569!-------------------------------------------------------------------------------
570
571 IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
572
573 CALL message(msg_info, &
574 "Correcting bending angle profile by LINEAR COMBINATION \n")
575 ro_data%bangle_method = TRIM(config%occ_method) // "(LC)"
576
577 CALL ropp_pp_linear_combination( out_ba%impact_L1, &
578 out_ba%bangle_L1, &
579 out_ba%impact_L2, &
580 out_ba%bangle_L2, &
581 out_ba%impact_opt, &
582 out_ba%bangle_opt )
583
584 imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0))))
585 imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin)))
586
587!-------------------------------------------------------------------------------
588! 18. Ionospheric correct bending angle profile by statistical optimization
589!-------------------------------------------------------------------------------
590
591 ELSE
592
593 CALL message(msg_info, &
594 "Correcting bending angle profile by STATISTICAL OPTIMISATION \n")
595 ro_data%bangle_method = TRIM(config%occ_method) // " (StatOpt)"
596
597! 18.1 Retrieve model bending angles on L1 impact parameter levels
598
599 CALL ropp_io_init(mod_ba, nbi)
600
601 IF (INDEX(config%method, "GMSIS" ) == 1 ) THEN
602 CALL ropp_pp_search_model_refraction( config%mfile, &
603 ro_data%dtocc%month, &
604 ro_data%georef%lat, &
605 ro_data%georef%lon, &
606 smt_ba%impact, &
607 smt_ba%bangle, &
608 out_ba%impact_L1, &
609 mod_ba%bangle, &
610 config )
611
612 ELSE IF (INDEX(config%method, "MSIS" ) == 1 ) THEN
613 CALL ropp_pp_model_refraction( config%mfile, &
614 ro_data%dtocc%month, &
615 ro_data%georef%lat, &
616 ro_data%georef%lon, &
617 out_ba%impact_L1, &
618 mod_ba%bangle, &
619 config )
620
621
622 ELSE IF (INDEX(config%method, "BG" ) == 1 ) THEN
623 CALL ropp_pp_bg_refraction( config%bfile, &
624 ro_data%dtocc%month, &
625 ro_data%georef%lat, &
626 ro_data%georef%lon, &
627 out_ba%impact_L1, &
628 mod_ba%bangle, &
629 config )
630
631 ELSE
632 config%obs_ok = .FALSE.
633 CALL message(msg_warn, &
634 "Statistical optimisation method " // config%method // &
635 " not supported")
636 EXIT
637 ENDIF
638
639! 18.2 Fit model bending angle with (smoothed) observed bending angles
640
641 CALL ropp_pp_fit_model_refraction( smt_ba%impact, &
642 smt_ba%bangle, &
643 out_ba%impact_L1, &
644 mod_ba%bangle, &
645 config )
646
647! 18.3 Perform ionospheric correction with statistical optimization
648
649 imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0))))
650 imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin)))
651
652 WLC(:) = 0.0 !HGL!
653
654 ! Merge observed bending angles with model above observation top
655 WHERE ( out_ba%impact_L1 > config%Pmax - 5000.0_wp )
656 out_ba%bangle_L1(:) = mod_ba%bangle(:)
657 out_ba%bangle_L2(:) = mod_ba%bangle(:)
658 out_ba%bangle_L1_sigma(:) = 0.0_wp
659 out_ba%bangle_L2_sigma(:) = 0.0_wp
660 END WHERE
661
662 IF (imin < imax) THEN
663 out_ba%impact_opt(:) = out_ba%impact_L1(:)
664 out_ba%bangle_opt(:) = out_ba%bangle_L1(:)
665 CALL ropp_pp_ionospheric_correction( &
666 out_ba%impact_L1(imin:imax), & ! L1 impact
667 out_ba%bangle_L1(imin:imax), & ! L1 bending
668 out_ba%impact_L2(imin:imax), & ! L2 impact
669 out_ba%bangle_L2(imin:imax), & ! L2 bending
670 out_ba%impact_L1(imin:imax), & ! Model impact
671 mod_ba%bangle(imin:imax), & ! Model bending
672 config, & ! Configuration
673 out_ba%impact_opt(imin:imax), & ! Opt impact
674 out_ba%bangle_opt(imin:imax), & ! Opt bending
675 diag, & ! Diagnostics !HGL!
676 WLC(imin:imax)) ! Weight of the linear combination of L1 and L2 !HGL!
677 diag%err_neut = diag%err_neut + &
678 out_ba%bangle_L1_sigma(imin:imax)**2
679 diag%sq = 100.0_wp * &
680 MAXVAL(SQRT(diag%err_neut(:))/out_ba%bangle_opt(imin:imax))
681 ELSE
682 config%obs_ok = .FALSE.
683 CALL message(msg_warn, &
684 "Cannot perform ionospheric correction: Imin >= Imax \n")
685 ENDIF
686
687 ENDIF
688
689! 18.4 Interpolate LC bending angles, lat_tp, lon_tp, azimuth_tp to output grid
690
691 out_ba%impact = out_ba%impact_opt
692 CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
693 bangle%bangle, out_ba%bangle)
694 CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
695 bangle%lat_tp, out_ba%lat_tp)
696 CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
697 bangle%lon_tp, out_ba%lon_tp)
698 CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
699 bangle%azimuth_tp, out_ba%azimuth_tp)
700
701 WRITE(nstr, '(i10)') imax-imin+1
702 WRITE(pstr2, '(f6.1)') out_ba%impact(imax)/1000.0_wp
703 WRITE(pstr1, '(f6.1)') out_ba%impact(imin)/1000.0_wp
704 CALL message(msg_info, &
705 TRIM(nstr) // " data points in output between " // &
706 TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n")
707
708 END IF
709
710!-------------------------------------------------------------------------------
711! 19. Perform inverse Abel transform to compute refractivity
712!-------------------------------------------------------------------------------
713
714 IF (config%obs_ok) THEN
715
716 ro_data%refrac_method = "ABEL (" // TRIM(config%abel) // ")"
717
718! 19.1 Abel inversion of corrected bending angle profile
719
720 IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
721
722 CALL message(msg_info, &
723 "Retrieving refractivity profile by " // TRIM(config%abel) // &
724 " ABEL TRANSFORM \n")
725
726 IF ( INDEX(config%abel, "LIN" ) == 1 ) THEN
727 CALL ropp_pp_invert_LIN(out_ba%impact_opt, out_ba%bangle_opt, &
728 out_ba%impact_opt, out_refrac%refrac)
729 ELSE IF ( INDEX(config%abel, "EXP" ) == 1 ) THEN
730 CALL ropp_pp_invert_EXP(out_ba%impact_opt, out_ba%bangle_opt, &
731 out_ba%impact_opt, out_refrac%refrac)
732 ELSE
733 config%obs_ok = .FALSE.
734 CALL message(msg_warn, &
735 "Abel integral method " // config%abel // " not supported")
736 ENDIF
737
738! 19.2 Abel inversion of corrected bending angle profile with stat opt
739
740 ELSE
741
742 CALL message(msg_info, &
743 "Retrieving refractivity profile by " // TRIM(config%abel) // &
744 " ABEL TRANSFORM with STAT OPT \n")
745
746 CALL ropp_pp_invert_refraction( config%mfile, &
747 ro_data%dtocc%month, &
748 ro_data%georef%lat, &
749 ro_data%georef%lon, &
750 out_ba%impact_opt, &
751 out_ba%bangle_opt, &
752 out_refrac%geop_refrac, &
753 out_refrac%refrac, &
754 config )
755
756 ENDIF
757
758!-------------------------------------------------------------------------------
759! 20. Compute output height scales
760!-------------------------------------------------------------------------------
761
762 ! Override default output height scales wrt geoid if requested
763 IF (earth_ellipsoid) ro_data%georef%undulation = ropp_MDFV
764
765 IF (ro_data%georef%undulation > ropp_MDTV) THEN
766 CALL message(msg_info, "Writing output altitude scales " // &
767 "with respect to EGM96 GEOID")
768 out_refrac%alt_refrac = &
769 ((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) &
770 - (ro_data%GEOref%roc + ro_data%GEOref%undulation))
771 ELSE
772 IF (.not. earth_ellipsoid) &
773 CALL message(msg_warn, "Invalid undulation calculated. " // &
774 "Check for valid EGM geoid coefficient and correction file.")
775 CALL message(msg_info, "Writing output altitude scales " // &
776 "with respect to WGS84 ELLIPSOID.")
777 out_refrac%alt_refrac = &
778 ((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) &
779 - ro_data%GEOref%roc)
780 ENDIF
781 out_refrac%geop_refrac = &
782 geometric2geopotential(ro_data%georef%lat, out_refrac%alt_refrac)
783
784!-------------------------------------------------------------------------------
785! 21. Copy retrieved profiles to RO structure
786! - only output data within observed range
787!-------------------------------------------------------------------------------
788
789 ! Override default output only within observed range if requested
790 IF (.not. output_full) THEN
791 out_ba%npoints = Imax - Imin + 1
792 out_refrac%npoints = Imax - Imin + 1
793 ENDIF
794
795 CALL ropp_io_roprof2roprof(out_ba, ro_data%lev1b)
796 CALL ropp_io_roprof2roprof(out_refrac, ro_data%lev2a)
797
798
799!-------------------------------------------------------------------------------
800! 22. Compute Tdry and Pdry
801!-------------------------------------------------------------------------------
802
803 IF (config%output_tdry) THEN
804 CALL message(msg_info, "Computing dry temperature and pressure \n")
805 CALL ropp_io_init(ro_data%lev2b, ro_data%lev2a%npoints)
806 ro_data%lev2b%geop = ro_data%lev2a%geop_refrac
807 ro_data%lev2b%shum = 0.0_wp
808 CALL ropp_pp_tdry(ro_data%georef%lat, ro_data%lev2a%alt_refrac, &
809 ro_data%lev2a%refrac, ro_data%lev2b%shum, ro_data%lev2b%temp, &
810 ro_data%lev2b%press, config%ztop_invert)
811 ro_data%lev2b%press = ro_data%lev2b%press / 100.0_wp ! convert to hPa
812 ENDIF
813 ENDIF
814
815!-------------------------------------------------------------------------------
816! 23. Write data
817!-------------------------------------------------------------------------------
818
819 CALL ropp_io_free(ro_data%vlist) !! interim, avoid writing lcf data
820 IF (config%output_diag) THEN
821 CALL message(msg_info, "Writing additional diagnostic output \n")
822 CALL ropp_pp_diag2roprof(diag, ro_data)
823 ENDIF
824
825 CALL ropp_io_addvar(ro_data, & !HGL!
826 name = "LC_weight", & !HGL!
827 long_name = "SO-weight of observed LC bending angle", & !HGL!
828 units = " ", & !HGL!
829 range = (/0.0_wp,1.0_wp/), & !HGL!
830 data = WLC(imin:imax)) !HGL!
831
832 CALL message(msg_info, "Writing to output file " // TRIM(ofile) // "\n")
833 CALL ropp_io_write(ro_data, ofile, append = .TRUE., ranchk = ranchk )
834
835!-------------------------------------------------------------------------------
836! 24. Clean up
837!-------------------------------------------------------------------------------
838
839 CALL ropp_io_free(ro_data)
840 CALL ropp_io_free(smt_ba)
841 CALL ropp_io_free(bangle)
842 CALL ropp_io_free(out_ba)
843 CALL ropp_io_free(out_refrac)
844 CALL ropp_io_free(mod_ba)
845 IF (ASSOCIATED(diag%ba_ion)) DEALLOCATE(diag%ba_ion)
846 IF (ASSOCIATED(diag%err_ion)) DEALLOCATE(diag%err_ion)
847 IF (ASSOCIATED(diag%err_neut)) DEALLOCATE(diag%err_neut)
848 IF (ASSOCIATED(diag%CTimpact)) DEALLOCATE(diag%CTimpact)
849 IF (ASSOCIATED(diag%CTamplitude)) DEALLOCATE(diag%CTamplitude)
850 IF (ASSOCIATED(diag%CTamplitude_smt)) DEALLOCATE(diag%CTamplitude_smt)
851 IF (ASSOCIATED(diag%CTimpactL2)) DEALLOCATE(diag%CTimpactL2)
852 IF (ASSOCIATED(diag%CTamplitudeL2)) DEALLOCATE(diag%CTamplitudeL2)
853 IF (ASSOCIATED(diag%CTamplitudeL2_smt)) DEALLOCATE(diag%CTamplitudeL2_smt)
854
855 IF (ALLOCATED(WLC)) DEALLOCATE(WLC) !HGL!
856
857 END DO
858 END DO
859 CONTAINS
860
861!-------------------------------------------------------------------------------
862! 25. Usage information
863!-------------------------------------------------------------------------------
864
865 SUBROUTINE usage()
866
867 PRINT *, 'ropp_pp_occ_tool - Pre-processor tool'
868 PRINT *, ' Calculate corrected bending angle and '
869 PRINT *, ' refractivity from L1 and L2 phase data'
870 PRINT *, 'ropp_pp_tool [<options>] <input_file(s)> -o <output_file>'
871 PRINT *, 'Valid options are:'
872 PRINT *, ' -h give (this) help.'
873 PRINT *, ' -o <output_file> name of netCDF/ropp output file.'
874 PRINT *, ' -c <config_file> name of configuration file.'
875 PRINT *, ' -m <method> ionospheric correction method '
876 PRINT *, ' [NONE,MSIS,GMSIS,BG], (default MSIS).'
877 PRINT *, ' -mfile <file> model refractivity coefficients file path'
878 PRINT *, ' (default local search).'
879 PRINT *, ' -bfile <file> background model atmospheric profile file path'
880 PRINT *, ' (if using BG ionospheric correction method).'
881 PRINT *, ' -navfile <file> external navigation bit *_txt file path'
882 PRINT *, ' (default internal correction).'
883 PRINT *, ' -occ <occ_method> processing method, WO or GO (default WO)'
884 PRINT *, ' -filter <filt_method> filtering method, slpoly or optest'
885 PRINT *, ' (default slpoly, sliding polynomial)'
886 PRINT *, ' -fit apply 2-parameter regression fit to model'
887 PRINT *, ' -ellipsoid output height with respect to WGS84 ellipsoid'
888 PRINT *, ' (default output wrt EGM96 geoid)'
889 PRINT *, ' -d output additional diagnostics.'
890 PRINT *, ' -no_ranchk no range check performed before writing output'
891 PRINT *, ' -V give some version information.'
892 PRINT *, ''
893
894 END SUBROUTINE usage
895
896!-------------------------------------------------------------------------------
897! 26. Version information
898!-------------------------------------------------------------------------------
899
900 SUBROUTINE version_info(version)
901 CHARACTER(len = *) :: version
902 PRINT *, 'ropp_pp_occ_tool - Pre-processor occ tool: '
903 PRINT *, ' Calculate corrected bending angles '
904 PRINT *, ' and refractivity from excess phase'
905 PRINT *, ''
906 PRINT *, 'This program is part of ROPP version ' // TRIM(version)
907 PRINT *, ''
908 END SUBROUTINE version_info
909
910END PROGRAM ropp_pp_occ_tool