var/da/da_obs_io/da_read_obs_chem_sfc.inc

   1 subroutine da_read_obs_chem_sfc (iv, filename, grid)
   2
   3    !-------------------------------------------------------------------------
   4    ! Purpose:        Read surface chem IC observation files
   5    !
   6    ! History: 03/2019  Creation (Wei Sun)
   7    !          06/2021  Add maximum thresholds for each species (Soyoung Ha)
   8    ! Caution: chem_cv_options == 10 or 20 assumes all six species in [ug/m3]
   9    !          chem_cv_options == 108 assumes pm25 and pm10 in [ug/m3], but
  10    !                                         so2, no2, o3, co in [ppmv]
  11    !-------------------------------------------------------------------------
  12
  13    implicit none
  14
  15    type (iv_type),    intent(inout) :: iv
  16    character(len=*),  intent(in)    :: filename
  17    type(domain),     intent(in)     :: grid     ! first guess state.
  18
  19    ! chem_cv_options == 108: Maximum thresholds for observation values - treat as missing above these values
  20    ! FIXME: We may want to put these maximum thresholds in namelist later.
  21    real, dimension(6)  :: max_pm = (/500.,999.,2.0, 2.0, 2.0, 50./)  ! [ug/m3, ug/m3, ppmv, ppmv, ppmv, ppmv]
  22    real                :: pm25, pm10, so2, no2, o3, co
  23
  24    character (len =   6)        :: SiteID
  25    character (len =  10)        :: fmt_name
  26
  27    character (len = 160)        :: info_string
  28
  29    integer                      :: i, j, iost, nlevels, old_nlevels, fm,iunit
  30    integer                      :: ichem, ipm
  31    real                         :: dattmp,lattmp,lontmp
  32
  33    type (singl_level_type)      :: platform
  34    logical                      :: outside
  35    logical                      :: outside_all
  36    integer                      :: surface_level
  37    real                         :: height_error, u_comp, v_comp
  38    integer                      :: nlocal(num_ob_indexes)
  39    integer                      :: ilocal(num_ob_indexes)
  40    integer                      :: ntotal(num_ob_indexes)
  41
  42    integer                      :: ndup, n, report, obs_index
  43
  44    real*8                       :: obs_time, analysis_time
  45    integer                      :: iyear, imonth, iday, ihour, imin, isite
  46    real                         :: tdiff, dlat_earth,dlon_earth,crit
  47    integer                      :: itt,itx,iout
  48    real, allocatable            :: in(:), out(:)
  49    logical                      :: found, iuse, thin_3d, is_surface
  50    integer                      :: i1,j1,k, levs
  51    real                         :: dx,dy,dxm,dym,zk
  52    real                         :: v_p(kms:kme),v_h(kms:kme)
  53
  54    if (trace_use) call da_trace_entry("da_read_obs_chem_sfc")
  55
  56  ! Initialize counts
  57    ntotal(:) = iv%info(:)%ptotal(iv%time-1)
  58    nlocal(:) = iv%info(:)%plocal(iv%time-1)
  59    ilocal    = nlocal
  60
  61    if ( thin_conv_ascii ) then
  62        do n = 1, num_ob_indexes
  63           if ( n == radar ) cycle
  64           call cleangrids_conv(n)
  65        end do
  66    end if
  67
  68    ! open file
  69    ! =========
  70
  71    call da_get_unit(iunit)
  72    open(unit   = iunit,     &
  73       FILE   = trim(filename), &
  74       FORM   = 'FORMATTED',  &
  75       ACCESS = 'SEQUENTIAL', &
  76       iostat =  iost,     &
  77       STATUS = 'OLD')
  78
  79    if (iost /= 0) then
  80       write(unit=message(1),fmt='(A,I5,A)') &
  81          "Error",iost," opening chem obs file "//trim(filename)
  82       call da_warning(__FILE__,__LINE__,message(1:1))
  83       call da_free_unit(iunit)
  84       if (trace_use) call da_trace_exit("da_read_obs_chem_sfc")
  85       return
  86    end if
  87
  88    allocate ( platform%chem (num_chemic_surf) )
  89    obs_index = chemic_surf
  90
  91
  92    imin=0     ! minutes
  93
  94    ! loop over records
  95    ! -----------------
  96
  97    report = 0 ! report number in file
  98
  99    reports: &
 100    do
 101
 102       report = report+1
 103
 104       ! Initialization
 105       ! =============================
 106       platform%chem(:)%inv = missing_r
 107       platform%chem(:)%qc  = missing
 108       platform%chem(:)%error = missing_r
 109
 110       ! read station general info
 111       ! =============================
 112
 113       if (chem_cv_options == 10) then !!!! ium data !!!!
 114             read (iunit, *, iostat =  iost)  &
 115             iyear, imonth, iday, ihour, imin, &
 116             platform%info%lat,                &
 117             platform%info%lon,                &
 118             platform%chem(PARAM_FIRST_SCALAR)%inv        !!! pm25 only !!!
 119             write (platform%info%date_char,'(i4,3(i2))')  iyear, imonth, iday, ihour
 120             platform%info%id = "333"
 121
 122          if (platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR)%inv>=500) then
 123             platform%chem(PARAM_FIRST_SCALAR)%inv=missing_r
 124          else
 125             platform%chem(PARAM_FIRST_SCALAR)%qc=333
 126          end if
 127
 128       else if (chem_cv_options == 20) then
 129          if (chemicda_opt == 1) then
 130             read (iunit, *, iostat =  iost)  &
 131             iyear, imonth, iday, ihour, imin, &
 132             platform%info%lat,                &
 133             platform%info%lon,                &
 134             platform%chem(PARAM_FIRST_SCALAR)%inv        !!! pm25 only !!!
 135          else if (chemicda_opt == 2) then
 136             read (iunit, *, iostat =  iost)  &
 137             iyear, imonth, iday, ihour, imin, &
 138             platform%info%lat,                &
 139             platform%info%lon,                &
 140             platform%chem(PARAM_FIRST_SCALAR)%inv,  &    !!! pm25 read !!!
 141             platform%chem(PARAM_FIRST_SCALAR)%inv        !!! pm10 only !!!
 142          else if (chemicda_opt == 3) then
 143             read (iunit, *, iostat =  iost)  &
 144             iyear, imonth, iday, ihour, imin, &
 145             platform%info%lat,                &
 146             platform%info%lon,                &
 147             platform%chem(PARAM_FIRST_SCALAR)%inv,  &
 148             platform%chem(PARAM_FIRST_SCALAR+1)%inv        !!! pm25 & pm10 !!!
 149          else if (chemicda_opt == 4) then
 150             read (iunit, *, iostat =  iost)  &
 151             iyear, imonth, iday, ihour, imin, &
 152             platform%info%lat,                &
 153             platform%info%lon,                &
 154             platform%chem(PARAM_FIRST_SCALAR)%inv,  &
 155             platform%chem(PARAM_FIRST_SCALAR)%inv,  &
 156             platform%chem(PARAM_FIRST_SCALAR)%inv,  &      !!! so2 !!!
 157             platform%chem(PARAM_FIRST_SCALAR+1)%inv,  &    !!! no2 !!!
 158             platform%chem(PARAM_FIRST_SCALAR+2)%inv,  &    !!! o3 !!!
 159             platform%chem(PARAM_FIRST_SCALAR+3)%inv        !!! co !!!
 160
 161             platform%chem(PARAM_FIRST_SCALAR+3)%inv = platform%chem(PARAM_FIRST_SCALAR+3)%inv*1000 !!! mg to ug !!!
 162
 163          else if (chemicda_opt == 5) then
 164             read (iunit, *, iostat =  iost)  &
 165             iyear, imonth, iday, ihour, imin, &
 166             platform%info%lat,                &
 167             platform%info%lon,                &
 168             platform%chem(PARAM_FIRST_SCALAR)%inv,  &      !!! pm25 !!!
 169             platform%chem(PARAM_FIRST_SCALAR+1)%inv,  &    !!! pm10 !!!
 170             platform%chem(PARAM_FIRST_SCALAR+2)%inv,  &    !!! so2 !!!
 171             platform%chem(PARAM_FIRST_SCALAR+3)%inv,  &    !!! no2 !!!
 172             platform%chem(PARAM_FIRST_SCALAR+4)%inv,  &    !!! o3 !!!
 173             platform%chem(PARAM_FIRST_SCALAR+5)%inv        !!! co !!!
 174
 175             platform%chem(PARAM_FIRST_SCALAR+5)%inv = platform%chem(PARAM_FIRST_SCALAR+5)%inv*1000 !!! mg to ug !!!
 176          end if
 177
 178          write (platform%info%date_char,'(i4,3(i2))')  iyear, imonth, iday, ihour
 179          platform%info%id = "333"
 180
 181          if (platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR)%inv>=500) then
 182             platform%chem(PARAM_FIRST_SCALAR)%inv=missing_r
 183          else
 184             platform%chem(PARAM_FIRST_SCALAR)%qc=333
 185          end if
 186
 187          if (chemicda_opt == 3) then
 188            if (platform%chem(PARAM_FIRST_SCALAR+1)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR+1)%inv>=500.or.platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR+1)%inv.lt.platform%chem(PARAM_FIRST_SCALAR)%inv) then
 189               platform%chem(PARAM_FIRST_SCALAR+1)%inv=missing_r
 190            else
 191               platform%chem(PARAM_FIRST_SCALAR+1)%inv=platform%chem(PARAM_FIRST_SCALAR+1)%inv-platform%chem(PARAM_FIRST_SCALAR)%inv
 192               platform%chem(PARAM_FIRST_SCALAR+1)%qc=333
 193            end if
 194
 195          else if (chemicda_opt == 4) then
 196            do ichem = PARAM_FIRST_SCALAR+1, num_chemic_surf  !!! PARAM_FIRST_SCALAR+1 to PARAM_FIRST_SCALAR+3 !!!
 197                 if (ichem.eq.num_chemic_surf) then
 198                     if (platform%chem(ichem)%inv==-999.or.platform%chem(ichem)%inv>=500*1000) then
 199                         platform%chem(ichem)%inv=missing_r
 200                     else
 201                         platform%chem(ichem)%qc=333
 202                     end if
 203                 else
 204                     if (platform%chem(ichem)%inv==-999.or.platform%chem(ichem)%inv>=500) then
 205                         platform%chem(ichem)%inv=missing_r
 206                     else
 207                         platform%chem(ichem)%qc=333
 208                     end if
 209                 end if
 210            end do
 211
 212          else if (chemicda_opt == 5) then
 213            do ichem = PARAM_FIRST_SCALAR+1, num_chemic_surf  !!! PARAM_FIRST_SCALAR+1 to PARAM_FIRST_SCALAR+5 !!!
 214                 if (ichem.eq.num_chemic_surf) then
 215                     if (platform%chem(ichem)%inv==-999.or.platform%chem(ichem)%inv>=500*1000) then
 216                         platform%chem(ichem)%inv=missing_r
 217                     else
 218                         platform%chem(ichem)%qc=333
 219                     end if
 220                 else if (ichem.eq.PARAM_FIRST_SCALAR+1) then
 221                     if (platform%chem(PARAM_FIRST_SCALAR+1)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR+1)%inv>=500.or.platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR+1)%inv.lt.platform%chem(PARAM_FIRST_SCALAR)%inv) then
 222                        platform%chem(PARAM_FIRST_SCALAR+1)%inv=missing_r
 223                     else
 224                        platform%chem(PARAM_FIRST_SCALAR+1)%inv=platform%chem(PARAM_FIRST_SCALAR+1)%inv-platform%chem(PARAM_FIRST_SCALAR)%inv
 225                        platform%chem(PARAM_FIRST_SCALAR+1)%qc=333
 226                     end if
 227                 else
 228                     if (platform%chem(ichem)%inv==-999.or.platform%chem(ichem)%inv>=500) then
 229                         platform%chem(ichem)%inv=missing_r
 230                     else
 231                         platform%chem(ichem)%qc=333
 232                     end if
 233                 end if
 234            end do
 235
 236          end if
 237
 238       else if (chem_cv_options == 108) then
 239
 240          read (iunit, fmt='(A6,F10.6,1X,F10.6,1X,I4,3I2,2F9.1,4F9.3)', iostat=iost)  &
 241                SiteID,                                   &
 242                platform%info%lat,                        &
 243                platform%info%lon,                        &
 244                iyear, imonth, iday, ihour,               &
 245                pm25,  pm10,   so2,  no2,   o3, co
 246
 247          platform%info%id = trim(SiteID)
 248          read(SiteID,'(I6)') isite
 249          write (platform%info%date_char,'(i4,3(i2.2))')  iyear, imonth, iday, ihour
 250
 251          !!!!!!!!!!!!!!!!!!!
 252          ! Caution: Observations are assumed to have the same unit of ppmv.
 253          ! Unit conversion (to ppmv) must have been done in the preprocessing.
 254          ! (ex. The original chinese data had ug/m3 for gas species.)
 255          !!!!!!!!!!!!!!!!!!!
 256
 257          select case ( chemicda_opt )
 258          case ( 1 )
 259             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 260             platform%info%name                      = "pm25"
 261          case ( 2 )
 262             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm10
 263             platform%info%name                      = "pm10"
 264          case ( 3 )
 265             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 266             platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10
 267             platform%info%name                      = "all_pm"
 268          case ( 4 )
 269             platform%chem(PARAM_FIRST_SCALAR  )%inv = so2
 270             platform%chem(PARAM_FIRST_SCALAR+1)%inv = no2
 271             platform%chem(PARAM_FIRST_SCALAR+2)%inv = o3
 272             platform%chem(PARAM_FIRST_SCALAR+3)%inv = co
 273             platform%info%name                      = "gas"
 274          case ( 5 )
 275             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 276             platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10
 277             platform%chem(PARAM_FIRST_SCALAR+2)%inv = so2
 278             platform%chem(PARAM_FIRST_SCALAR+3)%inv = no2
 279             platform%chem(PARAM_FIRST_SCALAR+4)%inv = o3
 280             platform%chem(PARAM_FIRST_SCALAR+5)%inv = co
 281             platform%info%name                      = "all"
 282          case default
 283             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 284             platform%info%name                      = ""
 285          end select
 286
 287          ! Sanity check
 288          !!!!!!!!!!!!!!!!!!!
 289          write(*,'(A,6f7.1)') &
 290          "da_read_obs_chem_sfc: pm25, pm10 [ug/m3], so2, no2, o3, co [ppmv] obs are rejected if >= ", max_pm
 291          do ichem = PARAM_FIRST_SCALAR, num_chemic_surf
 292             ipm = ichem - PARAM_FIRST_SCALAR + 1
 293             if (platform%chem(ichem)%inv<0.or.abs(platform%chem(ichem)%inv)>=max_pm(ipm)) then
 294                 platform%chem(ichem)%inv=missing_r
 295             else
 296                 platform%chem(ichem)%qc = 0
 297             endif
 298          end do
 299
 300          if ((chemicda_opt == 3) .or. (chemicda_opt == 5)) then
 301           if ((platform%chem(PARAM_FIRST_SCALAR+1)%inv.lt.platform%chem(PARAM_FIRST_SCALAR)%inv) .or. &
 302               (platform%chem(PARAM_FIRST_SCALAR+1)%inv<0).or.(platform%chem(PARAM_FIRST_SCALAR)%inv<0)) then ! Feb-27-2022
 303               platform%chem(PARAM_FIRST_SCALAR+1)%inv = missing_r
 304               platform%chem(PARAM_FIRST_SCALAR+1)%qc  = missing
 305           else
 306               ! pm10 <= pm10 - pm25 residual
 307               platform%chem(PARAM_FIRST_SCALAR+1)%inv = platform%chem(PARAM_FIRST_SCALAR+1)%inv - &
 308                                                         platform%chem(PARAM_FIRST_SCALAR)%inv
 309           endif
 310          end if ! ((chemicda_opt == 3) .or. (chemicda_opt == 5)) then
 311
 312       end if !  if (chem_cv_options == 108)
 313
 314       if (iost /= 0) then
 315          ! FIX? This is expected, but its unclear how we handle failure
 316          ! here without assuming the fortran2003 convention on
 317          ! error statuses
 318          exit reports
 319       end if
 320
 321       if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
 322       ! Fix funny wind direction at Poles
 323       if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 324          platform%info%lon = 0.0
 325       end if
 326
 327       ! read model location
 328       ! =========================
 329
 330       ! Check if outside of the time range:
 331
 332       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 333       if ( obs_time < time_slots(0) .or. &
 334            obs_time >= time_slots(num_fgat_time) ) then
 335          if (print_detail_obs) then
 336             write(unit=stdout, fmt='(a)') '*** Outside of the time range:'
 337             write(unit=stdout, fmt=fmt_info) &
 338             platform%info%platform,    &
 339             platform%info%date_char,   &
 340             platform%info%name,        &
 341             platform%info%levels,      &
 342             platform%info%lat,         &
 343             platform%info%lon,         &
 344             platform%info%elv,         &
 345             platform%info%id
 346          end if
 347          cycle
 348       endif
 349
 350       ! Restrict to a range of reports, useful for debugging
 351
 352       if (report < report_start) then
 353          cycle
 354       end if
 355
 356       if (report > report_end) then
 357          exit
 358       end if
 359
 360       call da_llxy (platform%info, platform%loc, outside, outside_all)
 361
 362       if (outside_all) then
 363          cycle reports
 364       end if
 365
 366       if (print_detail_obs) then
 367          ! Simplistic approach, could be improved to get it all done on PE 0
 368          if (.NOT. outside) then
 369              write(unit=stdout,fmt='(A,I5,A,A,3X,2F8.2,2I4,2F7.1,4F9.3)') &
 370              "Report",report," at ",trim(platform%info%id),platform%info%lat,platform%info%lon, &
 371              platform%loc%j,  platform%loc%i, pm25, pm10, so2, no2, o3, co
 372             !write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') &
 373             !   "Report",report," at",platform%info%lon,"E",platform%info%lat, &
 374             !   "N on processor", myproc,", position", platform%loc%x,platform%loc%y
 375          end if
 376       end if
 377
 378       call da_get_julian_time (iyear,imonth,iday,ihour,imin,analysis_time)
 379       tdiff = abs((obs_time - analysis_time)-0.02)
 380       dlat_earth = platform%info%lat
 381       dlon_earth = platform%info%lon
 382       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 383       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 384       dlat_earth = dlat_earth * deg2rad
 385       dlon_earth = dlon_earth * deg2rad
 386
 387       levs = 1
 388
 389       ! Loop over duplicating obs for global
 390       ndup = 1
 391       if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
 392
 393       ! It is possible that logic for counting obs is incorrect for the
 394       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
 395       ! TBH:  20050913
 396
 397       if (.not.outside) then
 398          if (print_detail_obs .and. ndup > 1) then
 399             write(unit=stdout, fmt = fmt_info) &
 400                platform%info%platform,    &
 401                platform%info%date_char,   &
 402                platform%info%name,        &
 403                platform%info%levels,      &
 404                platform%info%lat,         &
 405                platform%info%lon,         &
 406                platform%info%elv,         &
 407                platform%info%id
 408
 409             write(unit=stdout, fmt = '(a,2i5,4f10.3)') &
 410                ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
 411                platform%loc%i,  platform%loc%j,   &
 412                platform%loc%dx, platform%loc%dxm, &
 413               platform%loc%dy, platform%loc%dym
 414          end if
 415       end if
 416
 417       dup_loop: do n = 1, ndup
 418          is_surface=.true.
 419          if (use_chemic_surfobs) then
 420             if ( ntotal(chemic_surf) == max_synop_input ) cycle reports
 421             if (n==1) ntotal(chemic_surf) = ntotal(chemic_surf)+1
 422             if (outside) cycle reports
 423             if ( thin_conv_ascii ) then
 424                crit = tdiff
 425                call map2grids_conv(chemic_surf,dlat_earth,dlon_earth,crit,nlocal(chemic_surf),itx,1,itt,ilocal(chemic_surf),iuse)
 426                if ( .not. iuse ) cycle reports
 427             else
 428                nlocal(chemic_surf) = nlocal(chemic_surf) + 1
 429                ilocal(chemic_surf) = ilocal(chemic_surf) + 1
 430             end if
 431
 432             allocate ( iv%chemic_surf(ilocal(chemic_surf))%chem(num_chemic_surf))
 433
 434             do ichem = PARAM_FIRST_SCALAR, num_chemic_surf
 435                 iv%chemic_surf(ilocal(chemic_surf))%chem(ichem) = platform%chem(ichem)
 436             end do
 437
 438          else
 439
 440             write(unit=message(1), fmt='(a)') 'unsaved obs found:'
 441             write(unit=message(2), fmt='(2a)') &
 442                'platform%info%platform=', platform%info%platform
 443             write(unit=message(3), fmt='(a, i3)') &
 444                'platform%info%levels=', platform%info%levels
 445             call da_warning(__FILE__,__LINE__,message(1:3))
 446             cycle
 447
 448          end if
 449
 450          if( is_surface .or. (obs_index == gpspw) .or. (levs > 0 .and. .not. thin_conv_ascii) .or. &
 451             (levs > 0 .and. (thin_conv_ascii .and. (obs_index /=  airep .and. obs_index /= tamdar))) ) then
 452             iv%info(obs_index)%name(ilocal(obs_index))      = platform%info%name
 453             iv%info(obs_index)%platform(ilocal(obs_index))  = platform%info%platform
 454             iv%info(obs_index)%id(ilocal(obs_index))        = platform%info%id
 455             iv%info(obs_index)%date_char(ilocal(obs_index)) = platform%info%date_char
 456             ! nlevels adjusted for some obs types so use that
 457             iv%info(obs_index)%levels(ilocal(obs_index))    = min(iv%info(obs_index)%max_lev, levs)
 458             iv%info(obs_index)%lat(:,ilocal(obs_index))     = platform%info%lat
 459             iv%info(obs_index)%lon(:,ilocal(obs_index))     = platform%info%lon
 460             iv%info(obs_index)%elv(ilocal(obs_index))       = platform%info%elv
 461             iv%info(obs_index)%pstar(ilocal(obs_index))     = platform%info%pstar
 462
 463             iv%info(obs_index)%slp(ilocal(obs_index))           = platform%loc%slp
 464             iv%info(obs_index)%pw(ilocal(obs_index))            = platform%loc%pw
 465             iv%info(obs_index)%x(:,ilocal(obs_index))           = platform%loc%x
 466             iv%info(obs_index)%y(:,ilocal(obs_index))           = platform%loc%y
 467             iv%info(obs_index)%i(:,ilocal(obs_index))           = platform%loc%i
 468             iv%info(obs_index)%j(:,ilocal(obs_index))           = platform%loc%j
 469             iv%info(obs_index)%dx(:,ilocal(obs_index))          = platform%loc%dx
 470             iv%info(obs_index)%dxm(:,ilocal(obs_index))         = platform%loc%dxm
 471             iv%info(obs_index)%dy(:,ilocal(obs_index))          = platform%loc%dy
 472             iv%info(obs_index)%dym(:,ilocal(obs_index))         = platform%loc%dym
 473             iv%info(obs_index)%proc_domain(:,ilocal(obs_index)) = platform%loc%proc_domain
 474
 475             iv%info(obs_index)%obs_global_index(nlocal(obs_index)) = ntotal(obs_index)
 476
 477          end if  ! for thin_conv_ascii skipping obs_index for which thin_3d is true like airep and tamdir
 478
 479          if (global .and. n < 2) then
 480             if (test_transforms) exit dup_loop
 481             if (platform%loc % i >= ide) then
 482                platform%loc%i = platform%loc % i - ide
 483             else if (platform%loc % i < ids) then
 484                platform%loc%i = platform%loc % i + ide
 485             end if
 486
 487             platform%loc%proc_domain = .not. platform%loc%proc_domain
 488          end if
 489       end do dup_loop
 490
 491    end do reports
 492
 493    close(iunit)
 494
 495    call da_free_unit(iunit)
 496
 497    ! thinning check
 498    if ( thin_conv_ascii ) then
 499       print*, 'da_read_obs_chem_sfc: thin_conv_ascii ',thin_conv_ascii
 500       do n = 1, num_ob_indexes
 501           if (n==radar) cycle
 502           allocate ( in(thinning_grid_conv(n)%itxmax) )
 503           allocate (out(thinning_grid_conv(n)%itxmax) )
 504             do i = 1, thinning_grid_conv(n)%itxmax
 505                in(i) = thinning_grid_conv(n)%score_crit(i)
 506             end do
 507 #ifdef DM_PARALLEL
 508             ! Get minimum crit and associated processor index.
 509             call mpi_reduce(in, out, thinning_grid_conv(n)%itxmax, true_mpi_real, mpi_min, root, comm, ierr)
 510             call wrf_dm_bcast_real (out, thinning_grid_conv(n)%itxmax)
 511 #else
 512             out = in
 513 #endif
 514             do i = 1, thinning_grid_conv(n)%itxmax
 515               if( out(i) < 9.99e6) iv%info(n)%thin_ntotal=  iv%info(n)%thin_ntotal + 1
 516                if ( abs(out(i)-thinning_grid_conv(n)%score_crit(i)) > 1.0E-10 ) then
 517                   thinning_grid_conv(n)%ibest_obs(i) = 0
 518                end if
 519             end do
 520 !            do j = iv%info(n)%plocal(iv%time -1)+1 , iv%info(n)%plocal(iv%time -1)+nlocal(n)
 521             do j = iv%info(n)%plocal(iv%time -1)+1 , nlocal(n)
 522                found = .false.
 523                do i = 1, thinning_grid_conv(n)%itxmax
 524                   if ( thinning_grid_conv(n)%ibest_obs(i) == j .and.         &
 525                        thinning_grid_conv(n)%score_crit(i) < 9.99e6 ) then
 526                    iv%info(n)%thin_nlocal =  iv%info(n)%thin_nlocal + 1
 527                      found = .true.
 528
 529                      exit
 530                   end if
 531                end do
 532                if ( .not. found ) then
 533                   iv%info(n)%thinned(:,j) = .true.
 534                end if
 535             end do
 536          deallocate( in  )
 537          deallocate( out )
 538       end do ! loop over num_ob_indexes
 539    end if  ! thin_conv_ascii
 540
 541    if (trace_use) call da_trace_exit("da_read_obs_chem_sfc")
 542
 543 end subroutine da_read_obs_chem_sfc
 544
 545