Context navigation

Back to Ticket #488

Ticket #488: ropp_pp_occ_tool_24042017.f90

File ropp_pp_occ_tool_24042017.f90, 58.5 KB (added by Ian Culverwell, 8 years ago)
ropp_pp_occ_tool_24042017.f90

Line
1	! $Id: ropp_pp_occ_tool.f90 2021 2009-01-16 10:49:04Z frhl $
2
3	PROGRAM ropp_pp_occ_tool
4
5	!***p Programs/ropp_pp_occ_tool *
6	!
7	! NAME
8	! ropp_pp_occ_tool
9	!
10	! SYNOPSIS
11	! Pre-processing tool to calculate refractivity and dry temperature
12	! profile from L1 and L2 excess phase radio occultation data using
13	! geometric optics or wave optics, ionospheric correction, statistical
14	! optimization, Abel transform and hydrostatic equation.
15	!
16	! > ropp_pp_occ_tool [<options>] <infile(s)>
17	!
18	! ARGUMENTS
19	! <infile(s)> One (or more) input file names.
20	!
21	! OPTIONS
22	! -c <config_file> name of configuration file
23	! -o <output_file> name of ROPP netCDF output file
24	! -m <method> ionospheric correction method [NONE,MSIS,GMSIS,BG]
25	! (default GMSIS)
26	! -mfile <file> model refractivity coefficients file
27	! (default 'MSIS_coeff.nc')
28	! -bfile <file> background model atmospheric profile file path
29	! (if using BG ionospheric correction method)
30	! -navfile <file> external navigation bit *_txt file path
31	! (default internal correction)
32	! -occ <method> processing method, WO or GO (default WO)
33	! -filter <method> filtering method, slpoly or optest
34	! (default slpoly, sliding polynomial)
35	! -fit apply 2-parameter regression fit to model
36	! -ellipsoid output height with respect to WGS84 ellipsoid
37	! (default output wrt EGM96 geoid)
38	! -d output additional diagnostics
39	! -h help
40	! -v version information
41	!
42	! DESCRIPTION
43	! This program reads RO L1 and L2 excess phase data as a function of time
44	! from the input data files and calculates vertical profiles L1 and L2
45	! bending angles, ionospheric corrected bending angle, refractivity and
46	! dry temperature. The result is written to a ROPP output file.
47	!
48	! NOTES
49	! If the input file is a multifile, or more than one input files are
50	! specified, the output file is a multifile.
51	!
52	! Already existing output files will be overwritten.
53	!
54	! ERRORS
55	! Program (shell) return codes:
56	! 0 = OK
57	! 1 = At least one Warning occured
58	! 2 = At least one Error occured
59	! 3 = A Fatal error occured
60	!
61	! EXAMPLE
62	! To calculate bending angle, refractivity and dry temperature from one
63	! of the example (single-) files in the data directory:
64	!
65	! > ropp_pp_occ_tool ../data/input.nc -o example_01.nc
66	!
67	! To calculate bending angle, refractivity and dry temperature profiles
68	! from all singlefiles in the data directory:
69	!
70	! > ropp_pp_occ_tool ../data/*.nc -o example_02.nc
71	!
72	! Note that the resulting example_02.nc file contains processed data from
73	! all example profiles.
74	!
75	! To calculate bending angle, refractivity and dry temperature profiles
76	! from all profiles contained in the multifile multi.nc:
77	!
78	! > ropp_pp_occ_tool ../data/multi.nc -o example_03.nc
79	!
80	! Since the multi_* file was generated by concatenating the other files
81	! in the data directory, example_02.nc and example_03.nc should be identical
82	! apart from the file names.
83	!
84	! AUTHOR
85	! Met Office, Exeter, UK.
86	! Any comments on this software should be given via the ROM SAF
87	! Helpdesk at http://www.romsaf.org
88	!
89	! COPYRIGHT
90	! (c) EUMETSAT. All rights reserved.
91	! For further details please refer to the file COPYRIGHT
92	! which you should have received as part of this distribution.
93	!
94	!****
95
96	!-------------------------------------------------------------------------------
97	! 1. Declarations
98	!-------------------------------------------------------------------------------
99
100	USE typesizes, ONLY: wp => EightByteReal
101	USE ropp_utils
102	USE ropp_io
103	USE ropp_io_types, ONLY: ROprof, L1atype, L1btype, L2atype
104	USE ropp_pp
105	USE ropp_pp_preproc
106	USE ropp_pp_types, ONLY: PPConfig, PPDiag
107
108	IMPLICIT NONE
109
110	TYPE(ROprof) :: ro_data ! Input RO data
111	TYPE(L1atype) :: Lev1a ! Temporary Level1a structure for storage
112	TYPE(L1btype) :: bangle ! Output bending angle profiles
113	TYPE(L1btype) :: cma_bangle ! for CMA data
114	TYPE(L1btype) :: out_ba ! Corrected bending angles
115	TYPE(L1btype) :: smt_ba ! Smoothed bending angle
116	TYPE(L1btype) :: mod_ba ! Model bending angle
117	TYPE(L2atype) :: out_refrac ! Retrieved refractivity profile
118	TYPE(L2btype) :: dum_meteo ! Level2b structure for dummy meteo variables
119
120	TYPE(ppConfig) :: config ! Configuration options
121	TYPE(ppDiag) :: diag ! Diagnostic variables
122
123	INTEGER :: idummy, imin, imax
124	INTEGER :: i, iargc, argc, k, j,itemp
125	INTEGER :: n_files, n_profiles, nbi
126	INTEGER :: ws_go_smooth, ws_go_full, ws_wo, ws_low
127	REAL(wp) :: Pmin, Pmax
128
129	! useful to save these
130
131	REAL(wp),DIMENSION(:), ALLOCATABLE :: snr_l2,phase_l2
132
133	LOGICAL :: give_help
134	LOGICAL :: earth_ellipsoid = .FALSE.
135	LOGICAL :: output_full = .FALSE.
136	LOGICAL :: twofit = .FALSE.
137	LOGICAL :: ranchk = .TRUE.
138
139	CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
140	CHARACTER(len = 4096) :: ofile
141	CHARACTER(len = 4096) :: ofile2 ! just the name of directory
142	CHARACTER(len = 4096) :: cfg_file = " "
143	CHARACTER(len = 256) :: buffer
144	CHARACTER(len = 64) :: outstr
145	CHARACTER(len = 80) :: mfile = " "
146	!2016-1-17
147	!CHARACTER(len = 80) :: bfile = " "
148	CHARACTER(len = 80) :: bfile = " "
149	!2016-1-17
150	CHARACTER(len = 80) :: navfile = " "
151	CHARACTER(len = 10) :: method = " "
152	CHARACTER(len = 10) :: occmethod = " "
153	CHARACTER(len = 10) :: filtmethod = " "
154	CHARACTER(len = 4) :: istr
155	CHARACTER(len = 10) :: nstr
156	CHARACTER(len = 10) :: pstr1
157	CHARACTER(len = 10) :: pstr2
158	CHARACTER(len = 256):: infoprint
159
160	! for SLTA 2017-3-6
161	INTEGER :: iLine,iNUM
162	REAL(wp),DIMENSION(:,:), ALLOCATABLE :: footVertiList
163	REAL(wp),DIMENSION(:), ALLOCATABLE :: P_H
164
165	REAL(wp) :: Pleo(3)
166	REAL(wp) :: Pgps(3)
167	REAL(wp) :: Pver(3)
168
169	! for L2 problems
170
171	REAL(wp) :: min_slta
172
173	! for LC
174	INTEGER :: Neg_num, T_num
175	REAL(wp) :: R
176
177	! for saving file
178	REAL(wp) :: ps1, psN ! Start/end impact parameter
179	INTEGER :: n ! Number of data points
180	INTEGER :: ocd ! Occ direction
181	! REAL(wp) :: min_slta_l1,min_slta_l2,last_slta_l1
182	REAL(wp) :: max_slta_B,min_slta_G
183	! for ave and std
184	REAL(wp) :: noise_estimate
185	REAL(wp) :: sum_phsl1,sum_phsl2,sum_snrl1,sum_snrl2
186	REAL(wp) :: sumd_phsl1,sumd_phsl2,sumd_snrl1,sumd_snrl2
187	REAL(wp) :: ave_phsl1,ave_phsl2,ave_snrl1,ave_snrl2
188	REAL(wp) :: std_phsl1,std_phsl2,std_snrl1,std_snrl2
189	INTEGER :: nl1,nl2
190
191	! for saving txt files
192	INTEGER :: inb,ii
193	CHARACTER(len = 4) :: yr
194	CHARACTER(len = 2) :: mn,dy,hr,mi,se
195	CHARACTER(len = 256) :: fname,Thinfile
196
197	! for leo height
198	REAL(wp),DIMENSION(:), ALLOCATABLE :: P_leo
199	REAL(wp) :: dHleo
200	!
201	!-------------------------------------------------------------------------------
202	! 2. Default settings
203	!-------------------------------------------------------------------------------
204
205	give_help = .FALSE.
206	ofile = "ropp_pp_occ.nc"
207	ofile2 = ""
208
209	CALL message(msg_noin, '')
210	CALL message(msg_noin, &
211	'----------------------------------------------------------------------')
212	CALL message(msg_noin, &
213	' ROPP Occultation Pre-processor Tool')
214	CALL message(msg_noin, &
215	'----------------------------------------------------------------------')
216	CALL message(msg_noin, '')
217
218	!-------------------------------------------------------------------------------
219	! 3. Command line arguments
220	!-------------------------------------------------------------------------------
221
222	argc = iargc()
223	i = 1
224	n_files = 0
225	ALLOCATE(ifiles(argc))
226
227	DO WHILE(i <= argc)
228	CALL getarg(i, buffer)
229	SELECT CASE (buffer)
230	CASE('-o') ! Output file name (netCDF output)
231	CALL getarg(i+1, buffer)
232	ofile = buffer
233	i = i + 1
234	! 20131001 FY3
235	CASE('-out') ! Output file name (netCDF output)
236	CALL getarg(i+1, buffer)
237	ofile2 = buffer
238	i = i + 1
239	CASE('-c') ! Configuration file name
240	CALL getarg(i+1, buffer)
241	cfg_file = buffer
242	i = i + 1
243	CASE('-m') ! Ionospheric correction method
244	CALL getarg(i+1, buffer)
245	method = buffer
246	i = i + 1
247	CASE('-mfile') ! Model refractivity file
248	CALL getarg(i+1, buffer)
249	mfile = buffer
250	i = i + 1
251	CASE('-bfile') ! Background atmosphere profile file
252	CALL getarg(i+1, buffer)
253	bfile = buffer
254	i = i + 1
255	CASE('-navfile') ! External navigation bit file
256	CALL getarg(i+1, buffer)
257	navfile = buffer
258	i = i + 1
259	CASE('-occ', '--occ') ! Occultation processing method
260	CALL getarg(i+1, buffer)
261	occmethod = buffer
262	i = i + 1
263	CASE('-filter', '--filter') ! Filtering method
264	CALL getarg(i+1, buffer)
265	filtmethod = buffer
266	i = i + 1
267	CASE('-fit', '--fit') ! 2 parameter fitting method
268	twofit = .true.
269	CASE('-ellipsoid') ! Output height wrt reference ellipsoid
270	earth_ellipsoid = .TRUE.
271	CASE('-full') ! Output full height grid (IC plus MSIS)
272	output_full = .TRUE.
273	CASE('-d') ! 'Diagnostic' output mode
274	msg_MODE = VerboseMode
275	CASE ('-no-ranchk') ! Use no rangecheck on output
276	ranchk = .FALSE.
277	CASE('-h', '--help', '?') ! Give some help
278	give_help = .TRUE.
279	CASE('-v', '-V', '--version') ! Give some version information
280	CALL version_info()
281	CALL EXIT(msg_exit_ok)
282	CASE default ! Input file name
283	n_files = n_files + 1
284	ifiles(n_files) = buffer
285	END SELECT
286	i = i + 1
287	END DO
288
289	IF ( n_files == 0 .AND. .NOT. give_help ) THEN
290	CALL message ( msg_error, "No input file(s) specified" )
291	END IF
292
293	IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN
294	CALL usage()
295	CALL EXIT(msg_exit_status)
296	ENDIF
297
298	! 3.1 Read configuration file (if exists), preserving command-line options
299
300	IF (method /= " ") config%method = method
301	IF (mfile /= " ") config%mfile = mfile
302	IF (bfile /= " ") config%bfile = bfile
303	IF (navfile /= " ") config%navbit_file = navfile
304	IF (occmethod /= " ") config%occ_method = occmethod
305	IF (filtmethod /= " ") config%filter_method = filtmethod
306
307	IF (cfg_file /= " ") THEN
308	CALL message(msg_info, &
309	"Reading configuration file " // TRIM(cfg_file) // ".\n")
310
311	CALL ropp_pp_read_config(cfg_file, config)
312	IF (method /= " ") config%method = method
313	IF (mfile /= " ") config%mfile = mfile
314	IF (bfile /= " ") config%bfile = bfile
315	IF (navfile /= " ") config%navbit_file = navfile
316	IF (occmethod /= " ") config%occ_method = occmethod
317	IF (filtmethod /= " ") config%filter_method = filtmethod
318	IF (twofit) config%nparm_fit = 2
319	ENDIF
320
321	!-------------------------------------------------------------------------------
322	! 4. Remove pre-existing output file
323	!-------------------------------------------------------------------------------
324
325	CALL file_delete(ofile, idummy)
326
327	!-------------------------------------------------------------------------------
328	! 5. Loop over all input files and profiles
329	!-------------------------------------------------------------------------------
330
331	DO k = 1, n_files
332
333	n_profiles = ropp_io_nrec(ifiles(k))
334
335	DO i = 1, n_profiles
336
337	WRITE(istr, '(i4)') i
338	WRITE(nstr, '(i6)') n_profiles
339	CALL message(msg_info, "Processing profile " // istr // " of " // nstr )
340
341	!-------------------------------------------------------------------------------
342	! 6. Read data
343	!-------------------------------------------------------------------------------
344
345	! 6.1 Read input file
346
347	CALL message(msg_info, &
348	"Reading input data file " // TRIM(ifiles(k)) // ".\n")
349	CALL ropp_io_read(ro_data, ifiles(k), rec=i, ranchk=.TRUE.)
350	CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")
351
352	IF (ro_data%Lev1a%Npoints == 0) THEN
353	config%obs_ok = .FALSE.
354	CALL ropp_io_free(ro_data%vlist)
355	! 2016-08-16 if L1a excess phase got NaN values, do not produce the .nc file
356	! CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk )
357	CALL message(msg_fatal, "FAILURE: No Level1a data in file " // &
358	ifiles(k) // "\n")
359	ENDIF
360
361	! 6.2 Shrink ro_data to correct size (for multiple profiles)
362
363	CALL ropp_io_init(Lev1a, ro_data%lev1a%npoints)
364	Lev1a = ro_data%lev1a
365	ro_data%lev1a = Lev1a
366	CALL ropp_io_free(Lev1a)
367	CALL ropp_io_free(ro_data%Lev1b)
368
369	!-------------------------------------------------------------------------------
370	! 7. Check coordinate frames
371	!-------------------------------------------------------------------------------
372
373	CALL ropp_pp_set_coordinates(ro_data)
374
375	!-------------------------------------------------------------------------------
376	! 8. Determine occultation point georeference information
377	!-------------------------------------------------------------------------------
378
379	CALL occ_point( ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
380	ro_data%georef%lat, ro_data%georef%lon, &
381	ro_data%georef%r_coc, ro_data%georef%roc, &
382	ro_data%georef%azimuth, &
383	ro_data%georef%undulation, &
384	config%egm96, config%corr_egm96)
385
386	WRITE(outstr,'(A,1x,i2.2,a1,i2.2,a1,i4.4)') &
387	"Occultation date (dd/mm/yyyy): = ", &
388	ro_data%dtocc%day, "/", &
389	ro_data%dtocc%month, "/", &
390	ro_data%dtocc%year
391	CALL message(msg_diag, outstr)
392
393	WRITE(outstr,'(A,1x,i2.2,a1,i2.2,a1,i2.2)') &
394	"Occultation time (hh:mm:ss): = ", &
395	ro_data%dtocc%hour, ":", &
396	ro_data%dtocc%minute, ":", &
397	ro_data%dtocc%second
398	CALL message(msg_diag, outstr)
399
400	WRITE(outstr,'(A,1x,a1,f7.2,a2,f7.2,a2)') &
401	"Occultation point (lat,lon): = ", &
402	"(", ro_data%georef%lat, "N,", &
403	ro_data%georef%lon, "E)"
404	CALL message(msg_diag, outstr)
405
406	WRITE(outstr,'(A,1x,e15.5)') &
407	"Undulation (m) = ", ro_data%georef%undulation
408	CALL message(msg_diag, outstr)
409
410	!-------------------------------------------------------------------------------
411	! xxx. science code for GNOS QC ---- FEB.to APR.2017
412	!-------------------------------------------------------------------------------
413	! CALCULATE THE straight LINE HEIGHT
414	iNUM = size(ro_data%lev1a%dtime)
415
416	ALLOCATE(P_H(iNUM))
417	ALLOCATE(P_leo(iNUM))
418
419	max_slta_B = -1.0E04
420
421	min_slta_G = 1.0E20
422
423	DO iLine=1,iNUM
424
425	Pleo(1) = ro_data%lev1a%r_leo(iLine,1) - ro_data%georef%r_coc(1)
426	Pleo(2) = ro_data%lev1a%r_leo(iLine,2) - ro_data%georef%r_coc(2)
427	Pleo(3) = ro_data%lev1a%r_leo(iLine,3) - ro_data%georef%r_coc(3) !LEO XYZ
428	Pgps(1) = ro_data%lev1a%r_gns(iLine,1) - ro_data%georef%r_coc(1)
429	Pgps(2) = ro_data%lev1a%r_gns(iLine,2) - ro_data%georef%r_coc(2)
430	Pgps(3) = ro_data%lev1a%r_gns(iLine,3) - ro_data%georef%r_coc(3) !GPS XYZ
431
432	CALL vertical_point(Pleo,Pgps,Pver) !
433
434	P_H(iLine) = sqrt(Pver(1)2 + Pver(2)2 + Pver(3)**2)/1000.0_wp - ro_data%georef%roc/1000.0_wp
435	P_leo(iLine) = sqrt(Pleo(1)2 + Pleo(2)2 + Pleo(3)**2)/1000.0_wp - ro_data%georef%roc/1000.0_wp
436
437	! find the last min_slta of L2
438
439	IF (ABS(ro_data%lev1a%snr_L2P(iLine) - ropp_MDFV) < 1.0_wp .OR. & ! amplitude missing ?
440	ABS(ro_data%lev1a%phase_L2(iLine) - ropp_MDFV) < 1.0_wp .OR. & ! phase missing ?
441	ABS(ro_data%lev1a%phase_L2(iLine) + 999.0_wp ) < 0.001_wp) THEN ! CMA missing data ?
442
443	max_slta_B = MAX(max_slta_B,P_H(iLine ))
444
445	ELSE
446
447	min_slta_G = MIN(min_slta_G,P_H(iLine ))
448
449
450	ENDIF
451	PRINT*,'Height xxx', iLine, ro_data%lev1a%dtime(iLine),P_H(iLine),&
452	ro_data%lev1a%snr_L2P(iLine),ro_data%lev1a%snr_L1ca(iLine),&
453	ro_data%lev1a%phase_L2(iLine),ro_data%lev1a%phase_L1(iLine),&
454	P_leo(iLine)
455
456
457	ENDDO
458
459	! determine the height of leo
460	dHleo = maxval(P_leo) - minval(P_leo)
461	print*,'dhleo',dHleo
462	if(dHleo > 10.0_wp)config%obs_ok = .FALSE.
463
464	! save in m. Added 6 km to the min SLTA.
465	PRINT*,'MIN_SLTA_G',min_slta_G,max_slta_B
466	config%hmax_wo = 20000.0_wp !18000.0_wp !22000.0_wp !25000.0_wp !20000.0_wp
467	min_slta = MAX(1.0E3_wp(min_slta_G+6.0_wp),config%hmax_wo + 0.5_wpconfig%fw_go_full)
468
469
470	!
471	! calculate the mean and std SNR and phase at the range of 60 - 80 km
472	!
473	DO iLine=1,iNUM
474
475	IF ( P_H(iLine) .gt. 60.0_wp .and. P_H(iLine) .lt.80.0_wp ) THEN
476
477	IF (ro_data%lev1a%phase_L1(iLine) .gt. -99999.0_wp .and.&
478	ro_data%lev1a%snr_L1ca(iLine) .gt. 0.0_wp ) THEN
479	sum_phsl1 = sum_phsl1 + ro_data%lev1a%phase_L1(iLine)
480	sum_snrl1 = sum_snrl1 + ro_data%lev1a%snr_L1ca(iLine)
481	sumd_phsl1 = sumd_phsl1 + ro_data%lev1a%phase_L1(iLine)**2
482	sumd_snrl1 = sumd_snrl1 + ro_data%lev1a%snr_L1ca(iLine)**2
483	nl1 = nl1 + 1
484	ENDIF
485	IF (ro_data%lev1a%phase_L2(iLine) .gt. -99999.0_wp .and.&
486	ro_data%lev1a%snr_L2P(iLine) .gt. 0.0_wp ) THEN
487	sum_phsl2 = sum_phsl2 + ro_data%lev1a%phase_L2(iLine)
488	sum_snrl2 = sum_snrl2 + ro_data%lev1a%snr_L2P(iLine)
489	sumd_phsl2 = sumd_phsl2 + ro_data%lev1a%phase_L2(iLine)**2
490	sumd_snrl2 = sumd_snrl2 + ro_data%lev1a%snr_L2P(iLine)**2
491	nl2 = nl2 + 1
492	ENDIF
493
494	ENDIF
495
496	ENDDO
497
498	DEALLOCATE(P_H)
499	DEALLOCATE(P_leo)
500
501	IF (nl1 > 1) THEN
502	ave_phsl1 = sum_phsl1/nl1
503	ave_snrl1 = sum_snrl1/nl1
504	std_phsl1 = sqrt((sumd_phsl1 - sum_phsl1*sum_phsl1/nl1) /(nl1-1))
505	std_snrl1 = sqrt((sumd_snrl1 - sum_snrl1*sum_snrl1/nl1) /(nl1-1))
506	ELSE
507	ave_phsl1 = -10.0_wp
508	ave_snrl1 = -10.0_wp
509	std_phsl1 = -10.0_wp
510	std_snrl1 = -10.0_wp
511	ENDIF
512
513	IF (nl2 > 1) THEN
514	ave_phsl2 = sum_phsl2/nl2
515	ave_snrl2 = sum_snrl2/nl2
516	std_phsl2 = sqrt((sumd_phsl2 - sum_phsl2*sum_phsl2/nl2) /(nl2-1))
517	std_snrl2 = sqrt((sumd_snrl2 - sum_snrl2*sum_snrl2/nl2) /(nl2-1))
518	ELSE
519	ave_phsl2 = -10.0_wp
520	ave_snrl2 = -10.0_wp
521	std_phsl2 = -10.0_wp
522	std_snrl2 = -10.0_wp
523	ENDIF
524
525	! Determine the occ direction
526
527	n = ro_data%lev1a%npoints
528	ps1 = impact_parameter(ro_data%lev1a%r_leo(1,:) - ro_data%georef%r_coc(:), &
529	ro_data%lev1a%r_gns(1,:) - ro_data%georef%r_coc)
530	psN = impact_parameter(ro_data%lev1a%r_leo(n,:) - ro_data%georef%r_coc(:), &
531	ro_data%lev1a%r_gns(n,:) - ro_data%georef%r_coc)
532
533	ocd = NINT(SIGN(1.0_wp, psN - ps1))
534
535
536	!-------------------------------------------------------------------------------
537	! 9. Mission-specific pre-processing and input data cut-off
538	!-------------------------------------------------------------------------------
539
540	!
541	! save the L2 excess phase and SNR before messed up in preprocess
542	!
543	ALLOCATE(snr_l2(ro_data%lev1a%npoints))
544	ALLOCATE(phase_l2(ro_data%lev1a%npoints))
545
546	phase_l2 = ro_data%Lev1a%phase_L2
547	snr_l2 = ro_data%Lev1a%snr_L2p
548
549	IF (config%obs_ok) THEN
550	CALL ropp_pp_preprocess(ro_data, config, diag)
551	ENDIF
552
553
554	!-------------------------------------------------------------------------------
555	! 10. Initialise bending angle structures
556	!-------------------------------------------------------------------------------
557
558	IF (config%obs_ok) THEN
559
560	CALL ropp_io_init(smt_ba, ro_data%Lev1a%Npoints)
561	CALL ropp_io_init(bangle, ro_data%Lev1a%Npoints)
562
563	CALL ropp_io_init(cma_bangle, ro_data%Lev1a%Npoints)
564
565	CALL tangent_point(ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, &
566	bangle%lat_tp, bangle%lon_tp, bangle%azimuth_tp)
567
568	!-------------------------------------------------------------------------------
569	! 11. Geometric optics processing
570	!-------------------------------------------------------------------------------
571
572	CALL message(msg_info, &
573	"Retrieving bending angle profiles by GEOMETRIC OPTICS \n")
574
575
576	! 11.1 Calculate smoothed bending angle
577
578	DO j=1,ro_data%Lev1a%Npoints
579	smt_ba%impact(j) = &
580	impact_parameter(ro_data%Lev1a%r_leo(j,:)-ro_data%georef%r_coc, &
581	ro_data%Lev1a%r_gns(j,:)-ro_data%georef%r_coc )
582	ENDDO
583	Pmax = MAXVAL(smt_ba%impact)
584	Pmin = MAX(MINVAL(smt_ba%impact), ro_data%georef%roc+2000.0_wp)
585	ws_go_smooth = CEILING( config%fw_go_smooth * &
586	(ro_data%Lev1a%Npoints-1)/ABS(Pmax-Pmin))
587
588	WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc
589	WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc
590	CALL message(msg_diag, "Smoothed bending angle. ")
591	CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
592	WRITE(nstr, '(i8)') ws_go_smooth
593	CALL message(msg_diag, "ws_go_smooth " // nstr)
594
595	IF (Pmax > Pmin) THEN
596
597	CALL ropp_pp_bending_angle_go( ro_data%Lev1a%dtime, &
598	ro_data%Lev1a%r_leo, &
599	ro_data%Lev1a%r_gns, &
600	ro_data%georef%r_coc, &
601	ro_data%Lev1a%phase_L1, &
602	ro_data%Lev1a%phase_L2, &
603	ws_go_smooth, &
604	config%filter_method, &
605	smt_ba%impact_L1, &
606	smt_ba%bangle_L1, &
607	smt_ba%impact_L2, &
608	smt_ba%bangle_L2 )
609
610
611	CALL ropp_pp_linear_combination( smt_ba%impact_L1, &
612	smt_ba%bangle_L1, &
613	smt_ba%impact_L2, &
614	smt_ba%bangle_L2, &
615	smt_ba%impact, &
616	smt_ba%bangle )
617
618	ELSE
619
620	config%obs_ok = .FALSE.
621	CALL message(msg_warn, "Cannot process profile (smt): Pmin > Pmax")
622
623	ENDIF
624
625	! 11.2 Calculate bending angle
626
627	bangle%impact = smt_ba%impact_L1
628	Pmax = MAXVAL(smt_ba%impact_L1)
629	Pmin = MAX(MINVAL(smt_ba%impact_L1), ro_data%georef%roc+2000.0_wp)
630	ws_go_full = CEILING(config%fw_go_full*(ro_data%Lev1a%Npoints - 1)/ &
631	ABS(Pmax - Pmin))
632
633	WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc
634	WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc
635
636	CALL message(msg_diag, "Full resolution bending angle. ")
637	CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
638	WRITE(nstr, '(i8)') ws_go_full
639	CALL message(msg_diag, "ws_go_full " // nstr)
640
641
642
643
644	IF (Pmax > Pmin) THEN
645
646	CALL ropp_pp_bending_angle_go(ro_data%Lev1a%dtime, &
647	ro_data%Lev1a%r_leo, &
648	ro_data%Lev1a%r_gns, &
649	ro_data%georef%r_coc, &
650	ro_data%Lev1a%phase_L1, &
651	ro_data%Lev1a%phase_L2, &
652	ws_go_full, &
653	config%filter_method, &
654	bangle%impact_L1, &
655	bangle%bangle_L1, &
656	bangle%impact_L2, &
657	bangle%bangle_L2 )
658	! new
659
660
661	cma_bangle%impact_L1 = bangle%impact_L1
662	cma_bangle%impact_L2 = bangle%impact_L2
663
664	cma_bangle%bangle_L1 = bangle%bangle_L1
665	cma_bangle%bangle_L2 = bangle%bangle_L2
666
667
668	WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints
669	WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp
670	WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp
671	CALL message(msg_info, &
672	TRIM(nstr) // " data points in output between " // &
673	TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n")
674
675	ELSE
676
677	config%obs_ok = .FALSE.
678	CALL message(msg_warn, "Cannot process profile (full): Pmin > Pmax")
679
680	ENDIF
681
682	!-------------------------------------------------------------------------------
683	! 12. Wave optics processing
684	!-------------------------------------------------------------------------------
685
686	IF ( INDEX(config%occ_method, "WO" ) == 1 ) THEN
687
688	CALL message(msg_info, &
689	"Retrieving bending angle profiles by WAVE OPTICS \n")
690
691	ws_wo = CEILING(config%fw_wo * (bangle%npoints-1.0_wp) / &
692	ABS(Pmax - Pmin))
693	ws_low = CEILING(config%fw_low * (bangle%npoints-1.0_wp) / &
694	ABS(Pmax - Pmin))
695
696	WRITE(nstr,'(I5)') ws_go_smooth
697	CALL message(msg_diag, 'WM = ' // nstr)
698	WRITE(nstr,'(I5)') ws_wo
699	CALL message(msg_diag, 'W = ' // nstr)
700	WRITE(nstr,'(I5)') ws_low
701	CALL message(msg_diag, 'WL = ' // nstr)
702
703	IF (Pmax > Pmin) THEN
704
705	CALL ropp_pp_bending_angle_wo( ro_data%Lev1a%dtime, &
706	ro_data%Lev1a%r_leo, &
707	ro_data%Lev1a%r_gns, &
708	ro_data%georef%r_coc, &
709	ro_data%georef%roc, &
710	ro_data%Lev1a%phase_L1, &
711	ro_data%Lev1a%phase_L2, &
712	ro_data%Lev1a%snr_L1ca, &
713	ro_data%Lev1a%snr_L2p, &
714	ws_go_smooth, &
715	ws_wo, &
716	ws_low, &
717	config%hmax_wo, &
718	config%filter_method, &
719	config%opt_DL2, &
720	config%cff, &
721	config%dsh, &
722	bangle%impact_L1, &
723	bangle%bangle_L1, &
724	bangle%bangle_L1_sigma, &
725	bangle%impact_L2, &
726	bangle%bangle_L2, &
727	bangle%bangle_L2_sigma, &
728	diag )
729
730	WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints
731	WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp
732	WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp
733	CALL message(msg_info, &
734	TRIM(nstr) // " data points in output between " // &
735	TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n")
736
737	!
738	! overwrite the L2
739	!
740
741
742	call cma_merge(ro_data%lev1a%npoints, &
743	min_slta, &
744	ro_data%georef%roc, &
745	cma_bangle,bangle,noise_estimate)
746
747
748
749	ELSE
750
751	config%obs_ok = .FALSE.
752	CALL message(msg_warn, "Cannot process profile (wo): Pmin > Pmax")
753
754	ENDIF
755
756
757	ENDIF
758
759	!-------------------------------------------------------------------------------
760	! 13. Check that profiles are monotonously increasing in height
761	! (index 1 towards surface) and valid bending angles
762	!-------------------------------------------------------------------------------
763
764	CALL ropp_pp_monotonous(bangle%impact_L1, -1)
765	CALL ropp_pp_monotonous(bangle%impact_L2, -1)
766	CALL ropp_pp_monotonous(smt_ba%impact, -1)
767
768	IF ( ANY(bangle%bangle_L1 <= ropp_MDTV) .OR. &
769	ANY(bangle%bangle_L2 <= ropp_MDTV) ) THEN
770	config%obs_ok = .FALSE.
771	CALL message(msg_warn, "Cannot process profile: " // &
772	"Invalid bending angles computed \n")
773	ENDIF
774
775	ENDIF
776
777	!-------------------------------------------------------------------------------
778	! 14. Compute un-optimised ionospheric corrected bending angle profile (LC)
779	!-------------------------------------------------------------------------------
780
781	IF (config%obs_ok) THEN
782
783	CALL ropp_pp_linear_combination( bangle%impact_L1, bangle%bangle_L1, &
784	bangle%impact_L2, bangle%bangle_L2, &
785	bangle%impact, bangle%bangle )
786
787
788
789	!-------------------------------------------------------------------------------
790	! 15. Define configuration parameters
791	!-------------------------------------------------------------------------------
792
793	config%r_curve = ro_data%georef%roc
794	config%npoints = bangle%npoints
795
796	config%Pmax = MAXVAL(bangle%impact_L1(:))
797	config%Pmin = MINVAL(bangle%impact_L1(:))
798
799	!-------------------------------------------------------------------------------
800	! 16. Interpolate input data to standard grid
801	!-------------------------------------------------------------------------------
802
803	! 16.1 Calculate size of standard grid
804
805	IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
806	Pmax = MAXVAL(bangle%impact_L1(:)) - 5000.0_wp
807	ELSE
808	Pmax = config%ztop_invert + config%r_curve
809	ENDIF
810	Pmin = MINVAL(bangle%impact_L1(:))
811	nbi = 1 + CEILING((Pmax - Pmin)/config%dpi)
812
813	WRITE(pstr2, '(f10.3)') Pmax-config%r_curve
814	WRITE(pstr1, '(f10.3)') Pmin-config%r_curve
815	CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
816	WRITE(nstr, '(I6)') nbi
817	CALL message(msg_diag, "Standard grid size nbi = " // nstr)
818
819	! 16.2 Initialise standard grid data structures for output
820
821	CALL ropp_io_init(out_ba, nbi)
822	CALL ropp_io_init(out_refrac, nbi)
823
824	! 16.3 Interpolate input data onto standard grid
825
826	CALL ropp_pp_merge_profile(bangle%impact_L1, bangle%bangle_L1, &
827	bangle%impact_L2, bangle%bangle_L2, &
828	out_ba%impact_L1, out_ba%bangle_L1, &
829	out_ba%impact_L2, out_ba%bangle_L2, &
830	Pmin, Pmax)
831
832	CALL ropp_pp_merge_profile(bangle%impact_L1,bangle%bangle_L1_sigma, &
833	bangle%impact_L2,bangle%bangle_L2_sigma, &
834	out_ba%impact_L1,out_ba%bangle_L1_sigma, &
835	out_ba%impact_L2,out_ba%bangle_L2_sigma, &
836	Pmin, Pmax)
837
838	!-------------------------------------------------------------------------------
839	! 17. Ionospheric correction of bending angle profile by linear combination
840	!-------------------------------------------------------------------------------
841
842	IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
843
844	CALL message(msg_info, &
845	"Correcting bending angle profile by LINEAR COMBINATION \n")
846	ro_data%bangle_method = TRIM(config%occ_method) // "(LC)"
847
848	CALL ropp_pp_linear_combination( out_ba%impact_L1, &
849	out_ba%bangle_L1, &
850	out_ba%impact_L2, &
851	out_ba%bangle_L2, &
852	out_ba%impact_opt, &
853	out_ba%bangle_opt )
854
855	imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0))))
856	imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin)))
857
858
859	!-------------------------------------------------------------------------------
860	! 18. Ionospheric correct bending angle profile by statistical optimization
861	!-------------------------------------------------------------------------------
862
863	ELSE
864
865	CALL message(msg_info, &
866	"Correcting bending angle profile by STATISTICAL OPTIMISATION \n")
867	ro_data%bangle_method = TRIM(config%occ_method) // " (StatOpt)"
868
869	! 18.1 Retrieve model bending angles on L1 impact parameter levels
870
871	CALL ropp_io_init(mod_ba, nbi)
872
873	IF (INDEX(config%method, "GMSIS" ) == 1 ) THEN
874	CALL ropp_pp_search_model_refraction( config%mfile, &
875	ro_data%dtocc%month, &
876	ro_data%georef%lat, &
877	ro_data%georef%lon, &
878	smt_ba%impact, &
879	smt_ba%bangle, &
880	out_ba%impact_L1, &
881	mod_ba%bangle, &
882	config )
883
884	ELSE IF (INDEX(config%method, "MSIS" ) == 1 ) THEN
885	CALL ropp_pp_model_refraction( config%mfile, &
886	ro_data%dtocc%month, &
887	ro_data%georef%lat, &
888	ro_data%georef%lon, &
889	out_ba%impact_L1, &
890	mod_ba%bangle, &
891	config )
892
893
894	ELSE IF (INDEX(config%method, "BG" ) == 1 ) THEN
895
896	CALL ropp_pp_bg_refraction( config%bfile, &
897	ro_data%dtocc%month, &
898	ro_data%georef%lat, &
899	ro_data%georef%lon, &
900	out_ba%impact_L1, &
901	mod_ba%bangle, &
902	config )
903
904	ELSE
905	config%obs_ok = .FALSE.
906	CALL message(msg_warn, &
907	"Statistical optimisation method " // config%method // &
908	" not supported")
909	EXIT
910	ENDIF
911
912	! 18.2 Fit model bending angle with (smoothed) observed bending angles
913
914	CALL ropp_pp_fit_model_refraction( smt_ba%impact, &
915	smt_ba%bangle, &
916	out_ba%impact_L1, &
917	mod_ba%bangle, &
918	config )
919
920	! 18.3 Perform ionospheric correction with statistical optimization
921
922	imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0))))
923	imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin)))
924
925	! Merge observed bending angles with model above observation top
926	WHERE ( out_ba%impact_L1 > config%Pmax - 5000.0_wp )
927	out_ba%bangle_L1(:) = mod_ba%bangle(:)
928	out_ba%bangle_L2(:) = mod_ba%bangle(:)
929	out_ba%bangle_L1_sigma(:) = 0.0_wp
930	out_ba%bangle_L2_sigma(:) = 0.0_wp
931	END WHERE
932
933	IF (imin < imax) THEN
934	out_ba%impact_opt(:) = out_ba%impact_L1(:)
935	out_ba%bangle_opt(:) = out_ba%bangle_L1(:)
936	CALL ropp_pp_ionospheric_correction( &
937	out_ba%impact_L1(imin:imax), & ! L1 impact
938	out_ba%bangle_L1(imin:imax), & ! L1 bending
939	out_ba%impact_L2(imin:imax), & ! L2 impact
940	out_ba%bangle_L2(imin:imax), & ! L2 bending
941	out_ba%impact_L1(imin:imax), & ! Model impact
942	mod_ba%bangle(imin:imax), & ! Model bending
943	config, & ! Configuration
944	out_ba%impact_opt(imin:imax), & ! Opt impact
945	out_ba%bangle_opt(imin:imax), & ! Opt bending
946	diag ) ! Diagnostics
947	diag%err_neut = diag%err_neut + &
948	out_ba%bangle_L1_sigma(imin:imax)**2
949	diag%sq = 100.0_wp * &
950	MAXVAL(SQRT(diag%err_neut(:))/out_ba%bangle_opt(imin:imax))
951	ELSE
952	config%obs_ok = .FALSE.
953	CALL message(msg_warn, &
954	"Cannot perform ionospheric correction: Imin >= Imax \n")
955	ENDIF
956
957	ENDIF
958
959	! 18.4 Interpolate LC bending angles, lat_tp, lon_tp, azimuth_tp to output grid
960
961	out_ba%impact = out_ba%impact_opt
962	CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
963	bangle%bangle, out_ba%bangle)
964	CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
965	bangle%lat_tp, out_ba%lat_tp)
966	CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
967	bangle%lon_tp, out_ba%lon_tp)
968	CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
969	bangle%azimuth_tp, out_ba%azimuth_tp)
970
971	WRITE(nstr, '(i10)') imax-imin+1
972	WRITE(pstr2, '(f6.1)') out_ba%impact(imax)/1000.0_wp
973	WRITE(pstr1, '(f6.1)') out_ba%impact(imin)/1000.0_wp
974	CALL message(msg_info, &
975	TRIM(nstr) // " data points in output between " // &
976	TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n")
977
978	END IF
979
980	!-------------------------------------------------------------------------------
981	! 19. Perform inverse Abel transform to compute refractivity
982	!-------------------------------------------------------------------------------
983
984	IF (config%obs_ok) THEN
985
986	ro_data%refrac_method = "ABEL (" // TRIM(config%abel) // ")"
987
988	! 19.1 Abel inversion of corrected bending angle profile
989
990	IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
991
992	CALL message(msg_info, &
993	"Retrieving refractivity profile by " // TRIM(config%abel) // &
994	" ABEL TRANSFORM \n")
995
996	IF ( INDEX(config%abel, "LIN" ) == 1 ) THEN
997	CALL ropp_pp_invert_LIN(out_ba%impact_opt, out_ba%bangle_opt, &
998	out_ba%impact_opt, out_refrac%refrac)
999	ELSE IF ( INDEX(config%abel, "EXP" ) == 1 ) THEN
1000	CALL ropp_pp_invert_EXP(out_ba%impact_opt, out_ba%bangle_opt, &
1001	out_ba%impact_opt, out_refrac%refrac)
1002	ELSE
1003	config%obs_ok = .FALSE.
1004	CALL message(msg_warn, &
1005	"Abel integral method " // config%abel // " not supported")
1006	ENDIF
1007
1008	! 19.2 Abel inversion of corrected bending angle profile with stat opt
1009
1010	ELSE
1011
1012	CALL message(msg_info, &
1013	"Retrieving refractivity profile by " // TRIM(config%abel) // &
1014	" ABEL TRANSFORM with STAT OPT \n")
1015
1016	CALL ropp_pp_invert_refraction( config%mfile, &
1017	ro_data%dtocc%month, &
1018	ro_data%georef%lat, &
1019	ro_data%georef%lon, &
1020	out_ba%impact_opt, &
1021	out_ba%bangle_opt, &
1022	out_refrac%geop_refrac, &
1023	out_refrac%refrac, &
1024	config )
1025
1026	ENDIF
1027
1028	!-------------------------------------------------------------------------------
1029	! 20. Compute output height scales
1030	!-------------------------------------------------------------------------------
1031
1032	! Override default output height scales wrt geoid if requested
1033	IF (earth_ellipsoid) ro_data%georef%undulation = ropp_MDFV
1034
1035	IF (ro_data%georef%undulation > ropp_MDTV) THEN
1036	CALL message(msg_info, "Writing output altitude scales " // &
1037	"with respect to EGM96 GEOID")
1038	out_refrac%alt_refrac = &
1039	((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) &
1040	- (ro_data%GEOref%roc + ro_data%GEOref%undulation))
1041	ELSE
1042	IF (.not. earth_ellipsoid) &
1043	CALL message(msg_warn, "Invalid undulation calculated. " // &
1044	"Check for valid EGM geoid coefficient and correction file.")
1045	CALL message(msg_info, "Writing output altitude scales " // &
1046	"with respect to WGS84 ELLIPSOID.")
1047	out_refrac%alt_refrac = &
1048	((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) &
1049	- ro_data%GEOref%roc)
1050	ENDIF
1051	out_refrac%geop_refrac = &
1052	geometric2geopotential(ro_data%georef%lat, out_refrac%alt_refrac)
1053
1054	!-------------------------------------------------------------------------------
1055	! 21. Compute Tdry and Pdry - though only Tdry is written to ro_data
1056	!-------------------------------------------------------------------------------
1057
1058	IF (config%output_tdry) THEN
1059	CALL message(msg_info, "Computing dry temperature \n")
1060	CALL ropp_io_init(dum_meteo, nbi) ! easy way to get dummy storage
1061	dum_meteo%shum = 0.0_wp
1062	CALL ropp_pp_tdry(ro_data%georef%lat, out_refrac%alt_refrac, &
1063	out_refrac%refrac, dum_meteo%shum, out_refrac%dry_temp, &
1064	dum_meteo%press)
1065	CALL ropp_io_free(dum_meteo)
1066	ENDIF
1067
1068	!-------------------------------------------------------------------------------
1069	! 22. Copy retrieved profiles to RO structure
1070	! - only output data within observed range
1071	!-------------------------------------------------------------------------------
1072
1073	! Override default output only within observed range if requested
1074	IF (.not. output_full) THEN
1075	out_ba%npoints = Imax - Imin + 1
1076	out_refrac%npoints = Imax - Imin + 1
1077	ENDIF
1078
1079	CALL ropp_io_roprof2roprof(out_ba, ro_data%lev1b)
1080	CALL ropp_io_roprof2roprof(out_refrac, ro_data%lev2a)
1081
1082	ENDIF
1083
1084	!-------------------------------------------------------------------------------
1085	! 23. Write data
1086	!-------------------------------------------------------------------------------
1087	IF (config%obs_ok) THEN
1088	CALL ropp_io_free(ro_data%vlist) !! interim, avoid writing lcf data
1089	IF (config%output_diag) THEN
1090	CALL message(msg_info, "Writing additional diagnostic output \n")
1091	CALL ropp_pp_diag2roprof(diag, ro_data)
1092	ENDIF
1093
1094	CALL message(msg_info, "Writing to output file " // TRIM(ofile) // "\n")
1095
1096	!!! new for FY3 GNOS (below)
1097	CALL message(msg_info, "Writing to output FY3 LEV2 file path " // TRIM(ofile) // "\n")
1098
1099	!20141216 config%HL2_limit (must add a new HL2_limit in config)
1100	! if the tangent height HL2_limit(in default is 20km from config file) then QC= 50;
1101	IF ((ro_data%exL2Type == -999) .OR. (ro_data%lowestTphL2P == ropp_MDFV ) .OR. (ro_data%lowestTphL2C == ropp_MDFV )) THEN
1102	ro_data%QC = -999
1103	ELSE
1104	ro_data%QC = 100
1105	ENDIF
1106	! if refractivity gets negative value, then QC=60
1107	DO itemp = 1,ro_data%Lev2a%Npoints
1108	IF (ro_data%Lev2a%refrac(itemp) < 0.0 ) THEN
1109	ro_data%QC = 60
1110	!EXIT
1111	ENDIF
1112	END DO
1113	IF (ro_data%exL2Type == 0) THEN !L2 is p code
1114	IF (ro_data%lowestTphL2P > config%HL2_limit) THEN
1115	ro_data%QC = 50
1116	ELSE IF (ro_data%lowestTphL2P == -999) THEN ! L2 gets no signal at all
1117	ro_data%QC = 0
1118	END IF
1119	ELSE
1120	IF (ro_data%lowestTphL2C > config%HL2_limit) THEN
1121	ro_data%QC = 50
1122	ELSE IF (ro_data%lowestTphL2C == -999) THEN
1123	ro_data%QC = 0
1124	END IF
1125	ENDIF
1126	!!!!! new for GNOS (above)
1127
1128	! Save snr info, noise_estimate , slta and L2_badness score into one file
1129	noise_estimate=noise_estimate*1E06
1130	OPEN(99,file = '/home/rd/dam1/mi_liao/gnosrun/ropp-8.0/gnossh/atECM/multi/SNR_INFO_test.txt',access = 'append')
1131
1132	IF (noise_estimate > 20.0_wp .or. (ocd == 1 .and. abs(ave_phsl1) <150.0_wp .and. abs(ave_phsl2) <150.0_wp ))THEN
1133	WRITE(99,'(5i5,2f12.3,i8,11f12.3)') ro_data%DTocc%year,ro_data%DTocc%month,ro_data%DTocc%day,ro_data%DTocc%hour,&
1134	ro_data%DTocc%minute,ro_data%GEOref%lat,ro_data%GEOref%lon,ocd,&
1135	ave_snrl1,std_snrl1,ave_snrl2,std_snrl2,ave_phsl1,std_phsl1,ave_phsl2,std_phsl2,&
1136	noise_estimate,min_slta_G,diag%L2_badness
1137
1138	CLOSE(99)
1139	ELSE
1140	CALL ropp_io_write(ro_data, ofile, ofile2,append=.TRUE., ranchk=ranchk, output_precision='double' ) !add output_precision='double' 20160831
1141	ENDIF
1142
1143	ENDIF
1144	!
1145
1146
1147	! ! thin into 247 levels and write into a txt file---for the purpose of making the input files to BA fm
1148	!
1149	! ThinFile = '/home/rd/dam1/mi_liao/gnosrun/ropp-8.0/ropp_io/data/ropp_thin_eum-247.dat'
1150	!
1151	! CALL ropp_io_thin(ro_data, trim(ThinFile))
1152	!
1153	! inb = 247
1154	! print*, inb
1155	! ! save
1156	!
1157	! write(yr,'(i4)')ro_data%dtocc%year
1158	! write(mn,'(i2.2)')ro_data%dtocc%month
1159	! write(dy,'(i2.2)')ro_data%dtocc%day
1160	! write(hr,'(i2.2)')ro_data%dtocc%hour
1161	! write(mi,'(i2.2)')ro_data%dtocc%minute
1162	! write(se,'(i2.2)')ro_data%dtocc%second
1163	!
1164	! fname = '/home/rd/dam1/mi_liao/gnosrun/ropp-8.0/gnossh/atECM/multi/bangle_'//yr//mn//dy//'_'//hr//mi//se//'.txt'
1165	!
1166	! OPEN(77,file = trim(fname))
1167	! WRITE(77,'(6i6,2f12.3,2f16.3)')ro_data%DTocc%year,ro_data%DTocc%month,ro_data%DTocc%day,ro_data%DTocc%hour,&
1168	! ro_data%DTocc%minute,ro_data%dtocc%second,ro_data%GEOref%lat,ro_data%GEOref%lon,ro_data%GEOref%roc,&
1169	! ro_data%GEOref%undulation
1170	!
1171	! DO ii = 1, inb
1172	! WRITE(77,*)ii,(ro_data%lev1b%impact(ii)-ro_data%GEOref%roc),ro_data%lev1b%bangle(ii)
1173	! ENDDO
1174	!
1175	! CLOSE(77)
1176	!-------------------------------------------------------------------------------
1177	! 24. Clean up
1178	!-------------------------------------------------------------------------------
1179
1180
1181	! new
1182
1183	DEALLOCATE(snr_l2)
1184	DEALLOCATE(phase_l2)
1185
1186	CALL ropp_io_free(ro_data)
1187	CALL ropp_io_free(smt_ba)
1188	CALL ropp_io_free(bangle)
1189
1190	! new
1191
1192	CALL ropp_io_free(cma_bangle)
1193
1194	CALL ropp_io_free(out_ba)
1195	CALL ropp_io_free(out_refrac)
1196	CALL ropp_io_free(mod_ba)
1197	IF (ASSOCIATED(diag%ba_ion)) DEALLOCATE(diag%ba_ion)
1198	IF (ASSOCIATED(diag%err_ion)) DEALLOCATE(diag%err_ion)
1199	IF (ASSOCIATED(diag%err_neut)) DEALLOCATE(diag%err_neut)
1200	IF (ASSOCIATED(diag%wt_data)) DEALLOCATE(diag%wt_data)
1201	IF (ASSOCIATED(diag%CTimpact)) DEALLOCATE(diag%CTimpact)
1202	IF (ASSOCIATED(diag%CTamplitude)) DEALLOCATE(diag%CTamplitude)
1203	IF (ASSOCIATED(diag%CTamplitude_smt)) DEALLOCATE(diag%CTamplitude_smt)
1204	IF (ASSOCIATED(diag%CTimpactL2)) DEALLOCATE(diag%CTimpactL2)
1205	IF (ASSOCIATED(diag%CTamplitudeL2)) DEALLOCATE(diag%CTamplitudeL2)
1206	IF (ASSOCIATED(diag%CTamplitudeL2_smt)) DEALLOCATE(diag%CTamplitudeL2_smt)
1207
1208	END DO
1209	END DO
1210
1211	CALL EXIT(msg_exit_status)
1212
1213	CONTAINS
1214
1215	!-------------------------------------------------------------------------------
1216	! 25. Usage information
1217	!-------------------------------------------------------------------------------
1218
1219	SUBROUTINE usage()
1220	PRINT *, 'Purpose:'
1221	PRINT *, ' Calculate corrected bending angle, refractivity and '
1222	PRINT *, ' dry temperature from L1 and L2 phase data'
1223	PRINT *, 'Usage:'
1224	PRINT *, ' > ropp_pp_occ_tool [<options>] <input_file(s)>'
1225	PRINT *, 'Options'
1226	PRINT *, ' -c <config_file> name of configuration file'
1227	PRINT *, ' -o <output_file> name of ROPP netCDF output file'
1228	PRINT *, ' -m <method> ionospheric correction method'
1229	PRINT *, ' [NONE, MSIS, GMSIS, BG] (default MSIS)'
1230	PRINT *, ' -mfile <file> model refractivity coefficients file path'
1231	PRINT *, ' (default local search)'
1232	PRINT *, ' -bfile <file> background model atmospheric profile file path'
1233	PRINT *, ' (if using BG ionospheric correction method)'
1234	PRINT , ' -navfile <file> external navigation bit _txt file path'
1235	PRINT *, ' (default internal correction)'
1236	PRINT *, ' -occ <method> processing method, WO or GO (default WO)'
1237	PRINT *, ' -filter <method> filtering method, slpoly or optest'
1238	PRINT *, ' default slpoly, sliding polynomial)'
1239	PRINT *, ' -fit apply 2-parameter regression fit to model'
1240	PRINT *, ' -ellipsoid output height with respect to WGS84 ellipsoid'
1241	PRINT *, ' (default output wrt EGM96 geoid)'
1242	PRINT *, ' -d output additional diagnostics'
1243	PRINT *, ' -h this help'
1244	PRINT *, ' -v version information'
1245	PRINT *, ''
1246	END SUBROUTINE usage
1247
1248	!-------------------------------------------------------------------------------
1249	! 26. Version information
1250	!-------------------------------------------------------------------------------
1251
1252	SUBROUTINE version_info()
1253	CHARACTER (LEN=40) :: version
1254	version = ropp_pp_version()
1255	PRINT *, 'ropp_pp_occ_tool - Pre-processor occ tool: '
1256	PRINT *, ' Calculate corrected bending angles, '
1257	PRINT *, ' refractivity and dry temperature '
1258	PRINT *, ' from excess phase'
1259	PRINT *, ''
1260	PRINT *, 'This program is part of ROPP (PP) Release ' // TRIM(version)
1261	PRINT *, ''
1262	END SUBROUTINE version_info
1263
1264	SUBROUTINE cma_merge(npoints, min_slta, roc, cma_bangle, bangle, noise_estimate)
1265
1266	INTEGER, INTENT(IN) :: npoints ! Number of data points
1267
1268	REAL(wp), INTENT(IN) :: roc,min_slta
1269	TYPE(L1btype), INTENT(IN) :: cma_bangle ! for CMA data
1270	TYPE(L1btype), INTENT(INOUT) :: bangle ! Output bending angle profiles
1271	REAL(wp), INTENT(INOUT) :: noise_estimate
1272
1273	! local
1274
1275	REAL(wp) :: new_bangle_l1(npoints)
1276	REAL(wp), ALLOCATABLE :: new_wo_bangle_l1(:)
1277
1278	REAL(wp) :: x_opt
1279
1280	REAL(wp) :: lowest_height, impact_height, wt, r_ion, H, Hy, Hy_sum,H_sq_sum
1281	REAL(wp) :: j_pen
1282
1283	INTEGER :: i,inum,n_wopt,iused
1284
1285
1286	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1287	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1288	! This science code for testing the gnos data.
1289	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1290	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1291
1292	!
1293	! interpolate L1 onto L2 impact parameters
1294	!
1295
1296	CALL ropp_pp_interpol&
1297	(cma_bangle%impact_l1 , cma_bangle%impact_l2, cma_bangle%bangle_l1, new_bangle_l1)
1298
1299	!
1300	! for the least squares fit.
1301	!
1302
1303	H_sq_sum = 0.0_wp
1304	Hy_sum = 0.0_wp
1305	inum = 0
1306
1307	lowest_height = 1.0E20
1308
1309	r_ion = roc + 3.0E5_wp ! assume peak of ionosphere at 300 km
1310
1311	DO i = 1,npoints
1312
1313	! print*,'snr',i,snr_l2(i),cma_bangle%bangle_l2(i)
1314
1315	impact_height = cma_bangle%impact_l2(i)-roc
1316
1317	! screen out missing L2 data.
1318
1319	IF (impact_height < MIN(min_slta + 2.0E4_wp,7.0E4_wp) .and. impact_height >= min_slta) THEN ! 70 km is a guess
1320
1321	H = r_ion/(r_ion2 - cma_bangle%impact_l2(i)2)**1.5_wp ! See RSR 17, eq. 3.2
1322
1323	!print*,' HHHHH', H
1324
1325	IF (ABS(cma_bangle%bangle_l2(i) - new_bangle_l1(i)) < 1.0E-3_wp) THEN ! screen out very large departures
1326
1327	Hy = H*(cma_bangle%bangle_l2(i) - new_bangle_l1(i))
1328
1329	Hy_sum = Hy_sum + Hy
1330
1331	! PRINT*,' HYYYYYY HEIGHT', impact_height, cma_bangle%impact_l2(i),bangle%impact_l2(i)
1332
1333	H_sq_sum = H_sq_sum + H*H
1334
1335	inum = inum + 1
1336
1337	ENDIF
1338
1339	ENDIF
1340
1341	ENDDO
1342	! print *,'inum1', inum
1343	! The SLTA where the data is missing
1344
1345	lowest_height = min_slta
1346
1347	! compute the least square estimate
1348
1349	x_opt = Hy_sum/H_sq_sum
1350
1351	! print*,'new-routine',x_opt,lowest_height
1352
1353	!
1354	! compute cost function and then the noise estimate
1355	!
1356
1357	j_pen = 0.0_wp
1358
1359	DO i = 1,npoints
1360
1361	impact_height = cma_bangle%impact_l2(i)-roc
1362
1363	! screen out missing L2 data.
1364
1365	IF (impact_height < MIN(min_slta + 2.0E4_wp,7.0E4_wp) .and. impact_height >= min_slta) THEN ! 70 km is a guess
1366
1367	H = r_ion/(r_ion2 - cma_bangle%impact_l2(i)2)**1.5_wp ! See RSR 17, eq. 3.2
1368
1369	IF (ABS(cma_bangle%bangle_l2(i) - new_bangle_l1(i)) < 1.0E-3_wp) THEN ! screen out very large departures
1370
1371	j_pen = j_pen + (Hx_opt-(cma_bangle%bangle_l2(i) - new_bangle_l1(i)))*2
1372
1373	ENDIF
1374
1375	ENDIF
1376
1377	ENDDO
1378
1379	!
1380	! Now estimate the fitting noise from the data. A standard deviation = sigma
1381	!
1382	IF (inum > 1) THEN
1383	noise_estimate = SQRT(MAX(j_pen,1.0E-18_wp)/REAL(MAX(inum,1)))
1384	ELSE
1385	noise_estimate = 9.9E-05_wp
1386	ENDIF
1387
1388
1389	! We use x_opt to extrapolate the L2-L1 differences
1390	!
1391	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392
1393	IF (inum > 10) THEN ! simple check that we fitted something.
1394
1395	! the number of bending angles
1396
1397	n_wopt = SIZE(bangle%impact_l2)
1398
1399	! allocate array for interpolation
1400
1401	ALLOCATE(new_wo_bangle_l1(n_wopt))
1402
1403	! interpolate to L2
1404
1405	CALL ropp_pp_interpol&
1406	(bangle%impact_l1 , bangle%impact_l2, bangle%bangle_l1, new_wo_bangle_l1)
1407
1408	DO i = 1,n_wopt
1409
1410
1411	impact_height = bangle%impact_l2(i)-roc
1412
1413	H = r_ion/(r_ion2 - bangle%impact_l2(i)2)**1.5_wp
1414
1415	! Try to identify any gross errors in fitting interval
1416
1417
1418
1419	IF (impact_height < lowest_height) THEN
1420
1421
1422	! overwrite with extrapolated (L2-L1) correction
1423
1424	bangle%bangle_l2(i) = new_wo_bangle_l1(i) + H*x_opt
1425
1426
1427	ENDIF
1428
1429
1430	ENDDO
1431
1432	! tidy up
1433
1434	DEALLOCATE(new_wo_bangle_l1)
1435
1436	ENDIF
1437
1438	END SUBROUTINE
1439
1440
1441
1442	SUBROUTINE vertical_point(Pleo,Pgps,Pver)
1443
1444	IMPLICIT NONE
1445
1446	real(wp), dimension(:), intent(in) :: Pleo
1447	real(wp), dimension(:), intent(in) :: Pgps
1448	real(wp), dimension(:), intent(out) :: Pver
1449	REAL*8 t
1450
1451	t=-((Pgps(1)-Pleo(1))Pleo(1)+(Pgps(2)-Pleo(2))Pleo(2)+(Pgps(3)-Pleo(3))*Pleo(3)) &
1452	/((Pgps(1)-Pleo(1))2+(Pgps(2)-Pleo(2))2+(Pgps(3)-Pleo(3))**2)
1453
1454	Pver(1)=(Pgps(1)-Pleo(1))*t+Pleo(1)
1455	Pver(2)=(Pgps(2)-Pleo(2))*t+Pleo(2)
1456	Pver(3)=(Pgps(3)-Pleo(3))*t+Pleo(3)
1457
1458
1459	END SUBROUTINE
1460
1461	END PROGRAM ropp_pp_occ_tool

Download in other formats:

Original format