Context navigation

Back to Ticket #314

Ticket #314: ropp_fm_bg2ro_1d.f90_14042015

File ropp_fm_bg2ro_1d.f90_14042015, 34.4 KB (added by Ian Culverwell, 10 years ago)
ropp_fm_bg2ro_1d.f90_14042015

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	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	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	ENDIF
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	ENDIF
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	ENDIF
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	ENDIF
439
440	IF (state%state_ok .AND. do_adj .AND. (.NOT. direct_ion)) THEN
441	DEALLOCATE (gradient_bangle)
442	ENDIF
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, obs_refrac2
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=200
559	INTEGER :: iter, imin
560	CHARACTER(LEN=3) :: smax_iter
561	CHARACTER(LEN=256) :: routine
562
563	! 16b.2 Initialisations
564	! ---------------------
565
566	CALL message_get_routine(routine)
567	CALL message_set_routine('set_obs_levels_refrac_247')
568
569	impact_height = impact_height_eum_247 ! The particular 247L set
570
571	max_refrac = ro_data%lev2a%range%refrac(2)
572
573	min_geop = ro_data%lev2b%range%geop(1)
574
575	ALLOCATE(obs_refrac%refrac(n))
576	ALLOCATE(obs_refrac%geop(n))
577	ALLOCATE(obs_refrac%weights(n))
578
579	ALLOCATE(obs_refrac2%refrac(n))
580	ALLOCATE(obs_refrac2%geop(n))
581	ALLOCATE(obs_refrac2%weights(n))
582
583	! 16b.3 Vertical geopotential height levels
584	! -----------------------------------------
585
586	! We seek solutions h_i of
587	! (nr = ) (1 + 1e-6 N(Z(h_i)))*(h_i + RoC + und) = a_i + RoC
588	! where a_i are the given impact heights, RoC is the radius of curvature
589	! and und is the undulation. h_i is the geometric height;
590	! Z_i is the geopotential height. Non-linear because N depends on Z.
591	! When geop_refrac is below the lowest state refrac, the refracs are
592	! (now) missing data. In such cases, replace the missing refracs with ones
593	! calculated with the old refrac scheme with this nullification switched off.
594	! Because the same principle is applied in set_obs_levels_bangle, the
595	! resulting impact heights are still good.
596
597	! 16b.4 Generate radius of curvature
598	! ----------------------------------
599
600	IF ( ro_data%GEOref%roc > ropp_MDTV ) THEN
601	roc = ro_data%GEOref%roc
602	ELSE
603	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
604	'setting equal to effective radius')
605	roc = R_eff(ro_data%GEOref%lat) ! Maybe +21km? See RSR 14.
606	ENDIF
607
608	! 16b.5 Generate undulation
609	! -------------------------
610
611	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
612	und = ro_data%GEOref%undulation
613	ELSE
614	CALL message(msg_warn, "Undulation missing ... will assume to " // &
615	"be zero when calculating geopotential heights.")
616	und = 0.0_wp
617	ENDIF
618
619	! 16b.6 Generate latitude
620	! -----------------------
621
622	IF ( ro_data%GEOref%lat > ropp_MDTV ) THEN
623	lat = ro_data%GEOref%lat
624	ELSE
625	CALL message(msg_warn, "Latitude missing ... will assume to " // &
626	"be zero when calculating geopotential heights.")
627	lat = 0.0_wp
628	END IF
629
630	! 16b.7 Initialise
631	! ----------------
632
633	IF ( (state%temp(1) > ropp_MDTV) .AND. (state%pres(1) > ropp_MDTV) ) THEN
634	refrac_surface = kappa1 * state%pres(1) / state%temp(1) ! roughly
635	ELSE
636	refrac_surface = refrac_surface_default
637	END IF
638
639	obs_refrac%geop = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*refrac_surface)) - &
640	(roc + und) ! roughly
641
642	CALL ropp_fm_refrac_1d(state, obs_refrac, nonull=.TRUE.) ! Initialise with the old interp scheme sans missing data
643
644	! 16b.8 Iterate on GPH
645	! --------------------
646
647	iter = 0
648
649	DO
650
651	tmp = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*obs_refrac%refrac)) - &
652	(roc + und) ! geometric height
653
654	obs_refrac%geop = geometric2geopotential(lat, tmp)
655
656	IF ( state%new_ref_op ) THEN
657	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
658	ELSE
659	CALL ropp_fm_refrac_1d(state, obs_refrac)
660	END IF
661
662	IF ( COUNT(obs_refrac%refrac <= ropp_MDTV) > 0 ) THEN ! rederive approximately some zapped refracs
663
664	obs_refrac2 = obs_refrac
665
666	CALL ropp_fm_refrac_1d(state, obs_refrac2, nonull=.TRUE.) ! Recalculate sans missing data
667
668	WHERE (obs_refrac%refrac <= ropp_MDTV) &
669	obs_refrac%refrac = obs_refrac2%refrac
670
671	END IF
672
673	WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac
674
675	WHERE ( obs_refrac%geop < min_geop ) obs_refrac%geop = min_geop
676
677	tmp = (1.0_wp + 1.0e-6_wpobs_refrac%refrac)(tmp + roc + und) - &
678	(impact_height + roc) ! residual
679
680	IF ( MAXVAL(ABS(tmp)) < tol ) EXIT
681
682	iter = iter + 1
683
684	IF ( iter > max_iter ) THEN
685	WRITE (smax_iter, '(i3)') max_iter
686	CALL message(msg_fatal, 'GPH iteration failed to converge in ' // &
687	smax_iter // ' iterations')
688	ENDIF
689
690	ENDDO
691
692	WRITE (smax_iter, '(i3)') iter
693	CALL message(msg_diag, 'GPH iteration converged in ' // &
694	smax_iter // ' iterations')
695
696	! 16b.9 Set other elements of obs_refrac structure
697	! ------------------------------------------------
698
699	obs_refrac%refrac = 0.0_wp ! Refrac will be calculated in the main program
700	obs_refrac%weights = 1.0_wp
701	obs_refrac%lon = ro_data%GEOref%lon
702	obs_refrac%lat = ro_data%GEOref%lat
703
704	! 16b.10 Clean up
705	! ---------------
706
707	CALL ropp_fm_free(obs_refrac2)
708
709	CALL message_set_routine(routine)
710
711	END SUBROUTINE set_obs_levels_refrac_247
712
713
714	!-------------------------------------------------------------------------------
715	! 17. Calculate obs levels for bending angle (consistent with obs_refrac)
716	!-------------------------------------------------------------------------------
717
718	SUBROUTINE set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle)
719
720	! 17.1 Declarations
721	! -----------------
722
723	USE typesizes, ONLY: wp => EightByteReal
724	USE geodesy
725	USE ropp_io_types
726	USE ropp_fm
727
728	IMPLICIT NONE
729
730	TYPE(ROprof) :: ro_data
731	TYPE(Obs1dRefrac) :: obs_refrac, obs_refrac2
732	TYPE(Obs1dbangle) :: obs_bangle
733	TYPE(State1dFM) :: state
734
735	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp
736	REAL(wp) :: max_refrac
737
738	INTEGER :: n, imin
739
740	CHARACTER(LEN=256) :: routine
741
742	CALL message_get_routine(routine)
743	CALL message_set_routine('set_obs_levels_bangle')
744
745	! 17.2 Allocate arrays
746	! --------------------
747
748	max_refrac = ro_data%lev2a%range%refrac(2)
749
750	n = SIZE(obs_refrac%geop)
751
752	obs_bangle%n_L1 = n
753
754	ALLOCATE(tmp(n))
755
756	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
757	ALLOCATE(obs_bangle%bangle(n))
758	ALLOCATE(obs_bangle%impact(n))
759	ALLOCATE(obs_bangle%weights(n))
760	ELSE
761	ALLOCATE(obs_bangle%bangle(2*n))
762	ALLOCATE(obs_bangle%impact(2*n))
763	ALLOCATE(obs_bangle%weights(2*n))
764	ENDIF
765
766	ALLOCATE(obs_refrac2%refrac(n))
767	ALLOCATE(obs_refrac2%geop(n))
768	ALLOCATE(obs_refrac2%weights(n))
769
770	! 17.3 Calculate refractivities
771	! -----------------------------
772
773	IF ( state%new_ref_op ) THEN
774	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
775	ELSE
776	CALL ropp_fm_refrac_1d(state, obs_refrac)
777	END IF
778
779	IF ( COUNT(obs_refrac%refrac <= ropp_MDTV) > 0 ) THEN ! rederive approximately some zapped refracs
780
781	obs_refrac2 = obs_refrac
782
783	CALL ropp_fm_refrac_1d(state, obs_refrac2, nonull=.TRUE.) ! Estimate sans missing data
784
785	WHERE (obs_refrac%refrac <= ropp_MDTV) &
786	obs_refrac%refrac = obs_refrac2%refrac
787
788	END IF
789
790	WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac
791
792	! 17.4 Set scalar arguments of the observation vector
793	! ---------------------------------------------------
794
795	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
796	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
797	obs_bangle%undulation = ro_data%GEOref%undulation
798	obs_bangle%lon = ro_data%GEOref%lon
799	obs_bangle%lat = ro_data%GEOref%lat
800	obs_bangle%azimuth = ro_data%GEOref%azimuth
801
802	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
803	obs_bangle%r_curve = ro_data%GEOref%roc
804	ELSE
805	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
806	'setting equal to effective radius')
807	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
808	ENDIF
809
810	! 17.5 Calculate levels to coincide with the geopotential levels
811	! --------------------------------------------------------------
812
813	tmp = geopotential2geometric(ro_data%GEOref%lat, obs_refrac%geop) + &
814	obs_bangle%r_curve
815
816	IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
817	tmp = tmp + obs_bangle%undulation
818	ELSE
819	CALL message(msg_warn, "Undulation missing. Will assume to " // &
820	"be zero when calculating impact parameters.")
821	ENDIF
822
823	tmp = (1.0_wp + obs_refrac%refrac1.e-6_wp) tmp
824
825	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
826	obs_bangle%impact = tmp
827	ELSE
828	obs_bangle%impact(1:n) = tmp
829	obs_bangle%impact(n+1:2*n) = tmp
830	ENDIF
831
832	! 17.6 Fill other arrays
833	! ----------------------
834
835	obs_bangle%bangle(:) = 0.0_wp
836	obs_bangle%weights(:) = 1.0_wp
837
838	! 17.7 Clean up
839	! -------------
840
841	IF ( state%new_ref_op ) THEN
842	CALL ropp_fm_refrac_1d_new(state, obs_refrac)
843	ELSE
844	CALL ropp_fm_refrac_1d(state, obs_refrac) ! Regenerate avec missing data
845	END IF
846
847	DEALLOCATE(tmp)
848
849	CALL ropp_fm_free(obs_refrac2)
850
851	CALL message_set_routine(routine)
852
853	END SUBROUTINE set_obs_levels_bangle
854
855	!-------------------------------------------------------------------------------
856	! 17b. Calculate standard 247 observation levels for bending angle
857	!-------------------------------------------------------------------------------
858
859	SUBROUTINE set_obs_levels_bangle_247(ro_data, obs_bangle)
860
861	! 17b.1 Declarations
862	! ------------------
863
864	USE typesizes, ONLY: wp => EightByteReal
865	USE geodesy
866	USE ropp_io_types
867	USE ropp_fm
868
869	IMPLICIT NONE
870
871	TYPE(ROprof) :: ro_data
872	TYPE(Obs1dbangle) :: obs_bangle
873
874	INTEGER, PARAMETER :: n=247
875	REAL(wp), DIMENSION(n) :: impact_height
876
877	CHARACTER(LEN=256) :: routine
878
879	CALL message_get_routine(routine)
880	CALL message_set_routine('set_obs_levels_bangle_247')
881
882	! 17b.1 Allocate arrays
883	! ---------------------
884
885	obs_bangle%n_L1 = n
886
887	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
888	ALLOCATE(obs_bangle%bangle(n))
889	ALLOCATE(obs_bangle%impact(n))
890	ALLOCATE(obs_bangle%weights(n))
891	ELSE
892	ALLOCATE(obs_bangle%bangle(2*n))
893	ALLOCATE(obs_bangle%impact(2*n))
894	ALLOCATE(obs_bangle%weights(2*n))
895	ENDIF
896
897	! 17b.2 Define standard impact heights
898	! ------------------------------------
899
900	impact_height = impact_height_eum_247 ! The particular 247L set
901
902	! 17b.3 Set scalar arguments of the observation vector
903	! ----------------------------------------------------
904
905	obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
906	obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
907	obs_bangle%undulation = ro_data%GEOref%undulation
908	obs_bangle%lon = ro_data%GEOref%lon
909	obs_bangle%lat = ro_data%GEOref%lat
910	obs_bangle%azimuth = ro_data%GEOref%azimuth
911
912	IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
913	obs_bangle%r_curve = ro_data%GEOref%roc
914	ELSE
915	CALL message (msg_warn, 'RoC missing from data structure ... ' // &
916	'setting equal to effective radius')
917	obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
918	ENDIF
919
920	! 17b.4 Calculate impact parameters
921	! ---------------------------------
922
923	IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
924
925	obs_bangle%impact = impact_height + obs_bangle%r_curve
926
927	ELSE
928
929	obs_bangle%impact(1:n) = impact_height + obs_bangle%r_curve
930	obs_bangle%impact(n+1:2*n) = impact_height + obs_bangle%r_curve
931
932	END IF
933
934	! 17b.5 Fill other arrays
935	! -----------------------
936
937	obs_bangle%bangle(:) = 0.0_wp
938	obs_bangle%weights(:) = 1.0_wp
939
940	! 17b.6 Clean up
941	! --------------
942
943	CALL message_set_routine(routine)
944
945	END SUBROUTINE set_obs_levels_bangle_247
946
947	!-------------------------------------------------------------------------------
948	! 18. Usage information
949	!-------------------------------------------------------------------------------
950
951	SUBROUTINE usage()
952	PRINT *, 'Purpose:'
953	PRINT *, ' Bending angles and refractivity forward model'
954	PRINT *, 'Usage:'
955	PRINT *, ' > ropp_fm_bg2ro_1d [<options>] <input_file(s)>'
956	PRINT *, 'Options:'
957	PRINT *, ' -o <output_file> name of ROPP netCDF output file'
958	PRINT *, ' -l <levels_file> optional name of (observation) file'
959	PRINT *, ' (non-multi) containing output level information'
960	PRINT *, ' -f forward model only, no gradient calculation'
961	PRINT *, ' -use_logp use log(pressure) for forward model'
962	PRINT *, ' -use_logq use log(spec humidity) for forward model'
963	PRINT *, ' -comp include non ideal gas compressibility'
964	PRINT *, ' -check_qsat include check against saturation'
965	PRINT *, ' -new_op use alternative refrac and bangle interpolation'
966	PRINT *, ' -direct_ion model L1 and L2 directly'
967	PRINT *, ' -247L output refrac and bangle on standard 247L'
968	PRINT *, ' -zmin <geop_min> minimum refractivity geopotential (gpm)'
969	PRINT *, ' -zmax <geop_max> maximum refractivity geopotential (gpm)'
970	PRINT *, ' -nz <n_geop> number of uniformly spaced refractivity geopotentials (gpm)'
971	PRINT *, ' -d output additional diagnostics'
972	PRINT *, ' -h this help'
973	PRINT *, ' -v version information'
974	PRINT *, ''
975	END SUBROUTINE usage
976
977	!-------------------------------------------------------------------------------
978	! 19. Version information
979	!-------------------------------------------------------------------------------
980
981	SUBROUTINE version_info()
982	CHARACTER (LEN=40) :: version
983	version = ropp_fm_version()
984	PRINT *, 'ropp_fm_bg2ro_1d - Bending angles and refractivity forward model.'
985	PRINT *, ''
986	PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
987	PRINT *, ''
988	END SUBROUTINE version_info
989
990	END PROGRAM ropp_fm_bg2ro_1d

Download in other formats:

Original format