subroutine da_scan_obs_ssmi (iv, filename) 2,11
!---------------------------------------------------------------------------
! Purpose: Read the header of a SSMI 2.0 GTS observation file
!---------------------------------------------------------------------------
implicit none
type(iv_type), intent(inout) :: iv
character(len=*), intent(in) :: filename
character (len = 10) :: fmt_name
character (len = 160) :: fmt_info, fmt_loc
character (len = 120) :: char_ned
integer :: iost, fm,iunit
type (model_loc_type) :: loc
type (info_type) :: info
type (field_type) :: speed, tpw
type (field_type) :: tb19v, tb19h, tb22v
type (field_type) :: tb37v, tb37h, tb85v, tb85h
logical :: isfilter,ipass
logical :: outside
integer :: irain, iprecip
if (trace_use) call da_trace_entry
("da_scan_obs_ssmi")
isfilter = .true. ! filter out rain points
irain = 0
! open FILE
call da_get_unit
(iunit)
open(unit = iunit, &
FILE = trim(filename), &
FORM = 'FORMATTED', &
ACCESS = 'SEQUENTIAL', &
iostat = iost, &
STATUS = 'OLD')
if (iost /= 0) then
call da_warning
(__FILE__,__LINE__, (/"Cannot open SSMI file "//filename/))
call da_free_unit
(iunit)
return
end if
rewind (unit = iunit)
! 2. read header
! ===============
! 2.1 read first line
! ---------------
read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
if (iost /= 0) then
call da_error
(__FILE__,__LINE__, (/"Cannot read SSMI file "//filename/))
end if
! 2.3 read number of reports
do
read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
if (iost /= 0) then
call da_error
(__FILE__,__LINE__, (/"Cannot read SSMI file "//filename/))
end if
if (char_ned(1:6) == 'NESTIX') exit
end do
do
read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
if (char_ned(1:6) == 'INFO ') exit
end do
read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
! read formats
read (unit=iunit, fmt = '(A,1X,A)') &
fmt_name, fmt_info, &
fmt_name, fmt_loc
! skip 1 line
read (unit=iunit, fmt = '(A)') fmt_name
! loop over records
reports: do
! read station general info
read (unit=iunit, fmt = fmt_info, iostat = iost) &
info % platform, &
info % date_char, &
info % name, &
info % levels, &
info % lat, &
info % lon, &
info % elv, &
info % id
read(unit=info % platform (4:6),fmt='(I3)') fm
if (iost /= 0) exit reports
select case(fm)
case (125) ;
! read surface wind speed and precipitable water
read (unit=iunit, fmt = fmt_loc) &
speed%inv, speed%qc, speed%error, &
tpw%inv, tpw%qc, tpw%error
case (126) ;
read (unit=iunit, fmt = fmt_loc) &
tb19v%inv, tb19v%qc, tb19v%error, &
tb19h%inv, tb19h%qc, tb19h%error, &
tb22v%inv, tb22v%qc, tb22v%error, &
tb37v%inv, tb37v%qc, tb37v%error, &
tb37h%inv, tb37h%qc, tb37h%error, &
tb85v%inv, tb85v%qc, tb85v%error, &
tb85h%inv, tb85h%qc, tb85h%error
case default;
write(unit=message(1), fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
write(unit=message(2), fmt='(a, 2f12.6)') 'info%(lon,lat)=', info%lon, info%lat
call da_warning
(__FILE__,__LINE__,message(1:2))
end select
! check if obs is in horizontal domain
! Compute the model horizontal coordinate x, y
! Check if obs is wihin horizontal domain
call da_llxy
(info, loc, outside)
select case(fm)
case (125) ;
if (.not. use_ssmiretrievalobs .or. iv%info(ssmi_rv)%ntotal == max_ssmi_rv_input) cycle reports
! Check if at least one field is present
if ((tpw % qc == missing) .AND. (speed % qc == missing)) then
cycle reports
end if
iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
if (outside) cycle reports
iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
case (126) ;
if (.not. use_ssmitbobs .or. iv%info(ssmi_tb)%ntotal == max_ssmi_tb_input) cycle reports
! Check if at least one field is present
if ((tb19v % qc == missing) .AND. (tb19h % qc == missing) .AND. &
(tb22v % qc == missing) .AND. &
(tb37v % qc == missing) .AND. (tb37h % qc == missing) .AND. &
(tb85v % qc == missing) .AND. (tb85h % qc == missing)) then
cycle reports
end if
! filter rain pixels
! ====================================
if (isfilter) then
ipass = .false.
iprecip = 0
call filter
(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
info%lat)
if (iprecip == 1) then
irain = irain + 1
cycle reports
end if
end if
iv%info(ssmi_tb)%ntotal = iv%info(ssmi_tb)%ntotal + 1
if (outside) cycle reports
iv%info(ssmi_tb)%nlocal = iv%info(ssmi_tb)%nlocal + 1
case default;
! Do nothing.
end select
end do reports
iv%info(ssmt1)%max_lev = 1
iv%info(ssmt2)%max_lev = 1
iv%info(ssmi_tb)%max_lev = 1
iv%info(ssmi_rv)%max_lev = 1
close(unit=iunit)
call da_free_unit
(iunit)
if (irain > 0) then
write(unit=stdout, fmt='(/,5x,a,i6/)') 'Rejected rain contaminated ssmi_tb =', irain
write(unit=stdout, fmt='(A)') ' '
end if
if (trace_use) call da_trace_exit
("da_scan_obs_ssmi")
end subroutine da_scan_obs_ssmi