da_ominusb.f90
References to this file elsewhere.
1 program daprog_ominusb
2
3 !-----------------------------------------------------------
4 ! Abstract:
5 ! Purpose: Program for tuning observation errors
6 ! (Hollingsworh method)
7 ! Ref: Tellus (1986) 38, pp.111-161 (Part I & II)
8 !-----------------------------------------------------------
9
10 use da_control, only : filename_len
11
12 implicit none
13
14 character*5, parameter :: missing_c = '*****'
15 integer, parameter :: ounit = 10
16 integer, parameter :: num_bins_p = 10
17 integer, parameter :: obs_qc_pointer = 0
18 integer, parameter :: max_stations = 5000
19 integer, parameter :: max_times = 200
20 integer, parameter :: inunit = 35
21 integer, parameter :: min_obs = 30
22 integer, parameter :: num_bins = 100
23 integer, parameter :: missing_i = -88
24 real, parameter :: bottom_pressure = 1000.0
25 real, parameter :: bin_width_p = 100.0
26 real, parameter :: missing_r = -888888.0
27 real, parameter :: max_distance = 5000 ! km
28 real, parameter :: pi = 3.1415926535897932346
29 ! real, parameter :: earth_radius = 6378.1500 ! km
30 real, parameter :: earth_radius = 6370.00 ! km, consistant with WRF
31
32 type sub_type
33 integer :: qcflag(1:max_times)
34 real :: ominusb(1:max_times)
35 real :: sigmao(1:max_times)
36 real :: pressure(1:max_times)
37 end type sub_type
38
39 type obs_type
40 character*5 :: id(1:max_stations)
41 integer :: num_reject
42 integer :: num_keep
43 integer :: num_stations
44 integer :: num_obs(1:max_stations)
45 real :: mean_omb
46 real :: stdv_omb
47 real :: lat(1:max_stations)
48 real :: lon(1:max_stations)
49 real :: mean_ominusb(1:max_stations)
50 real :: dis(1:max_stations,1:max_stations)
51 real :: cov(1:num_bins)
52 type (sub_type) :: data(1:max_stations)
53 end type obs_type
54
55 character(len=filename_len) :: filename
56 character*5 :: station_chosen
57 character*5 :: station_id
58 integer :: times, n, n1, n2, b
59 integer :: bin
60 integer :: num_times
61 integer :: qc
62 integer :: percent_reject
63 real :: target_p, p_ob
64 real :: lati, long, iv, error
65 logical :: surface_obs
66 ! type (obs_type) :: obs(1:num_bins_p)
67 type (obs_type),allocatable :: obs(:)
68
69
70 real :: bin_start_p(1:num_bins_p) ! Number of pressure bins
71 real :: obs_err_used(1:num_bins_p) ! Obs error currently used
72
73
74 ! for temperature it is 1.0
75 ! data obs_err_used/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/
76 data obs_err_used/ 1.1, 1.1, 1.1, 1.4, 1.8, 2.3, 2.8, 3.3, 3.3, 2.7/
77 ! 1000 900 800 700 600 500 400 300 200 100
78 print*,' num_bins_p = ',num_bins_p
79 allocate ( obs( 1: num_bins_p) )
80 do b = 1, num_bins_p
81 bin_start_p(b) = bottom_pressure - real(b-1) * bin_width_p
82
83 ! Initialize obs structure:
84 call da_init_obs( obs(b) )
85 end do
86
87 !-----------------------------------------------------------------------
88 ! [1.0] Read in O-B data and fill structure:
89 !-----------------------------------------------------------------------
90
91 times = 0
92 ! station_chosen = "MOWV3"
93 ! filename = "daprog_ominusb_poamvu_"//station_chosen(1:5)//"_omb.out"
94 ! station_chosen = "91652"
95 ! filename = "daprog_ominusb_soundt_"//station_chosen(1:5)//"_omb.out"
96 filename = "ominusb.out"
97 open(ounit, file = filename, status = "unknown" )
98
99 do ! Loop over time:
100 times = times + 1
101
102 do ! Loop over obs in file:
103
104 read(inunit,'(a5,2f9.3,3f17.7,i8)')station_id, lati, long, &
105 p_ob, iv, error, qc
106
107 ! [1.1] Exit when come to end of file markers:
108 ! if ( station_id(1:5) == '*****' .or. station_id == '*end*' ) exit
109 if ( station_id(1:5) == '*****' .or. station_id == '*end*' ) then
110 print*,' Hit station_id as ',station_id
111 exit
112 end if
113
114 ! [1.2] Store data:
115
116 ! if ( station_id(1:5) == station_chosen ) then
117
118 p_ob = 0.01 * p_ob ! Convert Pa to hPa
119
120 bin = int( (bottom_pressure - p_ob)/bin_width_p ) + 1
121 if (bin <= 0) bin = 1
122 if (bin > num_bins_p) bin = num_bins_p
123
124 call da_store_ominusb( station_id, times, qc, lati, long, &
125 p_ob, iv, error, obs(bin) )
126 ! end if
127
128 end do
129
130 ! [1.5] A station_id = '*end*' indicates the complete end of the file:
131 if ( station_id == '*end*' ) then
132 exit
133 end if
134
135 end do
136
137 num_times = times
138 write(0,'(A,i8)')' Number of analysis times read in = ', num_times
139 write(0,*)
140
141 !-----------------------------------------------------------------------
142 ! [2.0]: Calculate O-B statistics over all stations:
143 !-----------------------------------------------------------------------
144
145 write(0,'(A)')' P(hPA) Stations Obs-Used Obs-Reject Mean(O-B) STDV(O-B)'
146 write(ounit,'(A)')' P(hPA) Stations Obs-Used Obs-Reject Mean(O-B) STDV(O-B)'
147 ! write(0,'(a,a)')' Station chosen = ', station_chosen
148 do b = 1, num_bins_p
149 call da_obs_stats( num_times, obs(b) )
150
151 percent_reject = 0
152 if ( obs(b) % num_stations > 0 ) then
153 percent_reject = nint( 100.0 * real( obs(b) % num_reject) / &
154 real( obs(b) % num_keep + obs(b) % num_reject ) )
155 end if
156
157 write(ounit,'(i6,i8,2x,2i8,i3,2f13.6)')int(bin_start_p(b)), &
158 obs(b) % num_stations, obs(b) % num_keep, &
159 obs(b) % num_reject, percent_reject, &
160 obs(b) % mean_omb, obs(b) % stdv_omb
161 end do
162 write(0,*)
163
164 !-----------------------------------------------------------------------
165 ! [3.0]: Calculate matrix of distances between points:
166 !-----------------------------------------------------------------------
167
168 do b = 1, num_bins_p
169 n = obs(b) % num_stations
170 call da_get_distance( n, obs(b) % lat(1:n), obs(b) % lon(1:n), obs(b) % dis(1:n,1:n) )
171 write(0,'(4i8,2f13.6)')int(bin_start_p(b)), obs(b) % num_stations, &
172 obs(b) % num_keep, obs(b) % num_reject, &
173 obs(b) % mean_omb, obs(b) % stdv_omb
174 end do
175
176 !-----------------------------------------------------------------------
177 ! [4.0]: Calculate O-B covariances:
178 !-----------------------------------------------------------------------
179
180 do b = 1, num_bins_p
181 ! call da_bin_covariance( obs(b) % num_stations, num_times, max_distance, obs(b) )
182 call da_bin_covariance( obs(b) % num_stations, num_times, max_distance, obs(b), obs_err_used(b),b)
183 end do
184
185 contains
186
187 subroutine da_init_obs( obs )
188
189 implicit none
190
191 type (obs_type), intent(out) :: obs
192
193 integer :: n
194
195 obs % num_reject = 0
196 obs % num_keep = 0
197 obs % num_stations = 0
198 obs % num_obs(:) = 0
199 obs % id(:) = missing_c
200 obs % lat(:) = missing_r
201 obs % lon(:) = missing_r
202 do n = 1, max_stations
203 obs % data(n) % qcflag(:) = missing_i
204 obs % data(n) % ominusb(:) = missing_r
205 obs % data(n) % sigmao(:) = missing_r
206 obs % data(n) % pressure(:) = missing_r
207 end do
208
209 end subroutine da_init_obs
210
211 subroutine da_store_ominusb( station_id, times, qc, lati, long, p_ob, &
212 iv, error, obs )
213
214 implicit none
215
216 character*5, intent(in) :: station_id
217 integer, intent(in) :: times
218 integer, intent(in) :: qc
219 real, intent(in) :: lati
220 real, intent(in) :: long
221 real, intent(in) :: p_ob
222 real, intent(in) :: iv
223 real, intent(in) :: error
224 type(obs_type), intent(inout) :: obs
225
226 integer :: this_station
227 integer :: n
228 real :: pi_over_180
229
230 pi_over_180 = pi / 180.0
231
232 ! [1.0] Count total number of obs rejected by QC:
233
234 if ( qc < obs_qc_pointer ) then
235 obs % num_reject = obs % num_reject + 1
236 else
237 obs % num_keep = obs % num_keep + 1
238 end if
239
240 ! [2.0] Check if this station seen before:
241
242 if ( ANY(obs % id(:) == station_id) ) then
243 do n = 1, obs % num_stations
244 if ( obs % id(n) == station_id ) then
245 this_station = n
246 end if
247 end do
248 else ! New station:
249 obs % num_stations = obs % num_stations + 1
250 this_station = obs % num_stations
251 obs % id(this_station) = station_id
252
253 if ( long >= 0.0 ) then
254 obs % lon(this_station) = pi_over_180 * long
255 else
256 obs % lon(this_station) = pi_over_180 * ( 180.0 - long )
257 end if
258 obs % lat(this_station) = pi_over_180 * lati
259
260 end if
261
262 ! [3.0] Check if getting too many stations or times for array sizes:
263
264 if ( obs % num_stations > max_stations ) then
265 write(0,'(A)')' Need to increase max_stations. Stopping'
266 stop
267 end if
268
269 ! [4.0] Read in qc, O-B, and sigma_o:
270
271 obs % data(this_station) % qcflag(times) = qc
272 obs % data(this_station) % ominusb(times) = iv
273 obs % data(this_station) % sigmao(times) = error
274 obs % data(this_station) % pressure(times) = p_ob
275
276 end subroutine da_store_ominusb
277
278 subroutine da_obs_stats( num_times, obs )
279
280 integer, intent(in) :: num_times
281 type (obs_type), intent(inout) :: obs
282
283 integer, parameter :: num_bins = 101
284
285 integer :: n, nobs, nobs_station, times, b, count
286 integer :: bin_count(1:num_bins), maxcount
287 real :: sumobs1, sumobs2, sum_station
288 real :: min_bin, max_bin, bin_width, omb, x, z
289 real :: bin_start(1:num_bins)
290 real :: dommean(1:num_times)
291
292 obs % mean_omb = 0.0
293 obs % stdv_omb = 0.0
294
295 nobs = 0
296 sumobs1 = 0.0
297 sumobs2 = 0.0
298
299 do n = 1, obs % num_stations
300 nobs_station = 0
301 sum_station = 0.0
302 obs % mean_ominusb(n) = 0.0
303 do times = 1, num_times
304 if ( obs % data(n) % qcflag(times) >= obs_qc_pointer ) then
305 nobs = nobs + 1
306 sumobs1 = sumobs1 + obs % data(n) % ominusb(times)
307 sumobs2 = sumobs2 + obs % data(n) % ominusb(times) **2
308
309 nobs_station = nobs_station + 1
310 sum_station = sum_station + obs % data(n) % ominusb(times)
311 end if
312 end do
313 obs % num_obs(n) = nobs_station
314 if ( nobs_station > 0 ) then
315 obs % mean_ominusb(n) = sum_station / real(nobs_station)
316 end if
317 end do
318
319 ! Calculate basic statistics:
320 if ( nobs > 0 ) then
321 obs % mean_omb = sumobs1 / real(nobs)
322 obs % stdv_omb = sumobs2 / real(nobs) ! Actually mean square here.
323 obs % stdv_omb = sqrt( obs % stdv_omb - obs % mean_omb**2 )
324 end if
325
326 ! Get data distribution:
327 max_bin = 5.0 * obs % stdv_omb
328 min_bin = -max_bin
329 bin_width = 2.0 * max_bin / real(num_bins)
330 bin_count(:) = 0
331 do b = 1, num_bins
332 bin_start(b) = min_bin + real(b-1) * bin_width
333 end do
334
335 do n = 1, obs % num_stations
336 do times = 1, num_times
337 omb = obs % data(n) % ominusb(times)
338
339 if ( obs % data(n) % qcflag(times) >= obs_qc_pointer ) then
340 b = int( (omb-min_bin) / bin_width) + 1
341 if (b >= 1 .and. b <= num_bins) bin_count(b) = bin_count(b) + 1
342 end if
343 end do
344 end do
345
346 maxcount = maxval(bin_count(:))
347 ! write(0,'(a,i8)')' Max count = ', maxcount
348 ! write(0,'(a)')' Bin x=O-B z=(x-xm)/sd Count exp(-0.5*z*z)'
349 do b = 1, num_bins
350 x = bin_start(b) + 0.5 * bin_width
351 z = ( x - obs % mean_omb ) / obs % stdv_omb
352 ! write(0,'(i4,4f12.5)')b, x, z, bin_count(b)/real(maxcount), exp(-0.5*z*z)
353 end do
354 ! write(0,*)
355
356 ! Get time series of mean error:
357 dommean(:) = 0.0
358 ! write(0,'(a)')' Time NumObs Domain Mean'
359 do times = 1, num_times
360 count = 0
361 do n = 1, obs % num_stations
362 if ( obs % data(n) % qcflag(times) >= obs_qc_pointer ) then
363 dommean(times) = dommean(times) + obs % data(n) % ominusb(times)
364 count = count + 1
365 end if
366 end do
367 if ( count > 0 ) dommean(times) = dommean(times) / real(count)
368 ! write(0,'(i4,i8,f12.5)')times, count, dommean(times)
369 end do
370 ! write(0,*)
371
372 ! Remove station mean from O-B values (removes instrumental/model bias):
373
374 do n = 1, obs % num_stations
375 do times = 1, num_times
376 if ( obs % data(n) % qcflag(times) >= obs_qc_pointer ) then
377 obs % data(n) % ominusb(times) = obs % data(n) % ominusb(times) - &
378 obs % mean_ominusb(n)
379 ! obs % mean_omb
380 end if
381 end do
382 end do
383
384 end subroutine da_obs_stats
385
386 subroutine da_get_distance( num_stations, lat, lon, dis )
387
388 implicit none
389
390 integer, intent(in) :: num_stations
391 real, intent(in) :: lat(1:num_stations)
392 real, intent(in) :: lon(1:num_stations)
393 real, intent(out) :: dis(1:num_stations,1:num_stations)
394
395 integer :: n1, n2, min_n1, min_n2
396 real :: pi_over_2, colat1, colat2, londiff
397 real :: dist
398
399 pi_over_2 = 0.5 * pi
400
401 dis(1:num_stations,1:num_stations) = 0.0
402 do n1 = 1, num_stations
403 colat1 = pi_over_2 - lat(n1)
404
405 do n2 = n1, num_stations
406 if ( n1 /= n2 ) then
407 colat2 = pi_over_2 - lat(n2)
408 londiff = abs(lon(n2) - lon(n1))
409
410 dist = acos( cos(colat1) * cos(colat2) + &
411 sin(colat1) * sin(colat2) * cos(londiff) ) !in radians
412 dis(n1,n2) = earth_radius * dist ! in km
413 end if
414
415 end do
416 end do
417
418 end subroutine da_get_distance
419
420 subroutine da_bin_covariance( num_stations, num_times, max_distance, obs, obs_err_used, lvl)
421
422 implicit none
423
424 integer, intent(in) :: num_stations
425 integer, intent(in) :: num_times
426 real, intent(in) :: max_distance
427 type (obs_type), intent(inout) :: obs
428 real, intent(in) :: obs_err_used
429 integer, intent(in) :: lvl
430
431 integer :: b, n1, n2, times, sum_b
432 real :: bin_width, bin_dis
433 real :: sum_dis, sum_covar, sum2covar
434 real :: covar, gaussian
435 real :: x, y, sum_x, sum_x2, sum_xy, sum_y
436 real :: gradient, lengthscale, intercept
437 real :: bk_err_var, bk_err_stdv, ob_err_stdv
438 real :: avg_dis(1:num_bins)
439 real :: coverr(1:num_bins)
440 integer :: sum_obs(1:num_bins)
441 logical :: bin_pass(1:num_bins)
442
443 avg_dis(:) = 0.0
444 coverr(:) = 0.0
445 obs % cov(:) = 0.0
446 sum_obs(:) = 0
447
448 bin_width = max_distance / real(num_bins)
449
450 ! Calculate covariance and sum for all good obs:
451
452 do b = 1, num_bins
453 bin_dis = real(b-1) * bin_width
454 sum_dis = 0.0
455 sum_covar = 0.0
456 sum2covar = 0.0
457
458 do n1 = 1, num_stations
459 do n2 = n1, num_stations
460 if ( obs % dis(n1,n2) >= bin_dis .and. &
461 obs % dis(n1,n2) < bin_dis + bin_width .and. &
462 obs % dis(n1,n2) /= 0.0 ) then
463
464 do times = 1, num_times
465 if ( obs % data(n2) % qcflag(times) >= obs_qc_pointer .and. &
466 obs % data(n1) % qcflag(times) >= obs_qc_pointer ) then
467
468 covar = obs % data(n2) % ominusb(times) * &
469 obs % data(n1) % ominusb(times)
470
471 sum_obs(b) = sum_obs(b) + 1
472 sum_dis = sum_dis + obs % dis(n1,n2)
473 sum_covar = covar + sum_covar
474 sum2covar = covar**2 + sum2covar
475
476 end if
477 end do
478
479 end if
480 end do
481 end do
482
483 ! write(6,'(2i8,2f15.5)')b, sum_obs(b), sum_covar, obs % cov(b)
484
485 ! Calculate average separation and covariance of obs:
486 if ( sum_obs(b) > 0 ) then
487 avg_dis(b) = sum_dis / real(sum_obs(b))
488 obs % cov(b) = sum_covar / real(sum_obs(b))
489 sum2covar = sum2covar / real(sum_obs(b))
490 end if
491
492 ! Calculate 95% error bar for obs % cov estimate = 1.96 s/sqrt(n-1):
493 ! Needs min_obs observations in bin to be considered in Gaussian calc.
494 bin_pass(b) = .false.
495
496 if ( sum_obs(b) > min_obs ) then
497 coverr(b) = 1.96 * sqrt( sum2covar - obs % cov(b)**2 ) / &
498 sqrt( real(sum_obs(b) - 1) )
499
500 if ( coverr(b) > 0.0 .and. coverr(b) < obs % cov(b) ) then
501 bin_pass(b) = .true.
502 end if
503 end if
504
505 end do
506
507 ! Fit good data to a Gaussian to get lengthscale:
508 ! B = B0 exp (-r**2 / 8s**2) => ln B = ln B0 - r**2 / 8s**2
509 ! y = gradient * x + intercept => gradient = -1/8s**2, x = r**2, intercept = lnB0
510
511 sum_b = 0
512 lengthscale = 0.0
513 gradient = 0.0
514
515 if ( any(bin_pass) ) then
516 sum_x = 0.0
517 sum_x2= 0.0
518 sum_xy= 0.0
519 sum_y = 0.0
520
521 do b = 1, num_bins
522 if ( bin_pass(b) ) then
523 x = avg_dis(b) * avg_dis(b)
524 y = log( obs % cov(b) )
525 sum_b = sum_b + 1
526 sum_x = sum_x + x
527 sum_x2 = sum_x2 + x * x
528 sum_xy = sum_xy + x * y
529 sum_y = sum_y + y
530 end if
531 end do
532
533 if ( sum_b > 1 ) then
534 gradient = ( real(sum_b) * sum_xy - sum_x * sum_y ) / &
535 ( real(sum_b) * sum_x2 - sum_x * sum_x )
536 intercept = ( sum_x2 * sum_y - sum_x * sum_xy ) / &
537 ( real(sum_b) * sum_x2 - sum_x * sum_x )
538 if ( gradient < 0.0 ) then
539 lengthscale = sqrt( -1.0 / ( 8.0 * gradient ) )
540 bk_err_var = exp(intercept)
541 end if
542 end if
543 end if
544
545 ! Print out covariance data:
546
547 if (bin_pass(b)) &
548 write(0,'(a)')' Bin NumObs Separation(km) Covariance CovErr GaussFit'
549 gaussian = 0.0
550 do b = 1, num_bins
551 if (bin_pass(b)) &
552 write(0,'(a)')' Bin NumObs Separation(km) Covariance CovErr GaussFit'
553 if ( gradient < 0.0 ) then
554 gaussian = bk_err_var*exp(-0.125*(avg_dis(b)/lengthscale)**2)
555 end if
556
557 if (obs % cov(b) /= 0.0 ) then
558 write(0,'(i5,i7,3e14.6,l1,e13.6)')b, sum_obs(b), avg_dis(b), &
559 obs % cov(b), &
560 coverr(b), bin_pass(b), gaussian
561 if (obs % cov(b) /= 0.0 .and. bin_pass(b) ) &
562 write(70+lvl,'(2i7,4e14.6)')b, sum_obs(b), avg_dis(b), obs % cov(b),&
563 coverr(b), gaussian
564 end if
565 end do
566 write(0,*)
567
568 if ( sum_b > 1 ) then
569 ob_err_stdv = sqrt( obs % stdv_omb**2 - bk_err_var )
570 bk_err_stdv = sqrt( bk_err_var )
571 write(0,'(a)')' The following are derived from Gaussian fit (warning):'
572 write(0,'(3(a,e14.6),2(a,i3),a)')' Scale =', lengthscale, &
573 ' km, ob_err_stdv = ', ob_err_stdv, &
574 ' , bk_err_stdv = ', bk_err_stdv, &
575 ' using', sum_b, ' /', num_bins, ' bins.'
576 write(30,'(4f15.3)') ob_err_stdv, bk_err_stdv, obs_err_used, lengthscale
577 end if
578
579 end subroutine da_bin_covariance
580
581 end program daprog_ominusb