Context navigation

Back to Ticket #314

Ticket #314: ropp_fm_bg2ro_1d.f90

File ropp_fm_bg2ro_1d.f90, 31.3 KB (added by Ian Culverwell, 10 years ago)
ropp_fm_bg2ro_1d.f90

Line
1	! $Id: ropp_fm_bg2ro_1d.f90 4452 2015-01-29 14:42:02Z idculv $
2
3	PROGRAM ropp_fm_bg2ro_1d
4
5	!***p Programs/ropp_fm_bg2ro_1d *
6	!
7	! NAME
8	! ropp_fm_bg2ro_1d
9	!
10	! SYNOPSIS
11	! Calculate radio occultation pseudo observation from background model
12	! data using 1d forward models
13	!
14	! > ropp_fm_bg2ro_1d <infile(s)> -o <outfile>
15	! [-l] <levels_file> [-f] [-use_logp] [-use_logq]
16	! [-comp] [-check_qsat] [-new_op] [-direct_ion] [-247L]
17	! -zmin <geop_min> -zmax <geop_max> -nz <n_geop>
18	! [-d] [-h] [-v]
19	!
20	! ARGUMENTS
21	! <infile(s)> One (or more) input file names.
22	!
23	! OPTIONS
24	! -o <outfile> name of output file (default: bg2ro.nc)
25	! -l <levels_file> name of file containing a set of altitudes and
26	! impact heights. This will typically be an
27	! observation file. Only works for single file.
28	! -f forward model only, no gradient calculation.
29	! -use_logp use log(pressure) for forward model
30	! -use_logq use log(spec humidity) for forward model
31	! -comp include non ideal gas compressibility
32	! -check_qsat check against saturation
33	! -new_op use alternative refrac and bangle interpolation
34	! -direct_ion forward model L1 and L2 directly, using model ionosphere
35	! -247L output refrac and bangle on 'standard' 247L
36	! -zmin <geop_min> minimum refractivity geopotential (gpm)'
37	! -zmax <geop_max> maximum refractivity geopotential (gpm)'
38	! -nz <n_geop> number of uniformly spaced refractivity geopotentials (gpm)'
39	! -d output additional diagnostics
40	! -h help
41	! -v version information
42	!
43	! DESCRIPTION
44	! This program reads model data on model levels from the input data files
45	! and calculates vertical profiles of bending angle and refractivity using
46	! the 1d forward operators. The result is written to an ROPP formatted
47	! output file.
48	!
49	! NOTES
50	! If the input file is a multifile, or more than one input files are
51	! specified, the output file is a multifile.
52	!
53	! Already existing output files will be overwritten.
54	!
55	! EXAMPLES
56	! 1) To calculate bending angle and refractivity from one of the example
57	! (single-) files in the data directory:
58	!
59	! > ropp_fm_bg2ro_1d ../data/bgr20090401_000329_M02_2030337800_N0007_XXXX.nc
60	!
61	! 2) To calculate bending angle and refractivity profiles from all singlefiles
62	! in the data directory:
63	!
64	! > ropp_fm_bg2ro_1d ../data/bgr20090401*_N0007_XXXX.nc -o eg_02.nc
65	!
66	! Note that the resulting eg_02.nc file contains forward modelled data from
67	! all example profiles.
68	!
69	! 3) To calculate forward modelled bending angle and refractivity profiles from
70	! all profiles contained in the multifile bgr20090401_multi.nc:
71	!
72	! > ropp_fm_bg2ro_1d ../data/bgr20090401_multi.nc -o eg_03.nc
73	!
74	! Since the ecmwf_multi_* file was generated by concatenating the other
75	! files in the data directory, eg_02.nc and eg_03.nc should be identical
76	! apart from the file names.
77	!
78	! 4) To calculate forward modelled L1 and L2 bending angles, assuming
79	! a model ionosphere, from one of the example (single-) files
80	! in the data directory:
81	!
82	! > ropp_fm_bg2ro_1d --direct_ion ../data/bgr20090401_000329_M02_2030337800_N0007_YYYY.nc
83	!
84	! The model ionosphere comprises a set of Chapman layers, each of
85	! whose {Ne_max, H_peak, H_width} parameters are contained in the input file.
86	! (The L2ctype and state1dFM structures have been suitably extended.)
87	!
88	! SEE ALSO
89	! ropp_fm_bangle_1d
90	! ropp_fm_refrac_1d
91	! ropp_fm_refrac_1d_new
92	!
93	! AUTHOR
94	! Met Office, Exeter, UK.
95	! Any comments on this software should be given via the ROM SAF
96	! Helpdesk at http://www.romsaf.org
97	!
98	! COPYRIGHT
99	! (c) EUMETSAT. All rights reserved.
100	! For further details please refer to the file COPYRIGHT
101	! which you should have received as part of this distribution.
102	!
103	!****
104
105	!-------------------------------------------------------------------------------
106	! 1. Declarations
107	!-------------------------------------------------------------------------------
108
109	USE typesizes, ONLY: wp => EightByteReal
110	USE ropp_utils
111	USE ropp_io
112	USE ropp_io_types, ONLY: ROprof
113	USE ropp_fm
114	USE ropp_fm_types
115	USE ropp_fm_copy
116	USE ropp_fm_iono
117	USE ropp_fm_levels
118
119	IMPLICIT NONE
120
121	TYPE(ROprof) :: ro_data, l_data
122	TYPE(State1dFM) :: state
123	TYPE(Obs1dRefrac) :: obs_refrac
124	TYPE(Obs1dBangle) :: obs_bangle
125
126	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gradient_bangle
127	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gradient_refrac
128	REAL(wp) :: geop_min
129	REAL(wp) :: geop_max
130	INTEGER :: n_geop
131
132	INTEGER :: idummy
133	INTEGER :: i, iargc, argc, k
134	INTEGER :: n_files, n_profiles
135
136	LOGICAL :: give_help = .FALSE.
137	LOGICAL :: do_adj = .TRUE.
138	LOGICAL :: lfile_exist = .FALSE.
139	LOGICAL :: ranchk = .TRUE.
140	LOGICAL :: compress = .FALSE.
141	LOGICAL :: checksat = .FALSE.
142	LOGICAL :: new_op = .FALSE.
143	LOGICAL :: direct_ion = .FALSE.
144	LOGICAL :: use_247L = .FALSE.
145
146	CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
147	CHARACTER(len = 4096) :: ofile = 'bg2ro.nc'
148	CHARACTER(len = 4096) :: lfile
149	CHARACTER(len = 256) :: buffer
150	CHARACTER(len = 4) :: istr
151	CHARACTER(len = 6) :: nstr
152
153	!-------------------------------------------------------------------------------
154	! 2. Default settings
155	!-------------------------------------------------------------------------------
156
157	CALL message_set_routine ( 'ropp_fm_bg2ro_1d' )
158
159	CALL message(msg_noin, '')
160	CALL message(msg_noin, &
161	'-----------------------------------------------------------------------')
162	CALL message(msg_noin, &
163	' ROPP Forward Model' )
164	CALL message(msg_noin, &
165	'-----------------------------------------------------------------------')
166	CALL message(msg_noin, '')
167
168	!-------------------------------------------------------------------------------
169	! 3. Command line arguments
170	!-------------------------------------------------------------------------------
171
172	argc = iargc()
173	i = 1
174	n_files = 0
175	ALLOCATE(ifiles(argc))
176
177	geop_min = 200.0_wp
178	geop_max = 60000.0_wp
179	n_geop = 300
180
181	DO WHILE(i <= argc)
182	CALL getarg(i, buffer)
183	SELECT CASE (buffer)
184	CASE('-o') ! Output file name (netCDF output)
185	CALL getarg(i+1, buffer)
186	ofile = buffer
187	i = i + 1
188	CASE ('-l') ! Use level structure defined here
189	CALL getarg(i+1, buffer)
190	lfile = buffer
191	i = i + 1
192	lfile_exist = .TRUE.
193	CASE ('-f') ! Perform forward model only
194	do_adj = .FALSE.
195	CASE ('-use_logp') ! Use log(pressure) for FM
196	state%use_logp = .TRUE.
197	CASE ('-use_logq') ! Use log(shum) for FM
198	state%use_logq = .TRUE.
199	CASE ('-comp') ! Non ideal gas switch
200	compress = .TRUE.
201	CASE ('-check_qsat') ! Check against saturation
202	checksat = .TRUE.
203	CASE ('-new_op') ! New interpolation
204	new_op = .TRUE.
205	CASE ('-direct_ion') ! Model L1 and L2 bangles directly
206	direct_ion = .TRUE.
207	CASE ('--no-ranchk') ! Disable rangecheck; for experts only
208	CALL message ( msg_warn, "Range checking is disabled" )
209	ranchk = .FALSE. ! Do not document as a user option.
210	CASE ('-247L', '-247l', '-247') ! Output refrac and bangle on standard 247L
211	use_247L = .TRUE.
212	CASE ('-zmin', '-Zmin', '-ZMIN') ! Min refrac geopotential
213	CALL getarg(i+1, buffer)
214	READ(buffer, *) geop_min
215	i = i + 1
216	CASE ('-zmax', '-Zmax', '-ZMAX') ! Max refrac geopotential
217	CALL getarg(i+1, buffer)
218	READ(buffer, *) geop_max
219	i = i + 1
220	CASE ('-nz', '-Nz', '-NZ') ! Number of refractivity levels
221	CALL getarg(i+1, buffer)
222	READ(buffer, *) n_geop
223	i = i + 1
224	CASE('-d') ! Additional diagnostic mode
225	msg_MODE = VerboseMode
226	CASE('-h', '--help', '?') ! Help
227	give_help = .TRUE.
228	CASE('-v', '-V', '--version') ! Version information
229	CALL version_info()
230	CALL EXIT(msg_exit_ok)
231	CASE default ! Input file name
232	IF ( buffer(1:1) /= '-' ) THEN
233	n_files = n_files + 1
234	ifiles(n_files) = buffer
235	END IF
236	END SELECT
237	i = i + 1
238	END DO
239
240	IF ( n_files == 0 .AND. .NOT. give_help ) &
241	CALL message ( msg_error, "No input file(s) specified" )
242
243	IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN
244	CALL usage()
245	CALL EXIT(msg_exit_status)
246	ENDIF
247
248	!-------------------------------------------------------------------------------
249	! 4. Remove pre-existing output file
250	!-------------------------------------------------------------------------------
251
252	CALL file_delete(ofile, idummy)
253
254	!-------------------------------------------------------------------------------
255	! 4. Set default units
256	!-------------------------------------------------------------------------------
257
258	CALL ropp_fm_set_units(ro_data)
259
260	!-------------------------------------------------------------------------------
261	! 5. Loop over all input files
262	!-------------------------------------------------------------------------------
263
264	DO k = 1, n_files
265
266	!-------------------------------------------------------------------------------
267	! 6. Loop over all profiles
268	!-------------------------------------------------------------------------------
269
270	n_profiles = ropp_io_nrec(ifiles(k))
271
272	DO i = 1, n_profiles
273
274	WRITE(istr, '(i4)') i
275	WRITE(nstr, '(i6)') n_profiles
276	CALL message(msg_noin, '')
277	CALL message(msg_info, "Processing profile " // istr // " of " // nstr )
278
279	!-------------------------------------------------------------------------------
280	! 7. Read data
281	!-------------------------------------------------------------------------------
282
283	CALL ropp_io_read(ro_data, ifiles(k), rec = i, ranchk = ranchk)
284
285	CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")
286
287	IF (lfile_exist) THEN
288
289	CALL ropp_io_read(l_data, lfile, ranchk=ranchk)
290
291	ro_data%dtocc = l_data%dtocc
292	ro_data%georef = l_data%georef
293	ro_data%lev1b = l_data%lev1b
294	ro_data%lev2a = l_data%lev2a
295
296	CALL ropp_io_free(l_data)
297
298	ENDIF
299
300	IF (direct_ion) THEN
301
302	ro_data%lev2c%direct_ion = .TRUE.
303
304	CALL ropp_fm_iono_set_default(ro_data)
305
306	END IF
307
308	!-------------------------------------------------------------------------------
309	! 8. Copy data in RO structure to state and refrac obs vectors
310	!-------------------------------------------------------------------------------
311
312	IF (checksat) state%check_qsat = .TRUE. ! Switch on check against saturation
313
314	CALL ropp_fm_roprof2state(ro_data, state)
315
316	IF (compress) state%non_ideal = .TRUE. ! Switch on non-ideal gas compressibility option
317
318	IF (new_op) THEN
319	state%new_ref_op = .TRUE. ! Use new refrac interpolation
320	state%new_bangle_op = .TRUE. ! Use new bangle interpolation
321	ENDIF
322
323	IF (ro_data%lev2a%Npoints > 0) THEN
324	CALL ropp_fm_roprof2obs(ro_data, obs_refrac) ! Set to obs levels
325	obs_refrac%refrac(:) = ropp_MDFV
326	ELSE
327	IF (use_247L) THEN
328	CALL set_obs_levels_refrac_247(ro_data, obs_refrac, state) ! Pre-defined levels
329	ELSE
330	CALL set_obs_levels_refrac(ro_data, obs_refrac, &
331	geop_min=geop_min, geop_max=geop_max, n_geop=n_geop) ! Uniformly spaced levels
332	ENDIF
333	ENDIF
334
335	!-------------------------------------------------------------------------------
336	! 9. Calculate refractivity and its gradient
337	!-------------------------------------------------------------------------------
338
339	IF (state%state_ok) THEN
340
341	IF (state%new_ref_op) THEN
342	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
343	ELSE
344	CALL ropp_fm_refrac_1d(state, obs_refrac)
345	END IF
346
347	IF (do_adj) THEN
348	ALLOCATE(gradient_refrac(SIZE(obs_refrac%refrac),SIZE(state%state)))
349	CALL ropp_fm_refrac_1d_grad(state, obs_refrac, gradient_refrac)
350	ENDIF
351
352	ENDIF
353
354	!-------------------------------------------------------------------------------
355	! 10. Copy data in RO and refrac structure to bending angle obs vector
356	!-------------------------------------------------------------------------------
357
358	IF (ro_data%lev1b%Npoints > 0) THEN
359	CALL ropp_fm_roprof2obs(ro_data, obs_bangle) !!set to obs levels
360	obs_bangle%bangle(:) = ropp_MDFV
361	ELSE
362	CALL set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle) !default
363	END IF
364
365	!-------------------------------------------------------------------------------
366	! 11. Calculate bending angle and its gradient
367	!-------------------------------------------------------------------------------
368
369	IF (state%state_ok) THEN
370
371	CALL ropp_fm_bangle_1d(state, obs_bangle)
372
373	IF (do_adj .AND. (.NOT. direct_ion)) THEN
374	ALLOCATE(gradient_bangle(SIZE(obs_bangle%bangle),SIZE(state%state)))
375	CALL ropp_fm_bangle_1d_grad(state, obs_bangle, gradient_bangle)
376	END IF
377
378	ENDIF
379
380	!-------------------------------------------------------------------------------
381	! 12. Copy simulated observations to RO structure and write data
382	!-------------------------------------------------------------------------------
383
384	CALL ropp_fm_obs2roprof(obs_refrac, ro_data)
385
386	CALL ropp_fm_obs2roprof(obs_bangle, ro_data)
387
388	IF (state%state_ok) THEN
389
390	IF (do_adj) THEN
391	CALL ropp_io_addvar(ro_data, &
392	name = "gradient_refrac", &
393	long_name = &
394	"Gradient of the refractivity forward model", &
395	units = "1", &
396	range = (/MINVAL(gradient_refrac), &
397	MAXVAL(gradient_refrac)/), &
398	DATA = gradient_refrac)
399	ENDIF
400
401	IF (do_adj .AND. (.NOT. direct_ion)) THEN
402	CALL ropp_io_addvar(ro_data, &
403	name = "gradient_bangle" , &
404	long_name = &
405	"Gradient of the bending angle forward model", &
406	units = "rad", &
407	range = (/MINVAL(gradient_bangle), &
408	MAXVAL(gradient_bangle)/), &
409	DATA = gradient_bangle)
410	ENDIF
411
412	!-------------------------------------------------------------------------------
413	! 13. Update RO structure with computed state vector variables
414	!-------------------------------------------------------------------------------
415
416	CALL ropp_fm_state2roprof(state, ro_data)
417
418	END IF
419
420	!-------------------------------------------------------------------------------
421	! 14. Write data
422	!-------------------------------------------------------------------------------
423
424	IF (direct_ion) CALL ropp_fm_iono_unpack(ro_data)
425
426	CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk)
427
428	!-------------------------------------------------------------------------------
429	! 15. Clean up
430	!-------------------------------------------------------------------------------
431
432	IF (state%state_ok .AND. do_adj) THEN
433	DEALLOCATE (gradient_refrac)
434	ENDIF
435
436	IF (state%state_ok .AND. do_adj .AND. (.NOT. direct_ion)) THEN
437	DEALLOCATE (gradient_bangle)
438	ENDIF
439
440	CALL ropp_fm_free(state)
441	CALL ropp_fm_free(obs_refrac)
442	CALL ropp_fm_free(obs_bangle)
443	CALL ropp_io_free(ro_data)
444
445	END DO
446
447	END DO
448
449	CALL EXIT(msg_exit_status)
450
451	CONTAINS
452
453	!-------------------------------------------------------------------------------
454	! 16a. Calculate uniformly spaced observation levels for refractivity
455	!-------------------------------------------------------------------------------
456
457	SUBROUTINE set_obs_levels_refrac(ro_data, obs_refrac, &
458	geop_min, geop_max, n_geop)
459	!
460	! n_geop uniformly spaced geopotential heights between geop_min and geop_max (gpm).
461	!
462
463	! 16a.1 Declarations
464	! -----------------
465
466	USE typesizes, ONLY: wp => EightByteReal
467	USE ropp_io_types
468	USE ropp_fm
469
470	IMPLICIT NONE
471
472	! Input/Output variables
473	TYPE(ROprof) :: ro_data
474	TYPE(Obs1dRefrac) :: obs_refrac
475	REAL(wp), OPTIONAL :: geop_min
476	REAL(wp), OPTIONAL :: geop_max
477	INTEGER, OPTIONAL :: n_geop
478
479	! Local variables
480	REAL(wp) :: zmin, zmax
481	INTEGER :: i, nz
482	CHARACTER(len=256) :: routine
483
484	CALL message_get_routine(routine)
485	CALL message_set_routine('set_obs_levels_refrac')
486
487	! 16a. Vertical geopotential height levels between geop_min and geop_max
488	! ----------------------------------------------------------------------
489
490	zmin = 200.0_wp
491	IF ( PRESENT(geop_min) ) zmin = geop_min
492
493	zmax = 60.0e3_wp
494	IF ( PRESENT(geop_max) ) zmax = geop_max
495
496	nz = 300
497	IF ( PRESENT(n_geop) ) nz = n_geop
498
499	IF ( (zmin >= zmax) .OR. (n_geop <= 1) ) THEN
500	CALL message(msg_warn, 'Invalid geop_min, geop_max or n_geop' // &
501	' ... resetting all to default values')
502	zmin = 200.0_wp
503	zmax = 60000.0_wp
504	nz = 300
505	END IF
506
507	ALLOCATE(obs_refrac%refrac(nz))
508	ALLOCATE(obs_refrac%geop(nz))
509	ALLOCATE(obs_refrac%weights(nz))
510
511	obs_refrac%refrac(:) = 0.0_wp
512	obs_refrac%geop(:) = (/ (zmin+(i-1)*(zmax-zmin)/REAL(nz-1,wp), i=1,nz) /)
513	obs_refrac%weights(:) = 1.0_wp
514	obs_refrac%lon = ro_data%GEOref%lon
515	obs_refrac%lat = ro_data%GEOref%lat
516
517	CALL message_set_routine(routine)
518
519	END SUBROUTINE set_obs_levels_refrac
520
521	!-------------------------------------------------------------------------------
522	! 16b. Calculate standard 247 observation levels for refractivity
523	!-------------------------------------------------------------------------------
524
525	SUBROUTINE set_obs_levels_refrac_247(ro_data, obs_refrac, state)
526	!
527	! 247 geopotential height levels, going up to about 60 km, chosen in such a
528	! way that the resulting impact heights match the 'standard' ones in
529	! impact_height_eum_247 in ropp_fm/common/ropp_fm_levels.f90
530	!
531	! 16b.1 Declarations
532	! -----------------
533
534	USE typesizes, ONLY: wp => EightByteReal
535	USE ropp_io_types
536	USE ropp_fm
537
538	IMPLICIT NONE
539
540	TYPE(ROprof) :: ro_data
541	TYPE(Obs1dRefrac) :: obs_refrac, obs_refrac2
542	TYPE(State1dFM) :: state
543
544	INTEGER, PARAMETER :: n=247
545	REAL(wp), DIMENSION(n) :: impact_height
546	REAL(wp), DIMENSION(n) :: tmp
547	REAL(wp), PARAMETER :: refrac_surface_default=300.0_wp
548	REAL(wp), PARAMETER :: max_refrac=400.0_wp
549	REAL(wp) :: refrac_surface
550	REAL(wp), PARAMETER :: tol=1.0e-3_wp ! Demand < 1mm diffs
551	REAL(wp) :: roc, und
552	INTEGER, PARAMETER :: max_iter=500
553	INTEGER :: iter, imin
554	CHARACTER(LEN=3) :: smax_iter
555	CHARACTER(LEN=256) :: routine
556
557	! 16b.2 Initialisations
558	! ---------------------
559
560	CALL message_get_routine(routine)
561	CALL message_set_routine('set_obs_levels_refrac_247')
562
563	impact_height = impact_height_eum_247 ! The particular 247L set
564
565	ALLOCATE(obs_refrac%refrac(n))
566	ALLOCATE(obs_refrac%geop(n))
567	ALLOCATE(obs_refrac%weights(n))
568
569	ALLOCATE(obs_refrac2%refrac(n))
570	ALLOCATE(obs_refrac2%geop(n))
571	ALLOCATE(obs_refrac2%weights(n))
572
573	! 16b.3 Vertical geopotential height levels
574	! -----------------------------------------
575
576	! We seek solutions h_i of
577	! (nr = ) (1 + 1e-6 N(Z(h_i)))*(h_i + RoC + und) = a_i + RoC
578	! where a_i are the given impact heights, RoC is the radius of curvature
579	! and und is the undulation. h_i is the geometric height;
580	! Z_i is the geopotential height. Non-linear because N depends on h.
581	! Simple iteration converges to a solution (slowly).
582	! When geop_refrac is below the lowest state refrac, the refracs are
583	! (now) missing data. In such cases, replace the missing refracs with ones
584	! calculated with the old refrac scheme with this nullification switched off.
585	! Because the same principle is applied in set_obs_levels_bangle, the
586	! resulting impact heights are still good.
587
588	! 16b.4 Generate radius of curvature
589	! ----------------------------------
590
591	IF ( ro_data%GEOref%roc > ropp_MDTV ) THEN
592	roc = ro_data%GEOref%roc
593	ELSE
594	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
595	'setting equal to effective radius')
596	roc = R_eff(ro_data%GEOref%lat) ! Maybe +21km? See RSR 14.
597	ENDIF
598
599	! 16b.5 Generate undulation
600	! -------------------------
601
602	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
603	und = ro_data%GEOref%undulation
604	ELSE
605	CALL message(msg_warn, "Undulation missing ... will assume to " // &
606	"be zero when calculating geopotential heights.")
607	und = 0.0_wp
608	ENDIF
609
610	! 16b.6 Initialise
611	! ----------------
612
613	IF ( (state%temp(1) > ropp_MDTV) .AND. (state%pres(1) > ropp_MDTV) ) THEN
614	refrac_surface = kappa1 * state%pres(1) / state%temp(1) ! roughly
615	ELSE
616	refrac_surface = refrac_surface_default
617	END IF
618
619	obs_refrac%geop = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*refrac_surface)) - &
620	(roc + und) ! roughly
621
622	CALL ropp_fm_refrac_1d(state, obs_refrac, nonull=.TRUE.) ! Initialise with the old interp scheme sans missing data
623
624	! 16b.7 Iterate on GPH
625	! --------------------
626
627	iter = 0
628
629	DO
630
631	tmp = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*obs_refrac%refrac)) - &
632	(roc + und) ! geometric height
633
634	obs_refrac%geop = geometric2geopotential(ro_data%GEOref%lat, tmp)
635
636	IF ( state%new_ref_op ) THEN
637	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
638	ELSE
639	CALL ropp_fm_refrac_1d(state, obs_refrac)
640	END IF
641
642	IF ( COUNT(obs_refrac%refrac <= ropp_MDTV) > 0 ) THEN ! rederive approximately some zapped refracs
643
644	obs_refrac2 = obs_refrac
645
646	CALL ropp_fm_refrac_1d(state, obs_refrac2, nonull=.TRUE.) ! Recalculate sans missing data
647
648	WHERE (obs_refrac%refrac <= ropp_MDTV) &
649	obs_refrac%refrac = obs_refrac2%refrac
650
651	END IF
652
653	WHERE (obs_refrac%refrac < ropp_ZERO) obs_refrac%refrac = ropp_ZERO
654
655	WHERE (obs_refrac%refrac > max_refrac) obs_refrac%refrac = max_refrac
656
657	tmp = (1.0_wp + 1.0e-6_wpobs_refrac%refrac)(tmp + roc + und) - &
658	(impact_height + roc) ! residual
659
660	IF ( MAXVAL(ABS(tmp)) < tol ) EXIT
661
662	iter = iter + 1
663
664	IF ( iter > max_iter ) THEN
665	WRITE (smax_iter, '(i3)') max_iter
666	CALL message(msg_fatal, 'GPH iteration failed to converge in ' // &
667	smax_iter // ' iterations')
668	END IF
669
670	ENDDO
671
672	WRITE (smax_iter, '(i3)') iter
673	CALL message(msg_diag, 'GPH iteration converged in ' // &
674	smax_iter // ' iterations')
675
676	! 16b.8 Set other elements of obs_refrac structure
677	! ------------------------------------------------
678
679	obs_refrac%refrac = 0.0_wp ! Refrac will be calculated in the main program
680	obs_refrac%weights = 1.0_wp
681	obs_refrac%lon = ro_data%GEOref%lon
682	obs_refrac%lat = ro_data%GEOref%lat
683
684	! 16b.9 Clean up
685	! --------------
686
687	CALL ropp_fm_free(obs_refrac2)
688
689	CALL message_set_routine(routine)
690
691	END SUBROUTINE set_obs_levels_refrac_247
692
693
694	!-------------------------------------------------------------------------------
695	! 17. Calculate obs levels for bending angle (consistent with obs_refrac)
696	!-------------------------------------------------------------------------------
697
698	SUBROUTINE set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle)
699
700	! 17.1 Declarations
701	! -----------------
702
703	USE typesizes, ONLY: wp => EightByteReal
704	USE geodesy
705	USE ropp_io_types
706	USE ropp_fm
707
708	IMPLICIT NONE
709
710	TYPE(ROprof) :: ro_data
711	TYPE(Obs1dRefrac) :: obs_refrac, obs_refrac2
712	TYPE(Obs1dbangle) :: obs_bangle
713	TYPE(State1dFM) :: state
714
715	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp
716
717	INTEGER :: n, imin
718	REAL(wp), PARAMETER :: max_refrac=400.0_wp
719
720	CHARACTER(LEN=256) :: routine
721
722	CALL message_get_routine(routine)
723	CALL message_set_routine('set_obs_levels_bangle')
724
725	! 17.2 Allocate arrays
726	! --------------------
727
728	n = SIZE(obs_refrac%geop)
729
730	obs_bangle%n_L1 = n
731
732	ALLOCATE(tmp(n))
733
734	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
735	ALLOCATE(obs_bangle%bangle(n))
736	ALLOCATE(obs_bangle%impact(n))
737	ALLOCATE(obs_bangle%weights(n))
738	ELSE
739	ALLOCATE(obs_bangle%bangle(2*n))
740	ALLOCATE(obs_bangle%impact(2*n))
741	ALLOCATE(obs_bangle%weights(2*n))
742	ENDIF
743
744	ALLOCATE(obs_refrac2%refrac(n))
745	ALLOCATE(obs_refrac2%geop(n))
746	ALLOCATE(obs_refrac2%weights(n))
747
748	! 17.3 Calculate refractivities
749	! -----------------------------
750
751	IF ( state%new_ref_op ) THEN
752	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
753	ELSE
754	CALL ropp_fm_refrac_1d(state, obs_refrac)
755	END IF
756
757	IF ( COUNT(obs_refrac%refrac <= ropp_MDTV) > 0 ) THEN ! rederive approximately some zapped refracs
758
759	obs_refrac2 = obs_refrac
760
761	CALL ropp_fm_refrac_1d(state, obs_refrac2, nonull=.TRUE.) ! Estimate sans missing data
762
763	WHERE (obs_refrac%refrac <= ropp_MDTV) &
764	obs_refrac%refrac = obs_refrac2%refrac
765
766	END IF
767
768	WHERE (obs_refrac%refrac < ropp_ZERO) obs_refrac%refrac = ropp_ZERO
769
770	WHERE (obs_refrac%refrac > max_refrac) obs_refrac%refrac = max_refrac
771
772	! 17.4 Set scalar arguments of the observation vector
773	! ---------------------------------------------------
774
775	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
776	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
777	obs_bangle%undulation = ro_data%GEOref%undulation
778	obs_bangle%lon = ro_data%GEOref%lon
779	obs_bangle%lat = ro_data%GEOref%lat
780	obs_bangle%azimuth = ro_data%GEOref%azimuth
781
782	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
783	obs_bangle%r_curve = ro_data%GEOref%roc
784	ELSE
785	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
786	'setting equal to effective radius')
787	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
788	ENDIF
789
790	! 17.5 Calculate levels to coincide with the geopotential levels
791	! --------------------------------------------------------------
792
793	tmp = geopotential2geometric(ro_data%GEOref%lat, obs_refrac%geop) + &
794	obs_bangle%r_curve
795
796	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
797	tmp = tmp + obs_bangle%undulation
798	ELSE
799	CALL message(msg_warn, "Undulation missing. Will assume to " // &
800	"be zero when calculating impact parameters.")
801	ENDIF
802
803	tmp = (1.0_wp + obs_refrac%refrac1.e-6_wp) tmp
804
805	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
806	obs_bangle%impact = tmp
807	ELSE
808	obs_bangle%impact(1:n) = tmp
809	obs_bangle%impact(n+1:2*n) = tmp
810	ENDIF
811
812	! 17.6 Fill other arrays
813	! ----------------------
814
815	obs_bangle%bangle(:) = 0.0_wp
816	obs_bangle%weights(:) = 1.0_wp
817
818	! 17.7 Clean up
819	! -------------
820
821	IF ( state%new_ref_op ) THEN
822	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
823	ELSE
824	CALL ropp_fm_refrac_1d(state, obs_refrac) ! Regenerate avec missing data
825	END IF
826
827	DEALLOCATE(tmp)
828
829	CALL ropp_fm_free(obs_refrac2)
830
831	CALL message_set_routine(routine)
832
833	END SUBROUTINE set_obs_levels_bangle
834
835	!-------------------------------------------------------------------------------
836	! 18. Usage information
837	!-------------------------------------------------------------------------------
838
839	SUBROUTINE usage()
840	PRINT *, 'Purpose:'
841	PRINT *, ' Bending angles and refractivity forward model'
842	PRINT *, 'Usage:'
843	PRINT *, ' > ropp_fm_bg2ro_1d [<options>] <input_file(s)>'
844	PRINT *, 'Options:'
845	PRINT *, ' -o <output_file> name of ROPP netCDF output file'
846	PRINT *, ' -l <levels_file> optional name of (observation) file'
847	PRINT *, ' (non-multi) containing output level information'
848	PRINT *, ' -f forward model only, no gradient calculation'
849	PRINT *, ' -use_logp use log(pressure) for forward model'
850	PRINT *, ' -use_logq use log(spec humidity) for forward model'
851	PRINT *, ' -comp include non ideal gas compressibility'
852	PRINT *, ' -check_qsat include check against saturation'
853	PRINT *, ' -new_op use alternative refrac and bangle interpolation'
854	PRINT *, ' -direct_ion model L1 and L2 directly'
855	PRINT *, ' -247L output refrac and bangle on standard 247L'
856	PRINT *, ' -zmin <geop_min> minimum refractivity geopotential (gpm)'
857	PRINT *, ' -zmax <geop_max> maximum refractivity geopotential (gpm)'
858	PRINT *, ' -nz <n_geop> number of uniformly spaced refractivity geopotentials (gpm)'
859	PRINT *, ' -d output additional diagnostics'
860	PRINT *, ' -h this help'
861	PRINT *, ' -v version information'
862	PRINT *, ''
863	END SUBROUTINE usage
864
865	!-------------------------------------------------------------------------------
866	! 19. Version information
867	!-------------------------------------------------------------------------------
868
869	SUBROUTINE version_info()
870	CHARACTER (LEN=40) :: version
871	version = ropp_fm_version()
872	PRINT *, 'ropp_fm_bg2ro_1d - Bending angles and refractivity forward model.'
873	PRINT *, ''
874	PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
875	PRINT *, ''
876	END SUBROUTINE version_info
877
878	END PROGRAM ropp_fm_bg2ro_1d

Download in other formats:

Original format