var/da/da_util/da_tune_obs_hollingsworth2.f90

   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 = 2526
  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