subroutine da_read_obs_ssmi (iv, filename) 2,11
!---------------------------------------------------------------------------
! Purpose: Read 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
character (len = 160) :: 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
type (count_obs_number_type) :: count_obs_ssmir
type (count_obs_number_type) :: count_obs_ssmit
logical :: isfilter,ipass
logical :: outside, outside_all
integer :: irain, iprecip
integer :: n, ndup
integer :: nlocal(num_ob_indexes)
integer :: ntotal(num_ob_indexes)
if (trace_use) call da_trace_entry
("da_read_obs_ssmi")
nlocal(:) = iv%info(:)%plocal(iv%time-1)
ntotal(:) = iv%info(:)%ptotal(iv%time-1)
count_obs_ssmir = count_obs_number_type(0, 0, 0, 0)
count_obs_ssmit = count_obs_number_type(0, 0, 0, 0)
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
tb19v % error = tb19v % error + 2.0
tb19h % error = tb19h % error + 2.0
tb22v % error = tb22v % error + 2.0
tb37h % error = tb37h % error + 2.0
tb37v % error = tb37v % error + 2.0
tb85h % error = tb85h % error + 2.0
tb85v % error = tb85v % error + 2.0
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, outside_all)
if (outside_all) cycle reports
loc % pw % inv = missing_r
loc % pw % qc = missing_data
loc % slp = loc % pw
! Loop over duplicating obs for global
ndup = 1
if (global .and. (loc%i < ids .or. loc%i >= ide)) ndup= 2
! It is possible that logic for counting obs is incorrect for the
! global case with >1 MPI tasks due to obs duplication, halo, etc.
! TBH: 20050913
if (.not.outside) then
if (print_detail_obs .and. ndup > 1) then
write(unit=stdout, fmt = fmt_info) &
info%platform, &
info%date_char, &
info%name, &
info%levels, &
info%lat, &
info%lon, &
info%elv, &
info%id
write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
loc%i, loc%j, loc%dx, loc%dxm, loc%dy, loc%dym
end if
end if
dup_loop: do n = 1, ndup
select case(fm)
case (125) ;
if (.not. use_ssmiretrievalobs) cycle reports
if (n==1) ntotal(ssmi_rv) = ntotal(ssmi_rv)
if (outside) cycle reports
! Check if at least one field is present
if ((tpw % qc == missing_data) .AND. (speed % qc == missing_data)) then
count_obs_ssmir % num_missing = count_obs_ssmir % num_missing + 1
cycle reports
end if
! fill permanent structure
! ========================
nlocal(ssmi_rv) = nlocal(ssmi_rv) + 1
! Track serial obs index for reassembly of obs during bit-for-bit
! tests with different numbers of MPI tasks.
loc%obs_global_index = ntotal(ssmi_rv)
! One more data used
count_obs_ssmir % num_used = count_obs_ssmir % num_used + 1
! Initialize other non read fields
iv%info(ssmi_rv)%name(nlocal(ssmi_rv)) = info%name
iv%info(ssmi_rv)%platform(nlocal(ssmi_rv)) = info%platform
iv%info(ssmi_rv)%id(nlocal(ssmi_rv)) = info%id
iv%info(ssmi_rv)%date_char(nlocal(ssmi_rv)) = info%date_char
iv%info(ssmi_rv)%levels(nlocal(ssmi_rv)) = 1
iv%info(ssmi_rv)%lat(:,nlocal(ssmi_rv)) = info%lat
iv%info(ssmi_rv)%lon(:,nlocal(ssmi_rv)) = info%lon
iv%info(ssmi_rv)%elv(nlocal(ssmi_rv)) = info%elv
iv%info(ssmi_rv)%pstar(nlocal(ssmi_rv)) = info%pstar
iv%info(ssmi_rv)%slp(nlocal(ssmi_rv)) = loc%slp
iv%info(ssmi_rv)%pw(nlocal(ssmi_rv)) = loc%pw
iv%info(ssmi_rv)%x(:,nlocal(ssmi_rv)) = loc%x
iv%info(ssmi_rv)%y(:,nlocal(ssmi_rv)) = loc%y
iv%info(ssmi_rv)%i(:,nlocal(ssmi_rv)) = loc%i
iv%info(ssmi_rv)%j(:,nlocal(ssmi_rv)) = loc%j
iv%info(ssmi_rv)%dx(:,nlocal(ssmi_rv)) = loc%dx
iv%info(ssmi_rv)%dxm(:,nlocal(ssmi_rv)) = loc%dxm
iv%info(ssmi_rv)%dy(:,nlocal(ssmi_rv)) = loc%dy
iv%info(ssmi_rv)%dym(:,nlocal(ssmi_rv)) = loc%dym
iv%info(ssmi_rv)%proc_domain(:,nlocal(ssmi_rv)) = loc%proc_domain
iv%info(ssmi_rv)%obs_global_index(nlocal(ssmi_rv)) = ntotal(ssmi_rv)
iv % ssmi_rv (nlocal(ssmi_rv)) % speed = speed
iv % ssmi_rv (nlocal(ssmi_rv)) % tpw = tpw
case (126) ;
if (.not. use_ssmitbobs) cycle reports
if (n==1) ntotal(ssmi_tb) = ntotal(ssmi_tb) + 1
if (outside) cycle reports
! Check if at least one field is present
if ((tb19v % qc == missing_data) .AND. (tb19h % qc == missing_data) .AND. &
(tb22v % qc == missing_data) .AND. &
(tb37v % qc == missing_data) .AND. (tb37h % qc == missing_data) .AND. &
(tb85v % qc == missing_data) .AND. (tb85h % qc == missing_data)) then
count_obs_ssmit % num_missing = &
count_obs_ssmit % num_missing + 1
! write (unit=stdout,fmt=*) 'missing data'
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 .eq. 1) then
tb19v % qc = -88.0
tb19h % qc = -88.0
tb22v % qc = -88.0
tb37h % qc = -88.0
tb37v % qc = -88.0
tb85h % qc = -88.0
tb85v % qc = -88.0
irain = irain + 1
cycle reports
end if
end if
! fill permanent structure
! ========================
! One more data read in
nlocal(ssmi_tb) = nlocal(ssmi_tb) + 1
! Track serial obs index for reassembly of obs during bit-for-bit
! tests with different numbers of MPI tasks.
loc%obs_global_index = ntotal(ssmi_tb)
! One more data used
count_obs_ssmit % num_used = count_obs_ssmit % num_used + 1
! Initialize other non read fields
iv%info(ssmi_tb)%name(nlocal(ssmi_tb)) = info%name
iv%info(ssmi_tb)%platform(nlocal(ssmi_tb)) = info%platform
iv%info(ssmi_tb)%id(nlocal(ssmi_tb)) = info%id
iv%info(ssmi_tb)%date_char(nlocal(ssmi_tb)) = info%date_char
iv%info(ssmi_tb)%levels(nlocal(ssmi_tb)) = 1
iv%info(ssmi_tb)%lat(:,nlocal(ssmi_tb)) = info%lat
iv%info(ssmi_tb)%lon(:,nlocal(ssmi_tb)) = info%lon
iv%info(ssmi_tb)%elv(nlocal(ssmi_tb)) = info%elv
iv%info(ssmi_tb)%pstar(nlocal(ssmi_tb)) = info%pstar
iv%info(ssmi_tb)%slp(nlocal(ssmi_tb)) = loc%slp
iv%info(ssmi_tb)%pw(nlocal(ssmi_tb)) = loc%pw
iv%info(ssmi_tb)%x(:,nlocal(ssmi_tb)) = loc%x
iv%info(ssmi_tb)%y(:,nlocal(ssmi_tb)) = loc%y
iv%info(ssmi_tb)%i(:,nlocal(ssmi_tb)) = loc%i
iv%info(ssmi_tb)%j(:,nlocal(ssmi_tb)) = loc%j
iv%info(ssmi_tb)%dx(:,nlocal(ssmi_tb)) = loc%dx
iv%info(ssmi_tb)%dxm(:,nlocal(ssmi_tb)) = loc%dxm
iv%info(ssmi_tb)%dy(:,nlocal(ssmi_tb)) = loc%dy
iv%info(ssmi_tb)%dym(:,nlocal(ssmi_tb)) = loc%dym
iv%info(ssmi_tb)%proc_domain(:,nlocal(ssmi_tb)) = loc%proc_domain
iv%info(ssmi_tb)%obs_global_index(nlocal(ssmi_tb)) = ntotal(ssmi_tb)
iv % ssmi_tb (nlocal(ssmi_tb)) % tb19v = tb19v
iv % ssmi_tb (nlocal(ssmi_tb)) % tb19h = tb19h
iv % ssmi_tb (nlocal(ssmi_tb)) % tb22v = tb22v
iv % ssmi_tb (nlocal(ssmi_tb)) % tb37v = tb37v
iv % ssmi_tb (nlocal(ssmi_tb)) % tb37h = tb37h
iv % ssmi_tb (nlocal(ssmi_tb)) % tb85v = tb85v
iv % ssmi_tb (nlocal(ssmi_tb)) % tb85h = tb85h
case default;
! Do nothing.
end select
end do dup_loop
end do reports
close(iunit)
call da_free_unit(iunit)
write(unit=stdout, fmt='(/,25x,a)') ' used outdomain max_er_chk missing'
write(unit=stdout, fmt='(4x,a,4i11)') 'ssmi_rv_diag: ', count_obs_ssmir
write(unit=stdout, fmt='(4x,a,4i11)') 'ssmi_tb_diag: ', count_obs_ssmit
if (irain > 0) then
write(unit=stdout, fmt='(/,5x,a,i6/)') '** Rain contaminated SSMI_Tb =', irain
end if
write(unit=stdout, fmt='(/,a)') ' '
if (trace_use) call da_trace_exit("da_read_obs_ssmi")
end subroutine da_read_obs_ssmi