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