Context navigation

Back to Ticket #314

Ticket #314: ropp_fm_bg2ro_1d.f90.14072015

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

Line
1	! $Id: ropp_fm_bg2ro_1d.f90 4432 2015-01-19 14:26:46Z 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	CHARACTER(len = 4) :: slev1b, slev2a
153
154	!-------------------------------------------------------------------------------
155	! 2. Default settings
156	!-------------------------------------------------------------------------------
157
158	CALL message_set_routine ( 'ropp_fm_bg2ro_1d' )
159
160	CALL message(msg_noin, '')
161	CALL message(msg_noin, &
162	'-----------------------------------------------------------------------')
163	CALL message(msg_noin, &
164	' ROPP Forward Model' )
165	CALL message(msg_noin, &
166	'-----------------------------------------------------------------------')
167	CALL message(msg_noin, '')
168
169	!-------------------------------------------------------------------------------
170	! 3. Command line arguments
171	!-------------------------------------------------------------------------------
172
173	argc = iargc()
174	i = 1
175	n_files = 0
176	ALLOCATE(ifiles(argc))
177
178	geop_min = 200.0_wp
179	geop_max = 60000.0_wp
180	n_geop = 300
181
182	DO WHILE(i <= argc)
183	CALL getarg(i, buffer)
184	SELECT CASE (buffer)
185	CASE('-o') ! Output file name (netCDF output)
186	CALL getarg(i+1, buffer)
187	ofile = buffer
188	i = i + 1
189	CASE ('-l') ! Use level structure defined here
190	CALL getarg(i+1, buffer)
191	lfile = buffer
192	i = i + 1
193	lfile_exist = .TRUE.
194	CASE ('-f') ! Perform forward model only
195	do_adj = .FALSE.
196	CASE ('-use_logp') ! Use log(pressure) for FM
197	state%use_logp = .TRUE.
198	CASE ('-use_logq') ! Use log(shum) for FM
199	state%use_logq = .TRUE.
200	CASE ('-comp') ! Non ideal gas switch
201	compress = .TRUE.
202	CASE ('-check_qsat') ! Check against saturation
203	checksat = .TRUE.
204	CASE ('-new_op') ! New interpolation
205	new_op = .TRUE.
206	CASE ('-direct_ion') ! Model L1 and L2 bangles directly
207	direct_ion = .TRUE.
208	CASE ('--no-ranchk') ! Disable rangecheck; for experts only
209	CALL message ( msg_warn, "Range checking is disabled" )
210	ranchk = .FALSE. ! Do not document as a user option.
211	CASE ('-247L', '-247l', '-247') ! Output refrac and bangle on standard 247L
212	use_247L = .TRUE.
213	CASE ('-zmin', '-Zmin', '-ZMIN') ! Min refrac geopotential
214	CALL getarg(i+1, buffer)
215	READ(buffer, *) geop_min
216	i = i + 1
217	CASE ('-zmax', '-Zmax', '-ZMAX') ! Max refrac geopotential
218	CALL getarg(i+1, buffer)
219	READ(buffer, *) geop_max
220	i = i + 1
221	CASE ('-nz', '-Nz', '-NZ') ! Number of refractivity levels
222	CALL getarg(i+1, buffer)
223	READ(buffer, *) n_geop
224	i = i + 1
225	CASE('-d') ! Additional diagnostic mode
226	msg_MODE = VerboseMode
227	CASE('-h', '--help', '?') ! Help
228	give_help = .TRUE.
229	CASE('-v', '-V', '--version') ! Version information
230	CALL version_info()
231	CALL EXIT(msg_exit_ok)
232	CASE default ! Input file name
233	IF ( buffer(1:1) /= '-' ) THEN
234	n_files = n_files + 1
235	ifiles(n_files) = buffer
236	END IF
237	END SELECT
238	i = i + 1
239	END DO
240
241	IF ( n_files == 0 .AND. .NOT. give_help ) &
242	CALL message ( msg_error, "No input file(s) specified" )
243
244	IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN
245	CALL usage()
246	CALL EXIT(msg_exit_status)
247	END IF
248
249	!-------------------------------------------------------------------------------
250	! 4. Remove pre-existing output file
251	!-------------------------------------------------------------------------------
252
253	CALL file_delete(ofile, idummy)
254
255	!-------------------------------------------------------------------------------
256	! 4. Set default units
257	!-------------------------------------------------------------------------------
258
259	CALL ropp_fm_set_units(ro_data)
260
261	!-------------------------------------------------------------------------------
262	! 5. Loop over all input files
263	!-------------------------------------------------------------------------------
264
265	DO k = 1, n_files
266
267	!-------------------------------------------------------------------------------
268	! 6. Loop over all profiles
269	!-------------------------------------------------------------------------------
270
271	n_profiles = ropp_io_nrec(ifiles(k))
272
273	DO i = 1, n_profiles
274
275	WRITE(istr, '(i4)') i
276	WRITE(nstr, '(i6)') n_profiles
277	CALL message(msg_noin, '')
278	CALL message(msg_info, "Processing profile " // istr // " of " // nstr )
279
280	!-------------------------------------------------------------------------------
281	! 7. Read data
282	!-------------------------------------------------------------------------------
283
284	CALL ropp_io_read(ro_data, ifiles(k), rec = i, ranchk = ranchk)
285
286	CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")
287
288	IF (lfile_exist) THEN
289
290	CALL message(msg_info, "Reading levels information from " // lfile // " \n")
291
292	CALL ropp_io_read(l_data, lfile, ranchk=ranchk)
293
294	ro_data%dtocc = l_data%dtocc
295	ro_data%georef = l_data%georef
296	ro_data%lev1b = l_data%lev1b
297	ro_data%lev2a = l_data%lev2a
298
299	WRITE (slev1b, '(i4)') ro_data%lev1b%Npoints
300	CALL message(msg_info, "There are " // slev1b // &
301	" level 1b points in " // lfile // " \n")
302
303	WRITE (slev2a, '(i4)') ro_data%lev2a%Npoints
304	CALL message(msg_info, "There are " // slev2a // &
305	" level 2a points in " // lfile // " \n")
306
307	CALL ropp_io_free(l_data)
308
309	END IF
310
311	IF (direct_ion) THEN
312
313	ro_data%lev2c%direct_ion = .TRUE.
314
315	CALL ropp_fm_iono_set_default(ro_data)
316
317	END IF
318
319	!-------------------------------------------------------------------------------
320	! 8. Copy data in RO structure to state and refrac obs vectors
321	!-------------------------------------------------------------------------------
322
323	IF (checksat) state%check_qsat = .TRUE. ! Switch on check against saturation
324
325	CALL ropp_fm_roprof2state(ro_data, state)
326
327	IF (compress) state%non_ideal = .TRUE. ! Switch on non-ideal gas compressibility option
328
329	IF (new_op) THEN
330	state%new_ref_op = .TRUE. ! Use new refrac interpolation
331	state%new_bangle_op = .TRUE. ! Use new bangle interpolation
332	END IF
333
334	IF (ro_data%lev2a%Npoints > 0) THEN
335	CALL ropp_fm_roprof2obs(ro_data, obs_refrac) ! Set to obs levels
336	obs_refrac%refrac(:) = ropp_MDFV
337	ELSE
338	IF (use_247L) THEN
339	CALL set_obs_levels_refrac_247(ro_data, obs_refrac, state) ! Pre-defined levels
340	ELSE
341	CALL set_obs_levels_refrac(ro_data, obs_refrac, &
342	geop_min=geop_min, geop_max=geop_max, n_geop=n_geop) ! Uniformly spaced levels
343	END IF
344	END IF
345
346	!-------------------------------------------------------------------------------
347	! 9. Calculate refractivity and its gradient
348	!-------------------------------------------------------------------------------
349
350	IF (state%state_ok) THEN
351
352	IF (state%new_ref_op) THEN
353	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
354	ELSE
355	CALL ropp_fm_refrac_1d(state, obs_refrac)
356	END IF
357
358	IF (do_adj) THEN
359	ALLOCATE(gradient_refrac(SIZE(obs_refrac%refrac),SIZE(state%state)))
360	CALL ropp_fm_refrac_1d_grad(state, obs_refrac, gradient_refrac)
361	END IF
362
363	END IF
364
365	!-------------------------------------------------------------------------------
366	! 10. Copy data in RO and refrac structure to bending angle obs vector
367	!-------------------------------------------------------------------------------
368
369	IF (ro_data%lev1b%Npoints > 0) THEN
370	CALL ropp_fm_roprof2obs(ro_data, obs_bangle) ! Set to obs levels
371	obs_bangle%bangle(:) = ropp_MDFV
372	ELSE
373	IF (use_247L) THEN
374	CALL set_obs_levels_bangle_247(ro_data, obs_bangle) ! Pre-defined levels
375	ELSE
376	CALL set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle) !default
377	END IF
378	END IF
379
380	!-------------------------------------------------------------------------------
381	! 11. Calculate bending angle and its gradient
382	!-------------------------------------------------------------------------------
383
384	IF (state%state_ok) THEN
385
386	CALL ropp_fm_bangle_1d(state, obs_bangle)
387
388	IF (do_adj .AND. (.NOT. direct_ion)) THEN
389	ALLOCATE(gradient_bangle(SIZE(obs_bangle%bangle),SIZE(state%state)))
390	CALL ropp_fm_bangle_1d_grad(state, obs_bangle, gradient_bangle)
391	END IF
392
393	END IF
394
395	!-------------------------------------------------------------------------------
396	! 12. Copy simulated observations to RO structure and write data
397	!-------------------------------------------------------------------------------
398
399	CALL ropp_fm_obs2roprof(obs_refrac, ro_data)
400
401	CALL ropp_fm_obs2roprof(obs_bangle, ro_data)
402
403	IF (state%state_ok) THEN
404
405	IF (do_adj) THEN
406	CALL ropp_io_addvar(ro_data, &
407	name = "gradient_refrac", &
408	long_name = &
409	"Gradient of the refractivity forward model", &
410	units = "1", &
411	range = (/MINVAL(gradient_refrac), &
412	MAXVAL(gradient_refrac)/), &
413	DATA = gradient_refrac)
414	END IF
415
416	IF (do_adj .AND. (.NOT. direct_ion)) THEN
417	CALL ropp_io_addvar(ro_data, &
418	name = "gradient_bangle" , &
419	long_name = &
420	"Gradient of the bending angle forward model", &
421	units = "rad", &
422	range = (/MINVAL(gradient_bangle), &
423	MAXVAL(gradient_bangle)/), &
424	DATA = gradient_bangle)
425	END IF
426
427	!-------------------------------------------------------------------------------
428	! 13. Update RO structure with computed state vector variables
429	!-------------------------------------------------------------------------------
430
431	CALL ropp_fm_state2roprof(state, ro_data)
432
433	END IF
434
435	!-------------------------------------------------------------------------------
436	! 14. Write data
437	!-------------------------------------------------------------------------------
438
439	IF (direct_ion) CALL ropp_fm_iono_unpack(ro_data)
440
441	CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk)
442
443	!-------------------------------------------------------------------------------
444	! 15. Clean up
445	!-------------------------------------------------------------------------------
446
447	IF (state%state_ok .AND. do_adj) THEN
448	DEALLOCATE (gradient_refrac)
449	END IF
450
451	IF (state%state_ok .AND. do_adj .AND. (.NOT. direct_ion)) THEN
452	DEALLOCATE (gradient_bangle)
453	END IF
454
455	CALL ropp_fm_free(state)
456	CALL ropp_fm_free(obs_refrac)
457	CALL ropp_fm_free(obs_bangle)
458	CALL ropp_io_free(ro_data)
459
460	END DO
461
462	END DO
463
464	CALL EXIT(msg_exit_status)
465
466	CONTAINS
467
468	!-------------------------------------------------------------------------------
469	! 16a. Calculate uniformly spaced observation levels for refractivity
470	!-------------------------------------------------------------------------------
471
472	SUBROUTINE set_obs_levels_refrac(ro_data, obs_refrac, &
473	geop_min, geop_max, n_geop)
474	!
475	! n_geop uniformly spaced geopotential heights between geop_min and geop_max (gpm).
476	!
477
478	! 16a.1 Declarations
479	! -----------------
480
481	USE typesizes, ONLY: wp => EightByteReal
482	USE ropp_io_types
483	USE ropp_fm
484
485	IMPLICIT NONE
486
487	! Input/Output variables
488	TYPE(ROprof) :: ro_data
489	TYPE(Obs1dRefrac) :: obs_refrac
490	REAL(wp), OPTIONAL :: geop_min
491	REAL(wp), OPTIONAL :: geop_max
492	INTEGER, OPTIONAL :: n_geop
493
494	! Local variables
495	REAL(wp) :: zmin, zmax
496	INTEGER :: i, nz
497	CHARACTER(len=256) :: routine
498
499	CALL message_get_routine(routine)
500	CALL message_set_routine('set_obs_levels_refrac')
501
502	! 16a. Vertical geopotential height levels between geop_min and geop_max
503	! ----------------------------------------------------------------------
504
505	zmin = 200.0_wp
506	IF ( PRESENT(geop_min) ) zmin = geop_min
507
508	zmax = 60.0e3_wp
509	IF ( PRESENT(geop_max) ) zmax = geop_max
510
511	nz = 300
512	IF ( PRESENT(n_geop) ) nz = n_geop
513
514	IF ( (zmin >= zmax) .OR. (n_geop <= 1) ) THEN
515	CALL message(msg_warn, 'Invalid geop_min, geop_max or n_geop' // &
516	' ... resetting all to default values')
517	zmin = 200.0_wp
518	zmax = 60000.0_wp
519	nz = 300
520	END IF
521
522	ALLOCATE(obs_refrac%refrac(nz))
523	ALLOCATE(obs_refrac%geop(nz))
524	ALLOCATE(obs_refrac%weights(nz))
525
526	obs_refrac%refrac(:) = 0.0_wp
527	obs_refrac%geop(:) = (/ (zmin+(i-1)*(zmax-zmin)/REAL(nz-1, wp), i=1,nz) /)
528	obs_refrac%weights(:) = 1.0_wp
529	obs_refrac%lon = ro_data%GEOref%lon
530	obs_refrac%lat = ro_data%GEOref%lat
531
532	CALL message_set_routine(routine)
533
534	END SUBROUTINE set_obs_levels_refrac
535
536
537	!-------------------------------------------------------------------------------
538	! 16b. Calculate standard 247 observation levels for refractivity
539	!-------------------------------------------------------------------------------
540
541	SUBROUTINE set_obs_levels_refrac_247(ro_data, obs_refrac, state)
542	!
543	! 247 geopotential height levels, going up to about 60 km, chosen in such a
544	! way that the resulting impact heights match the 'standard' ones in
545	! impact_height_eum_247 in ropp_fm/common/ropp_fm_levels.f90
546	!
547	! 16b.1 Declarations
548	! -----------------
549
550	USE typesizes, ONLY: wp => EightByteReal
551	USE ropp_io_types
552	USE ropp_fm
553
554	IMPLICIT NONE
555
556	TYPE(ROprof) :: ro_data
557	TYPE(Obs1dRefrac) :: obs_refrac
558	TYPE(State1dFM) :: state
559
560	INTEGER, PARAMETER :: n=247
561	REAL(wp), DIMENSION(n) :: impact_height
562	REAL(wp), DIMENSION(n) :: tmp
563	REAL(wp), PARAMETER :: refrac_surface_default=300.0_wp
564	REAL(wp) :: refrac_surface
565	REAL(wp), PARAMETER :: tol=1.0e-3_wp ! Demand < 1mm diffs
566	REAL(wp) :: max_refrac
567	REAL(wp) :: min_geop
568	REAL(wp) :: roc, und, lat
569	INTEGER, PARAMETER :: max_iter=1000
570	INTEGER :: iter, imiss
571	INTEGER :: lowest_refrac_index
572
573	CHARACTER(LEN=3) :: smax_iter
574	CHARACTER(LEN=256) :: routine
575
576	! 16b.2 Initialisations
577	! ---------------------
578
579	CALL message_get_routine(routine)
580	CALL message_set_routine('set_obs_levels_refrac_247')
581
582	impact_height = impact_height_eum_247 ! The particular 247L set
583
584	max_refrac = ro_data%lev2a%range%refrac(2)
585
586	min_geop = ro_data%lev2b%range%geop(1)
587
588	ALLOCATE(obs_refrac%refrac(n))
589	ALLOCATE(obs_refrac%geop(n))
590	ALLOCATE(obs_refrac%weights(n))
591
592	! 16b.3 Vertical geopotential height levels
593	! -----------------------------------------
594
595	! We seek solutions h_i of
596	! (nr = ) (1 + 1e-6 N(Z(h_i)))*(h_i + RoC + und) = a_i + RoC
597	! where a_i are the given impact heights, RoC is the radius of curvature
598	! and und is the undulation. h_i is the geometric height;
599	! Z_i is the geopotential height. Non-linear because N depends on Z.
600	! When geop_refrac is below the lowest state refrac, the refracs are
601	! (now) missing data. In such cases, replace the missing refracs with ones
602	! calculated with the old refrac scheme with this nullification switched off.
603	! Because the same principle is applied in set_obs_levels_bangle, the
604	! resulting impact heights are still good.
605
606	! 16b.4 Generate radius of curvature
607	! ----------------------------------
608
609	IF ( ro_data%GEOref%roc > ropp_MDTV ) THEN
610	roc = ro_data%GEOref%roc
611	ELSE
612	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
613	'setting equal to effective radius')
614	roc = R_eff(ro_data%GEOref%lat) ! Maybe +21km? See RSR 14.
615	END IF
616
617	! 16b.5 Generate undulation
618	! -------------------------
619
620	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
621	und = ro_data%GEOref%undulation
622	ELSE
623	CALL message(msg_warn, "Undulation missing ... will assume to " // &
624	"be zero when calculating geopotential heights.")
625	und = 0.0_wp
626	END IF
627
628	! 16b.6 Generate latitude
629	! -----------------------
630
631	IF ( ro_data%GEOref%lat > ropp_MDTV ) THEN
632	lat = ro_data%GEOref%lat
633	ELSE
634	CALL message(msg_warn, "Latitude missing ... will assume to " // &
635	"be zero when calculating geopotential heights.")
636	lat = 0.0_wp
637	END IF
638
639	! 16b.7 Initialise
640	! ----------------
641
642	IF ( (state%temp(1) > ropp_MDTV) .AND. (state%pres(1) > ropp_MDTV) ) THEN
643	refrac_surface = kappa1 * state%pres(1) / state%temp(1) ! roughly
644	ELSE
645	refrac_surface = refrac_surface_default
646	END IF
647
648	obs_refrac%geop = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*refrac_surface)) - &
649	(roc + und) ! roughly
650
651	CALL ropp_fm_refrac_1d(state, obs_refrac, nonull=.TRUE.) ! Initialise with the old interp scheme sans missing data
652
653	! 16b.8 Iterate on GPH
654	! --------------------
655
656	iter = 0
657
658	DO
659
660	tmp = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*obs_refrac%refrac)) - &
661	(roc + und) ! geometric height
662
663	obs_refrac%geop = geometric2geopotential(lat, tmp)
664
665	IF ( state%new_ref_op ) THEN
666	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
667	ELSE
668	CALL ropp_fm_refrac_1d(state, obs_refrac)
669	END IF
670
671	IF ( ANY(obs_refrac%refrac <= ropp_MDTV) ) THEN ! rederive approximately the zapped refracs
672
673	lowest_refrac_index = SUM( MINLOC( (/ (i, i=1,n) /), mask=obs_refrac%refrac>ropp_MDTV ) )
674
675	IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < n) ) THEN
676
677	DO imiss=lowest_refrac_index-1,1,-1
678	obs_refrac%refrac(imiss) = obs_refrac%refrac(imiss+1) + &
679	(obs_refrac%refrac(lowest_refrac_index) - obs_refrac%refrac(lowest_refrac_index+1))
680	END DO
681
682	tmp(1:lowest_refrac_index-1) = ((impact_height(1:lowest_refrac_index-1) + roc) / &
683	(1.0_wp + 1.0e-6_wp*obs_refrac%refrac(1:lowest_refrac_index-1))) - &
684	(roc + und) ! geometric height exactly right in this region
685	obs_refrac%geop(1:lowest_refrac_index-1) = geometric2geopotential(lat, tmp(1:lowest_refrac_index-1))
686
687	END IF
688
689	END IF
690
691	tmp = (1.0_wp + 1.0e-6_wpobs_refrac%refrac)(tmp + roc + und) - &
692	(impact_height + roc) ! residual
693
694	IF ( MAXVAL(ABS(tmp)) < tol ) EXIT
695
696	iter = iter + 1
697
698	IF ( iter > max_iter ) THEN
699	WRITE (smax_iter, '(i3)') max_iter
700	CALL message(msg_fatal, 'GPH iteration failed to converge in ' // &
701	smax_iter // ' iterations')
702	END IF
703
704	WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac
705
706	WHERE ( obs_refrac%geop < min_geop ) obs_refrac%geop = min_geop
707
708	END DO
709
710	WRITE (smax_iter, '(i3)') iter
711	CALL message(msg_diag, 'GPH iteration converged in ' // &
712	smax_iter // ' iterations')
713
714	! 16b.9 Set other elements of obs_refrac structure
715	! ------------------------------------------------
716
717	obs_refrac%refrac = 0.0_wp ! Refrac will be calculated in the main program
718	obs_refrac%weights = 1.0_wp
719	obs_refrac%lon = ro_data%GEOref%lon
720	obs_refrac%lat = ro_data%GEOref%lat
721
722	! 16b.10 Clean up
723	! ---------------
724
725	CALL message_set_routine(routine)
726
727	END SUBROUTINE set_obs_levels_refrac_247
728
729
730	!-------------------------------------------------------------------------------
731	! 17. Calculate obs levels for bending angle (consistent with obs_refrac)
732	!-------------------------------------------------------------------------------
733
734	SUBROUTINE set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle)
735
736	! 17.1 Declarations
737	! -----------------
738
739	USE typesizes, ONLY: wp => EightByteReal
740	USE geodesy
741	USE ropp_io_types
742	USE ropp_fm
743
744	IMPLICIT NONE
745
746	TYPE(ROprof) :: ro_data
747	TYPE(Obs1dRefrac) :: obs_refrac
748	TYPE(Obs1dbangle) :: obs_bangle
749	TYPE(State1dFM) :: state
750
751	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp
752	REAL(wp) :: max_refrac
753	REAL(wp) :: min_geop
754
755	INTEGER :: n, imiss
756	INTEGER :: lowest_refrac_index
757
758	CHARACTER(LEN=256) :: routine
759
760	CALL message_get_routine(routine)
761	CALL message_set_routine('set_obs_levels_bangle')
762
763	! 17.2 Allocate arrays
764	! --------------------
765
766	max_refrac = ro_data%lev2a%range%refrac(2)
767
768	min_geop = ro_data%lev2b%range%geop(1)
769
770	n = SIZE(obs_refrac%geop)
771
772	obs_bangle%n_L1 = n
773
774	ALLOCATE(tmp(n))
775
776	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
777	ALLOCATE(obs_bangle%bangle(n))
778	ALLOCATE(obs_bangle%impact(n))
779	ALLOCATE(obs_bangle%weights(n))
780	ELSE
781	ALLOCATE(obs_bangle%bangle(2*n))
782	ALLOCATE(obs_bangle%impact(2*n))
783	ALLOCATE(obs_bangle%weights(2*n))
784	END IF
785
786	! 17.3 Calculate refractivities
787	! -----------------------------
788
789	IF ( state%new_ref_op ) THEN
790	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
791	ELSE
792	CALL ropp_fm_refrac_1d(state, obs_refrac)
793	END IF
794
795	IF ( ANY(obs_refrac%refrac <= ropp_MDTV) ) THEN ! rederive approximately the zapped refracs
796
797	lowest_refrac_index = SUM( MINLOC( (/ (i, i=1,n) /), mask=obs_refrac%refrac>ropp_MDTV ) )
798
799	IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < n) ) THEN
800
801	DO imiss=lowest_refrac_index-1,1,-1
802	obs_refrac%refrac(imiss) = obs_refrac%refrac(imiss+1) + &
803	(obs_refrac%refrac(lowest_refrac_index) - obs_refrac%refrac(lowest_refrac_index+1))
804	END DO
805
806	END IF
807
808	END IF
809
810	WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac
811
812	WHERE ( obs_refrac%geop < min_geop ) obs_refrac%geop = min_geop
813
814	! 17.4 Set scalar arguments of the observation vector
815	! ---------------------------------------------------
816
817	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
818	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
819	obs_bangle%undulation = ro_data%GEOref%undulation
820	obs_bangle%lon = ro_data%GEOref%lon
821	obs_bangle%lat = ro_data%GEOref%lat
822	obs_bangle%azimuth = ro_data%GEOref%azimuth
823
824	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
825	obs_bangle%r_curve = ro_data%GEOref%roc
826	ELSE
827	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
828	'setting equal to effective radius')
829	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
830	END IF
831
832	! 17.5 Calculate levels to coincide with the geopotential levels
833	! --------------------------------------------------------------
834
835	tmp = geopotential2geometric(ro_data%GEOref%lat, obs_refrac%geop) + &
836	obs_bangle%r_curve
837
838	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
839	tmp = tmp + obs_bangle%undulation
840	ELSE
841	CALL message(msg_warn, "Undulation missing. Will assume to " // &
842	"be zero when calculating impact parameters.")
843	END IF
844
845	tmp = (1.0_wp + obs_refrac%refrac1.e-6_wp) tmp
846
847	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
848	obs_bangle%impact = tmp
849	ELSE
850	obs_bangle%impact(1:n) = tmp
851	obs_bangle%impact(n+1:2*n) = tmp
852	END IF
853
854	! 17.6 Fill other arrays
855	! ----------------------
856
857	obs_bangle%bangle(:) = 0.0_wp
858	obs_bangle%weights(:) = 1.0_wp
859
860	! 17.7 Clean up
861	! -------------
862
863	IF ( state%new_ref_op ) THEN ! Regenerate avec missing data
864	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
865	ELSE
866	CALL ropp_fm_refrac_1d(state, obs_refrac)
867	END IF
868
869	DEALLOCATE(tmp)
870
871	CALL message_set_routine(routine)
872
873	END SUBROUTINE set_obs_levels_bangle
874
875
876	!-------------------------------------------------------------------------------
877	! 17b. Calculate standard 247 observation levels for bending angle
878	!-------------------------------------------------------------------------------
879
880	SUBROUTINE set_obs_levels_bangle_247(ro_data, obs_bangle)
881
882	! 17b.1 Declarations
883	! ------------------
884
885	USE typesizes, ONLY: wp => EightByteReal
886	USE geodesy
887	USE ropp_io_types
888	USE ropp_fm
889
890	IMPLICIT NONE
891
892	TYPE(ROprof) :: ro_data
893	TYPE(Obs1dbangle) :: obs_bangle
894
895	INTEGER, PARAMETER :: n=247
896	REAL(wp), DIMENSION(n) :: impact_height
897
898	CHARACTER(LEN=256) :: routine
899
900	CALL message_get_routine(routine)
901	CALL message_set_routine('set_obs_levels_bangle_247')
902
903	! 17b.2 Allocate arrays
904	! ---------------------
905
906	obs_bangle%n_L1 = n
907
908	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
909	ALLOCATE(obs_bangle%bangle(n))
910	ALLOCATE(obs_bangle%impact(n))
911	ALLOCATE(obs_bangle%weights(n))
912	ELSE
913	ALLOCATE(obs_bangle%bangle(2*n))
914	ALLOCATE(obs_bangle%impact(2*n))
915	ALLOCATE(obs_bangle%weights(2*n))
916	END IF
917
918	! 17b.3 Define standard impact heights
919	! ------------------------------------
920
921	impact_height = impact_height_eum_247 ! The particular 247L set
922
923	! 17b.4 Set scalar arguments of the observation vector
924	! ----------------------------------------------------
925
926	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
927	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
928	obs_bangle%undulation = ro_data%GEOref%undulation
929	obs_bangle%lon = ro_data%GEOref%lon
930	obs_bangle%lat = ro_data%GEOref%lat
931	obs_bangle%azimuth = ro_data%GEOref%azimuth
932
933	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
934	obs_bangle%r_curve = ro_data%GEOref%roc
935	ELSE
936	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
937	'setting equal to effective radius')
938	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
939	END IF
940
941	! 17b.5 Calculate impact parameters
942	! ---------------------------------
943
944	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
945	obs_bangle%impact = impact_height + obs_bangle%r_curve
946	ELSE
947	obs_bangle%impact(1:n) = impact_height + obs_bangle%r_curve
948	obs_bangle%impact(n+1:2*n) = impact_height + obs_bangle%r_curve
949	END IF
950
951	! 17b.6 Fill other arrays
952	! -----------------------
953
954	obs_bangle%bangle(:) = 0.0_wp
955	obs_bangle%weights(:) = 1.0_wp
956
957	! 17b.7 Clean up
958	! --------------
959
960	CALL message_set_routine(routine)
961
962	END SUBROUTINE set_obs_levels_bangle_247
963
964
965	!-------------------------------------------------------------------------------
966	! 18. Usage information
967	!-------------------------------------------------------------------------------
968
969	SUBROUTINE usage()
970	PRINT *, 'Purpose:'
971	PRINT *, ' Bending angles and refractivity forward model'
972	PRINT *, 'Usage:'
973	PRINT *, ' > ropp_fm_bg2ro_1d [<options>] <input_file(s)>'
974	PRINT *, 'Options:'
975	PRINT *, ' -o <output_file> name of ROPP netCDF output file'
976	PRINT *, ' -l <levels_file> optional name of (observation) file'
977	PRINT *, ' (non-multi) containing output level information'
978	PRINT *, ' -f forward model only, no gradient calculation'
979	PRINT *, ' -use_logp use log(pressure) for forward model'
980	PRINT *, ' -use_logq use log(spec humidity) for forward model'
981	PRINT *, ' -comp include non ideal gas compressibility'
982	PRINT *, ' -check_qsat include check against saturation'
983	PRINT *, ' -new_op use alternative refrac and bangle interpolation'
984	PRINT *, ' -direct_ion model L1 and L2 directly'
985	PRINT *, ' -247L output refrac and bangle on standard 247L'
986	PRINT *, ' -zmin <geop_min> minimum refractivity geopotential (gpm)'
987	PRINT *, ' -zmax <geop_max> maximum refractivity geopotential (gpm)'
988	PRINT *, ' -nz <n_geop> number of uniformly spaced refractivity geopotentials (gpm)'
989	PRINT *, ' -d output additional diagnostics'
990	PRINT *, ' -h this help'
991	PRINT *, ' -v version information'
992	PRINT *, ''
993	END SUBROUTINE usage
994
995
996	!-------------------------------------------------------------------------------
997	! 19. Version information
998	!-------------------------------------------------------------------------------
999
1000	SUBROUTINE version_info()
1001	CHARACTER (LEN=40) :: version
1002	version = ropp_fm_version()
1003	PRINT *, 'ropp_fm_bg2ro_1d - Bending angles and refractivity forward model.'
1004	PRINT *, ''
1005	PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
1006	PRINT *, ''
1007	END SUBROUTINE version_info
1008
1009
1010	END PROGRAM ropp_fm_bg2ro_1d

Download in other formats:

Original format