Context navigation

Back to Ticket #314

Ticket #314: ropp_fm_bg2ro_1d.f90.16062015

File ropp_fm_bg2ro_1d.f90.16062015, 35.3 KB (added by Ian Culverwell, 9 years ago)
ropp_fm_bg2ro_1d.f90.16062015

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
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	END IF
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	END IF
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	END IF
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	END IF
333	END IF
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	END IF
351
352	END IF
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	IF (use_247L) THEN
363	CALL set_obs_levels_bangle_247(ro_data, obs_bangle) ! Pre-defined levels
364	ELSE
365	CALL set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle) !default
366	END IF
367	END IF
368
369	!-------------------------------------------------------------------------------
370	! 11. Calculate bending angle and its gradient
371	!-------------------------------------------------------------------------------
372
373	IF (state%state_ok) THEN
374
375	CALL ropp_fm_bangle_1d(state, obs_bangle)
376
377	IF (do_adj .AND. (.NOT. direct_ion)) THEN
378	ALLOCATE(gradient_bangle(SIZE(obs_bangle%bangle),SIZE(state%state)))
379	CALL ropp_fm_bangle_1d_grad(state, obs_bangle, gradient_bangle)
380	END IF
381
382	END IF
383
384	!-------------------------------------------------------------------------------
385	! 12. Copy simulated observations to RO structure and write data
386	!-------------------------------------------------------------------------------
387
388	CALL ropp_fm_obs2roprof(obs_refrac, ro_data)
389
390	CALL ropp_fm_obs2roprof(obs_bangle, ro_data)
391
392	IF (state%state_ok) THEN
393
394	IF (do_adj) THEN
395	CALL ropp_io_addvar(ro_data, &
396	name = "gradient_refrac", &
397	long_name = &
398	"Gradient of the refractivity forward model", &
399	units = "1", &
400	range = (/MINVAL(gradient_refrac), &
401	MAXVAL(gradient_refrac)/), &
402	DATA = gradient_refrac)
403	END IF
404
405	IF (do_adj .AND. (.NOT. direct_ion)) THEN
406	CALL ropp_io_addvar(ro_data, &
407	name = "gradient_bangle" , &
408	long_name = &
409	"Gradient of the bending angle forward model", &
410	units = "rad", &
411	range = (/MINVAL(gradient_bangle), &
412	MAXVAL(gradient_bangle)/), &
413	DATA = gradient_bangle)
414	END IF
415
416	!-------------------------------------------------------------------------------
417	! 13. Update RO structure with computed state vector variables
418	!-------------------------------------------------------------------------------
419
420	CALL ropp_fm_state2roprof(state, ro_data)
421
422	END IF
423
424	!-------------------------------------------------------------------------------
425	! 14. Write data
426	!-------------------------------------------------------------------------------
427
428	IF (direct_ion) CALL ropp_fm_iono_unpack(ro_data)
429
430	CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk)
431
432	!-------------------------------------------------------------------------------
433	! 15. Clean up
434	!-------------------------------------------------------------------------------
435
436	IF (state%state_ok .AND. do_adj) THEN
437	DEALLOCATE (gradient_refrac)
438	END IF
439
440	IF (state%state_ok .AND. do_adj .AND. (.NOT. direct_ion)) THEN
441	DEALLOCATE (gradient_bangle)
442	END IF
443
444	CALL ropp_fm_free(state)
445	CALL ropp_fm_free(obs_refrac)
446	CALL ropp_fm_free(obs_bangle)
447	CALL ropp_io_free(ro_data)
448
449	END DO
450
451	END DO
452
453	CALL EXIT(msg_exit_status)
454
455	CONTAINS
456
457	!-------------------------------------------------------------------------------
458	! 16a. Calculate uniformly spaced observation levels for refractivity
459	!-------------------------------------------------------------------------------
460
461	SUBROUTINE set_obs_levels_refrac(ro_data, obs_refrac, &
462	geop_min, geop_max, n_geop)
463	!
464	! n_geop uniformly spaced geopotential heights between geop_min and geop_max (gpm).
465	!
466
467	! 16a.1 Declarations
468	! -----------------
469
470	USE typesizes, ONLY: wp => EightByteReal
471	USE ropp_io_types
472	USE ropp_fm
473
474	IMPLICIT NONE
475
476	! Input/Output variables
477	TYPE(ROprof) :: ro_data
478	TYPE(Obs1dRefrac) :: obs_refrac
479	REAL(wp), OPTIONAL :: geop_min
480	REAL(wp), OPTIONAL :: geop_max
481	INTEGER, OPTIONAL :: n_geop
482
483	! Local variables
484	REAL(wp) :: zmin, zmax
485	INTEGER :: i, nz
486	CHARACTER(len=256) :: routine
487
488	CALL message_get_routine(routine)
489	CALL message_set_routine('set_obs_levels_refrac')
490
491	! 16a. Vertical geopotential height levels between geop_min and geop_max
492	! ----------------------------------------------------------------------
493
494	zmin = 200.0_wp
495	IF ( PRESENT(geop_min) ) zmin = geop_min
496
497	zmax = 60.0e3_wp
498	IF ( PRESENT(geop_max) ) zmax = geop_max
499
500	nz = 300
501	IF ( PRESENT(n_geop) ) nz = n_geop
502
503	IF ( (zmin >= zmax) .OR. (n_geop <= 1) ) THEN
504	CALL message(msg_warn, 'Invalid geop_min, geop_max or n_geop' // &
505	' ... resetting all to default values')
506	zmin = 200.0_wp
507	zmax = 60000.0_wp
508	nz = 300
509	END IF
510
511	ALLOCATE(obs_refrac%refrac(nz))
512	ALLOCATE(obs_refrac%geop(nz))
513	ALLOCATE(obs_refrac%weights(nz))
514
515	obs_refrac%refrac(:) = 0.0_wp
516	obs_refrac%geop(:) = (/ (zmin+(i-1)*(zmax-zmin)/REAL(nz-1, wp), i=1,nz) /)
517	obs_refrac%weights(:) = 1.0_wp
518	obs_refrac%lon = ro_data%GEOref%lon
519	obs_refrac%lat = ro_data%GEOref%lat
520
521	CALL message_set_routine(routine)
522
523	END SUBROUTINE set_obs_levels_refrac
524
525
526	!-------------------------------------------------------------------------------
527	! 16b. Calculate standard 247 observation levels for refractivity
528	!-------------------------------------------------------------------------------
529
530	SUBROUTINE set_obs_levels_refrac_247(ro_data, obs_refrac, state)
531	!
532	! 247 geopotential height levels, going up to about 60 km, chosen in such a
533	! way that the resulting impact heights match the 'standard' ones in
534	! impact_height_eum_247 in ropp_fm/common/ropp_fm_levels.f90
535	!
536	! 16b.1 Declarations
537	! -----------------
538
539	USE typesizes, ONLY: wp => EightByteReal
540	USE ropp_io_types
541	USE ropp_fm
542
543	IMPLICIT NONE
544
545	TYPE(ROprof) :: ro_data
546	TYPE(Obs1dRefrac) :: obs_refrac
547	TYPE(State1dFM) :: state
548
549	INTEGER, PARAMETER :: n=247
550	REAL(wp), DIMENSION(n) :: impact_height
551	REAL(wp), DIMENSION(n) :: tmp
552	REAL(wp), PARAMETER :: refrac_surface_default=300.0_wp
553	REAL(wp) :: refrac_surface
554	REAL(wp), PARAMETER :: tol=1.0e-3_wp ! Demand < 1mm diffs
555	REAL(wp) :: max_refrac
556	REAL(wp) :: min_geop
557	REAL(wp) :: roc, und, lat
558	INTEGER, PARAMETER :: max_iter=1000
559	INTEGER :: iter, imiss
560	INTEGER :: lowest_refrac_index
561
562	CHARACTER(LEN=3) :: smax_iter
563	CHARACTER(LEN=256) :: routine
564
565	! 16b.2 Initialisations
566	! ---------------------
567
568	CALL message_get_routine(routine)
569	CALL message_set_routine('set_obs_levels_refrac_247')
570
571	impact_height = impact_height_eum_247 ! The particular 247L set
572
573	max_refrac = ro_data%lev2a%range%refrac(2)
574
575	min_geop = ro_data%lev2b%range%geop(1)
576
577	ALLOCATE(obs_refrac%refrac(n))
578	ALLOCATE(obs_refrac%geop(n))
579	ALLOCATE(obs_refrac%weights(n))
580
581	! 16b.3 Vertical geopotential height levels
582	! -----------------------------------------
583
584	! We seek solutions h_i of
585	! (nr = ) (1 + 1e-6 N(Z(h_i)))*(h_i + RoC + und) = a_i + RoC
586	! where a_i are the given impact heights, RoC is the radius of curvature
587	! and und is the undulation. h_i is the geometric height;
588	! Z_i is the geopotential height. Non-linear because N depends on Z.
589	! When geop_refrac is below the lowest state refrac, the refracs are
590	! (now) missing data. In such cases, replace the missing refracs with ones
591	! calculated with the old refrac scheme with this nullification switched off.
592	! Because the same principle is applied in set_obs_levels_bangle, the
593	! resulting impact heights are still good.
594
595	! 16b.4 Generate radius of curvature
596	! ----------------------------------
597
598	IF ( ro_data%GEOref%roc > ropp_MDTV ) THEN
599	roc = ro_data%GEOref%roc
600	ELSE
601	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
602	'setting equal to effective radius')
603	roc = R_eff(ro_data%GEOref%lat) ! Maybe +21km? See RSR 14.
604	END IF
605
606	! 16b.5 Generate undulation
607	! -------------------------
608
609	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
610	und = ro_data%GEOref%undulation
611	ELSE
612	CALL message(msg_warn, "Undulation missing ... will assume to " // &
613	"be zero when calculating geopotential heights.")
614	und = 0.0_wp
615	END IF
616
617	! 16b.6 Generate latitude
618	! -----------------------
619
620	IF ( ro_data%GEOref%lat > ropp_MDTV ) THEN
621	lat = ro_data%GEOref%lat
622	ELSE
623	CALL message(msg_warn, "Latitude missing ... will assume to " // &
624	"be zero when calculating geopotential heights.")
625	lat = 0.0_wp
626	END IF
627
628	! 16b.7 Initialise
629	! ----------------
630
631	IF ( (state%temp(1) > ropp_MDTV) .AND. (state%pres(1) > ropp_MDTV) ) THEN
632	refrac_surface = kappa1 * state%pres(1) / state%temp(1) ! roughly
633	ELSE
634	refrac_surface = refrac_surface_default
635	END IF
636
637	obs_refrac%geop = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*refrac_surface)) - &
638	(roc + und) ! roughly
639
640	CALL ropp_fm_refrac_1d(state, obs_refrac, nonull=.TRUE.) ! Initialise with the old interp scheme sans missing data
641
642	! 16b.8 Iterate on GPH
643	! --------------------
644
645	iter = 0
646
647	DO
648
649	tmp = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*obs_refrac%refrac)) - &
650	(roc + und) ! geometric height
651
652	obs_refrac%geop = geometric2geopotential(lat, tmp)
653
654	IF ( state%new_ref_op ) THEN
655	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
656	ELSE
657	CALL ropp_fm_refrac_1d(state, obs_refrac)
658	END IF
659
660	IF ( ANY(obs_refrac%refrac <= ropp_MDTV) ) THEN ! rederive approximately the zapped refracs
661
662	lowest_refrac_index = SUM( MINLOC( (/ (i, i=1,n) /), mask=obs_refrac%refrac>ropp_MDTV ) )
663
664	IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < n) ) THEN
665
666	DO imiss=lowest_refrac_index-1,1,-1
667	obs_refrac%refrac(imiss) = obs_refrac%refrac(imiss+1) + &
668	(obs_refrac%refrac(lowest_refrac_index) - obs_refrac%refrac(lowest_refrac_index+1))
669	END DO
670
671	tmp(1:lowest_refrac_index-1) = ((impact_height(1:lowest_refrac_index-1) + roc) / &
672	(1.0_wp + 1.0e-6_wp*obs_refrac%refrac(1:lowest_refrac_index-1))) - &
673	(roc + und) ! geometric height exactly right in this region
674	obs_refrac%geop(1:lowest_refrac_index-1) = geometric2geopotential(lat, tmp(1:lowest_refrac_index-1))
675
676	END IF
677
678	END IF
679
680	tmp = (1.0_wp + 1.0e-6_wpobs_refrac%refrac)(tmp + roc + und) - &
681	(impact_height + roc) ! residual
682
683	print*,'iter, MAXVAL(ABS(res)) = ', iter, MAXVAL(ABS(tmp))
684	IF ( MAXVAL(ABS(tmp)) < tol ) EXIT
685
686	iter = iter + 1
687
688	IF ( iter > max_iter ) THEN
689	WRITE (smax_iter, '(i3)') max_iter
690	CALL message(msg_fatal, 'GPH iteration failed to converge in ' // &
691	smax_iter // ' iterations')
692	END IF
693
694	WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac
695
696	WHERE ( obs_refrac%geop < min_geop ) obs_refrac%geop = min_geop
697
698	END DO
699
700	WRITE (smax_iter, '(i3)') iter
701	CALL message(msg_diag, 'GPH iteration converged in ' // &
702	smax_iter // ' iterations')
703
704	! 16b.9 Set other elements of obs_refrac structure
705	! ------------------------------------------------
706
707	obs_refrac%refrac = 0.0_wp ! Refrac will be calculated in the main program
708	obs_refrac%weights = 1.0_wp
709	obs_refrac%lon = ro_data%GEOref%lon
710	obs_refrac%lat = ro_data%GEOref%lat
711
712	! 16b.10 Clean up
713	! ---------------
714
715	CALL message_set_routine(routine)
716
717	END SUBROUTINE set_obs_levels_refrac_247
718
719
720	!-------------------------------------------------------------------------------
721	! 17. Calculate obs levels for bending angle (consistent with obs_refrac)
722	!-------------------------------------------------------------------------------
723
724	SUBROUTINE set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle)
725
726	! 17.1 Declarations
727	! -----------------
728
729	USE typesizes, ONLY: wp => EightByteReal
730	USE geodesy
731	USE ropp_io_types
732	USE ropp_fm
733
734	IMPLICIT NONE
735
736	TYPE(ROprof) :: ro_data
737	TYPE(Obs1dRefrac) :: obs_refrac
738	TYPE(Obs1dbangle) :: obs_bangle
739	TYPE(State1dFM) :: state
740
741	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp
742	REAL(wp) :: max_refrac
743	REAL(wp) :: min_geop
744
745	INTEGER :: n, imiss
746	INTEGER :: lowest_refrac_index
747
748	CHARACTER(LEN=256) :: routine
749
750	CALL message_get_routine(routine)
751	CALL message_set_routine('set_obs_levels_bangle')
752
753	! 17.2 Allocate arrays
754	! --------------------
755
756	max_refrac = ro_data%lev2a%range%refrac(2)
757
758	min_geop = ro_data%lev2b%range%geop(1)
759
760	n = SIZE(obs_refrac%geop)
761
762	obs_bangle%n_L1 = n
763
764	ALLOCATE(tmp(n))
765
766	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
767	ALLOCATE(obs_bangle%bangle(n))
768	ALLOCATE(obs_bangle%impact(n))
769	ALLOCATE(obs_bangle%weights(n))
770	ELSE
771	ALLOCATE(obs_bangle%bangle(2*n))
772	ALLOCATE(obs_bangle%impact(2*n))
773	ALLOCATE(obs_bangle%weights(2*n))
774	END IF
775
776	! 17.3 Calculate refractivities
777	! -----------------------------
778
779	IF ( state%new_ref_op ) THEN
780	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
781	ELSE
782	CALL ropp_fm_refrac_1d(state, obs_refrac)
783	END IF
784
785	IF ( ANY(obs_refrac%refrac <= ropp_MDTV) ) THEN ! rederive approximately the zapped refracs
786
787	lowest_refrac_index = SUM( MINLOC( (/ (i, i=1,n) /), mask=obs_refrac%refrac>ropp_MDTV ) )
788
789	IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < n) ) THEN
790
791	DO imiss=lowest_refrac_index-1,1,-1
792	obs_refrac%refrac(imiss) = obs_refrac%refrac(imiss+1) + &
793	(obs_refrac%refrac(lowest_refrac_index) - obs_refrac%refrac(lowest_refrac_index+1))
794	END DO
795
796	END IF
797
798	END IF
799
800	WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac
801
802	WHERE ( obs_refrac%geop < min_geop ) obs_refrac%geop = min_geop
803
804	! 17.4 Set scalar arguments of the observation vector
805	! ---------------------------------------------------
806
807	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
808	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
809	obs_bangle%undulation = ro_data%GEOref%undulation
810	obs_bangle%lon = ro_data%GEOref%lon
811	obs_bangle%lat = ro_data%GEOref%lat
812	obs_bangle%azimuth = ro_data%GEOref%azimuth
813
814	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
815	obs_bangle%r_curve = ro_data%GEOref%roc
816	ELSE
817	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
818	'setting equal to effective radius')
819	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
820	END IF
821
822	! 17.5 Calculate levels to coincide with the geopotential levels
823	! --------------------------------------------------------------
824
825	tmp = geopotential2geometric(ro_data%GEOref%lat, obs_refrac%geop) + &
826	obs_bangle%r_curve
827
828	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
829	tmp = tmp + obs_bangle%undulation
830	ELSE
831	CALL message(msg_warn, "Undulation missing. Will assume to " // &
832	"be zero when calculating impact parameters.")
833	END IF
834
835	tmp = (1.0_wp + obs_refrac%refrac1.e-6_wp) tmp
836
837	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
838	obs_bangle%impact = tmp
839	ELSE
840	obs_bangle%impact(1:n) = tmp
841	obs_bangle%impact(n+1:2*n) = tmp
842	END IF
843
844	! 17.6 Fill other arrays
845	! ----------------------
846
847	obs_bangle%bangle(:) = 0.0_wp
848	obs_bangle%weights(:) = 1.0_wp
849
850	! 17.7 Clean up
851	! -------------
852
853	IF ( state%new_ref_op ) THEN ! Regenerate avec missing data
854	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
855	ELSE
856	CALL ropp_fm_refrac_1d(state, obs_refrac)
857	END IF
858
859	DEALLOCATE(tmp)
860
861	CALL message_set_routine(routine)
862
863	END SUBROUTINE set_obs_levels_bangle
864
865
866	!-------------------------------------------------------------------------------
867	! 17b. Calculate standard 247 observation levels for bending angle
868	!-------------------------------------------------------------------------------
869
870	SUBROUTINE set_obs_levels_bangle_247(ro_data, obs_bangle)
871
872	! 17b.1 Declarations
873	! ------------------
874
875	USE typesizes, ONLY: wp => EightByteReal
876	USE geodesy
877	USE ropp_io_types
878	USE ropp_fm
879
880	IMPLICIT NONE
881
882	TYPE(ROprof) :: ro_data
883	TYPE(Obs1dbangle) :: obs_bangle
884
885	INTEGER, PARAMETER :: n=247
886	REAL(wp), DIMENSION(n) :: impact_height
887
888	CHARACTER(LEN=256) :: routine
889
890	CALL message_get_routine(routine)
891	CALL message_set_routine('set_obs_levels_bangle_247')
892
893	! 17b.2 Allocate arrays
894	! ---------------------
895
896	obs_bangle%n_L1 = n
897
898	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
899	ALLOCATE(obs_bangle%bangle(n))
900	ALLOCATE(obs_bangle%impact(n))
901	ALLOCATE(obs_bangle%weights(n))
902	ELSE
903	ALLOCATE(obs_bangle%bangle(2*n))
904	ALLOCATE(obs_bangle%impact(2*n))
905	ALLOCATE(obs_bangle%weights(2*n))
906	END IF
907
908	! 17b.3 Define standard impact heights
909	! ------------------------------------
910
911	impact_height = impact_height_eum_247 ! The particular 247L set
912
913	! 17b.4 Set scalar arguments of the observation vector
914	! ----------------------------------------------------
915
916	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
917	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
918	obs_bangle%undulation = ro_data%GEOref%undulation
919	obs_bangle%lon = ro_data%GEOref%lon
920	obs_bangle%lat = ro_data%GEOref%lat
921	obs_bangle%azimuth = ro_data%GEOref%azimuth
922
923	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
924	obs_bangle%r_curve = ro_data%GEOref%roc
925	ELSE
926	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
927	'setting equal to effective radius')
928	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
929	END IF
930
931	! 17b.5 Calculate impact parameters
932	! ---------------------------------
933
934	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
935	obs_bangle%impact = impact_height + obs_bangle%r_curve
936	ELSE
937	obs_bangle%impact(1:n) = impact_height + obs_bangle%r_curve
938	obs_bangle%impact(n+1:2*n) = impact_height + obs_bangle%r_curve
939	END IF
940
941	! 17b.6 Fill other arrays
942	! -----------------------
943
944	obs_bangle%bangle(:) = 0.0_wp
945	obs_bangle%weights(:) = 1.0_wp
946
947	! 17b.7 Clean up
948	! --------------
949
950	CALL message_set_routine(routine)
951
952	END SUBROUTINE set_obs_levels_bangle_247
953
954
955	!-------------------------------------------------------------------------------
956	! 18. Usage information
957	!-------------------------------------------------------------------------------
958
959	SUBROUTINE usage()
960	PRINT *, 'Purpose:'
961	PRINT *, ' Bending angles and refractivity forward model'
962	PRINT *, 'Usage:'
963	PRINT *, ' > ropp_fm_bg2ro_1d [<options>] <input_file(s)>'
964	PRINT *, 'Options:'
965	PRINT *, ' -o <output_file> name of ROPP netCDF output file'
966	PRINT *, ' -l <levels_file> optional name of (observation) file'
967	PRINT *, ' (non-multi) containing output level information'
968	PRINT *, ' -f forward model only, no gradient calculation'
969	PRINT *, ' -use_logp use log(pressure) for forward model'
970	PRINT *, ' -use_logq use log(spec humidity) for forward model'
971	PRINT *, ' -comp include non ideal gas compressibility'
972	PRINT *, ' -check_qsat include check against saturation'
973	PRINT *, ' -new_op use alternative refrac and bangle interpolation'
974	PRINT *, ' -direct_ion model L1 and L2 directly'
975	PRINT *, ' -247L output refrac and bangle on standard 247L'
976	PRINT *, ' -zmin <geop_min> minimum refractivity geopotential (gpm)'
977	PRINT *, ' -zmax <geop_max> maximum refractivity geopotential (gpm)'
978	PRINT *, ' -nz <n_geop> number of uniformly spaced refractivity geopotentials (gpm)'
979	PRINT *, ' -d output additional diagnostics'
980	PRINT *, ' -h this help'
981	PRINT *, ' -v version information'
982	PRINT *, ''
983	END SUBROUTINE usage
984
985
986	!-------------------------------------------------------------------------------
987	! 19. Version information
988	!-------------------------------------------------------------------------------
989
990	SUBROUTINE version_info()
991	CHARACTER (LEN=40) :: version
992	version = ropp_fm_version()
993	PRINT *, 'ropp_fm_bg2ro_1d - Bending angles and refractivity forward model.'
994	PRINT *, ''
995	PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
996	PRINT *, ''
997	END SUBROUTINE version_info
998
999
1000	END PROGRAM ropp_fm_bg2ro_1d

Download in other formats:

Original format