var/da/da_obs_io/da_read_obs_ascii.inc

   1 subroutine da_read_obs_ascii (iv, filename, uvq_direct, grid)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Read a GTS observation file
   5    !
   6    ! Modifications:
   7    !      10/26/2004 - Syed RH Rizvi
   8    !      Fix Obs-Long as -180 if it is +180
   9    !
  10    !      03/04/2005 - Syed RH Rizvi
  11    !      (a)  AMVs from Geostationary and Polar orbiting satellite are
  12    !           seperated & used as profile. Currently there is a provision
  13    !           to use MODIS winds only.
  14    !      (b)  MODIS obs error are currently assigned through namelist
  15    !           parameter (modis_cmv_error)
  16    !
  17    !      03/30/2005 - Syed RH Rizvi
  18    !      For global option duplicate obs at East/West boundary
  19    !
  20    !      08/30/2005 - Syed RH Rizvi
  21    !      Writing original errors & thicknesses desired for
  22    !      outputting QC obs with NTMAX = 0
  23    !
  24    !      02/28/2008 - Y.-R. Guo
  25    !      Satem thickness error should not be divided by gravity because
  26    !      the unit from obsproc is already meter, not geopotential meter.
  27    !
  28    !      03/19/2009 - Y.-R. Guo
  29    !      Added the time range check when reading in observations.
  30    !
  31    !      02/21/2013 - Syed RH Rizvi
  32    !      Updated with observation-thinning option
  33    !
  34    !      03/07/2014 - Feng Gao
  35    !      Updated wind_sd option for wind obs. types
  36    !
  37    !      12-2017 - Jamie Bresch
  38    !      Add a call to da_gpseph_create_ob for gpseph
  39    !---------------------------------------------------------------------------
  40
  41    implicit none
  42
  43    type (iv_type),    intent(inout) :: iv
  44    character(len=*),  intent(in)    :: filename
  45    type(domain),     intent(in)     :: grid     ! first guess state.
  46
  47    character (len =  10)        :: fmt_name
  48
  49    character (len = 160)        :: info_string
  50
  51 !   character (len = 160)        :: fmt_info, fmt_srfc, fmt_each
  52
  53    integer                      :: i, j, iost, nlevels, old_nlevels, fm,iunit
  54
  55    type (multi_level_type)      :: platform
  56    logical                      :: outside
  57    logical                      :: outside_all
  58    logical                      :: uvq_direct
  59    integer                      :: surface_level
  60    real                         :: height_error, u_comp, v_comp
  61    integer                      :: nlocal(num_ob_indexes)
  62    integer                      :: ilocal(num_ob_indexes)
  63    integer                      :: ntotal(num_ob_indexes)
  64
  65    integer                      :: ndup, n, report, obs_index
  66
  67    real*8                       :: obs_time, analysis_time
  68    integer                      :: iyear, imonth, iday, ihour, imin
  69    real                         :: tdiff, dlat_earth,dlon_earth,crit
  70    integer                      :: itt,itx,iout
  71    real, allocatable            :: in(:), out(:)
  72    logical                      :: found, iuse, thin_3d, is_surface
  73    integer                      :: i1,j1,k, levs
  74    real                         :: dx,dy,dxm,dym,zk
  75    real                         :: v_p(kms:kme),v_h(kms:kme)
  76
  77    ! for gpseph
  78    type(ob_in_mean_h) :: pseudo_ob
  79    integer            :: lowest_level
  80
  81    if (trace_use) call da_trace_entry("da_read_obs_ascii")
  82  ! Initialize counts
  83    ntotal(:) = iv%info(:)%ptotal(iv%time-1)
  84    nlocal(:) = iv%info(:)%plocal(iv%time-1)
  85    ilocal    = nlocal
  86    if ( thin_conv_ascii ) then
  87        do n = 1, num_ob_indexes
  88           if ( .not. allocated(thinning_grid_conv(n)%icount) ) cycle
  89           if ( n == radar ) cycle
  90           call cleangrids_conv(n)
  91        end do
  92    end if
  93    ! open file
  94    ! =========
  95
  96    call da_get_unit(iunit)
  97    open(unit   = iunit,     &
  98       FILE   = trim(filename), &
  99       FORM   = 'FORMATTED',  &
 100       ACCESS = 'SEQUENTIAL', &
 101       iostat =  iost,     &
 102       STATUS = 'OLD')
 103
 104    if (iost /= 0) then
 105       write(unit=message(1),fmt='(A,I5,A)') &
 106          "Error",iost," opening gts obs file "//trim(filename)
 107       call da_warning(__FILE__,__LINE__,message(1:1))
 108       call da_free_unit(iunit)
 109       if (trace_use) call da_trace_exit("da_read_obs_ascii")
 110       return
 111    end if
 112
 113    if ( use_gpsephobs ) then
 114       allocate (pseudo_ob%ref(kds:kde))
 115       allocate (pseudo_ob%lat(kds:kde))
 116       allocate (pseudo_ob%lon(kds:kde))
 117       allocate (pseudo_ob%azim(kds:kde))
 118    end if
 119
 120    ! read header
 121    ! ===========
 122
 123    do
 124       read (unit = iunit, fmt = '(A)', iostat = iost) info_string
 125       if (iost /= 0) then
 126          call da_warning(__FILE__,__LINE__, &
 127             (/"Problem reading gts obs header, skipping file"/))
 128          if (trace_use) call da_trace_exit("da_read_obs_ascii")
 129          return
 130       end if
 131
 132       if (info_string(1:6) == 'EACH  ') exit
 133    end do
 134
 135    !  read formats
 136    !  ------------
 137
 138    read (iunit, fmt = '(A,1X,A)', iostat = iost) &
 139       fmt_name, fmt_info, &
 140       fmt_name, fmt_srfc,  &
 141       fmt_name, fmt_each
 142
 143    if (iost /= 0) then
 144       call da_warning(__FILE__,__LINE__, &
 145          (/"Problem reading gts obs file, skipping"/))
 146          if (trace_use) call da_trace_exit("da_read_obs_ascii")
 147       return
 148    end if
 149
 150    if (print_detail_obs) then
 151       write(unit=stdout, fmt='(2a)') &
 152          'fmt_info=', fmt_info, &
 153          'fmt_srfc =', fmt_srfc,  &
 154          'fmt_each=', fmt_each
 155    end if
 156
 157    ! skip line
 158    ! ----------
 159
 160    ! read (iunit, fmt = '(a)') fmt_name
 161    read (iunit, fmt = '(a)') info_string
 162    if (info_string(1:21) == 'MODIFIED FILTERED OBS') then
 163       uvq_direct=.true.
 164    else
 165       uvq_direct=.false.
 166    endif
 167
 168    ! loop over records
 169    ! -----------------
 170
 171    report = 0 ! report number in file
 172
 173    reports: &
 174    do
 175
 176       report = report+1
 177
 178       ! read station general info
 179       ! =============================
 180
 181       read (iunit, fmt = fmt_info, iostat = iost) &
 182          platform%info%platform,    &
 183          platform%info%date_char,   &
 184          platform%info%name,        &
 185          platform%info%levels,      &
 186          platform%info%lat,         &
 187          platform%info%lon,         &
 188          platform%info%elv,         &
 189          platform%info%id
 190
 191       if (iost /= 0) then
 192          ! FIX? This is expected, but its unclear how we handle failure
 193          ! here without assuming the fortran2003 convention on
 194          ! error statuses
 195          exit reports
 196       end if
 197
 198       if (print_detail_obs) then
 199          write(unit=stdout, fmt=fmt_info) &
 200             platform%info%platform,    &
 201             platform%info%date_char,   &
 202             platform%info%name,        &
 203             platform%info%levels,      &
 204             platform%info%lat,         &
 205             platform%info%lon,         &
 206             platform%info%elv,         &
 207             platform%info%id
 208       end if
 209
 210       if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
 211       ! Fix funny wind direction at Poles
 212       if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 213          platform%info%lon = 0.0
 214       end if
 215
 216       if (platform%info%platform(6:6) == ' ') then
 217          read(platform%info%platform(4:5), '(I2)') fm
 218       else
 219          read(platform%info%platform(4:6), '(I3)') fm
 220       end if
 221
 222       ! read model location
 223       ! =========================
 224
 225       read (iunit, fmt = fmt_srfc)  &
 226          platform%loc%slp%inv, platform%loc%slp%qc, platform%loc%slp%error, &
 227          platform%loc%pw%inv, platform%loc%pw%qc, platform%loc%pw%error
 228
 229       ! read each level
 230       ! ---------------
 231
 232       do i = 1, platform%info%levels
 233          if (i > max_ob_levels) then
 234             ! platform%each only has size "max_ob_levels", so if we exceed this limit
 235             ! we should just read past these lines and not assign them to platform%each
 236             read (unit = iunit, fmt = trim (fmt_each))
 237             cycle
 238          end if
 239          platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
 240             field_type(missing_r, missing, missing_r, missing, missing_r), & ! u
 241             field_type(missing_r, missing, missing_r, missing, missing_r), & ! v
 242             field_type(missing_r, missing, missing_r, missing, missing_r), & ! p
 243             field_type(missing_r, missing, missing_r, missing, missing_r), & ! t
 244             field_type(missing_r, missing, missing_r, missing, missing_r), & ! q
 245             field_type(missing_r, missing, missing_r, missing, missing_r), & ! rh
 246             field_type(missing_r, missing, missing_r, missing, missing_r), & ! td
 247             field_type(missing_r, missing, missing_r, missing, missing_r))  ! speed
 248
 249          read (unit = iunit, fmt = trim (fmt_each)) &
 250             platform%each(i)%p%inv, platform%each(i)%p%qc, platform%each(i)%p%error, &
 251             platform%each(i)%speed%inv, platform%each(i)%speed%qc, &
 252             platform%each(i)%speed%error, &
 253             ! Here the 'direction' is stored in platform%each (i)%v:
 254             platform%each(i)%v%inv, platform%each(i)%v%qc, platform%each(i)%v%error, &
 255             platform%each(i)%height,       &
 256             platform%each(i)%height_qc,    &
 257             height_error,                  &
 258             platform%each(i)%t%inv, platform%each(i)%t%qc, platform%each(i)%t%error, &
 259             platform%each(i)%td%inv, platform%each(i)%td%qc, platform%each(i)%td%error, &
 260             platform%each(i)%rh%inv, platform%each(i)%rh%qc, platform%each(i)%rh%error
 261
 262          if (uvq_direct) then
 263             platform%each (i)%u = platform%each (i)%speed
 264             if (platform%each(i)%rh%inv .gt. missing_r)  &
 265                platform%each(i)%rh%inv=platform%each(i)%rh%inv / 1e3    !convert back to kg/kg
 266             if (platform%each(i)%rh%error .gt. 0.0)  &
 267                platform%each(i)%rh%error=platform%each(i)%rh%error / 1e3  !convert back to kg/kg
 268          end if
 269
 270          if (print_detail_obs) then
 271             write(unit=stdout, fmt = '(a, i6)') 'i=', i
 272 ! because now the field_type included: inv, qc, error, sens, and imp (YRG, 02/23/2009):
 273             write(unit=stdout, fmt = trim (fmt_each)) &
 274                platform%each(i)%p%inv, platform%each(i)%p%qc, platform%each(i)%p%error, &
 275                platform%each(i)%speed%inv, platform%each(i)%speed%qc, &
 276                platform%each(i)%speed%error,  &
 277                platform%each(i)%v%inv, platform%each(i)%v%qc, platform%each(i)%v%error, &
 278                platform%each(i)%height,       &
 279                platform%each(i)%height_qc,    &
 280                height_error,                  &
 281                platform%each(i)%t%inv, platform%each(i)%t%qc, platform%each(i)%t%error, &
 282                platform%each(i)%td%inv, platform%each(i)%td%qc, platform%each(i)%td%error, &
 283                platform%each(i)%rh%inv, platform%each(i)%rh%qc, platform%each(i)%rh%error
 284          end if
 285
 286          ! Uncomment the following if errors in reported heights (test later):
 287          ! if (fm == 97 .or. fm == 96 .or. fm == 42) then
 288          !    platform%each (i)%height_qc = -5 ! Aircraft heights are not above surface.
 289          ! end if
 290
 291          ! To convert the wind speed and direction to
 292          !      the model wind components: u, v
 293
 294          if (.not. uvq_direct) then
 295          if (fm /= 116 .and. fm /= 118 ) then !not for gpsref, gpseph
 296          if (platform%each (i)%speed%qc /= missing_data .and. &
 297              platform%each (i)%v%qc /= missing_data) then
 298
 299              call da_ffdduv (platform%each (i)%speed%inv, platform%each (i)%v%inv,     &
 300                              U_comp, V_comp, platform%info%lon, convert_fd2uv)
 301              platform%each (i)%u = platform%each (i)%speed
 302              platform%each (i)%v = platform%each (i)%v
 303              platform%each (i)%u%inv = U_comp
 304              platform%each (i)%v%inv = V_comp
 305
 306          else
 307             platform%each (i)%u%inv   = missing_r
 308             platform%each (i)%v%inv   = missing_r
 309             platform%each (i)%u%error = missing_r
 310             platform%each (i)%v%error = missing_r
 311             platform%each (i)%u%qc    = missing_data
 312             platform%each (i)%v%qc    = missing_data
 313          end if
 314          end if !not for gpsref, gpseph
 315          end if
 316       end do
 317
 318       ! Check if outside of the time range:
 319
 320       read (platform%info%date_char,'(i4,4(1x,i2))') &
 321                                     iyear, imonth, iday, ihour, imin
 322       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 323       if ( obs_time < time_slots(0) .or. &
 324            obs_time >= time_slots(num_fgat_time) ) then
 325          if (print_detail_obs) then
 326            write(unit=stdout, fmt='(a)') '*** Outside of the time range:'
 327            write(unit=stdout, fmt=fmt_info) &
 328             platform%info%platform,    &
 329             platform%info%date_char,   &
 330             platform%info%name,        &
 331             platform%info%levels,      &
 332             platform%info%lat,         &
 333             platform%info%lon,         &
 334             platform%info%elv,         &
 335             platform%info%id
 336          end if
 337          cycle
 338       endif
 339
 340       ! Restrict to a range of reports, useful for debugging
 341
 342       if (report < report_start) then
 343          cycle
 344       end if
 345
 346       if (report > report_end) then
 347          exit
 348       end if
 349
 350       call da_llxy (platform%info, platform%loc, outside, outside_all)
 351
 352       if (outside_all) then
 353          cycle reports
 354       end if
 355
 356       if (print_detail_obs) then
 357          ! Simplistic approach, could be improved to get it all done on PE 0
 358          if (.NOT. outside) then
 359             write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') &
 360                "Report",report," at",platform%info%lon,"E",platform%info%lat, &
 361                "N on processor", myproc,", position", platform%loc%x,platform%loc%y
 362          end if
 363       end if
 364
 365       read (analysis_date,'(i4,4(1x,i2))') &
 366                                     iyear, imonth, iday, ihour, imin
 367       call da_get_julian_time (iyear,imonth,iday,ihour,imin,analysis_time)
 368       tdiff = abs((obs_time - analysis_time)-0.02)
 369       dlat_earth = platform%info%lat
 370       dlon_earth = platform%info%lon
 371       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 372       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 373       dlat_earth = dlat_earth * deg2rad
 374       dlon_earth = dlon_earth * deg2rad
 375
 376       nlevels = platform%info%levels
 377       if (nlevels > max_ob_levels) then
 378          nlevels = max_ob_levels
 379          platform%info%levels = nlevels
 380       else if (nlevels < 1) then
 381          ! Not GPSPW and GPSZD data:
 382          if (fm /= 111 .and. fm /= 114) then
 383             cycle reports
 384          end if
 385       end if
 386
 387       levs = nlevels
 388
 389       if (fm /= 86 .and. fm /= 118) call da_obs_proc_station(platform,fm,uvq_direct)
 390
 391       ! Loop over duplicating obs for global
 392       ndup = 1
 393       if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
 394
 395       ! It is possible that logic for counting obs is incorrect for the
 396       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
 397       ! TBH:  20050913
 398
 399       if (.not.outside) then
 400          if (print_detail_obs .and. ndup > 1) then
 401             write(unit=stdout, fmt = fmt_info) &
 402                platform%info%platform,    &
 403                platform%info%date_char,   &
 404                platform%info%name,        &
 405                platform%info%levels,      &
 406                platform%info%lat,         &
 407                platform%info%lon,         &
 408                platform%info%elv,         &
 409                platform%info%id
 410
 411             write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
 412                ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
 413                platform%loc%i,  platform%loc%j,   &
 414                platform%loc%dx, platform%loc%dxm, &
 415               platform%loc%dy, platform%loc%dym
 416          end if
 417       end if
 418
 419       obs_index = fm_index(fm)
 420
 421       IF ( wind_sd ) THEN
 422          wind_sd_buoy   = .true.
 423          wind_sd_synop  = .true.
 424          wind_sd_ships  = .true.
 425          wind_sd_metar  = .true.
 426          wind_sd_sound  = .true.
 427          wind_sd_pilot  = .true.
 428          wind_sd_airep  = .true.
 429          wind_sd_qscat  = .true.
 430          wind_sd_tamdar = .true.
 431          wind_sd_geoamv = .true.
 432          wind_sd_mtgirs = .true.
 433        wind_sd_polaramv = .true.
 434        wind_sd_profiler = .true.
 435       END IF
 436
 437       dup_loop: do n = 1, ndup
 438          is_surface=.false.
 439          select case(fm)
 440
 441          case (12) ;
 442             if (.not. use_synopobs .or. ntotal(synop) == max_synop_input  ) cycle reports
 443             if (n==1) ntotal(synop) = ntotal(synop)+1
 444             if (outside) cycle reports
 445             if ( thin_conv_opt(synop) > no_thin ) then
 446                crit = tdiff
 447                call map2grids_conv(synop,dlat_earth,dlon_earth,crit,nlocal(synop),itx,1,itt,ilocal(synop),iuse)
 448                if ( .not. iuse ) cycle reports
 449             else
 450                nlocal(synop) = nlocal(synop) + 1
 451                ilocal(synop) = ilocal(synop) + 1
 452             end if
 453
 454             if (.not. wind_sd_synop) platform%each(1)%v%error = platform%each(1)%u%error
 455
 456             iv%synop(ilocal(synop))%h = platform%each(1)%height
 457             iv%synop(ilocal(synop))%u = platform%each(1)%u
 458             iv%synop(ilocal(synop))%v = platform%each(1)%v
 459             iv%synop(ilocal(synop))%t = platform%each(1)%t
 460             iv%synop(ilocal(synop))%q = platform%each(1)%q
 461             iv%synop(ilocal(synop))%p = platform%each(1)%p
 462
 463             if (wind_sd_synop .and. &
 464                 platform%each(1)%u%qc /= missing_data .and. platform%each(1)%v%qc /= missing_data ) &
 465                 call da_ffdduv_model (iv%synop(ilocal(synop))%u%inv, iv%synop(ilocal(synop))%v%inv, &
 466                                       platform%each(1)%u%inv, platform%each(1)%v%inv, convert_uv2fd)
 467
 468          case (13, 17) ;                  ! ships
 469             if (.not. use_shipsobs .or. ntotal(ships) == max_ships_input  ) cycle reports
 470             if (n==1) ntotal(ships) = ntotal(ships) + 1
 471             if (outside) cycle reports
 472              if ( thin_conv_opt(ships) > no_thin ) then
 473                crit = tdiff
 474                call map2grids_conv(ships,dlat_earth,dlon_earth,crit,nlocal(ships),itx,1,itt,ilocal(ships),iuse)
 475                 if ( .not. iuse ) cycle reports
 476              else
 477                 nlocal(ships) = nlocal(ships) + 1
 478                 ilocal(ships) = ilocal(ships) + 1
 479              end if
 480
 481             if (.not. wind_sd_ships) platform%each(1)%v%error = platform%each(1)%u%error
 482
 483             iv%ships(ilocal(ships))%h = platform%each(1)%height
 484             iv%ships(ilocal(ships))%u = platform%each(1)%u
 485             iv%ships(ilocal(ships))%v = platform%each(1)%v
 486             iv%ships(ilocal(ships))%t = platform%each(1)%t
 487             iv%ships(ilocal(ships))%p = platform%each(1)%p
 488             iv%ships(ilocal(ships))%q = platform%each(1)%q
 489
 490             if (wind_sd_ships .and. &
 491                 platform%each(1)%u%qc /= missing_data .and. platform%each(1)%v%qc /= missing_data ) &
 492                 call da_ffdduv_model (iv%ships(ilocal(ships))%u%inv, iv%ships(ilocal(ships))%v%inv, &
 493                                       platform%each(1)%u%inv, platform%each(1)%v%inv, convert_uv2fd)
 494
 495          case (15:16) ;
 496             if (.not. use_metarobs .or. ntotal(metar) == max_metar_input  ) cycle reports
 497             if (n==1) ntotal(metar) = ntotal(metar) + 1
 498             if (outside) cycle reports
 499              if ( thin_conv_opt(metar) > no_thin ) then
 500                crit = tdiff
 501                call map2grids_conv(metar,dlat_earth,dlon_earth,crit,nlocal(metar),itx,1,itt,ilocal(metar),iuse)
 502                 if ( .not. iuse ) cycle reports
 503              else
 504                 nlocal(metar) = nlocal(metar) + 1
 505                 ilocal(metar) = ilocal(metar) + 1
 506              end if
 507
 508             if (.not. wind_sd_metar) platform%each(1)%v%error = platform%each(1)%u%error
 509
 510             iv%metar(ilocal(metar))%h = platform%each(1)%height
 511             iv%metar(ilocal(metar))%u = platform%each(1)%u
 512             iv%metar(ilocal(metar))%v = platform%each(1)%v
 513             iv%metar(ilocal(metar))%t = platform%each(1)%t
 514             iv%metar(ilocal(metar))%p = platform%each(1)%p
 515             iv%metar(ilocal(metar))%q = platform%each(1)%q
 516
 517             if (wind_sd_metar .and. &
 518                 platform%each(1)%u%qc /= missing_data .and. platform%each(1)%v%qc /= missing_data ) &
 519                 call da_ffdduv_model (iv%metar(ilocal(metar))%u%inv, iv%metar(ilocal(metar))%v%inv, &
 520                                       platform%each(1)%u%inv, platform%each(1)%v%inv, convert_uv2fd)
 521
 522          case (32:34) ;
 523             if (.not. use_pilotobs .or. ntotal(pilot) == max_pilot_input  ) cycle reports
 524             if (n==1) ntotal(pilot) = ntotal(pilot) + 1
 525             if (outside) cycle reports
 526             if ( thin_conv_opt(pilot) > no_thin ) then
 527                crit = tdiff
 528                call map2grids_conv(pilot,dlat_earth,dlon_earth,crit,nlocal(pilot),itx,1,itt,ilocal(pilot),iuse)
 529                if ( .not. iuse ) cycle reports
 530             else
 531                nlocal(pilot) = nlocal(pilot) + 1
 532                ilocal(pilot) = ilocal(pilot) + 1
 533             end if
 534
 535             if (nlocal(pilot) == ilocal(pilot)) then
 536                allocate (iv%pilot(ilocal(pilot))%h(1:iv%info(pilot)%max_lev))
 537                allocate (iv%pilot(ilocal(pilot))%p(1:iv%info(pilot)%max_lev))
 538                allocate (iv%pilot(ilocal(pilot))%u(1:iv%info(pilot)%max_lev))
 539                allocate (iv%pilot(ilocal(pilot))%v(1:iv%info(pilot)%max_lev))
 540             end if
 541             do i = 1, nlevels
 542
 543                if (.not. wind_sd_pilot) platform%each(i)%v%error = platform%each(i)%u%error
 544
 545                iv%pilot(ilocal(pilot))%h(i) = platform%each(i)%height
 546                iv%pilot(ilocal(pilot))%p(i) = platform%each(i)%p%inv
 547                iv%pilot(ilocal(pilot))%u(i) = platform%each(i)%u
 548                iv%pilot(ilocal(pilot))%v(i) = platform%each(i)%v
 549
 550                if ( .not. wind_sd_pilot ) platform%each(i)%v%error = platform%each(i)%u%error
 551
 552                if (wind_sd_pilot .and. &
 553                    platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) then
 554                    call da_ffdduv_model (iv%pilot(ilocal(pilot))%u(i)%inv,iv%pilot(ilocal(pilot))%v(i)%inv, &
 555                                          platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
 556                end if
 557             end do
 558
 559          case (35:38) ;
 560             if (.not. use_soundobs .or. ntotal(sound) == max_sound_input  ) cycle reports
 561             if (n==1) ntotal(sound) = ntotal(sound) + 1
 562             if (n==1) ntotal(sonde_sfc) = ntotal(sonde_sfc) + 1
 563             if (outside) cycle reports
 564
 565 ! determine if level corresponds to surface
 566
 567             levs = nlevels
 568             surface_level = 0
 569             do i = 1, nlevels
 570                ! if (elevation and height are the same, it is surface.)
 571                if (platform%info%elv.ne.missing_r.and.platform%each(i)%height.ne.missing_r)then
 572                   if (abs(platform%info%elv - platform%each(i)%height) < 0.1) then
 573                      is_surface = .true.
 574                      surface_level = i
 575                      levs = nlevels - 1
 576                      exit
 577                   end if
 578                end if
 579             end do
 580
 581             if ( thin_conv_opt(sound) > no_thin .or. &
 582                  thin_conv_opt(sonde_sfc) > no_thin ) then
 583                crit = tdiff
 584                call map2grids_conv(sound,dlat_earth,dlon_earth,crit,nlocal(sound),itx,1,itt,ilocal(sound),iuse)
 585                crit = tdiff
 586                call map2grids_conv(sonde_sfc,dlat_earth,dlon_earth,crit,nlocal(sonde_sfc),itx,1,itt,ilocal(sonde_sfc),iuse)
 587                if ( .not. iuse ) cycle reports
 588             else
 589                nlocal(sound) = nlocal(sound) + 1
 590                ilocal(sound) = ilocal(sound) + 1
 591                nlocal(sonde_sfc) = nlocal(sonde_sfc) + 1
 592                ilocal(sonde_sfc) = ilocal(sonde_sfc) + 1
 593             end if
 594
 595             ! Search to see if we have surface obs.
 596
 597             surface_level = 0
 598
 599             do i = 1, nlevels
 600                ! if (elevation and height are the same, it is surface.)
 601                if (platform%info%elv.ne.missing_r.and.platform%each(i)%height.ne.missing_r)then
 602                   if (abs(platform%info%elv - platform%each(i)%height) < 0.1) then
 603                      surface_level = i
 604
 605                      if (.not. wind_sd_sound) platform%each(i)%v%error = platform%each(i)%u%error
 606
 607                      ! Save surface pressure.
 608                      iv%sonde_sfc(ilocal(sonde_sfc))%h = platform%each(i)%height
 609                      iv%sonde_sfc(ilocal(sonde_sfc))%u = platform%each(i)%u
 610                      iv%sonde_sfc(ilocal(sonde_sfc))%v = platform%each(i)%v
 611                      iv%sonde_sfc(ilocal(sonde_sfc))%t = platform%each(i)%t
 612                      iv%sonde_sfc(ilocal(sonde_sfc))%q = platform%each(i)%q
 613                      iv%sonde_sfc(ilocal(sonde_sfc))%p = platform%each(i)%p
 614
 615                      if (wind_sd_sound .and. &
 616                          platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
 617                          call da_ffdduv_model (iv%sonde_sfc(ilocal(sonde_sfc))%u%inv,iv%sonde_sfc(ilocal(sonde_sfc))%v%inv, &
 618                                                platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
 619                      exit
 620                   end if
 621                end if
 622             end do
 623
 624             ! processing the sonde_sfc data:
 625
 626             old_nlevels = nlevels
 627
 628             if (surface_level > 0) then
 629                 nlevels = nlevels - 1
 630             else
 631                iv%sonde_sfc(ilocal(sonde_sfc))%h = missing_r
 632                iv%sonde_sfc(ilocal(sonde_sfc))%u%inv   = missing_r
 633                iv%sonde_sfc(ilocal(sonde_sfc))%u%qc    = missing
 634                iv%sonde_sfc(ilocal(sonde_sfc))%u%error = abs(missing_r)
 635                iv%sonde_sfc(ilocal(sonde_sfc))%v = iv%sonde_sfc(ilocal(sonde_sfc))%u
 636                iv%sonde_sfc(ilocal(sonde_sfc))%t = iv%sonde_sfc(ilocal(sonde_sfc))%u
 637                iv%sonde_sfc(ilocal(sonde_sfc))%p = iv%sonde_sfc(ilocal(sonde_sfc))%u
 638                iv%sonde_sfc(ilocal(sonde_sfc))%q = iv%sonde_sfc(ilocal(sonde_sfc))%u
 639             end if
 640
 641             if (nlevels > 0) then
 642                if( nlocal(sound) == ilocal(sound)) then
 643                   allocate (iv%sound(ilocal(sound))%h (1:iv%info(sound)%max_lev))
 644                   allocate (iv%sound(ilocal(sound))%p (1:iv%info(sound)%max_lev))
 645                   allocate (iv%sound(ilocal(sound))%u (1:iv%info(sound)%max_lev))
 646                   allocate (iv%sound(ilocal(sound))%v (1:iv%info(sound)%max_lev))
 647                   allocate (iv%sound(ilocal(sound))%t (1:iv%info(sound)%max_lev))
 648                   allocate (iv%sound(ilocal(sound))%q (1:iv%info(sound)%max_lev))
 649                end if
 650
 651                j = 0
 652                do i = 1, old_nlevels
 653                   if (i == surface_level) cycle
 654
 655                   j=j+1
 656
 657                   if (.not. wind_sd_sound) platform%each(i)%v%error = platform%each(i)%u%error
 658
 659                   iv%sound(ilocal(sound))%h(j) = platform%each(i)%height
 660                   iv%sound(ilocal(sound))%p(j) = platform%each(i)%p%inv
 661                   iv%sound(ilocal(sound))%u(j) = platform%each(i)%u
 662                   iv%sound(ilocal(sound))%v(j) = platform%each(i)%v
 663                   iv%sound(ilocal(sound))%t(j) = platform%each(i)%t
 664                   iv%sound(ilocal(sound))%q(j) = platform%each(i)%q
 665
 666                   if (wind_sd_sound .and. &
 667                       platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
 668                       call da_ffdduv_model (iv%sound(ilocal(sound))%u(j)%inv,iv%sound(ilocal(sound))%v(j)%inv, &
 669                                             platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
 670                end do
 671
 672            end if
 673
 674
 675          case (101) ;
 676             if (.not. use_tamdarobs .or. ntotal(tamdar) == max_tamdar_input  ) cycle reports
 677
 678 ! determine if level corresponds to surface
 679             levs = nlevels
 680             surface_level = 0
 681             do i = 1, nlevels
 682                ! if (elevation and height are the same, it is surface.)
 683                if (platform%info%elv.ne.missing_r.and.platform%each(i)%height.ne.missing_r)then
 684                   if (abs(platform%info%elv - platform%each(i)%height) < 0.1) then
 685                      is_surface = .true.
 686                      surface_level = i
 687                      levs = nlevels - 1
 688                     exit
 689                   end if
 690                end if
 691             end do
 692
 693           if( is_surface) then
 694             if (n==1) ntotal(tamdar_sfc) = ntotal(tamdar_sfc) + 1
 695             if (outside) cycle reports
 696             if ( thin_conv_opt(tamdar_sfc) > no_thin ) then
 697                crit = tdiff
 698                call map2grids_conv(tamdar_sfc,dlat_earth,dlon_earth,crit,nlocal(tamdar_sfc),itx,1,itt,ilocal(tamdar_sfc),iuse)
 699                 if ( .not. iuse ) cycle reports
 700              else
 701                 nlocal(tamdar_sfc) = nlocal(tamdar_sfc) + 1
 702                 ilocal(tamdar_sfc) = ilocal(tamdar_sfc) + 1
 703              end if
 704
 705             if (.not. wind_sd_tamdar) platform%each(surface_level)%v%error = platform%each(surface_level)%u%error
 706
 707             iv%tamdar_sfc(ilocal(tamdar_sfc))%h = platform%each(surface_level)%height
 708             iv%tamdar_sfc(ilocal(tamdar_sfc))%u = platform%each(surface_level)%u
 709             iv%tamdar_sfc(ilocal(tamdar_sfc))%v = platform%each(surface_level)%v
 710             iv%tamdar_sfc(ilocal(tamdar_sfc))%t = platform%each(surface_level)%t
 711             iv%tamdar_sfc(ilocal(tamdar_sfc))%q = platform%each(surface_level)%q
 712             iv%tamdar_sfc(ilocal(tamdar_sfc))%p = platform%each(surface_level)%p
 713
 714             if (wind_sd_tamdar .and. &
 715                 platform%each(surface_level)%u%qc /= missing_data .and. platform%each(surface_level)%v%qc /= missing_data ) &
 716                 call da_ffdduv_model (iv%tamdar_sfc(ilocal(tamdar_sfc))%u%inv,iv%tamdar_sfc(ilocal(tamdar_sfc))%v%inv, &
 717                                       platform%each(surface_level)%u%inv, platform%each(surface_level)%v%inv, convert_uv2fd)
 718
 719            else ! not is_surface
 720
 721             if ( levs > 0 .and. n==1) ntotal(tamdar) = ntotal(tamdar) + 1
 722             if (outside) cycle reports
 723             if( thin_conv_opt(tamdar) <= no_thin ) then
 724                    nlocal(tamdar) = nlocal(tamdar) + 1
 725                    ilocal(tamdar) = ilocal(tamdar) + 1
 726                 if( nlocal(tamdar) == ilocal(tamdar)) then
 727                 allocate (iv%tamdar(ilocal(tamdar))%h (1:iv%info(tamdar)%max_lev))
 728                 allocate (iv%tamdar(ilocal(tamdar))%p (1:iv%info(tamdar)%max_lev))
 729                 allocate (iv%tamdar(ilocal(tamdar))%u (1:iv%info(tamdar)%max_lev))
 730                 allocate (iv%tamdar(ilocal(tamdar))%v (1:iv%info(tamdar)%max_lev))
 731                 allocate (iv%tamdar(ilocal(tamdar))%t (1:iv%info(tamdar)%max_lev))
 732                 allocate (iv%tamdar(ilocal(tamdar))%q (1:iv%info(tamdar)%max_lev))
 733              end if
 734
 735              do i = 1, nlevels
 736
 737                 if (.not. wind_sd_tamdar) platform%each(i)%v%error = platform%each(i)%u%error
 738
 739                 iv%tamdar(ilocal(tamdar))%h(i) = platform%each(i)%height
 740                 iv%tamdar(ilocal(tamdar))%p(i) = platform%each(i)%p%inv
 741                 iv%tamdar(ilocal(tamdar))%u(i) = platform%each(i)%u
 742                 iv%tamdar(ilocal(tamdar))%v(i) = platform%each(i)%v
 743                 iv%tamdar(ilocal(tamdar))%t(i) = platform%each(i)%t
 744                 iv%tamdar(ilocal(tamdar))%q(i) = platform%each(i)%q
 745
 746                 if (wind_sd_tamdar .and. &
 747                     platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
 748                     call da_ffdduv_model (iv%tamdar(ilocal(tamdar))%u(i)%inv,iv%tamdar(ilocal(tamdar))%v(i)%inv, &
 749                                           platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
 750              end do
 751
 752             else ! if thin_conv_opt > no_thin
 753               thin_3d=.true.
 754               i1   = platform%loc%i
 755               j1   = platform%loc%j
 756               dx   = platform%loc%dx
 757               dy   = platform%loc%dy
 758               dxm  = platform%loc%dxm
 759               dym  = platform%loc%dym
 760
 761                  do k=kms,kme
 762                   v_p(k) = dym*(dxm*grid%xb%p(i1,j1,k)+dx*grid%xb%p(i1+1,j1,k)) + &
 763                             dy*(dxm*grid%xb%p(i1,j1+1,k)+dx*grid%xb%p(i1+1,j1+1,k))
 764                  end do
 765                  do k=kms,kme
 766                   v_h(k) = dym*(dxm*grid%xb%h(i1,j1,k)+dx*grid%xb%h(i1+1,j1,k)) + &
 767                             dy*(dxm*grid%xb%h(i1,j1+1,k)+dx*grid%xb%h(i1+1,j1+1,k))
 768                  end do
 769               do k= 1, nlevels
 770                zk = missing_r
 771                if( platform%each(k)%p%qc  >= 0 ) then
 772                  call da_to_zk(platform%each(k)%p%inv, v_p, v_interp_p, zk)
 773                else if( platform%each(k)%height_qc  >= 0 ) then
 774                  call da_to_zk(platform%each(k)%height , v_h, v_interp_h, zk)
 775                else
 776                  write(unit=message(1),fmt='(A)')' For tamdar: neither height nor pressure qc is good'
 777                  call da_error(__FILE__,__LINE__,message(1:1))
 778                end if
 779                if ( zk == missing_r ) cycle
 780                crit = tdiff
 781                call map2grids_conv(tamdar,dlat_earth,dlon_earth,crit,nlocal(tamdar),itx,1,itt,ilocal(tamdar),iuse,zk,thin_3d)
 782                if ( .not. iuse ) cycle
 783                iv%info(tamdar)%levels(ilocal(tamdar))    = 1
 784                iv%info(tamdar)%name(ilocal(tamdar))      = platform%info%name
 785                iv%info(tamdar)%platform(ilocal(tamdar))  = platform%info%platform
 786                iv%info(tamdar)%id(ilocal(tamdar))        = platform%info%id
 787                iv%info(tamdar)%date_char(ilocal(tamdar)) = platform%info%date_char
 788                iv%info(tamdar)%lat(:,ilocal(tamdar))     = platform%info%lat
 789                iv%info(tamdar)%lon(:,ilocal(tamdar))     = platform%info%lon
 790                iv%info(tamdar)%elv(ilocal(tamdar))       = platform%info%elv
 791                iv%info(tamdar)%pstar(ilocal(tamdar))     = platform%info%pstar
 792
 793                iv%info(tamdar)%slp(ilocal(tamdar))           = platform%loc%slp
 794                iv%info(tamdar)%pw(ilocal(tamdar))            = platform%loc%pw
 795                iv%info(tamdar)%x(:,ilocal(tamdar))           = platform%loc%x
 796                iv%info(tamdar)%y(:,ilocal(tamdar))           = platform%loc%y
 797                iv%info(tamdar)%i(:,ilocal(tamdar))           = platform%loc%i
 798                iv%info(tamdar)%j(:,ilocal(tamdar))           = platform%loc%j
 799                iv%info(tamdar)%dx(:,ilocal(tamdar))          = platform%loc%dx
 800                iv%info(tamdar)%dxm(:,ilocal(tamdar))         = platform%loc%dxm
 801                iv%info(tamdar)%dy(:,ilocal(tamdar))          = platform%loc%dy
 802                iv%info(tamdar)%dym(:,ilocal(tamdar))         = platform%loc%dym
 803                iv%info(tamdar)%proc_domain(:,ilocal(tamdar)) = platform%loc%proc_domain
 804
 805                iv%info(tamdar)%obs_global_index(ilocal(tamdar)) = ntotal(tamdar)
 806                if( nlocal(tamdar) == ilocal(tamdar)) then
 807                 allocate (iv%tamdar(ilocal(tamdar))%h(1:iv%info(tamdar)%max_lev))
 808                 allocate (iv%tamdar(ilocal(tamdar))%p(1:iv%info(tamdar)%max_lev))
 809                 allocate (iv%tamdar(ilocal(tamdar))%u(1:iv%info(tamdar)%max_lev))
 810                 allocate (iv%tamdar(ilocal(tamdar))%v(1:iv%info(tamdar)%max_lev))
 811                 allocate (iv%tamdar(ilocal(tamdar))%t(1:iv%info(tamdar)%max_lev))
 812                 allocate (iv%tamdar(ilocal(tamdar))%q(1:iv%info(tamdar)%max_lev))
 813                end if
 814
 815                do i = 1, 1
 816
 817                   if (.not. wind_sd_tamdar) platform%each(k)%v%error = platform%each(k)%u%error
 818
 819                   iv%tamdar(ilocal(tamdar))%h(i) = platform%each(k)%height
 820                   iv%tamdar(ilocal(tamdar))%p(i) = platform%each(k)%p%inv
 821                   iv%tamdar(ilocal(tamdar))%u(i) = platform%each(k)%u
 822                   iv%tamdar(ilocal(tamdar))%v(i) = platform%each(k)%v
 823                   iv%tamdar(ilocal(tamdar))%t(i) = platform%each(k)%t
 824                   iv%tamdar(ilocal(tamdar))%q(i) = platform%each(k)%q
 825
 826                   if (wind_sd_tamdar .and. &
 827                       platform%each(k)%u%qc /= missing_data .and. platform%each(k)%v%qc /= missing_data ) &
 828                       call da_ffdduv_model (iv%tamdar(ilocal(tamdar))%u(i)%inv,iv%tamdar(ilocal(tamdar))%v(i)%inv, &
 829                                             platform%each(k)%u%inv, platform%each(k)%v%inv, convert_uv2fd)
 830                end do
 831               end do ! loop over k levels
 832
 833             end if ! if over thin_conv_ascii
 834
 835            end if ! if  is_surface
 836
 837          case (161) ;
 838             if (.not. use_mtgirsobs .or. ntotal(mtgirs) == max_mtgirs_input ) cycle reports
 839             if (n==1) ntotal(mtgirs) = ntotal(mtgirs) + 1
 840             if (outside) cycle reports
 841             if (outside) cycle reports
 842              if ( thin_conv_opt(mtgirs) > no_thin ) then
 843                crit = tdiff
 844                call map2grids_conv(mtgirs,dlat_earth,dlon_earth,crit,nlocal(mtgirs),itx,1,itt,ilocal(mtgirs),iuse)
 845                 if ( .not. iuse ) cycle reports
 846              else
 847                 nlocal(mtgirs) = nlocal(mtgirs) + 1
 848                 ilocal(mtgirs) = ilocal(mtgirs) + 1
 849              end if
 850
 851             if (nlevels > 0) then
 852
 853                if( nlocal(mtgirs) == ilocal(mtgirs)) then
 854                   allocate (iv%mtgirs(ilocal(mtgirs))%h (1:iv%info(mtgirs)%max_lev))
 855                   allocate (iv%mtgirs(ilocal(mtgirs))%p (1:iv%info(mtgirs)%max_lev))
 856                   allocate (iv%mtgirs(ilocal(mtgirs))%u (1:iv%info(mtgirs)%max_lev))
 857                   allocate (iv%mtgirs(ilocal(mtgirs))%v (1:iv%info(mtgirs)%max_lev))
 858                   allocate (iv%mtgirs(ilocal(mtgirs))%t (1:iv%info(mtgirs)%max_lev))
 859                   allocate (iv%mtgirs(ilocal(mtgirs))%q (1:iv%info(mtgirs)%max_lev))
 860                end if
 861
 862                j = 0
 863                do i = 1, nlevels
 864
 865                   j=j+1
 866
 867                   if (.not. wind_sd_mtgirs) platform%each(i)%v%error = platform%each(i)%u%error
 868
 869                   iv%mtgirs(ilocal(mtgirs))%h(j) = platform%each(i)%height
 870                   iv%mtgirs(ilocal(mtgirs))%p(j) = platform%each(i)%p%inv
 871                   iv%mtgirs(ilocal(mtgirs))%u(j) = platform%each(i)%u
 872                   iv%mtgirs(ilocal(mtgirs))%v(j) = platform%each(i)%v
 873                   iv%mtgirs(ilocal(mtgirs))%t(j) = platform%each(i)%t
 874                   iv%mtgirs(ilocal(mtgirs))%q(j) = platform%each(i)%q
 875
 876                   if(iv%mtgirs(ilocal(mtgirs))%q(j)%inv.ne.missing_r)then
 877                      iv%mtgirs(ilocal(mtgirs))%q(j)%inv = iv%mtgirs(ilocal(mtgirs))%q(j)%inv/1000.0
 878                   endif
 879                   if(iv%mtgirs(ilocal(mtgirs))%q(j)%error.ne.missing_r)then
 880                      iv%mtgirs(ilocal(mtgirs))%q(j)%error = iv%mtgirs(ilocal(mtgirs))%q(j)%error/1000.0/100.0
 881                   endif
 882
 883                   if (wind_sd_mtgirs .and. &
 884                       platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
 885                       call da_ffdduv_model (iv%mtgirs(ilocal(mtgirs))%u(j)%inv,iv%mtgirs(ilocal(mtgirs))%v(j)%inv, &
 886                                             platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
 887                end do
 888             end if
 889
 890          case (86)    ;
 891             ! Reject cloudy satem obs.
 892             if (.not.use_satemobs .or. ntotal(satem) == max_satem_input) cycle reports
 893             if (platform%loc%pw%inv > 10.0) cycle reports
 894             if (n==1) ntotal(satem) = ntotal(satem) + 1
 895             if (outside) cycle reports
 896              if ( thin_conv_opt(satem) > no_thin ) then
 897                crit = tdiff
 898                call map2grids_conv(satem,dlat_earth,dlon_earth,crit,nlocal(satem),itx,1,itt,ilocal(satem),iuse)
 899                 if ( .not. iuse ) cycle reports
 900              else
 901                 nlocal(satem) = nlocal(satem) + 1
 902                 ilocal(satem) = ilocal(satem) + 1
 903              end if
 904
 905
 906             ! The satem ref_p is put in the SLP position in OBSPROC data stream.
 907
 908             iv%satem(nlocal(satem))%ref_p= platform%loc%slp%inv
 909
 910             if( nlocal(satem) == ilocal(satem)) then
 911              allocate (iv%satem(ilocal(satem))%p        (1:iv%info(satem)%max_lev))
 912              allocate (iv%satem(ilocal(satem))%thickness(1:iv%info(satem)%max_lev))
 913              allocate (iv%satem(ilocal(satem))%org_thickness(1:iv%info(satem)%max_lev))
 914             end if
 915
 916             iv%satem(ilocal(satem))%p(1) = platform%each(1)%p%inv
 917             iv%satem(ilocal(satem))%thickness(1) = platform%each(1)%t
 918             !  write original thickness errors for filtered obs
 919             if (anal_type_qcobs) then
 920                do i = 1, nlevels
 921                   iv%satem(ilocal(satem))%org_thickness(i) = platform%each(i)%t
 922                end do
 923             end if
 924
 925             ! Splitting the reported satem data into smaller layers.
 926
 927             do i = 2, nlevels
 928                iv%satem(ilocal(satem))%p(i) = platform%each(i)%p%inv
 929                iv%satem(ilocal(satem))%thickness(i) = platform%each(i)%t
 930                if (platform%each(i)%t%qc /= missing_data   .and. &
 931                   platform%each(i-1)%t%qc /= missing_data) then
 932                   iv%satem(ilocal(satem))%thickness(i)%inv =            &
 933                   platform%each(i)%t%inv - platform%each(i-1)%t%inv
 934                else
 935                   iv%satem(ilocal(satem))%thickness(i)%inv = missing_r
 936                   iv%satem(ilocal(satem))%thickness(i)%qc  = missing_data
 937                end if
 938             end do
 939
 940             ! Thickness error (m):
 941
 942             do i = nlevels, 2, -1
 943                iv%satem(ilocal(satem))%thickness(i)%error = &
 944                sqrt(iv%satem(ilocal(satem))%thickness(i )%error ** 2 + &
 945                   iv%satem(ilocal(satem))%thickness(i-1)%error ** 2)
 946 !                  iv%satem(ilocal(satem))%thickness(i-1)%error ** 2) / gravity
 947             end do
 948
 949             iv%satem(ilocal(satem))%thickness(1)%error = &
 950                            sqrt(iv%satem(ilocal(satem))%thickness(1)%error ** 2 + &
 951                            platform%loc%pw%error ** 2)
 952 !                           platform%loc%pw%error ** 2) / gravity
 953
 954
 955             ! Geostationary ot Polar orbitting Satellite AMVs:
 956          case (88)    ;
 957             if (index(platform%info%name, 'MODIS') > 0 .or. &
 958                 index(platform%info%name, 'modis') > 0 .or. &
 959                 index(platform%info%id, 'AVHRR') > 0)  then
 960                if (.not.use_polaramvobs .or. ntotal(polaramv) == max_polaramv_input) cycle reports
 961                if (n==1) ntotal(polaramv) = ntotal(polaramv) + 1
 962                if (outside) cycle reports
 963                if (n==1) ntotal(Polaramv) = ntotal(polaramv) + 1
 964                 if ( thin_conv_opt(polaramv) > no_thin ) then
 965                   crit = tdiff
 966                   call map2grids_conv(polaramv,dlat_earth,dlon_earth,crit,nlocal(polaramv),itx,1,itt,ilocal(polaramv),iuse)
 967                    if ( .not. iuse ) cycle reports
 968                 else
 969                    nlocal(polaramv) = nlocal(polaramv) + 1
 970                    ilocal(polaramv) = ilocal(polaramv) + 1
 971                 end if
 972
 973                if (nlocal(polaramv) == ilocal(polaramv) ) then
 974                 allocate (iv%polaramv(ilocal(polaramv))%p (1:iv%info(polaramv)%max_lev))
 975                 allocate (iv%polaramv(ilocal(polaramv))%u (1:iv%info(polaramv)%max_lev))
 976                 allocate (iv%polaramv(ilocal(polaramv))%v (1:iv%info(polaramv)%max_lev))
 977                end if
 978
 979                do i = 1, nlevels
 980
 981                   if (.not. wind_sd_polaramv) platform%each(i)%v%error = platform%each(i)%u%error
 982
 983                   iv%polaramv(ilocal(polaramv))%p(i) = platform%each(i)%p%inv
 984                   iv%polaramv(ilocal(polaramv))%u(i) = platform%each(i)%u
 985                   iv%polaramv(ilocal(polaramv))%v(i) = platform%each(i)%v
 986
 987                   if (wind_sd_polaramv .and. &
 988                       platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
 989                       call da_ffdduv_model (iv%polaramv(ilocal(polaramv))%u(i)%inv,iv%polaramv(ilocal(polaramv))%v(i)%inv, &
 990                                             platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
 991                end do
 992                obs_index = polaramv ! geoamv is the fm_index value for 88
 993             else
 994                if (.not.use_geoamvobs .or. ntotal(geoamv) == max_geoamv_input) cycle reports
 995                if (n==1) ntotal(geoamv) = ntotal(geoamv) + 1
 996                if (outside) cycle reports
 997                 if ( thin_conv_opt(geoamv) > no_thin ) then
 998                   crit = tdiff
 999                   call map2grids_conv(geoamv,dlat_earth,dlon_earth,crit,nlocal(geoamv),itx,1,itt,ilocal(geoamv),iuse)
1000                    if ( .not. iuse ) cycle reports
1001                 else
1002                    nlocal(geoamv) = nlocal(geoamv) + 1
1003                    ilocal(geoamv) = ilocal(geoamv) + 1
1004                 end if
1005
1006                if( nlocal(geoamv) == ilocal(geoamv) )then
1007                 allocate (iv%geoamv(ilocal(geoamv))%p (1:iv%info(geoamv)%max_lev))
1008                 allocate (iv%geoamv(ilocal(geoamv))%u (1:iv%info(geoamv)%max_lev))
1009                 allocate (iv%geoamv(ilocal(geoamv))%v (1:iv%info(geoamv)%max_lev))
1010                end if
1011
1012                do i = 1, nlevels
1013
1014                   if (.not. wind_sd_geoamv) platform%each(i)%v%error = platform%each(i)%u%error
1015
1016                   iv%geoamv(ilocal(geoamv))%p(i) = platform%each(i)%p%inv
1017                   iv%geoamv(ilocal(geoamv))%u(i) = platform%each(i)%u
1018                   iv%geoamv(ilocal(geoamv))%v(i) = platform%each(i)%v
1019
1020                   if (wind_sd_geoamv .and. &
1021                       platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
1022                       call da_ffdduv_model (iv%geoamv(ilocal(geoamv))%u(i)%inv,iv%geoamv(ilocal(geoamv))%v(i)%inv, &
1023                                             platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
1024                end do
1025
1026             end if
1027
1028          case (42,96:97) ;
1029
1030             if (.not.use_airepobs .or. ntotal(airep) == max_airep_input) cycle reports
1031             if (n==1) ntotal(airep) = ntotal(airep) + 1
1032             if (outside) cycle reports
1033
1034             if( thin_conv_opt(airep) <= no_thin ) then
1035                    nlocal(airep) = nlocal(airep) + 1
1036                    ilocal(airep) = ilocal(airep) + 1
1037               if( nlocal(airep) == ilocal(airep)) then
1038               allocate (iv%airep(ilocal(airep))%h        (1:iv%info(airep)%max_lev))
1039               allocate (iv%airep(ilocal(airep))%p        (1:iv%info(airep)%max_lev))
1040               allocate (iv%airep(ilocal(airep))%u        (1:iv%info(airep)%max_lev))
1041               allocate (iv%airep(ilocal(airep))%v        (1:iv%info(airep)%max_lev))
1042               allocate (iv%airep(ilocal(airep))%t        (1:iv%info(airep)%max_lev))
1043               allocate (iv%airep(ilocal(airep))%q        (1:iv%info(airep)%max_lev))
1044              end if
1045
1046              do i = 1, nlevels
1047
1048                 if (.not. wind_sd_airep) platform%each(i)%v%error = platform%each(i)%u%error
1049
1050                 iv%airep(ilocal(airep))%h(i) = platform%each(i)%height
1051                 iv%airep(ilocal(airep))%p(i) = platform%each(i)%p%inv
1052                 iv%airep(ilocal(airep))%u(i) = platform%each(i)%u
1053                 iv%airep(ilocal(airep))%v(i) = platform%each(i)%v
1054                 iv%airep(ilocal(airep))%t(i) = platform%each(i)%t
1055                 iv%airep(ilocal(airep))%q(i) = platform%each(i)%q
1056
1057                 if (wind_sd_airep .and. &
1058                     platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
1059                     call da_ffdduv_model (iv%airep(ilocal(airep))%u(i)%inv,iv%airep(ilocal(airep))%v(i)%inv, &
1060                                           platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
1061              end do
1062
1063             else ! if thin_conv_opt > no_thin
1064               thin_3d=.true.
1065               i1   = platform%loc%i
1066               j1   = platform%loc%j
1067               dx   = platform%loc%dx
1068               dy   = platform%loc%dy
1069               dxm  = platform%loc%dxm
1070               dym  = platform%loc%dym
1071
1072                  do k=kms,kme
1073                   v_p(k) = dym*(dxm*grid%xb%p(i1,j1,k)+dx*grid%xb%p(i1+1,j1,k)) + &
1074                             dy*(dxm*grid%xb%p(i1,j1+1,k)+dx*grid%xb%p(i1+1,j1+1,k))
1075                  end do
1076                  do k=kms,kme
1077                   v_h(k) = dym*(dxm*grid%xb%h(i1,j1,k)+dx*grid%xb%h(i1+1,j1,k)) + &
1078                             dy*(dxm*grid%xb%h(i1,j1+1,k)+dx*grid%xb%h(i1+1,j1+1,k))
1079                  end do
1080               do k= 1, nlevels
1081                zk = missing_r
1082                if( platform%each(k)%p%qc  >= 0 ) then
1083                  call da_to_zk(platform%each(k)%p%inv, v_p, v_interp_p, zk)
1084                else if( platform%each(k)%height_qc  >= 0 ) then
1085                  call da_to_zk(platform%each(k)%height , v_h, v_interp_h, zk)
1086                else
1087                  write(unit=message(1),fmt='(A)')' For airep: neither height nor pressure qc is good'
1088                  call da_error(__FILE__,__LINE__,message(1:1))
1089                end if
1090                if ( zk == missing_r ) cycle
1091                crit = tdiff
1092                call map2grids_conv(airep,dlat_earth,dlon_earth,crit,nlocal(airep),itx,1,itt,ilocal(airep),iuse,zk,thin_3d)
1093                if ( .not. iuse ) cycle
1094                iv%info(airep)%levels(ilocal(airep))    = 1
1095                iv%info(airep)%name(ilocal(airep))      = platform%info%name
1096                iv%info(airep)%platform(ilocal(airep))  = platform%info%platform
1097                iv%info(airep)%id(ilocal(airep))        = platform%info%id
1098                iv%info(airep)%date_char(ilocal(airep)) = platform%info%date_char
1099                iv%info(airep)%lat(:,ilocal(airep))     = platform%info%lat
1100                iv%info(airep)%lon(:,ilocal(airep))     = platform%info%lon
1101                iv%info(airep)%elv(ilocal(airep))       = platform%info%elv
1102                iv%info(airep)%pstar(ilocal(airep))     = platform%info%pstar
1103
1104                iv%info(airep)%slp(ilocal(airep))           = platform%loc%slp
1105                iv%info(airep)%pw(ilocal(airep))            = platform%loc%pw
1106                iv%info(airep)%x(:,ilocal(airep))           = platform%loc%x
1107                iv%info(airep)%y(:,ilocal(airep))           = platform%loc%y
1108                iv%info(airep)%i(:,ilocal(airep))           = platform%loc%i
1109                iv%info(airep)%j(:,ilocal(airep))           = platform%loc%j
1110                iv%info(airep)%dx(:,ilocal(airep))          = platform%loc%dx
1111                iv%info(airep)%dxm(:,ilocal(airep))         = platform%loc%dxm
1112                iv%info(airep)%dy(:,ilocal(airep))          = platform%loc%dy
1113                iv%info(airep)%dym(:,ilocal(airep))         = platform%loc%dym
1114                iv%info(airep)%proc_domain(:,ilocal(airep)) = platform%loc%proc_domain
1115
1116                iv%info(airep)%obs_global_index(ilocal(airep)) = ntotal(airep)
1117                if( nlocal(airep) == ilocal(airep)) then
1118                 allocate (iv%airep(ilocal(airep))%h        (1:iv%info(airep)%max_lev))
1119                 allocate (iv%airep(ilocal(airep))%p        (1:iv%info(airep)%max_lev))
1120                 allocate (iv%airep(ilocal(airep))%u        (1:iv%info(airep)%max_lev))
1121                 allocate (iv%airep(ilocal(airep))%v        (1:iv%info(airep)%max_lev))
1122                 allocate (iv%airep(ilocal(airep))%t        (1:iv%info(airep)%max_lev))
1123                 allocate (iv%airep(ilocal(airep))%q        (1:iv%info(airep)%max_lev))
1124                end if
1125
1126                do i = 1, 1
1127
1128                   if (.not. wind_sd_airep) platform%each(k)%v%error = platform%each(k)%u%error
1129
1130                   iv%airep(ilocal(airep))%h(i) = platform%each(k)%height
1131                   iv%airep(ilocal(airep))%p(i) = platform%each(k)%p%inv
1132                   iv%airep(ilocal(airep))%u(i) = platform%each(k)%u
1133                   iv%airep(ilocal(airep))%v(i) = platform%each(k)%v
1134                   iv%airep(ilocal(airep))%t(i) = platform%each(k)%t
1135                   iv%airep(ilocal(airep))%q(i) = platform%each(k)%q
1136
1137                   if (wind_sd_airep .and. &
1138                       platform%each(k)%u%qc /= missing_data .and. platform%each(k)%v%qc /= missing_data ) &
1139                       call da_ffdduv_model (iv%airep(ilocal(airep))%u(i)%inv,iv%airep(ilocal(airep))%v(i)%inv, &
1140                                             platform%each(k)%u%inv, platform%each(k)%v%inv, convert_uv2fd)
1141
1142                end do
1143               end do ! loop over k levels
1144
1145             end if ! if over thin_conv_ascii
1146
1147          case (111, 114) ;
1148             if((.not.use_gpspwobs  .and. fm == 111) .or. ntotal(gpspw) == max_gpspw_input ) cycle reports
1149             if((.not.use_gpsztdobs .and. fm == 114) .or. ntotal(gpspw) == max_gpspw_input ) cycle reports
1150             if (n==1) ntotal(gpspw) = ntotal(gpspw) + 1
1151             if (outside) cycle reports
1152                 if ( thin_conv_opt(gpspw) > no_thin ) then
1153                   crit = tdiff
1154                   call map2grids_conv(gpspw,dlat_earth,dlon_earth,crit,nlocal(gpspw),itx,1,itt,ilocal(gpspw),iuse)
1155                    if ( .not. iuse ) cycle reports
1156                 else
1157                    nlocal(gpspw) = nlocal(gpspw) + 1
1158                    ilocal(gpspw) = ilocal(gpspw) + 1
1159                 end if
1160
1161             iv%gpspw(ilocal(gpspw))%tpw  = platform%loc%pw
1162
1163          case (116) ;
1164             if(.not.use_gpsrefobs  .or. ntotal(gpsref) == max_gpsref_input ) cycle reports
1165             if ( ob_format_gpsro /= ob_format_ascii ) cycle reports
1166             if (n==1) ntotal(gpsref) = ntotal(gpsref) + 1
1167             if (outside) cycle reports
1168                 if ( thin_conv_opt(gpsref) > no_thin ) then
1169                   crit = tdiff
1170                   call map2grids_conv(gpsref,dlat_earth,dlon_earth,crit,nlocal(gpsref),itx,1,itt,ilocal(gpsref),iuse)
1171                    if ( .not. iuse ) cycle reports
1172                 else
1173                    nlocal(gpsref) = nlocal(gpsref) + 1
1174                    ilocal(gpsref) = ilocal(gpsref) + 1
1175                 end if
1176             if( nlocal(gpsref) == ilocal(gpsref)) then
1177                allocate (iv%gpsref(ilocal(gpsref))%h (1:iv%info(gpsref)%max_lev))
1178                allocate (iv%gpsref(ilocal(gpsref))%ref(1:iv%info(gpsref)%max_lev))
1179                allocate (iv%gpsref(ilocal(gpsref))%  p(1:iv%info(gpsref)%max_lev))
1180                allocate (iv%gpsref(ilocal(gpsref))%  t(1:iv%info(gpsref)%max_lev))
1181                allocate (iv%gpsref(ilocal(gpsref))%  q(1:iv%info(gpsref)%max_lev))
1182             end if
1183
1184             do i = 1, nlevels
1185                iv%gpsref(ilocal(gpsref))%h(i)   = platform%each(i)%height
1186
1187                ! In OBSPROC, use "td" field to store "gpsref"
1188                iv%gpsref(ilocal(gpsref))%ref(i) = platform%each(i)%td
1189
1190                ! check height, only keep data below and above certain height
1191                if ( iv%gpsref(ilocal(gpsref))%h(i) > top_km_gpsro*1000.0 .or. &
1192                     iv%gpsref(ilocal(gpsref))%h(i) < bot_km_gpsro*1000.0 ) then
1193                   iv%gpsref(ilocal(gpsref))%ref(i)%qc = -77
1194                end if
1195
1196                ! Keep the retrieved p and t (and q) as "field_type":
1197                iv%gpsref(ilocal(gpsref))%p(i)   = platform%each(i)%p
1198                iv%gpsref(ilocal(gpsref))%t(i)   = platform%each(i)%t
1199                iv%gpsref(ilocal(gpsref))%q(i)   = platform%each(i)%q
1200             end do
1201
1202          CASE (118) ;
1203
1204             if ( .not.use_GpsephObs .or. ntotal(gpseph) == max_gpseph_input ) cycle reports
1205             if ( ob_format_gpsro /= ob_format_ascii ) cycle reports
1206             if (n==1) ntotal(gpseph) = ntotal(gpseph) + 1
1207             if ( gpseph_loadbalance ) then
1208                if ( myproc /= mod( (ntotal(gpseph)-1), num_procs ) ) cycle reports
1209             else
1210                if (outside) cycle reports
1211             endif
1212
1213             if (gpsro_drift ==  0) then
1214                ! lat stored in speed, lon stored in v
1215                ! replacing all levels with one lat/lon for gpsro_drift=0
1216                platform%each(1:nlevels)%speed%inv = platform%info%lat
1217                platform%each(1:nlevels)%v%inv = platform%info%lon
1218             end if
1219
1220             !create pseudo_ob on grid mean altitude for gpseph
1221             call da_gpseph_create_ob(              &
1222                nlevels,                            &
1223                platform%each(1:nlevels)%height,    &
1224                platform%each(1:nlevels)%td%inv,    & ! ref stored in td
1225                platform%each(1:nlevels)%speed%inv, & ! lat stored in speed
1226                platform%each(1:nlevels)%v%inv,     & ! lon stored in v
1227                platform%each(1:nlevels)%rh%inv,    & ! azim stored in rh
1228                pseudo_ob, lowest_level)
1229
1230             ! cycle when no valid levels
1231             if ( lowest_level < 1 ) cycle reports
1232
1233             nlocal(gpseph) = nlocal(gpseph) + 1
1234             ilocal(gpseph) = ilocal(gpseph) + 1
1235
1236             iv%gpseph(ilocal(gpseph))%level1 = lowest_level+1
1237             iv%gpseph(ilocal(gpseph))%level2 = kde-1
1238             iv%gpseph(ilocal(gpseph))%rfict  = platform%loc%pw%inv
1239
1240             levs = kde - kds + 1
1241             allocate (iv%gpseph(ilocal(gpseph))%h  (kds:kde))
1242             allocate (iv%gpseph(ilocal(gpseph))%ref(kds:kde))
1243             allocate (iv%gpseph(ilocal(gpseph))%eph(kds:kde))
1244             ! iv%gpseph includes additional azim/lat/lon info of each level
1245             allocate (iv%gpseph(ilocal(gpseph))%lat(kds:kde))
1246             allocate (iv%gpseph(ilocal(gpseph))%lon(kds:kde))
1247             allocate (iv%gpseph(ilocal(gpseph))%azim(kds:kde))
1248
1249             ! initialize eph
1250             ! eph will be calculated in da_get_innov_vector_gpseph
1251             iv%gpseph(ilocal(gpseph))%eph(:)%inv   = missing_r
1252             iv%gpseph(ilocal(gpseph))%eph(:)%qc    = missing_data
1253             iv%gpseph(ilocal(gpseph))%eph(:)%error = xmiss
1254
1255             do k = kts, kte
1256                iv%gpseph(ilocal(gpseph))%h(k)         = global_h_mean(k)*1000.0 !km to m
1257                iv%gpseph(ilocal(gpseph))%lat(k)       = pseudo_ob%lat(k)
1258                iv%gpseph(ilocal(gpseph))%lon(k)       = pseudo_ob%lon(k)
1259                iv%gpseph(ilocal(gpseph))%azim(k)      = pseudo_ob%azim(k)
1260                iv%gpseph(ilocal(gpseph))%ref(k)%inv   = pseudo_ob%ref(k)
1261                iv%gpseph(ilocal(gpseph))%ref(k)%qc    = missing_data
1262                iv%gpseph(ilocal(gpseph))%ref(k)%error = xmiss
1263             end do
1264
1265          case (121) ;
1266             if(.not.use_ssmt1obs  .or. ntotal(ssmt1) == max_ssmt1_input ) cycle reports
1267             if (n==1) ntotal(ssmt1) = ntotal(ssmt1) + 1
1268             if (outside) cycle reports
1269             if (n==1) ntotal(ssmt1) = ntotal(ssmt1) + 1
1270                 if ( thin_conv_opt(ssmt1) > no_thin ) then
1271                   crit = tdiff
1272                   call map2grids_conv(ssmt1,dlat_earth,dlon_earth,crit,nlocal(ssmt1),itx,1,itt,ilocal(ssmt1),iuse)
1273                    if ( .not. iuse ) cycle reports
1274                 else
1275                    nlocal(ssmt1) = nlocal(ssmt1) + 1
1276                    ilocal(ssmt1) = ilocal(ssmt1) + 1
1277                 end if
1278
1279             if ( nlocal(ssmt1) == ilocal(ssmt1)) then
1280              allocate (iv%ssmt1(ilocal(ssmt1))%h (1:iv%info(ssmt1)%max_lev))
1281              allocate (iv%ssmt1(ilocal(ssmt1))%p (1:iv%info(ssmt1)%max_lev))
1282              allocate (iv%ssmt1(ilocal(ssmt1))%t (1:iv%info(ssmt1)%max_lev))
1283             end if
1284
1285             do i = 1, nlevels
1286               iv%ssmt1(ilocal(ssmt1))%h(i) = platform%each(i)%height
1287               iv%ssmt1(ilocal(ssmt1))%p(i) = platform%each(i)%p%inv
1288               iv%ssmt1(ilocal(ssmt1))%t(i) = platform%each(i)%t
1289             end do
1290
1291          case (122) ;
1292             if(.not.use_ssmt2obs  .or. ntotal(ssmt2) == max_ssmt2_input ) cycle reports
1293             if (n==1) ntotal(ssmt2) = ntotal(ssmt2) + 1
1294             if (outside) cycle reports
1295                 if ( thin_conv_opt(ssmt2) > no_thin ) then
1296                   crit = tdiff
1297                   call map2grids_conv(ssmt2,dlat_earth,dlon_earth,crit,nlocal(ssmt2),itx,1,itt,ilocal(ssmt2),iuse)
1298                    if ( .not. iuse ) cycle reports
1299                 else
1300                    nlocal(ssmt2) = nlocal(ssmt2) + 1
1301                    ilocal(ssmt2) = ilocal(ssmt2) + 1
1302                 end if
1303             if ( nlocal(ssmt2) == ilocal(ssmt2)) then
1304              allocate (iv%ssmt2(ilocal(ssmt2))%h (1:iv%info(ssmt2)%max_lev))
1305              allocate (iv%ssmt2(ilocal(ssmt2))%p (1:iv%info(ssmt2)%max_lev))
1306              allocate (iv%ssmt2(ilocal(ssmt2))%rh(1:iv%info(ssmt2)%max_lev))
1307             end if
1308
1309             do i = 1, nlevels
1310                iv%ssmt2(ilocal(ssmt2))% h(i) = platform%each(i)%height
1311                iv%ssmt2(ilocal(ssmt2))% p(i) = platform%each(i)%p%inv
1312                iv%ssmt2(ilocal(ssmt2))%rh(i) = platform%each(i)%rh
1313             end do
1314
1315          case (281)    ;
1316             if(.not.use_qscatobs  .or. ntotal(qscat) == max_qscat_input ) cycle reports
1317             if (n==1) ntotal(qscat) = ntotal(qscat) + 1
1318             if (outside) cycle reports
1319                 if ( thin_conv_opt(qscat) > no_thin ) then
1320                   crit = tdiff
1321                   call map2grids_conv(qscat,dlat_earth,dlon_earth,crit,nlocal(qscat),itx,1,itt,ilocal(qscat),iuse)
1322                    if ( .not. iuse ) cycle reports
1323                 else
1324                    nlocal(qscat) = nlocal(qscat) + 1
1325                    ilocal(qscat) = ilocal(qscat) + 1
1326                 end if
1327
1328 !            if (nlocal(qscat) == ilocal(qscat)) then
1329
1330                 if (.not. wind_sd_qscat) platform%each(1)%v%error = platform%each(1)%u%error
1331
1332                 iv%qscat(ilocal(qscat))%h = platform%each(1)%height
1333                 iv%qscat(ilocal(qscat))%u = platform%each(1)%u
1334                 iv%qscat(ilocal(qscat))%v = platform%each(1)%v
1335
1336                 if (wind_sd_qscat .and. &
1337                     platform%each(1)%u%qc /= missing_data .and. platform%each(1)%v%qc /= missing_data ) &
1338                     call da_ffdduv_model (iv%qscat(ilocal(qscat))%u%inv,iv%qscat(ilocal(qscat))%v%inv, &
1339                                           platform%each(1)%u%inv, platform%each(1)%v%inv, convert_uv2fd)
1340
1341 !            end if
1342
1343             ! Impose minimum observation error = 1.0m/s for Quikscat data:
1344             iv%qscat(ilocal(qscat))%u%error = max(platform%each(1)%u%error,1.0)
1345             iv%qscat(ilocal(qscat))%v%error = max(platform%each(1)%v%error,1.0)
1346
1347          case (132) ; ! profiler
1348             if (.not. use_profilerobs .or. ntotal(profiler) == max_profiler_input ) cycle reports
1349             if (n==1) ntotal(profiler) = ntotal(profiler) + 1
1350             if (outside) cycle reports
1351                 if ( thin_conv_opt(profiler) > no_thin ) then
1352                   crit = tdiff
1353                   call map2grids_conv(profiler,dlat_earth,dlon_earth,crit,nlocal(profiler),itx,1,itt,ilocal(profiler),iuse)
1354                    if ( .not. iuse ) cycle reports
1355                 else
1356                    nlocal(profiler) = nlocal(profiler) + 1
1357                    ilocal(profiler) = ilocal(profiler) + 1
1358                 end if
1359
1360             if (nlocal(profiler) == ilocal(profiler) ) then
1361                allocate (iv%profiler(ilocal(profiler))%h (1:iv%info(profiler)%max_lev))
1362                allocate (iv%profiler(ilocal(profiler))%p (1:iv%info(profiler)%max_lev))
1363                allocate (iv%profiler(ilocal(profiler))%u (1:iv%info(profiler)%max_lev))
1364                allocate (iv%profiler(ilocal(profiler))%v (1:iv%info(profiler)%max_lev))
1365             end if
1366
1367             do i = 1, nlevels
1368
1369                if (.not. wind_sd_profiler) platform%each(i)%v%error = platform%each(i)%u%error
1370
1371                iv%profiler(ilocal(profiler))%h(i) = platform%each(i)%height
1372                iv%profiler(ilocal(profiler))%p(i) = platform%each(i)%p%inv
1373                iv%profiler(ilocal(profiler))%u(i) = platform%each(i)%u
1374                iv%profiler(ilocal(profiler))%v(i) = platform%each(i)%v
1375
1376                if (wind_sd_profiler .and. &
1377                    platform%each(i)%u%qc /= missing_data .and. platform%each(i)%v%qc /= missing_data ) &
1378                    call da_ffdduv_model (iv%profiler(ilocal(profiler))%u(i)%inv,iv%profiler(ilocal(profiler))%v(i)%inv, &
1379                                          platform%each(i)%u%inv, platform%each(i)%v%inv, convert_uv2fd)
1380             end do
1381
1382          case (135) ; ! Bogus
1383             if (.not. use_bogusobs .or. ntotal(bogus) == max_bogus_input ) cycle reports
1384             if (n==1) ntotal(bogus) = ntotal(bogus) + 1
1385             if (outside) cycle reports
1386             if (n==1) ntotal(bogus) = ntotal(bogus) + 1
1387                 if ( thin_conv_opt(bogus) > no_thin ) then
1388                   crit = tdiff
1389                   call map2grids_conv(bogus,dlat_earth,dlon_earth,crit,nlocal(bogus),itx,1,itt,ilocal(bogus),iuse)
1390                    if ( .not. iuse ) cycle reports
1391                 else
1392                    nlocal(bogus) = nlocal(bogus) + 1
1393                    ilocal(bogus) = ilocal(bogus) + 1
1394                 end if
1395
1396             if (ilocal(bogus) > max_bogus_input) then
1397                write(unit=message(1),fmt='(A,I6,A,I6)') &
1398                   'Bogus #=', ilocal(bogus), ' > max_bogus_input=', max_bogus_input
1399                call da_error(__FILE__,__LINE__,message(1:1))
1400             end if
1401
1402             if (nlocal(bogus) == ilocal(bogus)) then
1403                allocate (iv%bogus(ilocal(bogus))%h (1:iv%info(bogus)%max_lev))
1404                allocate (iv%bogus(ilocal(bogus))%p (1:iv%info(bogus)%max_lev))
1405                allocate (iv%bogus(ilocal(bogus))%u (1:iv%info(bogus)%max_lev))
1406                allocate (iv%bogus(ilocal(bogus))%v (1:iv%info(bogus)%max_lev))
1407                allocate (iv%bogus(ilocal(bogus))%t (1:iv%info(bogus)%max_lev))
1408                allocate (iv%bogus(ilocal(bogus))%q (1:iv%info(bogus)%max_lev))
1409             end if
1410
1411             do i = 1, nlevels
1412                if (.not. wind_sd) platform%each(i)%v%error = platform%each(i)%u%error
1413                iv%bogus(ilocal(bogus))%h(i) = platform%each(i)%height
1414                iv%bogus(ilocal(bogus))%p(i) = platform%each(i)%p%inv
1415                iv%bogus(ilocal(bogus))%u(i) = platform%each(i)%u
1416                iv%bogus(ilocal(bogus))%v(i) = platform%each(i)%v
1417                iv%bogus(ilocal(bogus))%t(i) = platform%each(i)%t
1418                iv%bogus(ilocal(bogus))%q(i) = platform%each(i)%q
1419             end do
1420
1421             iv%bogus(ilocal(bogus))%slp    = platform%loc%slp
1422
1423             if (print_detail_obs) then
1424                write(unit=stdout,fmt=*) 'nlevels=', nlevels
1425                write(unit=stdout,fmt=*) 'iv%info(bogus)%nlocal,slp', ilocal(bogus),  &
1426                   iv % bogus (ilocal(bogus)) % slp
1427                do i=1,nlevels
1428                   write(unit=stdout,fmt=*) 'nlocal(bogus), i ', nlocal(bogus),i
1429                   write(unit=stdout,fmt=*) 'iv%bogus(nlocal(bogus))%u,v=',  &
1430                      iv%bogus(ilocal(bogus))%u(i),iv%bogus(ilocal(bogus))%v(i)
1431                   write(unit=stdout,fmt=*) 'iv%bogus(nlocal(bogus))%t,q=',  &
1432                      iv%bogus(ilocal(bogus))%t(i),iv%bogus(ilocal(bogus))%q(i)
1433                end do
1434                write(unit=stdout,fmt='(2(a,i4))') 'nlocal(bogus)=',ilocal(bogus), &
1435                   'nlevels=',nlevels
1436             end if
1437
1438          case (18,19) ;             ! buoy
1439             if (.not. use_buoyobs .or. ntotal(buoy) == max_buoy_input ) cycle reports
1440             if (n==1) ntotal(buoy) = ntotal(buoy) + 1
1441             if (outside) cycle reports
1442                 if ( thin_conv_opt(buoy) > no_thin ) then
1443                   crit = tdiff
1444                   call map2grids_conv(buoy,dlat_earth,dlon_earth,crit,nlocal(buoy),itx,1,itt,ilocal(buoy),iuse)
1445                    if ( .not. iuse ) cycle reports
1446                 else
1447                    nlocal(buoy) = nlocal(buoy)  + 1
1448                    ilocal(buoy) = ilocal(buoy)  + 1
1449                 end if
1450
1451             if (.not. wind_sd_buoy) platform%each(1)%v%error = platform%each(1)%u%error
1452
1453             iv%buoy(ilocal(buoy))%h = platform%each(1)%height
1454             iv%buoy(ilocal(buoy))%u = platform%each(1)%u
1455             iv%buoy(ilocal(buoy))%v = platform%each(1)%v
1456             iv%buoy(ilocal(buoy))%t = platform%each(1)%t
1457             iv%buoy(ilocal(buoy))%p = platform%each(1)%p
1458             iv%buoy(ilocal(buoy))%q = platform%each(1)%q
1459
1460             if (wind_sd_buoy .and. &
1461                 platform%each(1)%u%qc /= missing_data .and. platform%each(1)%v%qc /= missing_data ) &
1462                 call da_ffdduv_model (iv%buoy(ilocal(buoy))%u%inv,iv%buoy(ilocal(buoy))%v%inv, &
1463                                       platform%each(1)%u%inv, platform%each(1)%v%inv, convert_uv2fd)
1464
1465
1466          case (133)    ;         ! AIRS retrievals
1467             if (.not. use_airsretobs .or. ntotal(airsr) == max_airsr_input ) cycle reports
1468             if (n==1) ntotal(airsr) = ntotal(airsr) + 1
1469             if (outside) cycle reports
1470                 if ( thin_conv_opt(airsr) > no_thin ) then
1471                   crit = tdiff
1472                   call map2grids_conv(airsr,dlat_earth,dlon_earth,crit,nlocal(airsr),itx,1,itt,ilocal(airsr),iuse)
1473                    if ( .not. iuse ) cycle reports
1474                 else
1475                    nlocal(airsr) = nlocal(airsr)  + 1
1476                    ilocal(airsr) = ilocal(airsr)  + 1
1477                 end if
1478
1479             if (nlocal(airsr) == ilocal(airsr)) then
1480                allocate (iv%airsr(ilocal(airsr))%h (1:iv%info(airsr)%max_lev))
1481                allocate (iv%airsr(ilocal(airsr))%p (1:iv%info(airsr)%max_lev))
1482                allocate (iv%airsr(ilocal(airsr))%t (1:iv%info(airsr)%max_lev))
1483                allocate (iv%airsr(ilocal(airsr))%q (1:iv%info(airsr)%max_lev))
1484             end if
1485             do i = 1, nlevels
1486                iv%airsr(ilocal(airsr))%h(i) = platform%each(i)%height
1487                iv%airsr(ilocal(airsr))%p(i) = platform%each(i)%p%inv
1488                iv%airsr(ilocal(airsr))%t(i) = platform%each(i)%t
1489                iv%airsr(ilocal(airsr))%q(i) = platform%each(i)%q
1490             end do
1491
1492          case default;
1493
1494             write(unit=message(1), fmt='(a)') 'unsaved obs found:'
1495             write(unit=message(2), fmt='(2a)') &
1496                'platform%info%platform=', platform%info%platform
1497             write(unit=message(3), fmt='(a, i3)') &
1498                'platform%info%levels=', platform%info%levels
1499             call da_warning(__FILE__,__LINE__,message(1:3))
1500             cycle
1501          end select
1502          if( is_surface .or. (obs_index == gpspw) .or. (levs > 0 .and. .not. thin_conv_ascii) .or. &
1503             (levs > 0 .and. (thin_conv_ascii .and. (obs_index /=  airep .and. obs_index /= tamdar))) ) then
1504             iv%info(obs_index)%name(ilocal(obs_index))      = platform%info%name
1505             iv%info(obs_index)%platform(ilocal(obs_index))  = platform%info%platform
1506             iv%info(obs_index)%id(ilocal(obs_index))        = platform%info%id
1507             iv%info(obs_index)%date_char(ilocal(obs_index)) = platform%info%date_char
1508             ! nlevels adjusted for some obs types so use that
1509             iv%info(obs_index)%levels(ilocal(obs_index))    = min(iv%info(obs_index)%max_lev, levs)
1510             iv%info(obs_index)%lat(:,ilocal(obs_index))     = platform%info%lat
1511             iv%info(obs_index)%lon(:,ilocal(obs_index))     = platform%info%lon
1512             iv%info(obs_index)%elv(ilocal(obs_index))       = platform%info%elv
1513             iv%info(obs_index)%pstar(ilocal(obs_index))     = platform%info%pstar
1514
1515             iv%info(obs_index)%slp(ilocal(obs_index))           = platform%loc%slp
1516             iv%info(obs_index)%pw(ilocal(obs_index))            = platform%loc%pw
1517             iv%info(obs_index)%x(:,ilocal(obs_index))           = platform%loc%x
1518             iv%info(obs_index)%y(:,ilocal(obs_index))           = platform%loc%y
1519             iv%info(obs_index)%i(:,ilocal(obs_index))           = platform%loc%i
1520             iv%info(obs_index)%j(:,ilocal(obs_index))           = platform%loc%j
1521             iv%info(obs_index)%dx(:,ilocal(obs_index))          = platform%loc%dx
1522             iv%info(obs_index)%dxm(:,ilocal(obs_index))         = platform%loc%dxm
1523             iv%info(obs_index)%dy(:,ilocal(obs_index))          = platform%loc%dy
1524             iv%info(obs_index)%dym(:,ilocal(obs_index))         = platform%loc%dym
1525             iv%info(obs_index)%proc_domain(:,ilocal(obs_index)) = platform%loc%proc_domain
1526
1527             iv%info(obs_index)%obs_global_index(nlocal(obs_index)) = ntotal(obs_index)
1528
1529             ! special case for sonde_sfc, duplicate sound info
1530             if (obs_index == sound) then
1531                iv%info(sonde_sfc)%name(ilocal(sonde_sfc))      = platform%info%name
1532                iv%info(sonde_sfc)%platform(ilocal(sonde_sfc))  = platform%info%platform
1533                iv%info(sonde_sfc)%id(ilocal(sonde_sfc))        = platform%info%id
1534                iv%info(sonde_sfc)%date_char(ilocal(sonde_sfc)) = platform%info%date_char
1535                iv%info(sonde_sfc)%levels(ilocal(sonde_sfc))    = 1
1536                iv%info(sonde_sfc)%lat(:,ilocal(sonde_sfc))     = platform%info%lat
1537                iv%info(sonde_sfc)%lon(:,ilocal(sonde_sfc))     = platform%info%lon
1538                iv%info(sonde_sfc)%elv(ilocal(sonde_sfc))       = platform%info%elv
1539                iv%info(sonde_sfc)%pstar(ilocal(sonde_sfc))     = platform%info%pstar
1540
1541                iv%info(sonde_sfc)%slp(ilocal(sonde_sfc))           = platform%loc%slp
1542                iv%info(sonde_sfc)%pw(ilocal(sonde_sfc))            = platform%loc%pw
1543                iv%info(sonde_sfc)%x(:,ilocal(sonde_sfc))           = platform%loc%x
1544                iv%info(sonde_sfc)%y(:,ilocal(sonde_sfc))           = platform%loc%y
1545                iv%info(sonde_sfc)%i(:,ilocal(sonde_sfc))           = platform%loc%i
1546                iv%info(sonde_sfc)%j(:,ilocal(sonde_sfc))           = platform%loc%j
1547                iv%info(sonde_sfc)%dx(:,ilocal(sonde_sfc))          = platform%loc%dx
1548                iv%info(sonde_sfc)%dxm(:,ilocal(sonde_sfc))         = platform%loc%dxm
1549                iv%info(sonde_sfc)%dy(:,ilocal(sonde_sfc))          = platform%loc%dy
1550                iv%info(sonde_sfc)%dym(:,ilocal(sonde_sfc))         = platform%loc%dym
1551                iv%info(sonde_sfc)%proc_domain(:,ilocal(sonde_sfc)) = platform%loc%proc_domain
1552                iv%info(sonde_sfc)%obs_global_index(ilocal(sonde_sfc)) = ntotal(obs_index)
1553             end if
1554
1555             if (is_surface .and. obs_index == tamdar) then
1556                iv%info(tamdar_sfc)%name(ilocal(tamdar_sfc))      = platform%info%name
1557                iv%info(tamdar_sfc)%platform(ilocal(tamdar_sfc))  = platform%info%platform
1558                iv%info(tamdar_sfc)%id(ilocal(tamdar_sfc))        = platform%info%id
1559                iv%info(tamdar_sfc)%date_char(ilocal(tamdar_sfc)) = platform%info%date_char
1560                iv%info(tamdar_sfc)%levels(ilocal(tamdar_sfc))    = 1
1561                iv%info(tamdar_sfc)%lat(:,ilocal(tamdar_sfc))     = platform%info%lat
1562                iv%info(tamdar_sfc)%lon(:,ilocal(tamdar_sfc))     = platform%info%lon
1563                iv%info(tamdar_sfc)%elv(ilocal(tamdar_sfc))       = platform%info%elv
1564                iv%info(tamdar_sfc)%pstar(ilocal(tamdar_sfc))     = platform%info%pstar
1565
1566                iv%info(tamdar_sfc)%slp(ilocal(tamdar_sfc))           = platform%loc%slp
1567                iv%info(tamdar_sfc)%pw(ilocal(tamdar_sfc))            = platform%loc%pw
1568                iv%info(tamdar_sfc)%x(:,ilocal(tamdar_sfc))           = platform%loc%x
1569                iv%info(tamdar_sfc)%y(:,ilocal(tamdar_sfc))           = platform%loc%y
1570                iv%info(tamdar_sfc)%i(:,ilocal(tamdar_sfc))           = platform%loc%i
1571                iv%info(tamdar_sfc)%j(:,ilocal(tamdar_sfc))           = platform%loc%j
1572                iv%info(tamdar_sfc)%dx(:,ilocal(tamdar_sfc))          = platform%loc%dx
1573                iv%info(tamdar_sfc)%dxm(:,ilocal(tamdar_sfc))         = platform%loc%dxm
1574                iv%info(tamdar_sfc)%dy(:,ilocal(tamdar_sfc))          = platform%loc%dy
1575                iv%info(tamdar_sfc)%dym(:,ilocal(tamdar_sfc))         = platform%loc%dym
1576                iv%info(tamdar_sfc)%proc_domain(:,ilocal(tamdar_sfc)) = platform%loc%proc_domain
1577
1578                iv%info(tamdar_sfc)%obs_global_index(ilocal(tamdar_sfc)) = ntotal(tamdar_sfc)
1579             end if
1580          end if  ! for thin_conv_ascii skipping obs_index for which thin_3d is true like airep and tamdir
1581
1582          if (obs_index == gpseph) then
1583             iv%info(gpseph)%levels(ilocal(gpseph))   = kde - kds + 1
1584             iv%info(gpseph)%lat(:,ilocal(gpseph))    = iv%gpseph(ilocal(gpseph))%lat(:)
1585             iv%info(gpseph)%lon(:,ilocal(gpseph))    = iv%gpseph(ilocal(gpseph))%lon(:)
1586             if (gpsro_drift ==  0) then
1587                iv%info(gpseph)%lat(:,ilocal(gpseph)) = platform%info%lat
1588                iv%info(gpseph)%lon(:,ilocal(gpseph)) = platform%info%lon
1589             end if
1590          end if
1591
1592          if (global .and. n < 2) then
1593             if (test_transforms) exit dup_loop
1594             if (platform%loc % i >= ide) then
1595                platform%loc%i = platform%loc % i - ide
1596             else if (platform%loc % i < ids) then
1597                platform%loc%i = platform%loc % i + ide
1598             end if
1599
1600             platform%loc%proc_domain = .not. platform%loc%proc_domain
1601          end if
1602       end do dup_loop
1603    end do reports
1604
1605    if ( use_gpsephobs ) then
1606       deallocate (pseudo_ob%ref)
1607       deallocate (pseudo_ob%lat)
1608       deallocate (pseudo_ob%lon)
1609       deallocate (pseudo_ob%azim)
1610    end if
1611
1612    close(iunit)
1613
1614    call da_free_unit(iunit)
1615
1616    ! thinning check
1617    if ( thin_conv_ascii ) then
1618       do n = 1, num_ob_indexes
1619           if ( .not. allocated(thinning_grid_conv(n)%icount) ) cycle
1620           if (n==radar) cycle
1621           allocate ( in(thinning_grid_conv(n)%itxmax) )
1622           allocate (out(thinning_grid_conv(n)%itxmax) )
1623             do i = 1, thinning_grid_conv(n)%itxmax
1624                in(i) = thinning_grid_conv(n)%score_crit(i)
1625             end do
1626 #ifdef DM_PARALLEL
1627             ! Get minimum crit and associated processor index.
1628             call mpi_reduce(in, out, thinning_grid_conv(n)%itxmax, true_mpi_real, mpi_min, root, comm, ierr)
1629             call wrf_dm_bcast_real (out, thinning_grid_conv(n)%itxmax)
1630 #else
1631             out = in
1632 #endif
1633             do i = 1, thinning_grid_conv(n)%itxmax
1634               if( out(i) < 9.99e6) iv%info(n)%thin_ntotal=  iv%info(n)%thin_ntotal + 1
1635                if ( abs(out(i)-thinning_grid_conv(n)%score_crit(i)) > 1.0E-10 ) then
1636                   thinning_grid_conv(n)%ibest_obs(i) = 0
1637                end if
1638             end do
1639 !            do j = iv%info(n)%plocal(iv%time -1)+1 , iv%info(n)%plocal(iv%time -1)+nlocal(n)
1640             do j = iv%info(n)%plocal(iv%time -1)+1 , nlocal(n)
1641                found = .false.
1642                do i = 1, thinning_grid_conv(n)%itxmax
1643                   if ( thinning_grid_conv(n)%ibest_obs(i) == j .and.         &
1644                        thinning_grid_conv(n)%score_crit(i) < 9.99e6 ) then
1645                    iv%info(n)%thin_nlocal =  iv%info(n)%thin_nlocal + 1
1646                      found = .true.
1647
1648                      exit
1649                   end if
1650                end do
1651                if ( .not. found ) then
1652                   iv%info(n)%thinned(:,j) = .true.
1653                end if
1654             end do
1655          deallocate( in  )
1656          deallocate( out )
1657       end do ! loop over num_ob_indexes
1658    end if  ! thin_conv_ascii
1659    if (trace_use) call da_trace_exit("da_read_obs_ascii")
1660
1661 end subroutine da_read_obs_ascii
1662
1663