var/da/da_obs_io/da_scan_obs_chem_sfc.inc

   1 subroutine da_scan_obs_chem_sfc (iv, filename, grid)
   2
   3    !-------------------------------------------------------------------------
   4    ! Purpose:        Read surface chem IC observation files
   5    !
   6    ! History: 03/2019  Creation (Wei Sun)
   7    !
   8    !          06/2021  Add maximum thresholds for each species (Soyoung Ha)
   9    ! Caution: chem_cv_options == 10 or 20 assumes all six species in [ug/m3]
  10    !          chem_cv_options == 108 assumes pm25 and pm10 in [ug/m3], but
  11    !                                         so2, no2, o3, co in [ppmv]
  12    !-------------------------------------------------------------------------
  13
  14    implicit none
  15
  16    type (iv_type),    intent(inout) :: iv
  17    character(len=*),  intent(in)    :: filename
  18    type(domain),     intent(in)     :: grid     ! first guess state.
  19
  20    ! chem_cv_options == 108: Maximum thresholds for observation values - treat as missing above these values
  21    ! FIXME: We may want to put these maximum thresholds in namelist later.
  22    real, dimension(6)  :: max_pm = (/500.,999.,2.0, 2.0, 2.0, 50./)  ! [ug/m3, ug/m3, ppmv, ppmv, ppmv, ppmv]
  23    real                         :: pm25, pm10, so2, no2, o3, co
  24
  25    character (len =   6)        :: SiteID
  26    character (len =  10)        :: fmt_name
  27
  28    character (len = 160)        :: info_string
  29
  30    integer                      :: i, j, iost, nlevels, old_nlevels, fm,iunit
  31    integer                      :: ichem, ipm
  32    real                         :: dattmp,lattmp,lontmp
  33
  34    type (singl_level_type)      :: platform
  35    logical                      :: outside
  36    logical                      :: outside_all
  37    integer                      :: surface_level
  38    real                         :: height_error, u_comp, v_comp
  39    integer                      :: nlocal(num_ob_indexes)
  40    integer                      :: ilocal(num_ob_indexes)
  41    integer                      :: ntotal(num_ob_indexes)
  42
  43    integer                      :: ndup, n, report, obs_index
  44
  45    real*8                       :: obs_time, analysis_time
  46    integer                      :: iyear, imonth, iday, ihour, imin, isite
  47    real                         :: tdiff, dlat_earth,dlon_earth,crit
  48    integer                      :: itt,itx,iout
  49    real, allocatable            :: in(:), out(:)
  50    logical                      :: found, iuse, thin_3d, is_surface
  51    integer                      :: i1,j1,k, levs
  52    real                         :: dx,dy,dxm,dym,zk
  53    real                         :: v_p(kms:kme),v_h(kms:kme)
  54
  55    if (trace_use) call da_trace_entry("da_scan_obs_chem_sfc")
  56
  57  ! Initialize counts
  58
  59    if ( thin_conv_ascii ) then
  60        do n = 1, num_ob_indexes
  61           if ( n == radar ) cycle
  62           call cleangrids_conv(n)
  63        end do
  64    end if
  65
  66    ! open file
  67    ! =========
  68
  69    call da_get_unit(iunit)
  70    open(unit   = iunit,     &
  71       FILE   = trim(filename), &
  72       FORM   = 'FORMATTED',  &
  73       ACCESS = 'SEQUENTIAL', &
  74       iostat =  iost,     &
  75       STATUS = 'OLD')
  76
  77    if (iost /= 0) then
  78       write(unit=message(1),fmt='(A,I5,A)') &
  79          "Error",iost," opening chemic surface obs file "//trim(filename)
  80       call da_warning(__FILE__,__LINE__,message(1:1))
  81       call da_free_unit(iunit)
  82       if (trace_use) call da_trace_exit("da_scan_obs_chem_sfc")
  83       return
  84    end if
  85
  86    allocate ( platform%chem (num_chemic_surf) )
  87    obs_index = chemic_surf
  88
  89    imin=0     ! minutes
  90
  91    ! loop over records
  92    ! -----------------
  93
  94    report = 0 ! report number in file
  95
  96    reports: &
  97    do
  98
  99       report = report+1
 100
 101       ! Initialization
 102       ! =============================
 103       platform%chem(:)%inv = missing_r
 104       platform%chem(:)%qc  = missing
 105       platform%chem(:)%error = missing_r
 106
 107       ! read station general info
 108       ! =============================
 109
 110       if (chem_cv_options == 10) then !!!! ium data !!!!
 111
 112         write(*,*) 'read ium data'
 113             read (iunit, *, iostat =  iost)  &
 114             iyear, imonth, iday, ihour, imin, &
 115             platform%info%lat,                &
 116             platform%info%lon,                &
 117             platform%chem(PARAM_FIRST_SCALAR)%inv        !!! pm25 only !!!
 118             write (platform%info%date_char,'(i4,3(i2))')  iyear, imonth, iday, ihour
 119             platform%info%id = "333"
 120
 121       else if (chem_cv_options == 20) then
 122
 123          if (chemicda_opt == 1) then
 124             read (iunit, *, iostat =  iost)  &
 125             iyear, imonth, iday, ihour, imin, &
 126             platform%info%lat,                &
 127             platform%info%lon,                &
 128             platform%chem(PARAM_FIRST_SCALAR)%inv        !!! pm25 only !!!
 129          else if (chemicda_opt == 2) then
 130             read (iunit, *, iostat =  iost)  &
 131             iyear, imonth, iday, ihour, imin, &
 132             platform%info%lat,                &
 133             platform%info%lon,                &
 134             platform%chem(PARAM_FIRST_SCALAR)%inv,  &    !!! pm25 read !!!
 135             platform%chem(PARAM_FIRST_SCALAR)%inv        !!! pm10 only !!!
 136          else if (chemicda_opt == 3) then
 137             read (iunit, *, iostat =  iost)  &
 138             iyear, imonth, iday, ihour, imin, &
 139             platform%info%lat,                &
 140             platform%info%lon,                &
 141             platform%chem(PARAM_FIRST_SCALAR)%inv,  &
 142             platform%chem(PARAM_FIRST_SCALAR+1)%inv        !!! pm10 & pm25 !!!
 143          else if (chemicda_opt == 4) then
 144             read (iunit, *, iostat =  iost)  &
 145             iyear, imonth, iday, ihour, imin, &
 146             platform%info%lat,                &
 147             platform%info%lon,                &
 148             platform%chem(PARAM_FIRST_SCALAR)%inv,  &
 149             platform%chem(PARAM_FIRST_SCALAR)%inv,  &
 150             platform%chem(PARAM_FIRST_SCALAR)%inv,  &      !!! so2 !!!
 151             platform%chem(PARAM_FIRST_SCALAR+1)%inv,  &    !!! no2 !!!
 152             platform%chem(PARAM_FIRST_SCALAR+2)%inv        !!! o3 !!!
 153             write (platform%info%date_char,'(i4,3(i2))')  iyear, imonth, iday, ihour
 154          else if (chemicda_opt == 5) then
 155             read (iunit, *, iostat =  iost)  &
 156             iyear, imonth, iday, ihour, imin, &
 157             platform%info%lat,                &
 158             platform%info%lon,                &
 159             platform%chem(PARAM_FIRST_SCALAR)%inv,  &      !!! pm25 !!!
 160             platform%chem(PARAM_FIRST_SCALAR+1)%inv,  &    !!! pm10 !!!
 161             platform%chem(PARAM_FIRST_SCALAR+2)%inv,  &    !!! so2 !!!
 162             platform%chem(PARAM_FIRST_SCALAR+3)%inv,  &    !!! no2 !!!
 163             platform%chem(PARAM_FIRST_SCALAR+4)%inv,  &    !!! o3 !!!
 164             platform%chem(PARAM_FIRST_SCALAR+5)%inv        !!! co !!!
 165             write (platform%info%date_char,'(i4,3(i2))')  iyear, imonth, iday, ihour
 166             platform%info%id = "333"
 167          end if
 168
 169       else if (chem_cv_options == 108) then
 170
 171       read (iunit, fmt='(A6,F10.6,1X,F10.6,1X,I4,3I2,2F9.1,4F9.3)', iostat=iost)  &
 172             SiteID,                                   &
 173             platform%info%lat,                        &
 174             platform%info%lon,                        &
 175             iyear, imonth, iday, ihour,               &
 176             pm25,  pm10,   so2,  no2,   o3, co
 177
 178       platform%info%id = trim(SiteID)
 179       read(SiteID,'(I6)') isite
 180       write(platform%info%date_char,'(i4,3(i2.2))')  iyear, imonth, iday, ihour
 181
 182       select case ( chemicda_opt )
 183       case ( 1 )
 184             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 185             platform%info%name                      = "pm25"
 186       case ( 2 )
 187             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm10
 188             platform%info%name                      = "pm10"
 189       case ( 3 )
 190             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 191             platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10
 192             platform%info%name                      = "all_pm"
 193       case ( 4 )
 194             platform%chem(PARAM_FIRST_SCALAR  )%inv = so2
 195             platform%chem(PARAM_FIRST_SCALAR+1)%inv = no2
 196             platform%chem(PARAM_FIRST_SCALAR+2)%inv = o3
 197             platform%chem(PARAM_FIRST_SCALAR+3)%inv = co
 198             platform%info%name                      = "gas"
 199       case ( 5 )
 200             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 201             platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10
 202             platform%chem(PARAM_FIRST_SCALAR+2)%inv = so2
 203             platform%chem(PARAM_FIRST_SCALAR+3)%inv = no2
 204             platform%chem(PARAM_FIRST_SCALAR+4)%inv = o3
 205             platform%chem(PARAM_FIRST_SCALAR+5)%inv = co
 206             platform%info%name                      = "all"
 207       case default
 208             platform%chem(PARAM_FIRST_SCALAR  )%inv = pm25
 209             platform%info%name                      = ""
 210       end select
 211
 212       ! Sanity check
 213       !!!!!!!!!!!!!!!!!!!
 214       do ichem = PARAM_FIRST_SCALAR, num_chemic_surf
 215          ipm = ichem - PARAM_FIRST_SCALAR + 1
 216          if (platform%chem(ichem)%inv<=0.or.abs(platform%chem(ichem)%inv)>=max_pm(ipm)) then
 217              platform%chem(ichem)%inv=missing_r
 218          else
 219              platform%chem(ichem)%qc = 0
 220          endif
 221       end do
 222
 223       end if !(chem_cv_options == 108)
 224
 225       if (iost /= 0) then
 226          ! FIX? This is expected, but its unclear how we handle failure
 227          ! here without assuming the fortran2003 convention on
 228          ! error statuses
 229          exit reports
 230       end if
 231
 232       if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
 233       ! Fix funny wind direction at Poles
 234       if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 235          platform%info%lon = 0.0
 236       end if
 237
 238       ! read model location
 239       ! =========================
 240
 241       ! Check if outside of the time range:
 242
 243       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 244       if ( obs_time < time_slots(0) .or. &
 245            obs_time >= time_slots(num_fgat_time) ) then
 246          if (print_detail_obs) then
 247            write(unit=stdout, fmt='(a)') '*** Outside of the time range:'
 248            write(unit=stdout, fmt=fmt_info) &
 249             platform%info%platform,    &
 250             platform%info%date_char,   &
 251             platform%info%name,        &
 252             platform%info%levels,      &
 253             platform%info%lat,         &
 254             platform%info%lon,         &
 255             platform%info%elv,         &
 256             platform%info%id
 257          end if
 258          cycle
 259       endif
 260
 261       ! Restrict to a range of reports, useful for debugging
 262
 263       if (report < report_start) then
 264          cycle
 265       end if
 266
 267       if (report > report_end) then
 268          exit
 269       end if
 270
 271       call da_llxy (platform%info, platform%loc, outside, outside_all)
 272
 273       if (outside_all) then
 274          cycle reports
 275       end if
 276
 277       if (print_detail_obs) then
 278          ! Simplistic approach, could be improved to get it all done on PE 0
 279          if (.NOT. outside) then
 280             write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') &
 281                "Report",report," at",platform%info%lon,"E",platform%info%lat, &
 282                "N on proc", myproc," at x/y:", platform%loc%x,platform%loc%y
 283          end if
 284       end if
 285
 286       call da_get_julian_time (iyear,imonth,iday,ihour,imin,analysis_time)
 287       tdiff = abs((obs_time - analysis_time)-0.02)
 288       dlat_earth = platform%info%lat
 289       dlon_earth = platform%info%lon
 290       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 291       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 292       dlat_earth = dlat_earth * deg2rad
 293       dlon_earth = dlon_earth * deg2rad
 294
 295       levs = 1
 296
 297       ! Loop over duplicating obs for global
 298       ndup = 1
 299       if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
 300
 301       ! It is possible that logic for counting obs is incorrect for the
 302       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
 303       ! TBH:  20050913
 304
 305       if (.not.outside) then
 306          if (print_detail_obs .and. ndup > 1) then
 307             write(unit=stdout, fmt = fmt_info) &
 308                platform%info%platform,    &
 309                platform%info%date_char,   &
 310                platform%info%name,        &
 311                platform%info%levels,      &
 312                platform%info%lat,         &
 313                platform%info%lon,         &
 314                platform%info%elv,         &
 315                platform%info%id
 316
 317             write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
 318                ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
 319                platform%loc%i,  platform%loc%j,   &
 320                platform%loc%dx, platform%loc%dxm, &
 321               platform%loc%dy, platform%loc%dym
 322          end if
 323       end if
 324
 325       dup_loop: do n = 1, ndup
 326          is_surface=.true.
 327          if (use_chemic_surfobs) then
 328             if (.not. use_chemic_surfobs .or. ntotal(chemic_surf) == max_synop_input  ) cycle reports
 329             if (n==1) iv%info(chemic_surf)%ntotal = iv%info(chemic_surf)%ntotal + 1
 330             if (outside) cycle reports
 331             if ( thin_conv_ascii ) then
 332                crit = tdiff
 333 !!!               call map2grids_conv(chemic_surf,dlat_earth,dlon_earth,crit,nlocal(chemic_surf),itx,1,itt,ilocal(chemic_surf),iuse)
 334                call map2grids_conv(chemic_surf,dlat_earth,dlon_earth,crit,iv%info(chemic_surf)%nlocal,itx,1,itt,iout,iuse)
 335                if ( .not. iuse ) cycle reports
 336             else
 337                 iv%info(chemic_surf)%nlocal = iv%info(chemic_surf)%nlocal + 1
 338             end if
 339
 340          else
 341
 342             write(unit=message(1), fmt='(a)') 'unsaved obs found:'
 343             write(unit=message(2), fmt='(2a)') &
 344                'platform%info%platform=', platform%info%platform
 345             write(unit=message(3), fmt='(a, i3)') &
 346                'platform%info%levels=', platform%info%levels
 347             call da_warning(__FILE__,__LINE__,message(1:3))
 348             cycle
 349
 350          end if
 351
 352       end do dup_loop
 353
 354    end do reports
 355
 356    close(iunit)
 357
 358    call da_free_unit(iunit)
 359
 360    if (trace_use) call da_trace_exit("da_scan_obs_chem_sfc")
 361
 362 end subroutine da_scan_obs_chem_sfc
 363
 364