var/da/da_tools/da_buddy_qc.inc

   1 SUBROUTINE da_buddy_qc ( numobs, m_max, station_id, xob, yob, obs, qc_flag_small,    &
   2                          name, pressure, dx, buddy_weight , buddy_difference, &
   3                          iunit, print_buddy, xob_g, yob_g, obs_g, qc_flag_small_g, num_recv)
   4
   5 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   6 !   Imported from BUDDY in RAWINS
   7
   8 !   Notes:
   9 !      This routine performs error checks by comparing the difference
  10 !         between the first guess field
  11 !         and observations at the observing locations (xob, yob)
  12 !         with the averaged difference at the nearby stations within
  13 !         a given distance. If the difference value at the observation
  14 !         point differs from the average of the nearby difference values
  15 !         by more than a certain maximum (dependent on variable types,
  16 !         pressure levels, and BUDWGT parameter in namelist), the
  17 !         observation is flagged as suspect, and add the flag message
  18 !         to the observation records.
  19
  20 !   *** Note that x and y are in natural coordinates now.
  21 !         Need to check if the x and y are used correctly...
  22
  23 !   other variables:
  24 !    err:       difference between observations and first guess fields
  25 !
  26 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  27 !  This subroutine da_buddy_qc is adapted by Y.-R.Guo, NCAR/MMM, based on
  28 !  the WRFV3/OBSGRID/src/qc2_module.F90.
  29 !
  30 !  The algorithm for buddy check used here could be named as "Two Passes
  31 !  Innovation Buddy Check": the buddy members to a specific obs are
  32 !  determined by going through two times of checks for innovation deprture
  33 !  from the buddy member's mean values.
  34 !
  35 !  In OBSGRID, the buddy check is completed at each of  the pressure levels
  36 !  for all obs, i.e. all obs types are mixed to gether. In order to use this
  37 !  buddy check in wrfvar, first we should sort the obs into the different
  38 !  pressure bins for 3-D observations,such as SOUND, then apply the buddy
  39 !  check for each of the obs types in da_get_innov_vector_?????.inc. For the
  40 !  2-D observations, such as SYNOP, no binning is needed.
  41 !
  42 !  Most of the tolerances for each of variables at the different pressure
  43 !  levels are adapted from WRFV3/OBSGRID. The tolerances for specific humidity
  44 !  used here are converted from rh and specified temperature. The tolerance
  45 !  for PMSL is reduced to 350 Pa defined in da_control/da_control.f90.
  46 !
  47 !                                  Yong-Run Guo, 10/10/2008
  48 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  49 !
  50
  51   IMPLICIT NONE
  52
  53   INTEGER,                       intent(in)    :: numobs, m_max, iunit
  54   REAL,    DIMENSION ( m_max ), intent(in)    :: obs , xob , yob
  55   INTEGER, DIMENSION ( m_max ), intent(inout) :: qc_flag_small
  56   REAL,    DIMENSION ( m_max, 0:num_procs-1 ), intent(in), optional :: obs_g , xob_g , yob_g
  57   INTEGER, DIMENSION ( m_max, 0:num_procs-1 ), intent(inout), optional :: qc_flag_small_g
  58   INTEGER, DIMENSION ( 0:num_procs-1 ),        intent(in), optional :: num_recv
  59   REAL,                    intent(in)          :: buddy_weight,   &
  60                                                   buddy_difference , &
  61                                                   dx          , &
  62                                                   pressure
  63   CHARACTER(LEN = 5), DIMENSION(m_max), intent(in) :: station_id
  64   CHARACTER(LEN = 8),                    intent(in) :: name
  65   LOGICAL           ,                    intent(in) :: print_buddy
  66
  67   !  err:         difference between observations and first guess fields
  68
  69   INTEGER                                    :: num, num3, num4, numj, &
  70                                                 kob, n, ierr
  71   REAL                                       :: range,              &
  72                                                 distance,           &
  73                                                 sum,                &
  74                                                 average,            &
  75                                                 x, y,               &
  76                                                 diff_numj,          &
  77                                                 diff_check_1,       &
  78                                                 diff_check_2
  79   REAL                                          plevel,             &
  80                              pp, tt, rh, t_error, rh_error, q_error
  81   REAL  , DIMENSION ( numobs )               :: err,                &
  82                                                 difmax,             &
  83                                                 diff_at_obs
  84 #ifdef DM_PARALLEL
  85   REAL  , DIMENSION ( m_max*num_procs )      :: diff
  86   INTEGER, DIMENSION (3)                     :: i_dummy
  87 #else
  88   REAL  , DIMENSION ( numobs )               :: diff
  89 #endif
  90   INTEGER , DIMENSION ( numobs )             :: buddy_num,          &
  91                                                 buddy_num_final
  92   REAL                                       :: scale
  93   INTEGER, DIMENSION ( m_max )               :: qc_flag
  94
  95
  96   real :: f_qv_from_rh
  97   external f_qv_from_rh
  98
  99 ! [1.0] Store original qc flags:
 100   qc_flag = qc_flag_small
 101
 102   plevel = pressure
 103
 104 ! [2.0] Define distance within which buddy check is to find neighboring stations:
 105
 106  ! 2.1 Define the horizontal range for buddy check:
 107   IF ( plevel .LE. 201.0 ) THEN
 108      range = 650.
 109   ELSE IF ( plevel .LE. 1000.1 ) THEN
 110      range = 300.
 111   ELSE
 112      range = 500.
 113   END IF
 114
 115  ! 2.2 Define tolerance value for comparison:
 116
 117   error_allowance_by_field : SELECT CASE ( name )
 118
 119      CASE ( 'UU      ' , 'VV      ' )
 120         IF      ( plevel .LT. 401. ) THEN
 121            difmax = buddy_difference * 1.7
 122         ELSE IF ( plevel .LT. 701. ) THEN
 123            difmax = buddy_difference * 1.4
 124         ELSE IF ( plevel .LT. 851. ) THEN
 125            difmax = buddy_difference * 1.1
 126         ELSE
 127            difmax = buddy_difference
 128         END IF
 129
 130      CASE ( 'PMSL    ' )
 131         difmax = buddy_difference
 132
 133      CASE ( 'TT      ' )
 134         IF      ( plevel .LT. 401. ) THEN
 135            difmax = buddy_difference * 0.6
 136         ELSE IF ( plevel .LT. 701. ) THEN
 137            difmax = buddy_difference * 0.8
 138         ELSE
 139            difmax = buddy_difference
 140         END IF
 141
 142      CASE ( 'QQ      ' )
 143 !         difmax = buddy_difference
 144
 145     ! Convert the rh error to q error:
 146
 147         t_error = 2.0
 148         rh_error = buddy_difference
 149
 150         pp = pressure * 100.0
 151         if ( pp > 85000.0 ) then
 152           rh = 90.0
 153           tt = 300.0
 154         else if ( pp > 70000.0 ) then
 155           rh = 80.0
 156           tt = 290.0
 157         else if ( pp > 35000.0 ) then
 158           rh = 70.0
 159           tt = 280.0
 160         else if ( pp > 10000.0 ) then
 161           rh = 60.0
 162           tt = 260.0
 163         else
 164           rh = 50.0
 165           tt = 240.0
 166         endif
 167
 168         q_error = f_qv_from_rh( rh_error, t_error, rh, tt, pp)
 169         difmax = q_error
 170
 171 !         print '("p, t, rh, t_error, rh_error, difmax:",5f12.2,e15.5)', &
 172 !            pp, tt, rh, t_error, rh_error, q_error
 173
 174   END SELECT error_allowance_by_field
 175
 176   difmax = difmax * buddy_weight
 177
 178 ! [3.0] Compute the errors for buddy check:
 179
 180  ! 3.1 Loop through station locations (numobs) to compute buddy average
 181  !      at each observation locations.
 182
 183 !   station_loop_2 : DO num = numobs, 1, -1
 184   station_loop_2 : DO num = 1, numobs
 185
 186      average = 0.0
 187  ! No buddy check is applied to a bad data (qc < 00:
 188      if ( qc_flag_small(num) < 0 ) cycle  station_loop_2
 189
 190      !  Compute distance between all pairs of observations,
 191      !  find number of observations within a given distance: buddy_num,
 192      !  and compute the difference between first guess and obs.
 193
 194      buddy_num(num) = 0
 195 #ifdef DM_PARALLEL
 196      if (present(xob_g) .and. present(yob_g) .and. present(obs_g) .and. present(qc_flag_small_g)) then
 197      else
 198          call da_error(__FILE__,__LINE__, &
 199                       (/"Buddy check in parallel mode needs global observations"/))
 200      endif
 201      pe_loop : DO n = 0, num_procs-1
 202      station_loop_3 : DO num3 = 1, num_recv(n)
 203
 204         x = xob(num) - xob_g(num3,n)
 205         y = yob(num) - yob_g(num3,n)
 206         IF ( ABS(x) .GT. 1.E-20 .OR. ABS(y) .GT. 1.E-20 ) THEN
 207            distance = SQRT ( x*x + y*y ) * dx
 208            IF ( distance .LE. range .AND. distance .GE. 0.1 &
 209                 .and. qc_flag_small_g(num3,n) >= 0 ) THEN
 210 ! Only good data are included in buddy check group:
 211               buddy_num(num) = buddy_num(num) + 1
 212 ! innovation (O-B) ==> diff:
 213               diff(buddy_num(num)) = obs_g(num3,n)
 214            END IF
 215         END IF
 216
 217      END DO station_loop_3
 218      END DO pe_loop
 219 #else
 220 !      station_loop_3 : DO num3 = numobs, 1, -1
 221      station_loop_3 : DO num3 = 1, numobs, 1
 222
 223         IF ( num3 .NE. num ) THEN
 224            x = xob(num) - xob(num3)
 225            y = yob(num) - yob(num3)
 226            distance = SQRT ( x*x + y*y ) * dx
 227            IF ( distance .LE. range .AND. distance .GE. 0.1 &
 228                 .and. qc_flag_small(num3) >= 0 ) THEN
 229 ! Only good data are included in buddy check group:
 230               buddy_num(num) = buddy_num(num) + 1
 231 ! innovation (O-B) ==> diff:
 232               diff(buddy_num(num)) = obs(num3)
 233            END IF
 234         END IF
 235
 236      END DO station_loop_3
 237 #endif
 238
 239 ! innovation at the specific obs(num):
 240      diff_at_obs(num) = obs(num)
 241 !
 242 ! Summation of innovations over the surrounding obs:
 243      sum = 0.
 244      DO numj = 1, buddy_num(num)
 245         sum = sum + diff(numj)
 246      END DO
 247
 248      !  Check to see if there are any buddies, compute the mean innovation.
 249
 250      IF ( buddy_num(num) .GT. 0 ) average = sum / buddy_num(num)
 251
 252  ! 3.2 Check if there is any bad obs among the obs located within the
 253  !  the radius surrounding the test ob.
 254
 255      diff_check_1 = difmax (1) * 1.25
 256      diff_check_2 = difmax (1)
 257
 258      check_bad_ob : IF ( buddy_num(num) .GE. 2 ) THEN
 259
 260         kob = buddy_num(num)
 261         remove_bad_ob : DO numj = 1, buddy_num(num)
 262            diff_numj = ABS ( diff ( numj ) - average )
 263
 264            IF ( diff ( numj ) .GT. diff_check_1  &
 265                .AND. diff_numj .GT. diff_check_2 ) THEN
 266     ! Bad obs:                 Innovation itself: diff(numj) > diff_check_1
 267     !        The distance between the innovation and average > diff_check_2
 268               kob = kob - 1
 269               sum = sum - diff ( numj )
 270            END IF
 271         END DO remove_bad_ob
 272
 273 ! The final number of buddies:
 274         buddy_num_final(num) = kob
 275
 276         !  We may have removed too many observations.
 277
 278         IF ( kob .GE. 2 ) THEN
 279         !  Information for buddy check for specific obs(num)
 280            average = sum / kob
 281            err(num) = diff_at_obs(num) - average
 282         ELSE
 283         !  No buddy check completed for this specific obs(num):
 284            err(num)     = 0.
 285 ! Not change the flag (YRG, 02/10/2009)
 286 !            qc_flag_small(num) = qc_flag_small(num) + no_buddies
 287 #ifdef DM_PARALLEL
 288 ! Not change the flag (YRG, 02/10/2009)
 289 !            qc_flag_small_g(num,myproc) = qc_flag_small_g(num,myproc) + no_buddies
 290            i_dummy(1) = qc_flag_small_g(num,myproc)
 291            i_dummy(2) = num
 292            i_dummy(3) = myproc
 293            call mpi_bcast ( i_dummy, 3, mpi_integer, myproc, comm, ierr )
 294            qc_flag_small_g(i_dummy(2),i_dummy(3)) = i_dummy(1)
 295 #endif
 296         END IF
 297
 298      ELSE check_bad_ob
 299         ! Too less buddy numbers < 2, no buddy check completed:
 300         err(num)     = 0.
 301 ! Not change the flag (YRG, 02/10/20090
 302 !         qc_flag_small(num) = qc_flag_small(num) + no_buddies
 303 #ifdef DM_PARALLEL
 304 ! Not change the flag (YRG, 02/10/20090
 305 !         qc_flag_small_g(num,myproc) = qc_flag_small_g(num,myproc) + no_buddies
 306          i_dummy(1) = qc_flag_small_g(num,myproc)
 307          i_dummy(2) = num
 308          i_dummy(3) = myproc
 309          call mpi_bcast ( i_dummy, 3, mpi_integer, myproc, comm, ierr )
 310          qc_flag_small_g(i_dummy(2),i_dummy(3)) = i_dummy(1)
 311 #endif
 312
 313      END IF check_bad_ob
 314
 315  ! 3.3 If the buddy number is ONLY 2, increase the tolerance value:
 316      IF ( buddy_num_final(num) .EQ. 2 ) difmax(num) = 1.2 * difmax(num)
 317
 318   END DO station_loop_2
 319
 320 ! [4.0] Reset the qc flags according to the buddy check errors:
 321 !
 322 !       If an observation has been flagged as failing the buddy
 323 !       check criteria (error [err] larger than the allowable
 324 !       difference [difmax]), then qc_flag for this variable,
 325 !       level, time has to be modified.
 326
 327    station_loop_4 : do n = 1, numobs
 328    ! if the difference at station: n from the buddy-average > difmax,
 329    ! failed the buddy check:
 330      if ( abs(err(n)) > difmax(n) ) then
 331         qc_flag_small(n) = fails_buddy_check
 332      endif
 333    enddo station_loop_4
 334
 335 ! [5.0] This notifies the user which observations have been flagged as
 336 !       suspect by the buddy check routine.
 337
 338   IF ( print_buddy ) THEN
 339      n    = 0
 340      numj = 0
 341      num3 = 0
 342      num4 = 0
 343      station_loop_5 : DO num = 1 , numobs
 344
 345         IF ( name(1:8) .EQ. 'PMSL    ' ) THEN
 346            scale = 100.
 347         ELSE IF ( name(1:8) .EQ. 'QQ      ' ) THEN
 348            scale = 0.001
 349         ELSE
 350            scale = 1.0
 351         END IF
 352
 353         IF ( buddy_num_final(num) .GE. 2 ) THEN
 354           n = n + 1
 355           IF ( ABS ( err(num) ) .GT. difmax(num) ) THEN
 356              numj = numj + 1
 357              WRITE ( unit=iunit, FMT = '(3I8,&
 358                    &  " ID=",a5,&
 359                    &  ",NAME="  ,a8, &
 360                    &  ",BUDDY TOLERANCE=",f5.1,&
 361                    &  ",LOC=(" ,f6.1,",",f6.1,")",&
 362                    &  ",OBS INV="       ,f7.2,&
 363                    &  ",DIFF ERR="      ,f7.2,&
 364                    &  ",NOBS BUDDY="    ,I4,&
 365                    &  ",QC_FLAG=",I3,&
 366                    &  "  QC_FLAG_NEW=",I3,"  FAILED")' ) num, n, numj,&
 367                   station_id(num),name,difmax(num)/scale, &
 368                   xob(num),yob(num),obs(num)/scale,&
 369                   err(num)/scale, buddy_num_final(num), &
 370                   qc_flag(num), qc_flag_small(num)
 371           ELSE
 372              num3 = num3 + 1
 373 !              WRITE ( unit=iunit, FMT = '(3I8,&
 374 !                    &  " ID=",a5,&
 375 !                    &  ",NAME="  ,a8, &
 376 !                    &  ",BUDDY TOLERANCE=",f5.1,&
 377 !                    &  ",LOC=(" ,f6.1,",",f6.1,")",&
 378 !                    &  ",OBS INV="       ,f7.2,&
 379 !                    &  ",DIFF ERR="      ,f7.2,&
 380 !                    &  ",NOBS BUDDY="    ,I4,&
 381 !                    &  ",QC_FLAG=",I3,&
 382 !                    &  "  QC_FLAG_NEW=",I3,"  PASSED")' ) num, n, num3,&
 383 !                   station_id(num),name,difmax(num)/scale, &
 384 !                   xob(num),yob(num),obs(num)/scale,&
 385 !                   err(num)/scale, buddy_num_final(num), &
 386 !                   qc_flag(num), qc_flag_small(num)
 387           ENDIF
 388         ELSE
 389           num4 = num4 + 1
 390 !              WRITE ( unit=iunit, FMT = '(3I8,&
 391 !                    &  " ID=",a5,&
 392 !                    &  ",NAME="  ,a8, &
 393 !                    &  ",BUDDY TOLERANCE=",f5.1,&
 394 !                    &  ",LOC=(" ,f6.1,",",f6.1,")",&
 395 !                    &  ",OBS INV="       ,f7.2,&
 396 !                    &  ",DIFF ERR="      ,f7.2,&
 397 !                    &  ",NOBS BUDDY="    ,I4,&
 398 !                    &  ",QC_FLAG=",I3,&
 399 !                    &  "  QC_FLAG_NEW=",I3,"NO BUDDY")' ) num, n, num4,&
 400 !                   station_id(num),name,difmax(num)/scale, &
 401 !                   xob(num),yob(num),obs(num)/scale,&
 402 !                   err(num)/scale, buddy_num_final(num), &
 403 !                   qc_flag(num), qc_flag_small(num)
 404         END IF
 405      END DO station_loop_5
 406
 407      write(unit=iunit,fmt = '(5x,"NOB=",i6,&
 408           & "  Toltal N_buddy-checked =",i6,&
 409           & "  N_Passed =",i6,"  N_Failed=",i6,"  NO_Buddy=",i6, &
 410           & "  RANGE=",f10.2,"  DIFMAX=",f10.2)') &
 411           numobs, n, num3, numj, num4, range, difmax(1)/scale
 412   END IF
 413
 414 END SUBROUTINE da_buddy_qc
 415