Ticket #264: ropp_fm_bg2ro_1d.f90

File ropp_fm_bg2ro_1d.f90, 16.5 KB (added by kmk, 14 years ago)

1D forward-model tool with -l flag added for specifying output levels from a file

Line 
1! $Id: ropp_fm_bg2ro_1d.f90 3020 2011-09-14 08:56:19Z frdo $
2
3PROGRAM 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> [-f] [-use_logp]
15! [-use_logq] [-comp] [-d] [-h] [-v]
16!
17! ARGUMENTS
18! <infile(s)> One (or more) input file names.
19!
20! OPTIONS
21! -o <outfile> Name of output file (default: bg2ro.nc)
22! -f forward only, no adjoint calculation.
23! -use_logp use log(pressure) for forward model
24! -use_logq use log(spec humidity) for forward model
25! -comp include non ideal gas compressibility
26! -d output additional diagnostics
27! -h help
28! -v version information
29!
30! DESCRIPTION
31! This program reads model data on model levels from the input data files
32! and calculates vertical profiles of bending angle and refractivity using
33! the 1d forward operators. The result is written to an ROPP formatted
34! output file.
35!
36! NOTES
37! If the input file is a multifile, or more than one input files are
38! specified, the output file is a multifile.
39!
40! Already existing output files will be overwritten.
41!
42! EXAMPLE
43! To calculate bending angle and refractivity from one of the example
44! (single-) files in the data directory:
45!
46! > ropp_fm_bg2ro_1d ../data/bgr20090401_000329_M02_2030337800_N0007_XXXX.nc
47!
48! To calculate bending angle and refractivity profiles from all singlefiles
49! in the data directory:
50!
51! > ropp_fm_bg2ro_1d ../data/bgr20090401*_N0007_XXXX.nc -o eg_02.nc
52!
53! Note that the resulting eg_02.nc file contains forward modelled data from
54! all example profiles.
55!
56! To calculate forward modelled bending angle and refractivity profiles from
57! all profiles contained in the multifile bgr20090401_multi.nc:
58!
59! > ropp_fm_bg2ro_1d ../data/bgr20090401_multi.nc -o eg_03.nc
60!
61! Since the ecmwf_multi_* file was generated by concatenating the other
62! files in the data directory, eg_02.nc and eg_03.nc should be identical
63! apart from the file names.
64!
65! SEE ALSO
66! ropp_fm_bangle_1d
67! ropp_fm_refrac_1d
68!
69! AUTHOR
70! Met Office, Exeter, UK.
71! Any comments on this software should be given via the GRAS SAF
72! Helpdesk at http://www.grassaf.org
73!
74! COPYRIGHT
75! (c) EUMETSAT. All rights reserved.
76! For further details please refer to the file COPYRIGHT
77! which you should have received as part of this distribution.
78!
79!****
80
81!-------------------------------------------------------------------------------
82! 1. Declarations
83!-------------------------------------------------------------------------------
84
85 USE typesizes, ONLY: wp => EightByteReal
86 USE ropp_utils
87 USE ropp_io
88 USE ropp_io_types, ONLY: ROprof
89 USE ropp_fm
90 USE ropp_fm_types
91 USE ropp_fm_copy
92
93 IMPLICIT NONE
94
95 TYPE(ROprof) :: ro_data
96 TYPE(State1dFM) :: state
97 TYPE(Obs1dRefrac) :: obs_refrac
98 TYPE(Obs1dBangle) :: obs_bangle
99
100 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gradient_bangle
101 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gradient_refrac
102
103 INTEGER :: idummy
104 INTEGER :: i, iargc, argc, k
105 INTEGER :: n_files, n_profiles
106
107 LOGICAL :: give_help, do_adj
108 LOGICAL :: ranchk = .TRUE.
109 LOGICAL :: compress = .FALSE.
110
111 CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
112 CHARACTER(len = 4096) :: ofile
113 CHARACTER(len = 256) :: buffer
114 CHARACTER(len = 4) :: istr
115 CHARACTER(len = 6) :: nstr
116
117!-------------------------------------------------------------------------------
118! 2. Default settings
119!-------------------------------------------------------------------------------
120
121 do_adj = .TRUE.
122 give_help = .FALSE.
123 ofile = "bg2ro.nc"
124
125 CALL message(msg_noin, '')
126 CALL message(msg_noin, &
127 '-----------------------------------------------------------------------')
128 CALL message(msg_noin, &
129 ' ROPP Forward Model' )
130 CALL message(msg_noin, &
131 '-----------------------------------------------------------------------')
132 CALL message(msg_noin, '')
133
134!-------------------------------------------------------------------------------
135! 3. Command line arguments
136!-------------------------------------------------------------------------------
137
138 argc = iargc()
139 i = 1
140 n_files = 0
141 ALLOCATE(ifiles(argc))
142
143 DO WHILE(i <= argc)
144 CALL getarg(i, buffer)
145 SELECT CASE (buffer)
146 CASE('-o') ! Output file name (netCDF output)
147 CALL getarg(i+1, buffer)
148 ofile = buffer
149 i = i + 1
150 CASE ('-f') ! Perform forward model only
151 do_adj = .FALSE.
152 CASE ('-comp')
153 compress = .TRUE. ! non ideal gas switch
154 CASE ('-use_logp') ! Use log(pressure) for FM
155 state%use_logp = .TRUE.
156 CASE ('-use_logq') ! Use log(shum) for FM
157 state%use_logq = .TRUE.
158 CASE ('-no-ranchk') ! Use no rangecheck on output
159 ranchk = .FALSE.
160 CASE('-d') ! Additional diagnostic mode
161 msg_MODE = VerboseMode
162 CASE('-h', '--help', '?') ! help
163 give_help = .TRUE.
164 CASE('-v', '-V', '--version') ! version information
165 CALL version_info()
166 CALL EXIT(0)
167 CASE default ! Input file name
168 IF ( buffer(1:1) /= '-' ) THEN
169 n_files = n_files + 1
170 ifiles(n_files) = buffer
171 END IF
172 END SELECT
173 i = i + 1
174 END DO
175
176 IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN
177 CALL usage()
178 CALL EXIT(0)
179 ENDIF
180
181!-------------------------------------------------------------------------------
182! 4. Remove pre-existing output file
183!-------------------------------------------------------------------------------
184
185 CALL file_delete(ofile, idummy)
186
187!-------------------------------------------------------------------------------
188! 4. Set default units
189!-------------------------------------------------------------------------------
190
191 CALL ropp_fm_set_units(ro_data)
192
193!-------------------------------------------------------------------------------
194! 5. Loop over all input files
195!-------------------------------------------------------------------------------
196
197 DO k = 1, n_files
198
199!-------------------------------------------------------------------------------
200! 6. Loop over all profiles
201!-------------------------------------------------------------------------------
202
203 n_profiles = ropp_io_nrec(ifiles(k))
204
205 DO i = 1, n_profiles
206
207 WRITE(istr, '(i4)') i
208 WRITE(nstr, '(i6)') n_profiles
209 CALL message(msg_noin, '')
210 CALL message(msg_info, "Processing profile " // istr // " of " // nstr )
211
212!-------------------------------------------------------------------------------
213! 7. Read data
214!-------------------------------------------------------------------------------
215
216 CALL ropp_io_read(ro_data, ifiles(k), rec = i, ranchk = ranchk)
217 CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")
218
219!-------------------------------------------------------------------------------
220! 8. Copy data in RO structure to state and refrac obs vectors
221!-------------------------------------------------------------------------------
222
223 CALL ropp_fm_roprof2state(ro_data, state)
224
225! switch non-ideal gas on
226
227 IF (compress) state%non_ideal = .TRUE.
228
229 IF (ro_data%lev2a%Npoints > 0) THEN
230 CALL ropp_fm_roprof2obs(ro_data, obs_refrac) !!set to obs levels
231 obs_refrac%refrac(:) = ropp_MDFV
232 ELSE
233 CALL set_obs_levels_refrac(ro_data, obs_refrac) !!pre-defined levels
234 ENDIF
235
236!-------------------------------------------------------------------------------
237! 9. Calculate refractivity and its gradient
238!-------------------------------------------------------------------------------
239
240 IF (state%state_ok) THEN
241
242 CALL ropp_fm_refrac_1d(state, obs_refrac)
243
244 IF (do_adj) THEN
245 ALLOCATE(gradient_refrac(SIZE(obs_refrac%refrac),SIZE(state%state)))
246 CALL ropp_fm_refrac_1d_grad(state, obs_refrac, gradient_refrac)
247 ENDIF
248
249 ENDIF
250
251!-------------------------------------------------------------------------------
252! 10. Copy data in RO and refrac structure to bending angle obs vector
253!-------------------------------------------------------------------------------
254
255 IF (ro_data%lev1b%Npoints > 0) THEN
256 CALL ropp_fm_roprof2obs(ro_data, obs_bangle) !!set to obs levels
257 obs_bangle%bangle(:) = ropp_MDFV
258 ELSE
259 CALL set_obs_levels_bangle(ro_data, obs_refrac, obs_bangle) !default
260 ENDIF
261
262!-------------------------------------------------------------------------------
263! 11. Calculate bending angle and its gradient
264!-------------------------------------------------------------------------------
265
266 IF (state%state_ok) THEN
267
268 CALL ropp_fm_bangle_1d(state, obs_bangle)
269 IF (do_adj) THEN
270 ALLOCATE(gradient_bangle(SIZE(obs_bangle%bangle),SIZE(state%state)))
271 CALL ropp_fm_bangle_1d_grad(state, obs_bangle, gradient_bangle)
272 ENDIF
273
274 ENDIF
275
276!-------------------------------------------------------------------------------
277! 12. Copy simulated observations to RO structure and write data
278!-------------------------------------------------------------------------------
279
280 CALL ropp_fm_obs2roprof(obs_refrac, ro_data)
281 CALL ropp_fm_obs2roprof(obs_bangle, ro_data)
282
283 IF (state%state_ok) THEN
284
285 IF (do_adj) THEN
286 CALL ropp_io_addvar(ro_data, &
287 name = "gradient_refrac", &
288 long_name = &
289 "Gradient of the refractivity forward model", &
290 units = "1", &
291 range = (/MINVAL(gradient_refrac), &
292 MAXVAL(gradient_refrac)/), &
293 DATA = gradient_refrac)
294
295 CALL ropp_io_addvar(ro_data, &
296 name = "gradient_bangle" , &
297 long_name = &
298 "Gradient of the bending angle forward model", &
299 units = "rad", &
300 range = (/MINVAL(gradient_bangle), &
301 MAXVAL(gradient_bangle)/), &
302 DATA = gradient_bangle)
303 ENDIF
304
305!-------------------------------------------------------------------------------
306! 13. Update RO structure with computed state vector variables
307!-------------------------------------------------------------------------------
308
309 CALL ropp_fm_state2roprof(state, ro_data)
310
311 ENDIF
312
313!-------------------------------------------------------------------------------
314! 14. Write data
315!-------------------------------------------------------------------------------
316
317 CALL ropp_io_write(ro_data, ofile, append = .TRUE., ranchk = .TRUE. )
318
319!-------------------------------------------------------------------------------
320! 15. Clean up
321!-------------------------------------------------------------------------------
322
323 IF (state%state_ok .AND. do_adj) THEN
324 DEALLOCATE (gradient_refrac)
325 DEALLOCATE (gradient_bangle)
326 ENDIF
327
328 CALL ropp_fm_free(state)
329 CALL ropp_fm_free(obs_refrac)
330 CALL ropp_fm_free(obs_bangle)
331 CALL ropp_io_free(ro_data)
332
333 END DO
334 END DO
335
336CONTAINS
337
338!-------------------------------------------------------------------------------
339! 16. Calculate observation levels for refractivity
340!-------------------------------------------------------------------------------
341
342 SUBROUTINE set_obs_levels_refrac(ro_data, obs_refrac)
343
344! x.1 Declarations
345! ----------------
346
347 USE typesizes, ONLY: wp => EightByteReal
348 USE ropp_io_types
349 USE ropp_fm
350
351 IMPLICIT NONE
352
353 TYPE(ROprof) :: ro_data
354 TYPE(Obs1dRefrac) :: obs_refrac
355
356 INTEGER :: i, n, dummy
357
358 dummy = ro_data%Lev2a%Npoints !! fix nag ro_data 'unused dummy variable'
359
360! x.2 Vertical geopotential height levels between 200 and 60000 gpm
361! -----------------------------------------------------------------
362
363 n = INT(60.0_wp / 0.2_wp)
364
365 ALLOCATE(obs_refrac%refrac(n))
366 ALLOCATE(obs_refrac%geop(n))
367 ALLOCATE(obs_refrac%weights(n))
368
369 obs_refrac%refrac(:) = 0.0_wp
370 obs_refrac%geop(:) = (/ (i*200.0_wp, i = 1,n) /)
371 obs_refrac%weights(:) = 1.0_wp
372
373 END SUBROUTINE set_obs_levels_refrac
374
375
376!-------------------------------------------------------------------------------
377! 17. Calculate obs levels for bending angle (consistent with obs_refrac)
378!-------------------------------------------------------------------------------
379
380 SUBROUTINE set_obs_levels_bangle(ro_data, obs_refrac, obs_bangle)
381
382! x.1 Declarations
383! ----------------
384
385 USE typesizes, ONLY: wp => EightByteReal
386 USE geodesy
387 USE ropp_io_types
388 USE ropp_fm
389
390 IMPLICIT NONE
391
392 TYPE(ROprof) :: ro_data
393 TYPE(Obs1dRefrac) :: obs_refrac
394 TYPE(Obs1dbangle) :: obs_bangle
395
396 REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp
397
398 INTEGER :: n
399
400! x.2 Allocate arrays
401! -------------------
402
403 n = SIZE(obs_refrac%geop)
404
405 ALLOCATE(obs_bangle%bangle(n))
406 ALLOCATE(obs_bangle%impact(n))
407 ALLOCATE(obs_bangle%weights(n))
408
409 ALLOCATE(tmp(n))
410
411! x.3 Set scalar arguments of the observation vector
412! --------------------------------------------------
413
414 obs_bangle%g_sfc = gravity(ro_data%GEOref%lat)
415 obs_bangle%r_earth = R_eff(ro_data%GEOref%lat)
416 obs_bangle%undulation = ro_data%GEOref%undulation
417 IF (ro_data%GEOref%roc > 0.0_wp) THEN
418 obs_bangle%r_curve = ro_data%GEOref%roc
419 ELSE
420 obs_bangle%r_curve = obs_bangle%r_earth
421 ENDIF
422
423! x.4 Calculate levels to coincide with the geopotential levels
424! -------------------------------------------------------------
425
426 tmp = geopotential2geometric(ro_data%GEOref%lat, obs_refrac%geop) &
427 + obs_bangle%r_curve + obs_bangle%undulation
428
429 obs_bangle%impact(:) = (1.0_wp + obs_refrac%refrac*1.e-6_wp) * tmp
430
431! x.5 Fill other arrays
432! ---------------------
433
434 obs_bangle%bangle(:) = 0.0_wp
435 obs_bangle%weights(:) = 1.0_wp
436
437! x.6 Clean up
438! ------------
439
440 DEALLOCATE(tmp)
441
442 END SUBROUTINE set_obs_levels_bangle
443
444!-------------------------------------------------------------------------------
445! 18. Usage information
446!-------------------------------------------------------------------------------
447
448 SUBROUTINE usage()
449 PRINT *, 'Purpose:'
450 PRINT *, ' Bending angles and refractivity forward model'
451 PRINT *, 'Usage:'
452 PRINT *, ' > ropp_fm_bg2ro [<options>] <input_file(s)>'
453 PRINT *, 'Options:'
454 PRINT *, ' -o <output_file> name of ROPP netCDF output file'
455 PRINT *, ' -f forward only, no adjoint calculation.'
456 PRINT *, ' -use_logp use log(pressure) for forward model'
457 PRINT *, ' -use_logq use log(spec humidity) for forward model'
458 PRINT *, ' -comp include non ideal gas compressibility'
459 PRINT *, ' -d output additional diagnostics'
460 PRINT *, ' -h this help'
461 PRINT *, ' -v version information'
462 PRINT *, ''
463 END SUBROUTINE usage
464
465!-------------------------------------------------------------------------------
466! 19. Version information
467!-------------------------------------------------------------------------------
468
469 SUBROUTINE version_info()
470 CHARACTER (LEN=40) :: version
471 version = ropp_fm_version()
472 PRINT *, 'ropp_fm_bg2ro_1d - Bending angles and refractivity forward model.'
473 PRINT *, ''
474 PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
475 PRINT *, ''
476 END SUBROUTINE version_info
477
478END PROGRAM ropp_fm_bg2ro_1d