var/da/da_ssmi/da_read_obs_ssmi.inc

   1 subroutine da_read_obs_ssmi (iv, filename)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Read a SSMI 2.0 GTS observation file
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type(iv_type),    intent(inout) :: iv
  10    character(len=*), intent(in)    :: filename
  11
  12    character (len =  10)        :: fmt_name
  13
  14    character (len = 160)        :: fmt_info
  15    character (len = 160)        :: fmt_loc
  16    character (len = 120)        :: char_ned
  17
  18    integer                      :: iost, fm,iunit
  19
  20    type (model_loc_type)        :: loc
  21    type (info_type)             :: info
  22    type (field_type)            :: speed, tpw
  23
  24    type (field_type)            :: tb19v, tb19h, tb22v
  25    type (field_type)            :: tb37v, tb37h, tb85v, tb85h
  26
  27    type (count_obs_number_type) :: count_obs_ssmir
  28    type (count_obs_number_type) :: count_obs_ssmit
  29
  30    logical                      :: isfilter,ipass
  31    logical                      :: outside, outside_all
  32    integer                      :: irain, iprecip
  33    integer                      :: n, ndup
  34    integer                      :: nlocal(num_ob_indexes)
  35    integer                      :: ntotal(num_ob_indexes)
  36
  37    if (trace_use) call da_trace_entry("da_read_obs_ssmi")
  38
  39    nlocal(:) = iv%info(:)%plocal(iv%time-1)
  40    ntotal(:) = iv%info(:)%ptotal(iv%time-1)
  41
  42    count_obs_ssmir = count_obs_number_type(0, 0, 0, 0)
  43    count_obs_ssmit = count_obs_number_type(0, 0, 0, 0)
  44
  45    isfilter = .true. ! filter out rain points
  46    irain = 0
  47
  48    ! open file
  49
  50    call da_get_unit(iunit)
  51    open(unit   = iunit,     &
  52       FILE   = trim(filename), &
  53       FORM   = 'FORMATTED',  &
  54       ACCESS = 'SEQUENTIAL', &
  55       iostat =  iost,     &
  56       STATUS = 'OLD')
  57
  58    if (iost /= 0) then
  59       call da_warning(__FILE__,__LINE__, (/"Cannot open SSMI file "//filename/))
  60       call da_free_unit(iunit)
  61       return
  62    end if
  63
  64    rewind (unit = iunit)
  65
  66    ! 2.  read header
  67    ! ===============
  68
  69    ! 2.1 read first line
  70    !     ---------------
  71
  72    read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  73    if (iost /= 0) then
  74       call da_error(__FILE__,__LINE__, (/"Cannot read SSMI file"//filename/))
  75    end if
  76
  77    ! 2.3 read number of reports
  78    !     ---------------------
  79
  80    do
  81       read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  82       if (iost /= 0) then
  83          call da_error(__FILE__,__LINE__, (/"Cannot read SSMI file"//filename/))
  84       end if
  85       if (char_ned(1:6) == 'NESTIX') exit
  86    end do
  87
  88    do
  89      read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  90      if (char_ned(1:6) == 'INFO  ') exit
  91    end do
  92
  93    read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  94
  95    ! read formats
  96    ! ------------
  97
  98    read (unit=iunit, fmt = '(A,1X,A)') fmt_name, fmt_info, fmt_name, fmt_loc
  99
 100    !  skip 1 line
 101    !  -----------
 102
 103    read (unit=iunit, fmt = '(A)') fmt_name
 104
 105    !  loop over records
 106    !  -----------------
 107
 108    reports: do
 109       ! read station general info
 110       ! =========================
 111
 112       read (unit=iunit, fmt = fmt_info, iostat = iost) &
 113          info % platform,    &
 114          info % date_char,   &
 115          info % name,        &
 116          info % levels,      &
 117          info % lat,         &
 118          info % lon,         &
 119          info % elv,         &
 120          info % id
 121
 122       read(unit=info % platform (4:6),fmt='(I3)') fm
 123       if (iost /= 0) exit reports
 124
 125       select case(fm)
 126          case (125)    ;
 127             ! read surface wind speed and precipitable water
 128             read (unit=iunit, fmt = fmt_loc) speed%inv, speed%qc, speed%error, &
 129                                              tpw%inv, tpw%qc, tpw%error
 130          case (126)    ;
 131             read (unit=iunit, fmt = fmt_loc) &
 132                tb19v%inv, tb19v%qc, tb19v%error, &
 133                tb19h%inv, tb19h%qc, tb19h%error, &
 134                tb22v%inv, tb22v%qc, tb22v%error, &
 135                tb37v%inv, tb37v%qc, tb37v%error, &
 136                tb37h%inv, tb37h%qc, tb37h%error, &
 137                tb85v%inv, tb85v%qc, tb85v%error, &
 138                tb85h%inv, tb85h%qc, tb85h%error
 139
 140                tb19v % error = tb19v % error + 2.0
 141                tb19h % error = tb19h % error + 2.0
 142                tb22v % error = tb22v % error + 2.0
 143                tb37h % error = tb37h % error + 2.0
 144                tb37v % error = tb37v % error + 2.0
 145                tb85h % error = tb85h % error + 2.0
 146                tb85v % error = tb85v % error + 2.0
 147
 148          case default;
 149             write(unit=message(1), fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
 150             write(unit=message(2), fmt='(a, 2f12.6)') &
 151                'info%(lon,lat)=', info%lon, info%lat
 152             call da_warning(__FILE__,__LINE__,message(1:2))
 153       end select
 154
 155       ! check if obs is in horizontal domain
 156       ! ====================================
 157
 158       ! Compute the model horizontal coordinate x, y
 159       ! Check if obs is wihin horizontal domain
 160
 161       call da_llxy (info, loc, outside, outside_all)
 162
 163       if (outside_all) cycle reports
 164
 165       loc % pw  % inv = missing_r
 166       loc % pw  % qc  = missing_data
 167       loc % slp       = loc % pw
 168
 169       ! Loop over duplicating obs for global
 170       ndup = 1
 171       if (global .and. (loc%i < ids .or. loc%i >= ide)) ndup= 2
 172
 173       ! It is possible that logic for counting obs is incorrect for the
 174       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
 175       ! TBH:  20050913
 176
 177       if (.not.outside) then
 178          if (print_detail_obs .and. ndup > 1) then
 179             write(unit=stdout, fmt = fmt_info) &
 180                info%platform,    &
 181                info%date_char,   &
 182                info%name,        &
 183                info%levels,      &
 184                info%lat,         &
 185                info%lon,         &
 186                info%elv,         &
 187                info%id
 188
 189             write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
 190                ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
 191                loc%i, loc%j, loc%dx, loc%dxm, loc%dy, loc%dym
 192          end if
 193       end if
 194
 195       dup_loop: do n = 1, ndup
 196
 197       select case(fm)
 198          case (125) ;
 199             if (.not. use_ssmiretrievalobs) cycle reports
 200
 201             if (n==1) ntotal(ssmi_rv) = ntotal(ssmi_rv)
 202             if (outside) cycle reports
 203
 204             ! Check if at least one field is present
 205             if ((tpw % qc == missing_data) .AND. (speed % qc == missing_data)) then
 206                count_obs_ssmir % num_missing = count_obs_ssmir % num_missing + 1
 207                cycle reports
 208             end if
 209
 210             ! fill permanent structure
 211             ! ========================
 212
 213             nlocal(ssmi_rv) = nlocal(ssmi_rv) + 1
 214             ! Track serial obs index for reassembly of obs during bit-for-bit
 215             ! tests with different numbers of MPI tasks.
 216             loc%obs_global_index = ntotal(ssmi_rv)
 217
 218             !  One more data used
 219
 220             count_obs_ssmir % num_used = count_obs_ssmir % num_used + 1
 221
 222             !  Initialize other non read fields
 223
 224             iv%info(ssmi_rv)%name(nlocal(ssmi_rv))      = info%name
 225             iv%info(ssmi_rv)%platform(nlocal(ssmi_rv))  = info%platform
 226             iv%info(ssmi_rv)%id(nlocal(ssmi_rv))        = info%id
 227             iv%info(ssmi_rv)%date_char(nlocal(ssmi_rv)) = info%date_char
 228             iv%info(ssmi_rv)%levels(nlocal(ssmi_rv))    = 1
 229             iv%info(ssmi_rv)%lat(:,nlocal(ssmi_rv))     = info%lat
 230             iv%info(ssmi_rv)%lon(:,nlocal(ssmi_rv))     = info%lon
 231             iv%info(ssmi_rv)%elv(nlocal(ssmi_rv))       = info%elv
 232             iv%info(ssmi_rv)%pstar(nlocal(ssmi_rv))     = info%pstar
 233
 234             iv%info(ssmi_rv)%slp(nlocal(ssmi_rv))           = loc%slp
 235             iv%info(ssmi_rv)%pw(nlocal(ssmi_rv))            = loc%pw
 236             iv%info(ssmi_rv)%x(:,nlocal(ssmi_rv))           = loc%x
 237             iv%info(ssmi_rv)%y(:,nlocal(ssmi_rv))           = loc%y
 238             iv%info(ssmi_rv)%i(:,nlocal(ssmi_rv))           = loc%i
 239             iv%info(ssmi_rv)%j(:,nlocal(ssmi_rv))           = loc%j
 240             iv%info(ssmi_rv)%dx(:,nlocal(ssmi_rv))          = loc%dx
 241             iv%info(ssmi_rv)%dxm(:,nlocal(ssmi_rv))         = loc%dxm
 242             iv%info(ssmi_rv)%dy(:,nlocal(ssmi_rv))          = loc%dy
 243             iv%info(ssmi_rv)%dym(:,nlocal(ssmi_rv))         = loc%dym
 244             iv%info(ssmi_rv)%proc_domain(:,nlocal(ssmi_rv)) = loc%proc_domain
 245
 246             iv%info(ssmi_rv)%obs_global_index(nlocal(ssmi_rv)) = ntotal(ssmi_rv)
 247
 248             iv % ssmi_rv (nlocal(ssmi_rv)) % speed = speed
 249             iv % ssmi_rv (nlocal(ssmi_rv)) % tpw   = tpw
 250
 251          case (126) ;
 252             if (.not. use_ssmitbobs) cycle reports
 253
 254             if (n==1) ntotal(ssmi_tb) = ntotal(ssmi_tb) + 1
 255             if (outside) cycle reports
 256
 257             ! Check if at least one field is present
 258
 259             if ((tb19v % qc == missing_data) .AND. (tb19h % qc == missing_data)  .AND. &
 260                 (tb22v % qc == missing_data)                                .AND. &
 261                 (tb37v % qc == missing_data) .AND. (tb37h % qc == missing_data)  .AND. &
 262                 (tb85v % qc == missing_data) .AND. (tb85h % qc == missing_data)) then
 263                count_obs_ssmit % num_missing = &
 264                count_obs_ssmit % num_missing + 1
 265                ! write (unit=stdout,fmt=*) 'missing data'
 266                cycle reports
 267             end if
 268
 269             ! filter rain pixels
 270             !  ====================================
 271
 272             if (isfilter) then
 273                ipass = .false.
 274                iprecip = 0
 275                call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
 276                   tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, info%lat)
 277                if (iprecip .eq. 1) then
 278                   tb19v % qc    = -88.0
 279                   tb19h % qc    = -88.0
 280                   tb22v % qc    = -88.0
 281                   tb37h % qc    = -88.0
 282                   tb37v % qc    = -88.0
 283                   tb85h % qc    = -88.0
 284                   tb85v % qc    = -88.0
 285                   irain = irain + 1
 286                   cycle reports
 287                end if
 288             end if
 289
 290             ! fill permanent structure
 291             ! ========================
 292
 293             ! One more data read in
 294
 295             nlocal(ssmi_tb) = nlocal(ssmi_tb) + 1
 296             ! Track serial obs index for reassembly of obs during bit-for-bit
 297             ! tests with different numbers of MPI tasks.
 298             loc%obs_global_index = ntotal(ssmi_tb)
 299
 300             !  One more data used
 301
 302             count_obs_ssmit % num_used = count_obs_ssmit % num_used + 1
 303
 304             !  Initialize other non read fields
 305
 306             iv%info(ssmi_tb)%name(nlocal(ssmi_tb))      = info%name
 307             iv%info(ssmi_tb)%platform(nlocal(ssmi_tb))  = info%platform
 308             iv%info(ssmi_tb)%id(nlocal(ssmi_tb))        = info%id
 309             iv%info(ssmi_tb)%date_char(nlocal(ssmi_tb)) = info%date_char
 310             iv%info(ssmi_tb)%levels(nlocal(ssmi_tb))    = 1
 311             iv%info(ssmi_tb)%lat(:,nlocal(ssmi_tb))     = info%lat
 312             iv%info(ssmi_tb)%lon(:,nlocal(ssmi_tb))     = info%lon
 313             iv%info(ssmi_tb)%elv(nlocal(ssmi_tb))       = info%elv
 314             iv%info(ssmi_tb)%pstar(nlocal(ssmi_tb))     = info%pstar
 315
 316             iv%info(ssmi_tb)%slp(nlocal(ssmi_tb))           = loc%slp
 317             iv%info(ssmi_tb)%pw(nlocal(ssmi_tb))            = loc%pw
 318             iv%info(ssmi_tb)%x(:,nlocal(ssmi_tb))           = loc%x
 319             iv%info(ssmi_tb)%y(:,nlocal(ssmi_tb))           = loc%y
 320             iv%info(ssmi_tb)%i(:,nlocal(ssmi_tb))           = loc%i
 321             iv%info(ssmi_tb)%j(:,nlocal(ssmi_tb))           = loc%j
 322             iv%info(ssmi_tb)%dx(:,nlocal(ssmi_tb))          = loc%dx
 323             iv%info(ssmi_tb)%dxm(:,nlocal(ssmi_tb))         = loc%dxm
 324             iv%info(ssmi_tb)%dy(:,nlocal(ssmi_tb))          = loc%dy
 325             iv%info(ssmi_tb)%dym(:,nlocal(ssmi_tb))         = loc%dym
 326             iv%info(ssmi_tb)%proc_domain(:,nlocal(ssmi_tb)) = loc%proc_domain
 327
 328             iv%info(ssmi_tb)%obs_global_index(nlocal(ssmi_tb)) = ntotal(ssmi_tb)
 329
 330             iv % ssmi_tb (nlocal(ssmi_tb)) % tb19v = tb19v
 331             iv % ssmi_tb (nlocal(ssmi_tb)) % tb19h = tb19h
 332             iv % ssmi_tb (nlocal(ssmi_tb)) % tb22v = tb22v
 333             iv % ssmi_tb (nlocal(ssmi_tb)) % tb37v = tb37v
 334             iv % ssmi_tb (nlocal(ssmi_tb)) % tb37h = tb37h
 335             iv % ssmi_tb (nlocal(ssmi_tb)) % tb85v = tb85v
 336             iv % ssmi_tb (nlocal(ssmi_tb)) % tb85h = tb85h
 337
 338          case default;
 339             ! Do nothing.
 340       end select
 341       end do dup_loop
 342    end do reports
 343
 344    close(iunit)
 345    call da_free_unit(iunit)
 346
 347    write(unit=stdout, fmt='(/,25x,a)')   '     used   outdomain  max_er_chk   missing'
 348    write(unit=stdout, fmt='(4x,a,4i11)') 'ssmi_rv_diag:        ', count_obs_ssmir
 349    write(unit=stdout, fmt='(4x,a,4i11)') 'ssmi_tb_diag:        ', count_obs_ssmit
 350
 351    if (irain > 0) then
 352       write(unit=stdout, fmt='(/,5x,a,i6/)') '** Rain contaminated SSMI_Tb =', irain
 353    end if
 354
 355    write(unit=stdout, fmt='(/,a)') ' '
 356
 357    if (trace_use) call da_trace_exit("da_read_obs_ssmi")
 358
 359 end subroutine da_read_obs_ssmi
 360
 361