Ticket #182: ropp_fm_bg2ro_1d_AXEL0109.f90

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