SUBROUTINE da_buddy_qc ( numobs, m_max, station_id, xob, yob, obs, qc_flag_small, & 13,2
name, pressure, dx, buddy_weight , buddy_difference, &
iunit, print_buddy, xob_g, yob_g, obs_g, qc_flag_small_g, num_recv)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Imported from BUDDY in RAWINS
! Notes:
! This routine performs error checks by comparing the difference
! between the first guess field
! and observations at the observing locations (xob, yob)
! with the averaged difference at the nearby stations within
! a given distance. If the difference value at the observation
! point differs from the average of the nearby difference values
! by more than a certain maximum (dependent on variable types,
! pressure levels, and BUDWGT parameter in namelist), the
! observation is flagged as suspect, and add the flag message
! to the observation records.
! *** Note that x and y are in natural coordinates now.
! Need to check if the x and y are used correctly...
! other variables:
! err: difference between observations and first guess fields
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! This subroutine da_buddy_qc is adapted by Y.-R.Guo, NCAR/MMM, based on
! the WRFV3/OBSGRID/src/qc2_module.F90.
!
! The algorithm for buddy check used here could be named as "Two Passes
! Innovation Buddy Check": the buddy members to a specific obs are
! determined by going through two times of checks for innovation deprture
! from the buddy member's mean values.
!
! In OBSGRID, the buddy check is completed at each of the pressure levels
! for all obs, i.e. all obs types are mixed to gether. In order to use this
! buddy check in wrfvar, first we should sort the obs into the different
! pressure bins for 3-D observations,such as SOUND, then apply the buddy
! check for each of the obs types in da_get_innov_vector_?????.inc. For the
! 2-D observations, such as SYNOP, no binning is needed.
!
! Most of the tolerances for each of variables at the different pressure
! levels are adapted from WRFV3/OBSGRID. The tolerances for specific humidity
! used here are converted from rh and specified temperature. The tolerance
! for PMSL is reduced to 350 Pa defined in da_control/da_control.f90.
!
! Yong-Run Guo, 10/10/2008
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
IMPLICIT NONE
INTEGER, intent(in) :: numobs, m_max, iunit
REAL, DIMENSION ( m_max ), intent(in) :: obs , xob , yob
INTEGER, DIMENSION ( m_max ), intent(inout) :: qc_flag_small
REAL, DIMENSION ( m_max, 0:num_procs-1 ), intent(in), optional :: obs_g , xob_g , yob_g
INTEGER, DIMENSION ( m_max, 0:num_procs-1 ), intent(inout), optional :: qc_flag_small_g
INTEGER, DIMENSION ( 0:num_procs-1 ), intent(in), optional :: num_recv
REAL, intent(in) :: buddy_weight, &
buddy_difference , &
dx , &
pressure
CHARACTER(LEN = 5), DIMENSION(m_max), intent(in) :: station_id
CHARACTER(LEN = 8), intent(in) :: name
LOGICAL , intent(in) :: print_buddy
! err: difference between observations and first guess fields
INTEGER :: num, num3, num4, numj, &
kob, n, ierr
REAL :: range, &
distance, &
sum, &
average, &
x, y, &
diff_numj, &
diff_check_1, &
diff_check_2
REAL plevel, &
pp, tt, rh, t_error, rh_error, q_error
REAL , DIMENSION ( numobs ) :: err, &
difmax, &
diff_at_obs
#ifdef DM_PARALLEL
REAL , DIMENSION ( m_max*num_procs ) :: diff
INTEGER, DIMENSION (3) :: i_dummy
#else
REAL , DIMENSION ( numobs ) :: diff
#endif
INTEGER , DIMENSION ( numobs ) :: buddy_num, &
buddy_num_final
REAL :: scale
INTEGER, DIMENSION ( m_max ) :: qc_flag
real :: f_qv_from_rh
external f_qv_from_rh
! [1.0] Store original qc flags:
qc_flag = qc_flag_small
plevel = pressure
! [2.0] Define distance within which buddy check is to find neighboring stations:
! 2.1 Define the horizontal range for buddy check:
IF ( plevel .LE. 201.0 ) THEN
range = 650.
ELSE IF ( plevel .LE. 1000.1 ) THEN
range = 300.
ELSE
range = 500.
END IF
! 2.2 Define tolerance value for comparison:
error_allowance_by_field : SELECT CASE ( name )
CASE ( 'UU ' , 'VV ' )
IF ( plevel .LT. 401. ) THEN
difmax = buddy_difference * 1.7
ELSE IF ( plevel .LT. 701. ) THEN
difmax = buddy_difference * 1.4
ELSE IF ( plevel .LT. 851. ) THEN
difmax = buddy_difference * 1.1
ELSE
difmax = buddy_difference
END IF
CASE ( 'PMSL ' )
difmax = buddy_difference
CASE ( 'TT ' )
IF ( plevel .LT. 401. ) THEN
difmax = buddy_difference * 0.6
ELSE IF ( plevel .LT. 701. ) THEN
difmax = buddy_difference * 0.8
ELSE
difmax = buddy_difference
END IF
CASE ( 'QQ ' )
! difmax = buddy_difference
! Convert the rh error to q error:
t_error = 2.0
rh_error = buddy_difference
pp = pressure * 100.0
if ( pp > 85000.0 ) then
rh = 90.0
tt = 300.0
else if ( pp > 70000.0 ) then
rh = 80.0
tt = 290.0
else if ( pp > 35000.0 ) then
rh = 70.0
tt = 280.0
else if ( pp > 10000.0 ) then
rh = 60.0
tt = 260.0
else
rh = 50.0
tt = 240.0
endif
q_error = f_qv_from_rh
( rh_error, t_error, rh, tt, pp)
difmax = q_error
! print '("p, t, rh, t_error, rh_error, difmax:",5f12.2,e15.5)', &
! pp, tt, rh, t_error, rh_error, q_error
END SELECT error_allowance_by_field
difmax = difmax * buddy_weight
! [3.0] Compute the errors for buddy check:
! 3.1 Loop through station locations (numobs) to compute buddy average
! at each observation locations.
! station_loop_2 : DO num = numobs, 1, -1
station_loop_2 : DO num = 1, numobs
average = 0.0
! No buddy check is applied to a bad data (qc < 00:
if ( qc_flag_small(num) < 0 ) cycle station_loop_2
! Compute distance between all pairs of observations,
! find number of observations within a given distance: buddy_num,
! and compute the difference between first guess and obs.
buddy_num(num) = 0
#ifdef DM_PARALLEL
if (present(xob_g) .and. present(yob_g) .and. present(obs_g) .and. present(qc_flag_small_g)) then
else
call da_error(__FILE__,__LINE__, &
(/"Buddy check in parallel mode needs global observations"/))
endif
pe_loop : DO n = 0, num_procs-1
station_loop_3 : DO num3 = 1, num_recv(n)
x = xob(num) - xob_g(num3,n)
y = yob(num) - yob_g(num3,n)
IF ( ABS(x) .GT. 1.E-20 .OR. ABS(y) .GT. 1.E-20 ) THEN
distance = SQRT ( x*x + y*y ) * dx
IF ( distance .LE. range .AND. distance .GE. 0.1 &
.and. qc_flag_small_g(num3,n) >= 0 ) THEN
! Only good data are included in buddy check group:
buddy_num(num) = buddy_num(num) + 1
! innovation (O-B) ==> diff:
diff(buddy_num(num)) = obs_g(num3,n)
END IF
END IF
END DO station_loop_3
END DO pe_loop
#else
! station_loop_3 : DO num3 = numobs, 1, -1
station_loop_3 : DO num3 = 1, numobs, 1
IF ( num3 .NE. num ) THEN
x = xob(num) - xob(num3)
y = yob(num) - yob(num3)
distance = SQRT ( x*x + y*y ) * dx
IF ( distance .LE. range .AND. distance .GE. 0.1 &
.and. qc_flag_small(num3) >= 0 ) THEN
! Only good data are included in buddy check group:
buddy_num(num) = buddy_num(num) + 1
! innovation (O-B) ==> diff:
diff(buddy_num(num)) = obs(num3)
END IF
END IF
END DO station_loop_3
#endif
! innovation at the specific obs(num):
diff_at_obs(num) = obs(num)
!
! Summation of innovations over the surrounding obs:
sum = 0.
DO numj = 1, buddy_num(num)
sum = sum + diff(numj)
END DO
! Check to see if there are any buddies, compute the mean innovation.
IF ( buddy_num(num) .GT. 0 ) average = sum / buddy_num(num)
! 3.2 Check if there is any bad obs among the obs located within the
! the radius surrounding the test ob.
diff_check_1 = difmax (1) * 1.25
diff_check_2 = difmax (1)
check_bad_ob : IF ( buddy_num(num) .GE. 2 ) THEN
kob = buddy_num(num)
remove_bad_ob : DO numj = 1, buddy_num(num)
diff_numj = ABS ( diff ( numj ) - average )
IF ( diff ( numj ) .GT. diff_check_1 &
.AND. diff_numj .GT. diff_check_2 ) THEN
! Bad obs: Innovation itself: diff(numj) > diff_check_1
! The distance between the innovation and average > diff_check_2
kob = kob - 1
sum = sum - diff ( numj )
END IF
END DO remove_bad_ob
! The final number of buddies:
buddy_num_final(num) = kob
! We may have removed too many observations.
IF ( kob .GE. 2 ) THEN
! Information for buddy check for specific obs(num)
average = sum / kob
err(num) = diff_at_obs(num) - average
ELSE
! No buddy check completed for this specific obs(num):
err(num) = 0.
! Not change the flag (YRG, 02/10/2009)
! qc_flag_small(num) = qc_flag_small(num) + no_buddies
#ifdef DM_PARALLEL
! Not change the flag (YRG, 02/10/2009)
! qc_flag_small_g(num,myproc) = qc_flag_small_g(num,myproc) + no_buddies
i_dummy(1) = qc_flag_small_g(num,myproc)
i_dummy(2) = num
i_dummy(3) = myproc
call mpi_bcast ( i_dummy, 3, mpi_integer, myproc, comm, ierr )
qc_flag_small_g(i_dummy(2),i_dummy(3)) = i_dummy(1)
#endif
END IF
ELSE check_bad_ob
! Too less buddy numbers < 2, no buddy check completed:
err(num) = 0.
! Not change the flag (YRG, 02/10/20090
! qc_flag_small(num) = qc_flag_small(num) + no_buddies
#ifdef DM_PARALLEL
! Not change the flag (YRG, 02/10/20090
! qc_flag_small_g(num,myproc) = qc_flag_small_g(num,myproc) + no_buddies
i_dummy(1) = qc_flag_small_g(num,myproc)
i_dummy(2) = num
i_dummy(3) = myproc
call mpi_bcast ( i_dummy, 3, mpi_integer, myproc, comm, ierr )
qc_flag_small_g(i_dummy(2),i_dummy(3)) = i_dummy(1)
#endif
END IF check_bad_ob
! 3.3 If the buddy number is ONLY 2, increase the tolerance value:
IF ( buddy_num_final(num) .EQ. 2 ) difmax(num) = 1.2 * difmax(num)
END DO station_loop_2
! [4.0] Reset the qc flags according to the buddy check errors:
!
! If an observation has been flagged as failing the buddy
! check criteria (error [err] larger than the allowable
! difference [difmax]), then qc_flag for this variable,
! level, time has to be modified.
station_loop_4 : do n = 1, numobs
! if the difference at station: n from the buddy-average > difmax,
! failed the buddy check:
if ( abs(err(n)) > difmax(n) ) then
qc_flag_small(n) = fails_buddy_check
endif
enddo station_loop_4
! [5.0] This notifies the user which observations have been flagged as
! suspect by the buddy check routine.
IF ( print_buddy ) THEN
n = 0
numj = 0
num3 = 0
num4 = 0
station_loop_5 : DO num = 1 , numobs
IF ( name(1:8) .EQ. 'PMSL ' ) THEN
scale = 100.
ELSE IF ( name(1:8) .EQ. 'QQ ' ) THEN
scale = 0.001
ELSE
scale = 1.0
END IF
IF ( buddy_num_final(num) .GE. 2 ) THEN
n = n + 1
IF ( ABS ( err(num) ) .GT. difmax(num) ) THEN
numj = numj + 1
WRITE ( unit=iunit, FMT = '(3I8,&
& " ID=",a5,&
& ",NAME=" ,a8, &
& ",BUDDY TOLERANCE=",f5.1,&
& ",LOC=(" ,f6.1,",",f6.1,")",&
& ",OBS INV=" ,f7.2,&
& ",DIFF ERR=" ,f7.2,&
& ",NOBS BUDDY=" ,I4,&
& ",QC_FLAG=",I3,&
& " QC_FLAG_NEW=",I3," FAILED")' ) num, n, numj,&
station_id(num),name,difmax(num)/scale, &
xob(num),yob(num),obs(num)/scale,&
err(num)/scale, buddy_num_final(num), &
qc_flag(num), qc_flag_small(num)
ELSE
num3 = num3 + 1
! WRITE ( unit=iunit, FMT = '(3I8,&
! & " ID=",a5,&
! & ",NAME=" ,a8, &
! & ",BUDDY TOLERANCE=",f5.1,&
! & ",LOC=(" ,f6.1,",",f6.1,")",&
! & ",OBS INV=" ,f7.2,&
! & ",DIFF ERR=" ,f7.2,&
! & ",NOBS BUDDY=" ,I4,&
! & ",QC_FLAG=",I3,&
! & " QC_FLAG_NEW=",I3," PASSED")' ) num, n, num3,&
! station_id(num),name,difmax(num)/scale, &
! xob(num),yob(num),obs(num)/scale,&
! err(num)/scale, buddy_num_final(num), &
! qc_flag(num), qc_flag_small(num)
ENDIF
ELSE
num4 = num4 + 1
! WRITE ( unit=iunit, FMT = '(3I8,&
! & " ID=",a5,&
! & ",NAME=" ,a8, &
! & ",BUDDY TOLERANCE=",f5.1,&
! & ",LOC=(" ,f6.1,",",f6.1,")",&
! & ",OBS INV=" ,f7.2,&
! & ",DIFF ERR=" ,f7.2,&
! & ",NOBS BUDDY=" ,I4,&
! & ",QC_FLAG=",I3,&
! & " QC_FLAG_NEW=",I3,"NO BUDDY")' ) num, n, num4,&
! station_id(num),name,difmax(num)/scale, &
! xob(num),yob(num),obs(num)/scale,&
! err(num)/scale, buddy_num_final(num), &
! qc_flag(num), qc_flag_small(num)
END IF
END DO station_loop_5
write(unit=iunit,fmt = '(5x,"NOB=",i6,&
& " Toltal N_buddy-checked =",i6,&
& " N_Passed =",i6," N_Failed=",i6," NO_Buddy=",i6, &
& " RANGE=",f10.2," DIFMAX=",f10.2)') &
numobs, n, num3, numj, num4, range, difmax(1)/scale
END IF
END SUBROUTINE da_buddy_qc