var/da/da_gpsref/da_get_innov_vector_gpsref_deprecated.inc

   1 subroutine da_get_innov_vector_gpsref(it, num_qcstat_conv, grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !    Updated for Analysis on Arakawa-C grid
   6    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   7    !-----------------------------------------------------------------------
   8
   9    implicit none
  10
  11    integer,          intent(in)    :: it       ! External iteration.
  12    type(domain),     intent(in)    :: grid     ! first guess state.
  13    type(y_type),     intent(inout) :: ob       ! Observation structure.
  14    type(iv_type),    intent(inout) :: iv       ! O-B structure.
  15    integer,          intent(inout) :: num_qcstat_conv(:,:,:,:)
  16
  17    integer :: n        ! Loop counter.
  18    integer :: i, j, k, kk  ! Index dimension.
  19    real    :: dx, dxm, dz, dzm  ! Interpolation weights.
  20    real    :: dy, dym  ! Interpolation weights.
  21    real,allocatable :: model_ref(:,:) !Model gpsref at ob loc
  22    real    :: v_h(kms:kme), v_p(kms:kme)     ! Model value h at ob
  23
  24    integer :: Iu_ref, l
  25    integer,allocatable :: used_lev(:,:) ! for obs. data thinning
  26                                         ! record the closest level with model
  27    integer,allocatable :: qc(:)         ! record iv%gpsref(n)%ref(k)%qc
  28                                ! hor. location.
  29    real    :: distance_h       ! cal. h-difference between obs and model
  30    real,allocatable :: min_dis(:)   ! minimum difference
  31                                                ! hor. location.
  32 ! t_iwabuchi 20111216 T for QC
  33    real, allocatable :: int_t(:,:,:)     ! for T
  34
  35    ! For quality control
  36
  37    real   , parameter :: h_1 = 7000.0, h_2 = 25000.0
  38    ! Lidia Cucurull values:
  39    real   , parameter :: pcnt1 = 0.05, pcnt2 = 0.04, pcnt3 = 0.10
  40    ! testing values:
  41    ! real   , parameter :: pcnt1 = 0.02, pcnt2 = 0.01, pcnt3 = 0.03
  42    integer, parameter :: qc_below = -31, qc_middle = -32, qc_above = -33
  43
  44    integer, parameter :: qc_step1 = -34, qc_step2  = -35  ! refer to Poli et al. (2009)
  45    integer :: top_level
  46    real, allocatable :: dndz_obs(:),dndz_mod(:)
  47    real, allocatable :: dndz2_obs(:),dndz2_mod(:)
  48
  49    integer :: nn, na, nb, ntotal, nqc0, nqc1, nqc2, nqc3, n_rej
  50    real    :: percnt
  51    real    :: height_below(5000)
  52    character(len=40) :: name_qc(5000)
  53    character(len=2) :: c_it
  54    character(len=7) :: c_nt
  55
  56 ! t_iwabuchii 20111216 GSI regional QC scheme
  57    real :: g_height, g_lat, cutoff, stddev     ! QC cutoff GSI
  58    integer, parameter :: qc_cutoff = -36
  59    real, allocatable :: model_t(:,:)  ! Model value t at ob location
  60
  61
  62    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_gpsref")
  63
  64    n_rej = 0
  65 !
  66    if (iv%info(gpsref)%nlocal < 1) then
  67       if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
  68       return
  69    end if
  70
  71    ntotal = 0
  72
  73    allocate (model_ref(iv%info(gpsref)%max_lev,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
  74
  75    model_ref(:,:) = 0.0
  76
  77 ! t_iwabuchi 20111216 (hmn) allocate model_t for GSI regional QC
  78    allocate (model_t(iv%info(gpsref)%max_lev,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
  79    model_t(:,:) = 0.0
  80
  81
  82    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
  83
  84       do k=1, iv%info(gpsref)%levels(n)
  85 !sychen        if( iv%gpsref(n)%ref(k)%qc == fails_error_max .and. it > 1 ) &
  86 !sychen            iv%gpsref(n)%ref(k)%qc = 0
  87             if ( ( iv%gpsref(n)%ref(k)%qc == fails_error_max ) .or. &
  88                  ( iv%gpsref(n)%ref(k)%qc == qc_below ) .or. &
  89                  ( iv%gpsref(n)%ref(k)%qc == qc_middle ) .or. &
  90                  ( iv%gpsref(n)%ref(k)%qc == qc_above ) .or. &
  91                  ( iv%gpsref(n)%ref(k)%qc == qc_step1 ) .or. &
  92                  ( iv%gpsref(n)%ref(k)%qc == qc_step2 ) .or. &
  93                  ( iv%gpsref(n)%ref(k)%qc == qc_cutoff ) ) then  !hcl-202006
  94                if( it > 1 ) iv%gpsref(n)%ref(k)%qc = 0
  95             endif
  96       end do
  97
  98       ! Get cross pt. horizontal interpolation weights:
  99
 100       i   = iv%info(gpsref)%i(1,n)
 101       j   = iv%info(gpsref)%j(1,n)
 102       dx  = iv%info(gpsref)%dx(1,n)
 103       dy  = iv%info(gpsref)%dy(1,n)
 104       dxm = iv%info(gpsref)%dxm(1,n)
 105       dym = iv%info(gpsref)%dym(1,n)
 106
 107       if ( .not. pseudo_ref ) then
 108
 109          ! Get the zk from gpsref%h:
 110
 111          do k=kts,kte
 112             v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
 113                    + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
 114 !
 115             v_p(k) = dym*(dxm*grid%xb%p(i,j  ,k) + dx*grid%xb%p(i+1,j  ,k)) &
 116                    + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
 117          end do
 118          do k=1, iv%info(gpsref)%levels(n)
 119             if (iv%gpsref(n)%h(k) > 0.0) &
 120                call da_to_zk(iv%gpsref(n)%h(k), v_h, v_interp_h, iv%info(gpsref)%zk(k,n))
 121             if (iv%info(gpsref)%zk(k,n) < 0.0 .and. .not. anal_type_verify) then
 122                iv%gpsref(n)%ref(k)%qc = missing_data
 123             end if
 124          end do
 125       else
 126          iv%info(gpsref)%zk(:,n) = pseudo_z
 127       end if
 128
 129 !
 130 ! To assign the retrieved pressure to GPSREF data (YRG, 06/15/2011)............
 131    do k=1, iv%info(gpsref)%levels(n)
 132       kk = int (iv%info(gpsref)%zk(k,n))
 133       if (kk >= kts .and. kk+1 <= kte) then
 134         dz  = iv%info(gpsref)%zk(k,n) - real(kk)
 135         dzm = 1.0 - dz
 136         iv%gpsref(n)%p(k)%inv = v_p(kk  ) * dzm + v_p(kk+1) * dz
 137         ob%gpsref(n)%p(k) = iv%gpsref(n)%p(k)%inv
 138         iv%gpsref(n)%p(k)%qc = -5
 139       else
 140         n_rej = n_rej + 1
 141       endif
 142    end do
 143 !..............................................................................
 144
 145    end do
 146
 147    call da_convert_zk (iv%info(gpsref))
 148
 149 ! t_iwabuchi 20111216 (hmn's update) linear interpolation of log(n)
 150    call da_interp_lin_3d (grid%xb%reflog, iv%info(gpsref), model_ref)
 151    allocate (int_t(ims:ime, jms:jme, kms:kme))
 152    int_t = grid%xb%t
 153    call da_interp_lin_3d (int_t, iv%info(gpsref), model_t)
 154    deallocate(int_t)
 155 ! t_iwabuchi END
 156
 157
 158    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 159
 160 ! t_iwabuchi 20111216 compute exp of log (N)
 161    do k = 1, iv%info(gpsref)%levels(n)
 162      model_ref(k,n)=exp(model_ref(k,n))
 163 !     model_t(k,n)=exp(model_t(k,n))
 164    end do
 165 ! t_iwabuchi END
 166
 167       if ( ( .not. pseudo_ref ) .or. it > 1 ) then
 168          do k = 1, iv%info(gpsref)%levels(n)
 169             iv%gpsref(n)%ref(k)%inv = 0.0
 170
 171             if (ob%gpsref(n)%ref(k) > missing_r .AND. &
 172                  !iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer) then
 173                  iv%gpsref(n)%ref(k)%qc /= missing_data) then    !hcl-202006
 174                  iv%gpsref(n)%ref(k)%inv = ob%gpsref(n)%ref(k) - model_ref(k,n)
 175             end if
 176          end do
 177       else
 178          ob % gpsref(1)%ref(1) = model_ref(1,n) + iv %gpsref(1)%ref(1)%inv
 179       end if
 180    end do
 181
 182    if ( pseudo_ref ) then
 183       ! Done for pseudo_ref after getting iv%gpsref(n)%ref(k)%inv
 184       deallocate (model_ref)
 185       deallocate (model_t)
 186       if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
 187       return
 188    end if
 189
 190    ! Quality check 1: Gross error(departure from the background) check
 191     if ( check_max_iv ) &
 192        call da_check_max_iv_gpsref(iv, it, num_qcstat_conv, 1)
 193
 194 ! refer to Poli et al. (2009) ------------------------------------------------
 195 ! flag if fit in each of these two qc steps for both of obs. and model
 196 ! qc_step1: dN/dz < -50 km^-1
 197 ! qc_step2: abs(d^2N/dz^2) > 100 km^-2
 198 !                                  Shu-Ya Chen (2009-07-29)
 199    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 200       if (iv%info(gpsref)%levels(n) <= 2) cycle
 201       do k=1,iv%info(gpsref)%levels(n)
 202          if (model_ref(k,n) > 0.0) top_level=k
 203       end do
 204       allocate(dndz_obs(top_level))
 205       allocate(dndz_mod(top_level))
 206       allocate(dndz2_obs(top_level))
 207       allocate(dndz2_mod(top_level))
 208
 209       ! QC_STEP1
 210
 211         if (.not. anal_type_verify) then
 212
 213           ! check for bottom boundary (Forward Difference)
 214                dndz_obs(1)=(ob%gpsref(n)%ref(2)-ob%gpsref(n)%ref(1))/  &
 215                            ((iv%gpsref(n)%h(2)-iv%gpsref(n)%h(1))/1000.)
 216                dndz_mod(1)=(model_ref(2,n)-model_ref(1,n))/  &
 217                            ((iv%gpsref(n)%h(2)-iv%gpsref(n)%h(1))/1000.)
 218           ! check for upper boundary (Backward Difference)
 219                dndz_obs(top_level)= &
 220                        (ob%gpsref(n)%ref(top_level)-ob%gpsref(n)%ref(top_level-1))/  &
 221                        ((iv%gpsref(n)%h(top_level)-iv%gpsref(n)%h(top_level-1))/1000.)
 222                dndz_mod(top_level)= &
 223                        (model_ref(top_level,n)-model_ref(top_level-1,n))/  &
 224                        ((iv%gpsref(n)%h(top_level)-iv%gpsref(n)%h(top_level-1))/1000.)
 225           ! check for middle levels (Central Difference)
 226             do k=2, top_level-1
 227                dndz_obs(k)=(ob%gpsref(n)%ref(k+1)-ob%gpsref(n)%ref(k-1))/  &
 228                            ((iv%gpsref(n)%h(k+1)-iv%gpsref(n)%h(k-1))/1000.)
 229                dndz_mod(k)=(model_ref(k+1,n)-model_ref(k-1,n))/  &
 230                            ((iv%gpsref(n)%h(k+1)-iv%gpsref(n)%h(k-1))/1000.)
 231             end do
 232             do k=1, top_level
 233 ! hmn 20111206
 234                if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 235                if ((dndz_obs(k) < -50.) .or. (dndz_mod(k) < -50.)) then
 236                iv%gpsref(n)%ref(k)%qc = qc_step1
 237                end if
 238 ! hmn 20111206
 239                end if
 240             end do
 241
 242       ! QC_STEP2
 243
 244           ! check for bottom boundary
 245                dndz2_obs(1)=(dndz_obs(2)-dndz_obs(1))/  &
 246                             ((iv%gpsref(n)%h(2)-iv%gpsref(n)%h(1))/1000.)
 247                dndz2_mod(1)=(dndz_mod(2)-dndz_mod(1))/   &
 248                             ((iv%gpsref(n)%h(2)-iv%gpsref(n)%h(1))/1000.)
 249           ! check for upper boundary
 250                dndz2_obs(top_level)=(dndz_obs(top_level)-dndz_obs(top_level-1))/  &
 251                             ((iv%gpsref(n)%h(top_level)-iv%gpsref(n)%h(top_level-1))/1000.)
 252                dndz2_mod(top_level)=(dndz_mod(top_level)-dndz_mod(top_level-1))/   &
 253                             ((iv%gpsref(n)%h(top_level)-iv%gpsref(n)%h(top_level-1))/1000.)
 254           ! check for middle levels
 255             do k=2, top_level-1
 256                dndz2_obs(k)=(dndz_obs(k+1)-dndz_obs(k-1))/  &
 257                             ((iv%gpsref(n)%h(k+1)-iv%gpsref(n)%h(k-1))/1000.)
 258                dndz2_mod(k)=(dndz_mod(k+1)-dndz_mod(k-1))/   &
 259                             ((iv%gpsref(n)%h(k+1)-iv%gpsref(n)%h(k-1))/1000.)
 260             end do
 261             do k=1, top_level
 262 ! hmn 20111206
 263                if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 264                if ((abs(dndz2_obs(k)) > 100.) .or. (abs(dndz2_mod(k)) > 100.)) then
 265                iv%gpsref(n)%ref(k)%qc = qc_step2
 266                end if
 267                end if
 268
 269             end do
 270
 271         end if  ! end of if verify check
 272    deallocate(dndz_obs,dndz_mod)
 273    deallocate(dndz2_obs,dndz2_mod)
 274    end do  ! end of do iv%info(gpsref)%n1~n2
 275 !
 276 ! End of Poli's check. (2009) -------------------------------------------------
 277 !
 278
 279
 280 ! t_iwabuchi 20111216 GSI's regional QC
 281 ! START GSI-QC REGIONAL TI-GSI -------------------------------------------------
 282 ! ------------------------------------------------------------------------------
 283 ! GSI-QC Implementation, Ted Iwabuchi
 284 !   First release: 2011-03-04
 285 !   Review 3.3.1 : 2011-12-16
 286 !
 287 ! HISTORY:
 288 !   2011-05-04 : bug in formula fixed
 289 !   2011-05-11 : add abs, and validated
 290 !   2011-06-06 : reviewd, added comments
 291 !   2011-12-16 : implemented for 3.3.1
 292 !
 293 ! SUMMARY:
 294 ! Modified version of Cucurull, 2010 for NCEP GSI-Global DA
 295 ! NCEP-GSI regional DA  setupref.f90
 296 !
 297 !   O-B cutoff depending on assigned std dev for each height and latitude
 298 !
 299 !  > 30                                           all observation is rejected
 300 !  11 - 30 km    0.25 + 0.5 cos (lambda)   [c1]   3
 301 !   9 - 11 km    (11-h)/2 x c2 + (h-9)/2 x c1     3
 302 !   6 -  9 km    0.5 (if T<= 240)          [c2]   3
 303 !                0.01 x T^2 - 0.455 T + 52.075 (if T>240)
 304 !   4 -  6 km    (6-h)/2 x c3 + (h-4)/2 x c2    3
 305 !   0 -  4 km    1+2.5 cos (lambda)        [c3]   1
 306 !
 307 !   h, lambda, T
 308 ! ------------------------------------------------------------------------------
 309 !
 310 ! REFERENCES:
 311 !   GSI code: comGSI_v3Beta/src/main/setupref.f90
 312 !
 313 !
 314 ! NOTES:
 315 !  See TI-GSI blocks, and don't forget to declare type of variables
 316 !  used in this block:
 317 !
 318 ! 1  real :: g_height, g_lat, cutoff, stddev     ! QC cutoff GSI
 319 ! 2  real, allocatable :: model_t(:,:)  ! Model value t at ob location
 320 ! 3  integer, parameter :: qc_cutoff = -36
 321 !
 322
 323 !  process refracticity data n from n1 to n2
 324 !
 325 if (.not. anal_type_verify) then
 326    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 327        do k=1,iv%info(gpsref)%levels(n)
 328          if (model_ref(k,n) > 0.0) top_level=k
 329
 330          g_height = iv%gpsref(n)% h(k)
 331          g_lat    = cos ( radian * iv%info(gpsref)%lat(k,n) )
 332
 333          if (     g_height >= 0.0     .and. g_height < 4000.0) then
 334            cutoff = 1.0
 335            stddev = 1.0 + 2.5 * g_lat
 336
 337          else if (g_height >= 4000.0  .and. g_height < 6000.0) then
 338            cutoff = 3.0
 339            if (model_t(k,n) <= 240.0) then
 340              stddev = 0.5
 341            else
 342              stddev = 0.001 * model_t(k,n) * model_t(k,n) - 0.455 * model_t(k,n) + 52.075
 343            endif
 344            stddev = (6000.0 - g_height)/2000.0 * (1.0 + 2.5 * g_lat)  &
 345                     + (g_height - 4000.0)/2000.0 * stddev
 346
 347          else if (g_height >= 6000.0  .and. g_height < 9000.0) then
 348            cutoff = 3.0
 349            if (model_t(k,n) <= 240.0) then
 350              stddev = 0.5
 351            else
 352 !            0.01 x T^2 - 0.455 T + 52.075 (if T>240)
 353              stddev = 0.001 * model_t(k,n) * model_t(k,n) - 0.455 * model_t(k,n) + 52.075
 354            end if
 355
 356          else if (g_height >= 9000.0  .and. g_height < 11000.0) then
 357            cutoff = 3.0
 358            if (model_t(k,n) <= 240.0) then
 359              stddev = 0.5
 360            else
 361              stddev = 0.001 * model_t(k,n) * model_t(k,n) - 0.455 * model_t(k,n) + 52.075
 362            endif
 363            stddev = (11000.0 - g_height)/2000.0 * stddev  &
 364                    + (g_height - 9000)/2000.0 * (0.25 + 0.5 * g_lat)
 365
 366
 367          else if (g_height >= 11000.0 .and. g_height < 30000.0) then
 368            cutoff = 3.0
 369            stddev = 0.25 + 0.5 * g_lat
 370
 371          else if (g_height >= 30000.0 ) then
 372            cutoff = 0.0
 373
 374          endif
 375 ! hmn 20111206
 376                if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 377 ! Check innovation, stddev, and cutoff
 378          if (abs(iv%gpsref(n)%ref(k)%inv) >= stddev * cutoff) then
 379            iv%gpsref(n)%ref(k)%qc = qc_cutoff
 380 ! hmn 20111202
 381 !           qc_36(k,n)=1
 382 ! qc flag
 383          endif
 384 ! hmn 20111206
 385                end if
 386
 387        end do
 388    end do
 389 end if
 390
 391    deallocate (model_t)
 392 !
 393 ! End   GSI-QC REGIONAL TI-GSI ------------------------------------------------
 394 !
 395 ! t_iwabuchi END
 396
 397
 398
 399
 400    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 401       ! Quality check 2: Error percentage check.
 402
 403       if (.not. anal_type_verify) then
 404             do k=1, iv%info(gpsref)%levels(n)
 405
 406                ! incremetal refractivity or the relative error:
 407                !   abs[(O-B)/{(O+B)/2}]              (Lidia Cucurull 2005)
 408
 409                ntotal = ntotal + 1
 410                percnt = 2.0 * abs(iv%gpsref(n)%ref(k)%inv / &
 411                  (ob%gpsref(n)%ref(k) + model_ref(k,n)))
 412
 413 ! hmn 20111206
 414 !              if (iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer) then
 415                if (iv%gpsref(n)%ref(k)%qc /= missing_data) then
 416
 417                   if (iv%gpsref(n)%h(k) < h_1) then
 418                      if (percnt > pcnt1) iv%gpsref(n)%ref(k)%qc = qc_below
 419                   else if (iv%gpsref(n)%h(k) > h_2) then
 420                      if (percnt > pcnt3) iv%gpsref(n)%ref(k)%qc = qc_above
 421                   else
 422                      if (percnt > pcnt2) iv%gpsref(n)%ref(k)%qc = qc_middle
 423                   end if
 424                end if
 425             end do
 426       end if  ! end of if verify check
 427    end do
 428
 429    ! Quality check 3: Low levels quality control
 430
 431    if (.not. anal_type_verify) then
 432          ! Search for the GPS RO's name with the 'qc_below':
 433
 434        if ( maxval(iv%info(gpsref)%levels(:)) > 1 ) then ! gpsref in profiles
 435          nn = 0
 436          height_below = 0.0
 437          name_qc      = '                                       '
 438
 439          do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 440              nn = nn + 1
 441              iv%info(gpsref)%levels(n) = iv%info(gpsref)%levels(n)
 442              do k=1, iv%info(gpsref)%levels(n)
 443                 if (iv%gpsref(n)%ref(k)%qc == qc_below) then
 444                    name_qc(nn) = iv%info(gpsref)%name(n)
 445                    height_below(nn) = max(iv%gpsref(n)%h(k),height_below(nn))
 446                 end if
 447              end do
 448              if (height_below(nn) == 0.0) nn = nn - 1
 449          end do
 450
 451          ! Set the flag qc_below to the levels below percnt < pcnt1::
 452
 453          ntotal = 0
 454          nqc0   = 0
 455          nqc1   = 0
 456          nqc2   = 0
 457          nqc3   = 0
 458
 459          do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 460             do na = 1,nn
 461                if (iv%info(gpsref)%name(n) == name_qc(na)) then
 462                   do k=1, iv%info(gpsref)%levels(n)
 463 ! hmn 20111202
 464                      if (iv%gpsref(n)%h(k) < height_below(na) .and. &
 465 !                        iv%gpsref(n)%ref(k)%qc >= 0) then
 466 ! hmn 20111206
 467                          iv%gpsref(n)%ref(k)%qc /= missing_data) then
 468                        iv%gpsref(n)%ref(k)%qc = qc_below
 469                      end if
 470
 471                   end do
 472                   exit
 473                end if
 474             end do
 475
 476             do k=1, iv%info(gpsref)%levels(n)
 477                ntotal = ntotal + 1
 478                if (iv%gpsref(n)%ref(k)%qc == fails_error_max) nqc0 = nqc0 + 1
 479                if (iv%gpsref(n)%ref(k)%qc == qc_middle) nqc1 = nqc1 + 1
 480                if (iv%gpsref(n)%ref(k)%qc == qc_below) nqc2 = nqc2 + 1
 481                if (iv%gpsref(n)%ref(k)%qc == qc_above) nqc3 = nqc3 + 1
 482             end do
 483          end do
 484        else    ! gpsref not in profiles
 485           do na = iv%info(gpsref)%n1, iv%info(gpsref)%n2
 486              if ( iv%gpsref(na)%ref(1)%qc == qc_below) then
 487                 do nb = iv%info(gpsref)%n1, iv%info(gpsref)%n2
 488                    if ( iv%info(gpsref)%id(nb) == iv%info(gpsref)%id(na) .and. &
 489                         iv%info(gpsref)%name(nb) == iv%info(gpsref)%name(na) .and. &
 490 ! hmn 20111206
 491 !                        iv%gpsref(nb)%ref(1)%qc >= obs_qc_pointer .and.        &
 492                         iv%gpsref(nb)%ref(1)%qc /= missing_data .and.            &
 493
 494                         iv%gpsref(nb)%h(1) <= iv%gpsref(na)%h(1)   ) then
 495                       iv%gpsref(nb)%ref(1)%qc = qc_below
 496                    end if
 497                 end do
 498              end if
 499           end do
 500        end if  ! end of if gpsref profiles
 501    end if  ! end of if verify check
 502 12221 continue
 503
 504 ! print out the amounts of obs. rejection
 505
 506     if ( check_max_iv ) &
 507        call da_check_max_iv_gpsref(iv, it, num_qcstat_conv, 2)
 508
 509 ! ------------------------------------------------------------------------------
 510 !   GPSRO thinning  (Shu-Ya Chen 20090701)
 511    if (.not. anal_type_verify) then
 512     IF ( gpsref_thinning ) THEN
 513       allocate (used_lev(kms:kme,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
 514       used_lev(:,:) = missing_data
 515
 516        DO n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 517        allocate(min_dis(kms:kme))
 518        allocate(qc(iv%info(gpsref)%levels(n)))
 519        i   = iv%info(gpsref)%i(1,n)
 520        j   = iv%info(gpsref)%j(1,n)
 521        dx  = iv%info(gpsref)%dx(1,n)
 522        dy  = iv%info(gpsref)%dy(1,n)
 523        dxm = iv%info(gpsref)%dxm(1,n)
 524        dym = iv%info(gpsref)%dym(1,n)
 525
 526        ! Get the zk from gpsref%h:
 527           do k=kts,kte
 528              v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
 529                     + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
 530           end do
 531           do k=kts,kte
 532           min_dis(k)=1.0E10
 533              do l=1, iv%info(gpsref)%levels(n)
 534                 if ( iv%gpsref(n)%ref(l)%qc >= obs_qc_pointer ) then
 535                 distance_h=abs(iv%gpsref(n)%h(l)-v_h(k))
 536                 min_dis(k)=min(min_dis(k),distance_h)
 537                 if ( min_dis(k) == distance_h ) used_lev(k,n)=l
 538                 end if
 539              end do
 540           end do
 541
 542           write(533,*) 'obs_qc_pointer=',obs_qc_pointer,'missing_data=',missing_data
 543           do k=kts,kte
 544           write(533,*) n,k,'used_lev=',used_lev(k,n)
 545           enddo
 546
 547           do l=1, iv%info(gpsref)%levels(n)
 548           write(533,*) n,l,iv%gpsref(n)%ref(l)%qc
 549           enddo
 550
 551           do l=1, iv%info(gpsref)%levels(n)
 552           qc(l)=iv%gpsref(n)%ref(l)%qc
 553           end do
 554           do k=kts,kte
 555            qc(used_lev(k,n))=1   ! which level is closest to model level
 556           end do
 557        !  data thinning (set thinned levels to be -99)
 558           do l=1, iv%info(gpsref)%levels(n)
 559           if ( iv%gpsref(n)%ref(l)%qc >= obs_qc_pointer &
 560                .and. qc(l) /= 1 ) then
 561           iv%gpsref(n)%ref(l)%qc = -99
 562           end if
 563           end do
 564        deallocate(min_dis)
 565        deallocate(qc)
 566        END DO
 567        deallocate (used_lev)
 568     END IF
 569
 570 !    goto 12345
 571
 572 ! Write out GPS Ref data:
 573
 574   if ( write_iv_gpsref ) then
 575      DO n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 576      Iu_ref = 336
 577      write(c_it,'(I2.2)') it
 578      c_nt=iv%info(gpsref)%name(n)(8:11)//iv%info(gpsref)%name(n)(28:30)
 579      open (unit=Iu_ref, file='RO_Innov_'//iv%info(gpsref)%date_char(n)//'_'//c_nt//'.'//c_it, &
 580            form='formatted')
 581            write(unit=Iu_ref, fmt='(/i5,2x,a,2x,a,2x,4f10.3,i5)') n, &
 582               iv%info(gpsref)%date_char(n), iv%info(gpsref)%id(n),  &
 583               iv%info(gpsref)%lat(1,n)  , iv%info(gpsref)%lon(1,n), &
 584               iv%info(gpsref)%x(1,n)  , iv%info(gpsref)%y(1,n), &
 585               iv%info(gpsref)%levels(n)
 586            write(unit=Iu_ref, fmt='(a5,3x,6a14)') 'level','     height   ', &
 587                        '    Obs_ref   ','  model_ref   ','  Innov_ref   ', &
 588                        '  error_ref   ',' qc_ref       '
 589            do k = 1, iv%info(gpsref)%levels(n)
 590 !             if ( gpsref_thinning ) then
 591 !               if ( iv%gpsref(n)%ref(l)%qc >= obs_qc_pointer ) then
 592 !                  write(unit=Iu_ref, fmt='(i3,1x,5f14.3,i10)')  k, &
 593 !                  iv%gpsref(n)%h(k), ob%gpsref(n)%ref(k),   &
 594 !                  model_ref(k,n), iv%gpsref(n)%ref(k)%inv,   &
 595 !                  iv%gpsref(n)%ref(k)%error, iv%gpsref(n)%ref(k)%qc
 596 !               end if
 597 !             else
 598                write(unit=Iu_ref, fmt='(i3,1x,5f14.3,i10)')  k, &
 599                   iv%gpsref(n)%h(k), ob%gpsref(n)%ref(k),   &
 600                   model_ref(k,n), iv%gpsref(n)%ref(k)%inv,   &
 601                   iv%gpsref(n)%ref(k)%error, iv%gpsref(n)%ref(k)%qc
 602 !             end if
 603            end do
 604      close(Iu_ref)
 605      END DO
 606   end if ! write_iv_gpsref
 607 !12345 continue
 608    ! .........................................................................
 609    end if  ! end of verify check
 610
 611    deallocate (model_ref)
 612
 613    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
 614
 615 end subroutine da_get_innov_vector_gpsref