var/da/da_gpsref/da_qc_gpsref.inc

   1 subroutine da_qc_gpsref(it, grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: quality control for gpsro refractivity
   5    !          this subroutine is only called when anal_type_verify=false
   6    ! History:
   7    !   2020-06: QC code previously in da_get_innov_vector_gpsref
   8    !-----------------------------------------------------------------------
   9
  10    implicit none
  11
  12    integer,          intent(in)    :: it       ! External iteration
  13    type(domain),     intent(in)    :: grid     ! first guess state.
  14    type(y_type),     intent(inout) :: ob       ! Observation structure.
  15    type(iv_type),    intent(inout) :: iv       ! O-B structure.
  16
  17    integer :: n        ! Loop counter
  18    integer :: i, j, k  ! Index dimension
  19    real, allocatable :: model_ref(:,:) ! model gpsref at ob loc
  20    real, allocatable :: zdiff(:)       ! dz for dndz qc
  21    real, allocatable :: model_t(:,:)   ! model value t at ob location
  22    real, allocatable :: int_t(:,:,:)   ! xb%t for interpolating t for gsi qc
  23
  24    real :: max_iv_thresh
  25
  26    integer :: top_level
  27    real, allocatable :: dndz_obs(:),dndz_mod(:)
  28    real, allocatable :: dndz2_obs(:),dndz2_mod(:)
  29
  30    real :: h_1, h_2
  31    real :: pcnt1, pcnt2, pcnt3
  32    integer :: nn, na, nb
  33    real    :: percnt
  34    real    :: height_below(5000)
  35    character(len=40) :: name_qc(5000)
  36
  37    ! for QC cutoff GSI
  38    real :: g_height, g_lat, cutoff, stddev
  39
  40    ! for QC statistics
  41    integer :: nrej_maxiv, nrej_dndz, nrej_dndz2, nrej_cutoff
  42    integer :: nrej_pcnt_below, nrej_pcnt_above, nrej_pcnt_middle
  43    integer, parameter   :: nhh = 13
  44    real, dimension(nhh) :: hh = (/  500.0, 1000.0, 1500.0, 2000.0, 3000.0,  &
  45                                    5000.0, 6000.0, 7000.0, 8000.0, 10000.0, &
  46                                   12000.0, 15000.0, 18000.0 /)
  47    integer, dimension(nhh+1) :: nrej_allqc, nrej_height, ngood, ntotal
  48    integer :: ihh, khh
  49    integer :: iunit, ios
  50    character(len=30) :: filename
  51
  52    if (trace_use_dull) call da_trace_entry("da_qc_gpsref")
  53
  54 !   if (iv%info(gpsref)%nlocal < 1) then
  55 !      if (trace_use_dull) call da_trace_exit("da_qc_gpsref")
  56 !      return
  57 !   end if
  58
  59    write(unit=message(1),fmt='(A)') '   calling da_qc_gpsref'
  60    call da_message(message(1:1))
  61
  62    nrej_maxiv       = 0
  63    nrej_dndz        = 0
  64    nrej_dndz2       = 0
  65    nrej_pcnt_below  = 0
  66    nrej_pcnt_above  = 0
  67    nrej_pcnt_middle = 0
  68    nrej_cutoff      = 0
  69    nrej_allqc(:)    = 0
  70    nrej_height(:)   = 0
  71    ngood(:)         = 0
  72    ntotal(:)        = 0
  73
  74 if ( iv%info(gpsref)%n2 >= iv%info(gpsref)%n1 ) then
  75    allocate (model_ref(iv%info(gpsref)%max_lev,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
  76    model_ref(:,:) = 0.0
  77
  78    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
  79
  80       do k = 1, iv%info(gpsref)%levels(n)
  81          if ( iv%gpsref(n)%ref(k)%qc /= missing_data ) then
  82             model_ref(k,n) = ob%gpsref(n)%ref(k) - iv%gpsref(n)%ref(k)%inv
  83          end if
  84       end do
  85
  86    end do
  87
  88 !==========================================
  89 ! Quality check 1: Gross error check
  90 !==========================================
  91    ! Quality check : Gross error(departure from the background) check
  92    if ( check_max_iv ) then
  93       do n = iv%info(gpsref)%n1, iv%info(gpsref)%n2
  94          do k = 1, iv%info(gpsref)%levels(n)
  95             max_iv_thresh = iv%gpsref(n)%ref(k)%error * max_error_ref
  96             if ( iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer ) then
  97                if ( abs(iv%gpsref(n)%ref(k)%inv) > max_iv_thresh ) then
  98                   iv%gpsref(n)%ref(k)%qc = fails_error_max
  99                   if ( iv%info(gpsref)%proc_domain(1,n) ) nrej_maxiv = nrej_maxiv + 1
 100                end if
 101             end if
 102          end do
 103       end do
 104    end if
 105
 106 !==========================================
 107 !  Quality check 2.1: dNdz
 108 !  Quality check 2.2: d2N_dz2
 109 !==========================================
 110    if ( gpsref_qc_dndz_opt==1 .or. gpsref_qc_dndz2_opt==1 ) then
 111
 112 ! refer to Poli et al. (2009) ------------------------------------------------
 113 ! flag if fit in each of these two qc steps for both of obs. and model
 114 ! qc_dndz: dN/dz < -50 km^-1
 115 ! qc_dndz2: abs(d^2N/dz^2) > 100 km^-2
 116 !                                  Shu-Ya Chen (2009-07-29)
 117    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 118       if (iv%info(gpsref)%levels(n) <= 2) cycle
 119       top_level = 0 ! initialize
 120       do k=1,iv%info(gpsref)%levels(n)
 121          if (model_ref(k,n) > 0.0) top_level=k
 122       end do
 123       if ( top_level <= 2 ) cycle
 124
 125       ! zdiff values are used for both gpsref_qc_dndz_opt and gpsref_qc_dndz2_opt
 126       allocate (zdiff(1:top_level))
 127       zdiff(1)         = iv%gpsref(n)%h(2) - iv%gpsref(n)%h(1)
 128       zdiff(top_level) = iv%gpsref(n)%h(top_level) - iv%gpsref(n)%h(top_level-1)
 129       do k = 2, top_level-1
 130          zdiff(k) = iv%gpsref(n)%h(k+1) - iv%gpsref(n)%h(k-1)
 131       end do
 132       zdiff(:) = zdiff(:) * 0.001  ! meter to km
 133
 134       ! dndz values are used for both gpsref_qc_dndz_opt and gpsref_qc_dndz2_opt calculation
 135       allocate(dndz_obs(top_level))
 136       allocate(dndz_mod(top_level))
 137
 138       ! check for bottom boundary (Forward Difference)
 139       dndz_obs(1) = (ob%gpsref(n)%ref(2)-ob%gpsref(n)%ref(1)) / zdiff(1)
 140       dndz_mod(1) = (model_ref(2,n)-model_ref(1,n)) / zdiff(1)
 141       ! check for upper boundary (Backward Difference)
 142       dndz_obs(top_level) = &
 143          (ob%gpsref(n)%ref(top_level)-ob%gpsref(n)%ref(top_level-1)) / zdiff(top_level)
 144       dndz_mod(top_level) = &
 145          (model_ref(top_level,n)-model_ref(top_level-1,n))/ zdiff(top_level)
 146       ! check for middle levels (Central Difference)
 147       do k = 2, top_level-1
 148          dndz_obs(k) = (ob%gpsref(n)%ref(k+1)-ob%gpsref(n)%ref(k-1))/zdiff(k)
 149          dndz_mod(k) = (model_ref(k+1,n)-model_ref(k-1,n))/zdiff(k)
 150       end do
 151
 152       if ( gpsref_qc_dndz_opt==1 ) then
 153          ! setting qc flags according to dndz values
 154          do k = 1, top_level
 155             if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 156                ! default gpsref_qc_dndz_thresh = -50.0
 157                if ( (dndz_obs(k) < gpsref_qc_dndz_thresh) .or. &
 158                     (dndz_mod(k) < gpsref_qc_dndz_thresh) ) then
 159                   iv%gpsref(n)%ref(k)%qc = qcflag_dndz
 160                   if ( iv%info(gpsref)%proc_domain(1,n) ) nrej_dndz = nrej_dndz + 1
 161                end if
 162             end if
 163          end do
 164       end if ! gpsref_qc_dndz_opt
 165
 166       if ( gpsref_qc_dndz2_opt==1 ) then
 167
 168          allocate(dndz2_obs(top_level))
 169          allocate(dndz2_mod(top_level))
 170
 171          ! check for bottom boundary
 172          dndz2_obs(1) = (dndz_obs(2)-dndz_obs(1)) / zdiff(1)
 173          dndz2_mod(1) = (dndz_mod(2)-dndz_mod(1)) / zdiff(1)
 174          ! check for upper boundary
 175          dndz2_obs(top_level) = (dndz_obs(top_level)-dndz_obs(top_level-1)) / zdiff(top_level)
 176          dndz2_mod(top_level) = (dndz_mod(top_level)-dndz_mod(top_level-1)) / zdiff(top_level)
 177          ! check for middle levels
 178          do k = 2, top_level-1
 179             dndz2_obs(k) = (dndz_obs(k+1)-dndz_obs(k-1)) / zdiff(k)
 180             dndz2_mod(k) = (dndz_mod(k+1)-dndz_mod(k-1)) / zdiff(k)
 181          end do
 182
 183          ! setting qc flags according to dndz2 values
 184          do k=1, top_level
 185             if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 186                ! default gpsref_qc_dndz_thresh = -100.0
 187                if ( (abs(dndz2_obs(k)) > gpsref_qc_dndz2_thresh) .or. &
 188                     (abs(dndz2_mod(k)) > gpsref_qc_dndz2_thresh) ) then
 189                   iv%gpsref(n)%ref(k)%qc = qcflag_dndz2
 190                   if ( iv%info(gpsref)%proc_domain(1,n) ) nrej_dndz2 = nrej_dndz2 + 1
 191                end if
 192             end if
 193          end do
 194          deallocate(dndz2_obs,dndz2_mod)
 195       end if ! gpsref_qc_dndz2_opt
 196
 197       deallocate(dndz_obs,dndz_mod)
 198       deallocate (zdiff)
 199    end do  ! end of do iv%info(gpsref)%n1~n2
 200    end if ! gpsref_qc_dndz_opt or gpsref_qc_dndz2_opt
 201 !
 202 ! End of Poli's check. (2009) -------------------------------------------------
 203 !
 204
 205 !==========================================
 206 ! Quality check 3: GSI's regional QC
 207 !==========================================
 208    if ( gpsref_qc_gsi_opt == 1 ) then
 209
 210 ! t_iwabuchi 20111216 GSI's regional QC
 211 ! START GSI-QC REGIONAL TI-GSI -------------------------------------------------
 212 ! ------------------------------------------------------------------------------
 213 ! GSI-QC Implementation, Ted Iwabuchi
 214 !   First release: 2011-03-04
 215 !   Review 3.3.1 : 2011-12-16
 216 !
 217 ! HISTORY:
 218 !   2011-05-04 : bug in formula fixed
 219 !   2011-05-11 : add abs, and validated
 220 !   2011-06-06 : reviewd, added comments
 221 !   2011-12-16 : implemented for 3.3.1
 222 !
 223 ! SUMMARY:
 224 ! Modified version of Cucurull, 2010 for NCEP GSI-Global DA
 225 ! NCEP-GSI regional DA  setupref.f90
 226 !
 227 !   O-B cutoff depending on assigned std dev for each height and latitude
 228 !
 229 !  > 30                                           all observation is rejected
 230 !  11 - 30 km    0.25 + 0.5 cos (lambda)   [c1]   3
 231 !   9 - 11 km    (11-h)/2 x c2 + (h-9)/2 x c1     3
 232 !   6 -  9 km    0.5 (if T<= 240)          [c2]   3
 233 !                0.01 x T^2 - 0.455 T + 52.075 (if T>240)
 234 !   4 -  6 km    (6-h)/2 x c3 + (h-4)/2 x c2    3
 235 !   0 -  4 km    1+2.5 cos (lambda)        [c3]   1
 236 !
 237 !   h, lambda, T
 238 ! ------------------------------------------------------------------------------
 239 !
 240 ! REFERENCES:
 241 !   GSI code: comGSI_v3Beta/src/main/setupref.f90
 242 !
 243 !
 244 ! NOTES:
 245 !  See TI-GSI blocks, and don't forget to declare type of variables
 246 !  used in this block:
 247 !
 248 ! 1  real :: g_height, g_lat, cutoff, stddev     ! QC cutoff GSI
 249 ! 2  real, allocatable :: model_t(:,:)  ! Model value t at ob location
 250 ! 3  integer, parameter :: qc_cutoff = -36
 251 !
 252
 253 ! t_iwabuchi 20111216 (hmn) allocate model_t for GSI regional QC
 254    allocate (model_t(iv%info(gpsref)%max_lev,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
 255    model_t(:,:) = 0.0
 256
 257    allocate (int_t(ims:ime, jms:jme, kms:kme))
 258    int_t = grid%xb%t
 259    call da_interp_lin_3d (int_t, iv%info(gpsref), model_t)
 260    deallocate(int_t)
 261
 262 !  process refracticity data n from n1 to n2
 263 !
 264    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 265        do k=1,iv%info(gpsref)%levels(n)
 266
 267          g_height = iv%gpsref(n)% h(k)
 268          g_lat    = cos ( radian * iv%info(gpsref)%lat(k,n) )
 269
 270          if (     g_height >= 0.0     .and. g_height < 4000.0) then
 271            cutoff = 1.0
 272            stddev = 1.0 + 2.5 * g_lat
 273
 274          else if (g_height >= 4000.0  .and. g_height < 6000.0) then
 275            cutoff = 3.0
 276            if (model_t(k,n) <= 240.0) then
 277              stddev = 0.5
 278            else
 279              stddev = 0.001 * model_t(k,n) * model_t(k,n) - 0.455 * model_t(k,n) + 52.075
 280            endif
 281            stddev = (6000.0 - g_height)/2000.0 * (1.0 + 2.5 * g_lat)  &
 282                     + (g_height - 4000.0)/2000.0 * stddev
 283
 284          else if (g_height >= 6000.0  .and. g_height < 9000.0) then
 285            cutoff = 3.0
 286            if (model_t(k,n) <= 240.0) then
 287              stddev = 0.5
 288            else
 289 !            0.01 x T^2 - 0.455 T + 52.075 (if T>240)
 290              stddev = 0.001 * model_t(k,n) * model_t(k,n) - 0.455 * model_t(k,n) + 52.075
 291            end if
 292
 293          else if (g_height >= 9000.0  .and. g_height < 11000.0) then
 294            cutoff = 3.0
 295            if (model_t(k,n) <= 240.0) then
 296              stddev = 0.5
 297            else
 298              stddev = 0.001 * model_t(k,n) * model_t(k,n) - 0.455 * model_t(k,n) + 52.075
 299            endif
 300            stddev = (11000.0 - g_height)/2000.0 * stddev  &
 301                    + (g_height - 9000)/2000.0 * (0.25 + 0.5 * g_lat)
 302
 303
 304          else if (g_height >= 11000.0 .and. g_height < 30000.0) then
 305            cutoff = 3.0
 306            stddev = 0.25 + 0.5 * g_lat
 307
 308          else if (g_height >= 30000.0 ) then
 309            cutoff = 0.0
 310
 311          endif
 312
 313          if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 314             ! Check innovation, stddev, and cutoff
 315             if (abs(iv%gpsref(n)%ref(k)%inv) >= stddev * cutoff) then
 316               iv%gpsref(n)%ref(k)%qc = qcflag_cutoff
 317               if ( iv%info(gpsref)%proc_domain(1,n) ) nrej_cutoff = nrej_cutoff + 1
 318             endif
 319          end if
 320
 321       end do ! level
 322    end do ! n
 323
 324    deallocate (model_t)
 325
 326    end if ! gpsref_qc_gsi_opt
 327 !
 328 ! End   GSI-QC REGIONAL TI-GSI ------------------------------------------------
 329 !
 330 ! t_iwabuchi END
 331
 332 !==========================================
 333 ! Quality check 4: Error percentage check and Low levels quality control
 334 !==========================================
 335    if ( gpsref_qc_pcnt_opt == 1 ) then
 336
 337       ! Quality check : Error percentage check.
 338
 339       ! assign namelist settings to local variables
 340       h_1 = gpsref_qc_pcnt_h1        ! default 7000.0
 341       h_2 = gpsref_qc_pcnt_h2        ! default 25000.0
 342       pcnt1 = gpsref_qc_pcnt_below   ! default 0.05  ! test 0.02
 343       pcnt2 = gpsref_qc_pcnt_above   ! default 0.04  ! test 0.01
 344       pcnt3 = gpsref_qc_pcnt_middle  ! default 0.10  ! test 0.03
 345
 346       do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 347
 348          do k=1, iv%info(gpsref)%levels(n)
 349
 350             ! incremetal refractivity or the relative error:
 351             !   abs[(O-B)/{(O+B)/2}]              (Lidia Cucurull 2005)
 352
 353             percnt = 2.0 * abs(iv%gpsref(n)%ref(k)%inv / &
 354                      (ob%gpsref(n)%ref(k) + model_ref(k,n)))
 355
 356             if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 357
 358                if (iv%gpsref(n)%h(k) < h_1) then
 359                   if (percnt > pcnt1) iv%gpsref(n)%ref(k)%qc = qcflag_pcnt_below
 360                else if (iv%gpsref(n)%h(k) > h_2) then
 361                   if (percnt > pcnt3) iv%gpsref(n)%ref(k)%qc = qcflag_pcnt_above
 362                else
 363                   if (percnt > pcnt2) iv%gpsref(n)%ref(k)%qc = qcflag_pcnt_middle
 364                end if
 365             end if
 366          end do
 367       end do
 368
 369       ! Quality check : Low levels quality control
 370
 371       ! Search for the GPS RO's name with the 'qcflag_pcnt_below':
 372
 373       if ( maxval(iv%info(gpsref)%levels(:)) > 1 ) then ! gpsref in profiles
 374          nn = 0
 375          height_below = 0.0
 376          name_qc      = '                                       '
 377
 378          do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 379              nn = nn + 1
 380              iv%info(gpsref)%levels(n) = iv%info(gpsref)%levels(n)
 381              do k=1, iv%info(gpsref)%levels(n)
 382                 if (iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_below) then
 383                    name_qc(nn) = iv%info(gpsref)%name(n)
 384                    height_below(nn) = max(iv%gpsref(n)%h(k),height_below(nn))
 385                 end if
 386              end do
 387              if (height_below(nn) == 0.0) nn = nn - 1
 388          end do
 389
 390          ! Set the flag qcflag_pcnt_below to the levels below percnt < pcnt1::
 391
 392          do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 393             do na = 1,nn
 394                if (iv%info(gpsref)%name(n) == name_qc(na)) then
 395                   do k=1, iv%info(gpsref)%levels(n)
 396                      if (iv%gpsref(n)%h(k) < height_below(na) .and. &
 397                          iv%gpsref(n)%ref(k)%qc /= missing_data) then
 398                        iv%gpsref(n)%ref(k)%qc = qcflag_pcnt_below
 399                      end if
 400                   end do
 401                   exit
 402                end if
 403             end do
 404
 405          end do
 406       else    ! gpsref not in profiles
 407          do na = iv%info(gpsref)%n1, iv%info(gpsref)%n2
 408             if ( iv%gpsref(na)%ref(1)%qc == qcflag_pcnt_below) then
 409                do nb = iv%info(gpsref)%n1, iv%info(gpsref)%n2
 410                   if ( iv%info(gpsref)%id(nb) == iv%info(gpsref)%id(na)     .and. &
 411                        iv%info(gpsref)%name(nb) == iv%info(gpsref)%name(na) .and. &
 412                        iv%gpsref(nb)%ref(1)%qc /= missing_data              .and. &
 413                        iv%gpsref(nb)%h(1) <= iv%gpsref(na)%h(1)   ) then
 414                      iv%gpsref(nb)%ref(1)%qc = qcflag_pcnt_below
 415                   end if
 416                end do
 417             end if
 418          end do
 419       end if  ! end of if gpsref profiles
 420
 421       do n = iv%info(gpsref)%n1, iv%info(gpsref)%n2
 422          if ( .not. iv%info(gpsref)%proc_domain(1,n) ) cycle
 423          do k = 1, iv%info(gpsref)%levels(n)
 424             if (iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_middle) &
 425                nrej_pcnt_middle = nrej_pcnt_middle + 1
 426             if (iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_below)  &
 427                nrej_pcnt_below = nrej_pcnt_below + 1
 428             if (iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_above)  &
 429                nrej_pcnt_above = nrej_pcnt_above + 1
 430          end do
 431       end do
 432
 433    end if ! gpsref_qc_pcnt_opt
 434 ! ------------------------------------------------------------------------------
 435
 436    deallocate (model_ref)
 437
 438
 439    ! hcl-todo: add other qc counts and prints
 440
 441    nlocal_loop: do n = iv%info(gpsref)%n1, iv%info(gpsref)%n2
 442       if ( .not. iv%info(gpsref)%proc_domain(1,n) ) cycle nlocal_loop
 443       do k = 1, iv%info(gpsref)%levels(n)
 444          ! find the index of height bins for grouping statistics
 445          ihh = 0 ! initialize
 446          if ( iv%gpsref(n)%h(k) <= hh(1) ) then
 447             ihh = 1
 448          else if ( iv%gpsref(n)%h(k) > hh(nhh) ) then
 449             ihh = nhh+1
 450          else
 451             nhh_loop: do khh = 2, nhh
 452                if ( iv%gpsref(n)%h(k) >  hh(khh-1) .and. &
 453                     iv%gpsref(n)%h(k) <= hh(khh) ) then
 454                   ihh = khh
 455                   exit nhh_loop
 456                end if
 457             end do nhh_loop
 458          end if
 459          if ( ihh < 1 ) cycle nlocal_loop ! no valid h index found
 460          if ( iv%gpsref(n)%ref(k)%qc /= missing_data ) then
 461             ntotal(ihh) = ntotal(ihh) + 1
 462             if ( iv%gpsref(n)%ref(k)%qc < obs_qc_pointer .and. &
 463                  iv%gpsref(n)%ref(k)%qc /= qcflag_height ) then
 464                nrej_allqc(ihh) = nrej_allqc(ihh) + 1
 465             end if
 466             if ( iv%gpsref(n)%ref(k)%qc == qcflag_height ) then
 467                nrej_height(ihh) = nrej_height(ihh) + 1
 468             end if
 469             if ( iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer ) then
 470                ngood(ihh) = ngood(ihh) + 1
 471             end if
 472          end if
 473       end do
 474    end do nlocal_loop
 475
 476 end if ! n2 >= n1
 477
 478    call da_proc_sum_int (nrej_maxiv)
 479    call da_proc_sum_int (nrej_dndz)
 480    call da_proc_sum_int (nrej_dndz2)
 481    call da_proc_sum_int (nrej_pcnt_below)
 482    call da_proc_sum_int (nrej_pcnt_above)
 483    call da_proc_sum_int (nrej_pcnt_middle)
 484    call da_proc_sum_int (nrej_cutoff)
 485    call da_proc_sum_ints (nrej_allqc(:))
 486    call da_proc_sum_ints (nrej_height(:))
 487    call da_proc_sum_ints (ntotal(:))
 488    call da_proc_sum_ints (ngood(:))
 489
 490    if ( rootproc ) then
 491       if (num_fgat_time > 1) then
 492          write(filename,'(a,i2.2,a,i2.2)') 'qcstat_gpsref_', it, '_t', iv%time
 493       else
 494          write(filename,'(a,i2.2)') 'qcstat_gpsref_', it
 495       end if
 496       call da_get_unit(iunit)
 497       open(iunit,file=trim(filename),form='formatted',iostat=ios)
 498       if (ios /= 0) then
 499          write(unit=message(1),fmt='(a,a)') 'error opening the output file ', filename
 500          call da_error(__FILE__,__LINE__,message(1:1))
 501       end if
 502
 503       write(iunit, fmt='(/a/)') &
 504          'Quality Control Statistics for GPSRO Refractivity'   ! exclude missing_data
 505       write(iunit,'(20x,a,13i7,a,i5,a)') '  TOTAL', nint(hh(:)), ' >', nint(hh(nhh)), ' meter'
 506       write(iunit,'(a20,20i7)') ' ntotal           = ', sum(ntotal), ntotal(:)
 507       write(iunit,'(a20,20i7)') ' ngood            = ', sum(ngood), ngood(:)
 508       write(iunit,'(a20,20i7)') ' nrej_height      = ', sum(nrej_height), nrej_height(:)
 509       write(iunit,'(a20,20i7)') ' nrej_allqc       = ', sum(nrej_allqc), nrej_allqc(:)
 510       write(iunit,'(a)')
 511       write(iunit,'(a20,i7)') ' nrej_maxiv       = ', nrej_maxiv
 512       write(iunit,'(a20,i7)') ' nrej_dndz        = ', nrej_dndz
 513       write(iunit,'(a20,i7)') ' nrej_dndz2       = ', nrej_dndz2
 514       write(iunit,'(a20,i7)') ' nrej_pcnt_below  = ', nrej_pcnt_below
 515       write(iunit,'(a20,i7)') ' nrej_pcnt_middle = ', nrej_pcnt_middle
 516       write(iunit,'(a20,i7)') ' nrej_pcnt_above  = ', nrej_pcnt_above
 517       write(iunit,'(a20,i7)') ' nrej_cutoff      = ', nrej_cutoff
 518
 519       close(iunit)
 520       call da_free_unit(iunit)
 521    end if
 522
 523    if (trace_use_dull) call da_trace_exit("da_qc_gpsref")
 524
 525 end subroutine da_qc_gpsref