var/da/da_obs_io/da_scan_obs_ascii.inc

   1 subroutine da_scan_obs_ascii (iv, filename, grid)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Scan WRFVAR GTS observation file
   5    ! Updates:
   6    !       Date: 03/19/2009 -        Y.-R. Guo
   7    !           Added the time range check when reading in observations.
   8    !       Syed RH Rizvi NCAR/NESL/MMM/DAS Date:  02/21/2013
   9    !          Updated with thinning option
  10    !
  11    !      12-2017 - Jamie Bresch
  12    !      Add a call to da_gpseph_create_ob for gpseph
  13    !---------------------------------------------------------------------------
  14
  15    implicit none
  16
  17    type (iv_type),    intent(inout) :: iv
  18    character(*),      intent(in)    :: filename
  19    type(domain),     intent(in)     :: grid     ! first guess state.
  20
  21
  22    character (len =  10)   :: fmt_name
  23    character (len = 160)   :: info_string
  24 !   character (len = 160)   :: fmt_info
  25 !   character (len = 160)   :: fmt_srfc
  26 !   character (len = 160)   :: fmt_each
  27
  28    integer                 :: i, iost, fm, report, iunit
  29    type (multi_level_type) :: platform
  30    logical                 :: outside, outside_all
  31    real                    :: height_error
  32    integer                 :: ndup, n, obs_index
  33
  34    real*8                :: obs_time, analysis_time
  35    integer               :: iyear, imonth, iday, ihour, imin
  36    real                  :: tdiff, dlat_earth,dlon_earth,crit
  37    integer               :: itt,itx,iout
  38    logical               :: iuse, thin_3d, is_surface
  39    integer               :: i1,j1,k, nlevels, levs
  40    real                  :: dx,dy,dxm,dym,zk
  41    real                  :: v_p(kms:kme),v_h(kms:kme)
  42
  43    ! for gpseph
  44    type(ob_in_mean_h) :: pseudo_ob
  45    integer            :: lowest_level
  46
  47    if (trace_use) call da_trace_entry("da_scan_obs_ascii")
  48 ! Initialize
  49       if ( thin_conv_ascii ) then
  50           do n = 1, num_ob_indexes
  51              if ( .not. allocated(thinning_grid_conv(n)%icount) ) cycle
  52              if ( n == radar ) cycle
  53              call cleangrids_conv(n)
  54           end do
  55       end if
  56    ! open file
  57    ! ---------
  58    call da_get_unit(iunit)
  59    open(unit   = iunit,     &
  60       FILE   = trim(filename), &
  61       FORM   = 'FORMATTED',  &
  62       ACCESS = 'SEQUENTIAL', &
  63       iostat =  iost,     &
  64       STATUS = 'OLD')
  65
  66    if (iost /= 0) then
  67       write(unit=message(1),fmt='(A,I5,A)') &
  68          "Error",iost," opening gts obs file "//trim(filename)
  69       call da_warning(__FILE__,__LINE__,message(1:1))
  70       call da_free_unit(iunit)
  71       if (trace_use) call da_trace_exit("da_scan_obs_ascii")
  72       return
  73    end if
  74
  75    if ( use_gpsephobs ) then
  76       allocate (pseudo_ob%ref(kds:kde))
  77       allocate (pseudo_ob%lat(kds:kde))
  78       allocate (pseudo_ob%lon(kds:kde))
  79       allocate (pseudo_ob%azim(kds:kde))
  80    end if
  81
  82    ! read header
  83
  84    head_info: do
  85       read (unit = iunit, fmt = '(A)', iostat = iost) info_string
  86       if (iost /= 0) then
  87          write(unit=message(1),fmt='(A,I3,A,I3)') &
  88             "Error",iost,"reading gts obs header on unit",iunit
  89          call da_warning(__FILE__,__LINE__,message(1:1))
  90       if (trace_use) call da_trace_exit("da_scan_obs_ascii")
  91          return
  92       end if
  93       if (info_string(1:6) == 'EACH  ') exit
  94    end do head_info
  95
  96    ! read formats
  97
  98    read (iunit, fmt = '(A,1X,A)', iostat = iost) &
  99        fmt_name, fmt_info, &
 100        fmt_name, fmt_srfc,  &
 101        fmt_name, fmt_each
 102
 103    if (iost /= 0) then
 104       write(unit=message(1),fmt='(A,I3,A,I3)') &
 105          "Error",iost,"reading gts obs formats on unit",iunit
 106          call da_warning(__FILE__,__LINE__,message(1:1))
 107       if (trace_use) call da_trace_exit("da_scan_obs_ascii")
 108       return
 109    end if
 110
 111    ! skip units
 112    read (iunit, fmt = '(A)') fmt_name
 113
 114    ! loop over records
 115
 116    report = 0 ! report number in file
 117
 118    reports: do
 119       report = report+1
 120
 121       ! read station general info
 122
 123       read (iunit, fmt = fmt_info, iostat = iost) &
 124          platform%info%platform,    &
 125          platform%info%date_char,   &
 126          platform%info%name,        &
 127          platform%info%levels,      &
 128          platform%info%lat,         &
 129          platform%info%lon,         &
 130          platform%info%elv,         &
 131          platform%info%id
 132       if (iost /= 0) then
 133          ! FIX? This is expected, but its unclear how we handle failure
 134          ! here without assuming the fortran2003 convention on
 135          ! error statuses
 136          exit reports
 137       end if
 138
 139       if (print_detail_obs) then
 140          write(unit=stdout, fmt = fmt_info) &
 141             platform%info%platform,    &
 142             platform%info%date_char,   &
 143             platform%info%name,        &
 144             platform%info%levels,      &
 145             platform%info%lat,         &
 146             platform%info%lon,         &
 147             platform%info%elv,         &
 148             platform%info%id
 149       end if
 150
 151       if (platform%info%lon == 180.0) platform%info%lon =-180.000
 152       ! WHY?
 153       ! Fix funny wind direction at South Poles
 154       ! if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 155       !    platform%info%lon = 0.0
 156       ! end if
 157
 158       read (platform%info%platform(4:6), '(I3)') fm
 159
 160       ! read model location
 161       read (iunit, fmt = fmt_srfc)  &
 162          platform%loc%slp%inv, platform%loc%slp%qc, platform%loc%slp%error, &
 163          platform%loc%pw%inv, platform%loc%pw%qc, platform%loc%pw%error
 164
 165       ! read each level
 166
 167       nlevels= platform%info%levels
 168       if (nlevels > max_ob_levels) then
 169          write(unit=message(1), fmt='(4(2a,2x),a,2f8.2,2x,2(a,f9.2,2x),2(a,i0,2x))') &
 170             'Station: ', trim(platform%info%name), &
 171             'ID: ', trim(platform%info%id), &
 172             'Platform: ', trim(platform%info%platform), &
 173             'Date: ', trim(platform%info%date_char), &
 174             'At Location (lat, lon): ', platform%info%lat, platform%info%lon, &
 175             'At elevation: ', platform%info%elv, &
 176             'with pstar: ', platform%info%pstar, &
 177             'Has ', platform%info%levels, &
 178             'levels, which is greater than max_ob_levels: ', max_ob_levels
 179
 180          write (unit=message(2), fmt = '(A,1X,A,1X,A,1X,I4,2f9.3,f9.1,1X,A)') &
 181             platform%info%name,        &
 182             platform%info%platform,    &
 183             platform%info%date_char,   &
 184             platform%info%levels,      &
 185             platform%info%lat,         &
 186             platform%info%lon,         &
 187             platform%info%elv,         &
 188             platform%info%id
 189
 190          write (unit=message(3), fmt = '(A,I4,A)')'Only the first ',max_ob_levels,' levels will be used!'
 191
 192          call da_warning(__FILE__,__LINE__,message(1:3))
 193       end if
 194
 195       do i = 1, platform%info%levels
 196          if (i > max_ob_levels) then
 197             ! platform%each only has size "max_ob_levels", so if we exceed this limit
 198             ! we should just read past these lines and not assign them to platform%each
 199             read (unit = iunit, fmt = trim (fmt_each))
 200             cycle
 201          end if
 202
 203          platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
 204             field_type(missing_r, missing, missing_r, missing, missing_r), & ! u
 205             field_type(missing_r, missing, missing_r, missing, missing_r), & ! v
 206             field_type(missing_r, missing, missing_r, missing, missing_r), & ! p
 207             field_type(missing_r, missing, missing_r, missing, missing_r), & ! t
 208             field_type(missing_r, missing, missing_r, missing, missing_r), & ! q
 209             field_type(missing_r, missing, missing_r, missing, missing_r), & ! rh
 210             field_type(missing_r, missing, missing_r, missing, missing_r), & ! td
 211             field_type(missing_r, missing, missing_r, missing, missing_r))  ! speed
 212
 213          read (unit = iunit, fmt = trim (fmt_each)) &
 214             platform%each(i)%p%inv, platform%each(i)%p%qc, platform%each(i)%p%error, &
 215             platform%each(i)%speed%inv, platform%each(i)%speed%qc, &
 216             platform%each(i)%speed%error, &
 217             ! Here the 'direction' is stored in platform%each (i)%v:
 218             platform%each(i)%v%inv, platform%each(i)%v%qc, platform%each(i)%v%error, &
 219             platform%each(i)%height,       &
 220             platform%each(i)%height_qc,    &
 221             height_error,                   &
 222             platform%each(i)%t%inv, platform%each(i)%t%qc, platform%each(i)%t%error, &
 223             platform%each(i)%td%inv, platform%each(i)%td%qc, platform%each(i)%td%error, &
 224             platform%each(i)%rh%inv, platform%each(i)%rh%qc, platform%each(i)%rh%error
 225       end do
 226
 227       ! Check if outside of the time range:
 228
 229       read (platform%info%date_char,'(i4,4(1x,i2))') &
 230                                     iyear, imonth, iday, ihour, imin
 231       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 232       if ( obs_time < time_slots(0) .or. &
 233            obs_time >= time_slots(num_fgat_time) ) then
 234            cycle
 235       endif
 236
 237       ! Restrict to a range of reports, useful for debugging
 238       if (report < report_start) cycle
 239       if (report > report_end) exit
 240
 241       call da_llxy (platform%info, platform%loc, outside, outside_all)
 242       if( outside_all ) cycle reports
 243
 244       if (platform%info%levels < 1) then
 245          if (fm /= 111 .and. fm /= 114) then
 246             cycle reports
 247          end if
 248       end if
 249
 250       read (analysis_date,'(i4,4(1x,i2))') &
 251                                     iyear, imonth, iday, ihour, imin
 252       call da_get_julian_time (iyear,imonth,iday,ihour,imin,analysis_time)
 253       tdiff = abs((obs_time - analysis_time)-0.02)
 254       dlat_earth = platform%info%lat
 255       dlon_earth = platform%info%lon
 256       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 257       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 258       dlat_earth = dlat_earth * deg2rad
 259       dlon_earth = dlon_earth * deg2rad
 260
 261
 262       ! Loop over duplicating obs for global
 263       ndup = 1
 264       if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
 265
 266       if (test_transforms) ndup = 1
 267       obs_index = fm_index(fm)
 268       do n = 1, ndup
 269          select case(fm)
 270
 271          case (12) ;
 272             if (.not.use_synopobs .or. iv%info(synop)%ntotal == max_synop_input) cycle reports
 273             if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
 274             if (outside) cycle reports
 275              if ( thin_conv_opt(synop) > no_thin ) then
 276                crit = tdiff
 277                call map2grids_conv(synop,dlat_earth,dlon_earth,crit,iv%info(synop)%nlocal,itx,1,itt,iout,iuse)
 278                 if ( .not. iuse ) cycle reports
 279              else
 280                 iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
 281              end if
 282
 283          case (13, 17) ;
 284             if (.not.use_shipsobs .or. iv%info(ships)%ntotal == max_ships_input) cycle reports
 285             if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
 286             if (outside) cycle reports
 287              if ( thin_conv_opt(ships) > no_thin ) then
 288                crit = tdiff
 289                call map2grids_conv(ships,dlat_earth,dlon_earth,crit,iv%info(ships)%nlocal,itx,1,itt,iout,iuse)
 290                 if ( .not. iuse ) cycle reports
 291              else
 292                 iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
 293              end if
 294
 295
 296          case (15:16) ;
 297             if (.not.use_metarobs .or. iv%info(metar)%ntotal == max_metar_input) cycle reports
 298             if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
 299             if (outside) cycle reports
 300              if ( thin_conv_opt(metar) > no_thin ) then
 301                crit = tdiff
 302                call map2grids_conv(metar,dlat_earth,dlon_earth,crit,iv%info(metar)%nlocal,itx,1,itt,iout,iuse)
 303                 if ( .not. iuse ) cycle reports
 304              else
 305                 iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
 306              end if
 307
 308          case (32:34) ;
 309             if (.not.use_pilotobs .or. iv%info(pilot)%ntotal == max_pilot_input) cycle reports
 310             if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
 311             if (outside) cycle reports
 312              if ( thin_conv_opt(pilot) > no_thin ) then
 313                crit = tdiff
 314                call map2grids_conv(pilot,dlat_earth,dlon_earth,crit,iv%info(pilot)%nlocal,itx,1,itt,iout,iuse)
 315                 if ( .not. iuse ) cycle reports
 316              else
 317                 iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
 318              end if
 319
 320          case (35:38) ;
 321             if (.not.use_soundobs .or. iv%info(sound)%ntotal == max_sound_input) cycle reports
 322             if (n==1) iv%info(sound)%ntotal     = iv%info(sound)%ntotal + 1
 323             if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
 324             if (outside) cycle reports
 325              if ( thin_conv_opt(sound) > no_thin .or. &
 326                   thin_conv_opt(sonde_sfc) > no_thin ) then
 327                 crit = tdiff
 328                 call map2grids_conv(sound,dlat_earth,dlon_earth,crit,iv%info(sound)%nlocal,itx,1,itt,iout,iuse)
 329                 crit = tdiff
 330                 call map2grids_conv(sonde_sfc,dlat_earth,dlon_earth,crit,iv%info(sonde_sfc)%nlocal,itx,1,itt,iout,iuse)
 331                 if ( .not. iuse ) cycle reports
 332              else
 333                 iv%info(sound)%nlocal     = iv%info(sound)%nlocal + 1
 334                 iv%info(sonde_sfc)%nlocal = iv%info(sonde_sfc)%nlocal + 1
 335              end if
 336
 337          case (101) ;
 338             if (.not.use_tamdarobs .or. iv%info(tamdar)%ntotal == max_tamdar_input) cycle reports
 339
 340 ! determine if level corresponds to surface
 341             is_surface=.false.
 342             levs = nlevels
 343             do i = 1, nlevels
 344                ! if (elevation and height are the same, it is surface.)
 345                if (platform%info%elv.ne.missing_r.and.platform%each(i)%height.ne.missing_r)then
 346                   if (abs(platform%info%elv - platform%each(i)%height) < 0.1) then
 347                      is_surface = .true.
 348                      levs = nlevels - 1
 349                     exit
 350                   end if
 351                end if
 352             end do
 353
 354            if( is_surface) then
 355
 356              if (n==1) iv%info(tamdar_sfc)%ntotal = iv%info(tamdar_sfc)%ntotal + 1
 357              if (outside) cycle reports
 358              if ( thin_conv_opt(tamdar_sfc) > no_thin ) then
 359                 crit = tdiff
 360                 call map2grids_conv(tamdar_sfc,dlat_earth,dlon_earth,crit,iv%info(tamdar_sfc)%nlocal,itx,1,itt,iout,iuse)
 361                 if ( .not. iuse ) cycle reports
 362              else
 363                 iv%info(tamdar_sfc)%nlocal = iv%info(tamdar_sfc)%nlocal + 1
 364              end if
 365            else ! not is_surface
 366              if ( levs > 0 .and. n==1) iv%info(tamdar)%ntotal = iv%info(tamdar)%ntotal + 1
 367              if (outside) cycle reports
 368              if( thin_conv_opt(tamdar) <= no_thin ) then
 369                 iv%info(tamdar)%nlocal         = iv%info(tamdar)%nlocal + 1
 370              else  ! if thin_conv_opt > no_thin
 371                 thin_3d=.true.
 372                 i1   = platform%loc%i
 373                 j1   = platform%loc%j
 374                 dx   = platform%loc%dx
 375                 dy   = platform%loc%dy
 376                 dxm  = platform%loc%dxm
 377                 dym  = platform%loc%dym
 378                 do k=kms,kme
 379                    v_p(k) = dym*(dxm*grid%xb%p(i1,j1,k)+dx*grid%xb%p(i1+1,j1,k)) + &
 380                            dy*(dxm*grid%xb%p(i1,j1+1,k)+dx*grid%xb%p(i1+1,j1+1,k))
 381                 end do
 382                 do k=kms,kme
 383                    v_h(k) = dym*(dxm*grid%xb%h(i1,j1,k)+dx*grid%xb%h(i1+1,j1,k)) + &
 384                            dy*(dxm*grid%xb%h(i1,j1+1,k)+dx*grid%xb%h(i1+1,j1+1,k))
 385                 end do
 386                 do k=1,nlevels
 387                    zk = missing_r
 388                    if( platform%each(k)%p%qc  >= 0 ) then
 389                       call da_to_zk(platform%each(k)%p%inv, v_p, v_interp_p, zk)
 390                    else if( platform%each(k)%height_qc  >= 0 ) then
 391                       call da_to_zk(platform%each(k)%height, v_h, v_interp_h, zk)
 392                    else
 393                       write(unit=message(1),fmt='(A)')' For tamdar: neither height nor pressure qc is good'
 394                       call da_error(__FILE__,__LINE__,message(1:1))
 395                    end if
 396                    if ( zk == missing_r ) cycle
 397                    crit = tdiff
 398                    call map2grids_conv(tamdar,dlat_earth,dlon_earth,crit,iv%info(tamdar)%nlocal,itx,1,itt,iout,iuse,zk,thin_3d)
 399                    if ( .not. iuse ) cycle
 400                 end do ! loop over k levels
 401              end if ! if over thin_conv_ascii
 402
 403            end if ! if is_surface
 404
 405          case (161) ;
 406             if (.not.use_mtgirsobs .or. iv%info(mtgirs)%ntotal == max_mtgirs_input) cycle reports
 407             if (n==1) iv%info(mtgirs)%ntotal     = iv%info(mtgirs)%ntotal + 1
 408             if (outside) cycle reports
 409             if ( thin_conv_opt(mtgirs) > no_thin ) then
 410                crit = tdiff
 411                call map2grids_conv(mtgirs,dlat_earth,dlon_earth,crit,iv%info(mtgirs)%nlocal,itx,1,itt,iout,iuse)
 412                 if ( .not. iuse ) cycle reports
 413             else
 414                 iv%info(mtgirs)%nlocal     = iv%info(mtgirs)%nlocal + 1
 415             end if
 416
 417          case (86) ;
 418             if (.not.use_satemobs .or. iv%info(satem)%ntotal == max_satem_input) cycle reports
 419             ! Reject cloudy satem obs.
 420             if (platform%loc%pw%inv > 10.0) cycle reports
 421             if (n==1) iv%info(satem)%ntotal = iv%info(satem)%ntotal + 1
 422             if (outside) cycle reports
 423             if ( thin_conv_opt(satem) > no_thin ) then
 424                crit = tdiff
 425                call map2grids_conv(satem,dlat_earth,dlon_earth,crit,iv%info(satem)%nlocal,itx,1,itt,iout,iuse)
 426                 if ( .not. iuse ) cycle reports
 427             else
 428                 iv%info(satem)%nlocal = iv%info(satem)%nlocal + 1
 429             end if
 430
 431          case (88)    ;
 432             ! Geostationary or Polar orbitting Satellite AMVs:
 433             if (index(platform%info%name, 'MODIS') > 0 .or. &
 434                 index(platform%info%name, 'modis') > 0 .or. &
 435                 index(platform%info%id, 'AVHRR') > 0 )  then
 436                if (.not.use_polaramvobs .or. iv%info(polaramv)%ntotal == max_polaramv_input) cycle reports
 437                if (n==1) iv%info(polaramv)%ntotal = iv%info(polaramv)%ntotal + 1
 438                if (outside) cycle reports
 439                if ( thin_conv_opt(polaramv) > no_thin ) then
 440                   crit = tdiff
 441                   call map2grids_conv(polaramv,dlat_earth,dlon_earth,crit,iv%info(polaramv)%nlocal,itx,1,itt,iout,iuse)
 442                   if ( .not. iuse ) cycle reports
 443                else
 444                    iv%info(polaramv)%nlocal = iv%info(polaramv)%nlocal + 1
 445                end if
 446                obs_index = polaramv ! geoamv is the fm_index value for 88
 447             else
 448                if (.not.use_geoamvobs .or. iv%info(geoamv)%ntotal == max_geoamv_input) cycle reports
 449                if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
 450                if (outside) cycle reports
 451                if ( thin_conv_opt(geoamv) > no_thin ) then
 452                   crit = tdiff
 453                   call map2grids_conv(geoamv,dlat_earth,dlon_earth,crit,iv%info(geoamv)%nlocal,itx,1,itt,iout,iuse)
 454                    if ( .not. iuse ) cycle reports
 455                else
 456                    iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
 457                end if
 458             end if
 459
 460          case (42,96:97) ;
 461             if (.not.use_airepobs .or. iv%info(airep)%ntotal == max_airep_input) cycle reports
 462             if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
 463             if (outside) cycle reports
 464
 465             if( thin_conv_opt(airep) <= no_thin ) then
 466                iv%info(airep)%nlocal         = iv%info(airep)%nlocal + 1
 467             else  ! if thin_conv_opt > no_thin
 468                thin_3d=.true.
 469                i1   = platform%loc%i
 470                j1   = platform%loc%j
 471                dx   = platform%loc%dx
 472                dy   = platform%loc%dy
 473                dxm  = platform%loc%dxm
 474                dym  = platform%loc%dym
 475                do k=kms,kme
 476                   v_p(k) = dym*(dxm*grid%xb%p(i1,j1,k)+dx*grid%xb%p(i1+1,j1,k)) + &
 477                            dy*(dxm*grid%xb%p(i1,j1+1,k)+dx*grid%xb%p(i1+1,j1+1,k))
 478                end do
 479                do k=kms,kme
 480                   v_h(k) = dym*(dxm*grid%xb%h(i1,j1,k)+dx*grid%xb%h(i1+1,j1,k)) + &
 481                            dy*(dxm*grid%xb%h(i1,j1+1,k)+dx*grid%xb%h(i1+1,j1+1,k))
 482                end do
 483                do k=1,nlevels
 484                   zk = missing_r
 485                   if( platform%each(k)%p%qc  >= 0 ) then
 486                      call da_to_zk(platform%each(k)%p%inv, v_p, v_interp_p, zk)
 487                   else if( platform%each(k)%height_qc  >= 0 ) then
 488                      call da_to_zk(platform%each(k)%height, v_h, v_interp_h, zk)
 489                   else
 490                      write(unit=message(1),fmt='(A)')' For airep: neither height nor pressure qc is good'
 491                      call da_error(__FILE__,__LINE__,message(1:1))
 492                   end if
 493                   if ( zk == missing_r ) cycle
 494                   crit = tdiff
 495                   call map2grids_conv(airep,dlat_earth,dlon_earth,crit,iv%info(airep)%nlocal,itx,1,itt,iout,iuse,zk,thin_3d)
 496                   if ( .not. iuse ) cycle
 497                end do ! loop over k levels
 498
 499             end if ! if over thin_conv_ascii
 500
 501          case (111, 114) ;
 502             if ( (.not.use_gpspwobs  .and. fm == 111) .or. &
 503                   iv%info(gpspw)%ntotal == max_gpspw_input) cycle reports
 504             if ( (.not.use_gpsztdobs  .and. fm == 114) .or. &
 505                   iv%info(gpspw)%ntotal == max_gpspw_input) cycle reports
 506             if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
 507             if (outside) cycle reports
 508             if ( thin_conv_opt(gpspw) > no_thin ) then
 509                crit = tdiff
 510                call map2grids_conv(gpspw,dlat_earth,dlon_earth,crit,iv%info(gpspw)%nlocal,itx,1,itt,iout,iuse)
 511                if ( .not. iuse ) cycle reports
 512             else
 513                iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
 514             end if
 515
 516          case (116) ;
 517             if (.not.use_gpsrefobs .or. iv%info(gpsref)%ntotal == max_gpsref_input) cycle reports
 518             if ( ob_format_gpsro /= ob_format_ascii ) cycle reports
 519             if (n==1) iv%info(gpsref)%ntotal = iv%info(gpsref)%ntotal + 1
 520             if (outside) cycle reports
 521             if ( thin_conv_opt(gpsref) > no_thin ) then
 522                crit = tdiff
 523                call map2grids_conv(gpsref,dlat_earth,dlon_earth,crit,iv%info(gpsref)%nlocal,itx,1,itt,iout,iuse)
 524                if ( .not. iuse ) cycle reports
 525             else
 526                iv%info(gpsref)%nlocal = iv%info(gpsref)%nlocal + 1
 527             end if
 528
 529          case (118) ;
 530             if (.not.use_gpsephobs .or. iv%info(gpseph)%ntotal == max_gpseph_input) cycle reports
 531             if (n==1) iv%info(gpseph)%ntotal = iv%info(gpseph)%ntotal + 1
 532             if ( gpseph_loadbalance ) then
 533                if ( myproc /= mod( (iv%info(gpseph)%ntotal-1), num_procs ) ) cycle reports
 534             else
 535                if (outside) cycle reports
 536             endif
 537             if (gpsro_drift ==  0) then
 538                ! lat stored in speed, lon stored in v
 539                ! replacing all levels with one lat/lon for gpsro_drift=0
 540                platform%each(1:nlevels)%speed%inv = platform%info%lat
 541                platform%each(1:nlevels)%v%inv = platform%info%lon
 542             end if
 543             !create pseudo_ob on grid mean altitude for gpseph
 544             call da_gpseph_create_ob(              &
 545                nlevels,                            &
 546                platform%each(1:nlevels)%height,    &
 547                platform%each(1:nlevels)%td%inv,    & ! ref stored in td
 548                platform%each(1:nlevels)%speed%inv, & ! lat stored in speed
 549                platform%each(1:nlevels)%v%inv,     & ! lon stored in v
 550                platform%each(1:nlevels)%rh%inv,    & ! azim stored in rh
 551                pseudo_ob, lowest_level)
 552             ! cycle when no valid levels
 553             if ( lowest_level < 1 ) cycle reports
 554             iv%info(gpseph)%nlocal = iv%info(gpseph)%nlocal + 1
 555
 556           case (121) ;
 557             ! SSM/T1 temperatures
 558             if (.not.use_ssmt1obs .or. iv%info(ssmt1)%ntotal == max_ssmt1_input) cycle reports
 559             if (n==1) iv%info(ssmt1)%ntotal = iv%info(ssmt1)%ntotal + 1
 560             if (outside) cycle reports
 561             if ( thin_conv_opt(ssmt1) > no_thin ) then
 562                crit = tdiff
 563                call map2grids_conv(ssmt1,dlat_earth,dlon_earth,crit,iv%info(ssmt1)%nlocal,itx,1,itt,iout,iuse)
 564                if ( .not. iuse ) cycle reports
 565             else
 566                iv%info(ssmt1)%nlocal = iv%info(ssmt1)%nlocal + 1
 567             end if
 568
 569
 570          case (122) ;
 571             ! SSM/T2 relative humidities:
 572             if (.not.use_ssmt2obs .or. iv%info(ssmt2)%ntotal == max_ssmt2_input) cycle reports
 573             if (n==1) iv%info(ssmt2)%ntotal = iv%info(ssmt2)%ntotal + 1
 574             if (outside) cycle reports
 575             if ( thin_conv_opt(ssmt2) > no_thin ) then
 576                crit = tdiff
 577                call map2grids_conv(ssmt2,dlat_earth,dlon_earth,crit,iv%info(ssmt2)%nlocal,itx,1,itt,iout,iuse)
 578                if ( .not. iuse ) cycle reports
 579             else
 580                iv%info(ssmt2)%nlocal = iv%info(ssmt2)%nlocal + 1
 581             end if
 582
 583          case (281)    ;
 584             ! Scatterometer:
 585             if (.not.use_qscatobs .or. iv%info(qscat)%ntotal == max_qscat_input) cycle reports
 586             if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
 587             if (outside) cycle reports
 588             if ( thin_conv_opt(qscat) > no_thin ) then
 589                crit = tdiff
 590                call map2grids_conv(qscat,dlat_earth,dlon_earth,crit,iv%info(qscat)%nlocal,itx,1,itt,iout,iuse)
 591                if ( .not. iuse ) cycle reports
 592             else
 593                iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
 594             end if
 595          case (132) ;
 596             if (.not.use_profilerobs .or. iv%info(profiler)%ntotal == max_profiler_input) cycle reports
 597             if (n==1) iv%info(profiler)%ntotal = iv%info(profiler)%ntotal + 1
 598             if (outside) cycle reports
 599             if ( thin_conv_opt(profiler) > no_thin ) then
 600                crit = tdiff
 601                call map2grids_conv(profiler,dlat_earth,dlon_earth,crit,iv%info(profiler)%nlocal,itx,1,itt,iout,iuse)
 602                if ( .not. iuse ) cycle reports
 603             else
 604                iv%info(profiler)%nlocal = iv%info(profiler)%nlocal + 1
 605             end if
 606
 607          case (135) ;
 608             if (.not.use_bogusobs .or. iv%info(bogus)%ntotal == max_bogus_input) cycle reports
 609             if (n==1) iv%info(bogus)%ntotal = iv%info(bogus)%ntotal + 1
 610             if (outside) cycle reports
 611             if ( thin_conv_opt(bogus) > no_thin ) then
 612                crit = tdiff
 613                call map2grids_conv(bogus,dlat_earth,dlon_earth,crit,iv%info(bogus)%nlocal,itx,1,itt,iout,iuse)
 614                if ( .not. iuse ) cycle reports
 615             else
 616               iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
 617             end if
 618
 619          case (18,19) ;             ! buoy
 620             if (.not.use_buoyobs .or. iv%info(buoy)%ntotal == max_buoy_input) cycle reports
 621             if (n==1) iv%info(buoy)%ntotal = iv%info(buoy)%ntotal + 1
 622             if (outside) cycle reports
 623             if ( thin_conv_opt(buoy) > no_thin ) then
 624                crit = tdiff
 625                call map2grids_conv(buoy,dlat_earth,dlon_earth,crit,iv%info(buoy)%nlocal,itx,1,itt,iout,iuse)
 626                if ( .not. iuse ) cycle reports
 627             else
 628                iv%info(buoy)%nlocal = iv%info(buoy)%nlocal + 1
 629             end if
 630
 631          case (133) ;               ! AIRS retrievals
 632             if (.not.use_airsretobs .or. iv%info(airsr)%ntotal == max_airsr_input) cycle reports
 633             if (n==1) iv%info(airsr)%ntotal = iv%info(airsr)%ntotal + 1
 634             if (outside) cycle reports
 635             if ( thin_conv_opt(airsr) > no_thin ) then
 636                crit = tdiff
 637                call map2grids_conv(airsr,dlat_earth,dlon_earth,crit,iv%info(airsr)%nlocal,itx,1,itt,iout,iuse)
 638                if ( .not. iuse ) cycle reports
 639             else
 640                iv%info(airsr)%nlocal = iv%info(airsr)%nlocal + 1
 641             end if
 642
 643          case default;
 644             write(unit=message(1), fmt='(a)') 'unsaved obs found:'
 645             write(unit=message(2), fmt='(2a)') &
 646                'platform%info%platform=', platform%info%platform
 647             write(unit=message(3), fmt='(a, i3)') &
 648                  'platform%info%levels=', platform%info%levels
 649             call da_warning(__FILE__,__LINE__,message(1:3))
 650             cycle
 651          end select
 652          iv%info(obs_index)%max_lev = max(iv%info(obs_index)%max_lev, platform%info%levels)
 653       end do        !  loop over duplicate
 654    end do reports
 655
 656    if ( use_gpsephobs ) then
 657       deallocate (pseudo_ob%ref)
 658       deallocate (pseudo_ob%lat)
 659       deallocate (pseudo_ob%lon)
 660       deallocate (pseudo_ob%azim)
 661    end if
 662
 663    iv%info(sonde_sfc)%max_lev=1
 664    iv%info(tamdar_sfc)%max_lev=1
 665    iv%info(synop)%max_lev=1   ! To prevent some bad observations with more than 1 level, should I add more?
 666    iv%info(ships)%max_lev=1
 667    iv%info(qscat)%max_lev=1
 668    iv%info(metar)%max_lev=1
 669
 670    iv%info(gpseph)%max_lev= kde-kds+1
 671
 672    close(iunit)
 673
 674    call da_free_unit(iunit)
 675
 676    if (trace_use) call da_trace_exit("da_scan_obs_ascii")
 677
 678 end subroutine da_scan_obs_ascii
 679
 680