subroutine da_qc_amsr2 (it, i, nchan, ob, iv) 1,13
!---------------------------------------------------------------------------
! Purpose: perform quality control for amsr2 data.
!---------------------------------------------------------------------------
implicit none
integer, intent(in) :: it ! outer loop count
integer, intent(in) :: i ! sensor index.
integer, intent(in) :: nchan ! number of channel
type (y_type), intent(in) :: ob ! Observation structure.
type (iv_type), intent(inout) :: iv ! O-B structure.
real, parameter :: clwcutofx(1:14) = &
(/0.350, 0.350, 0.350, 0.350, 0.350, 0.350, 0.300, &
0.300, 0.250, 0.250, 0.100, 0.100, 0.020, 0.020 /)
! local variables
integer :: n,k,isflg,ios,fgat_rad_unit
real :: si, coscon, sincon, bearaz, sun_zenith, sun_azimuth, sun_glint
integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
nrej_omb_std(nchan), &
nrej_clw(nchan), num_proc_domain, &
nrej_land, nrej_glint
character(len=30) :: filename
if (trace_use) call da_trace_entry
("da_qc_amsr2")
ngood(:) = 0
nrej(:) = 0
nrej_omb_abs(:) = 0
nrej_omb_std(:) = 0
nrej_clw(:) = 0
nrej_land = 0
nrej_glint = 0
num_proc_domain = 0
coscon = cos( (90.0 - 55.0)*deg_to_rad)
sincon = sin( (90.0 - 55.0)*deg_to_rad)
do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
if (iv%instid(i)%info%proc_domain(1,n)) &
num_proc_domain = num_proc_domain + 1
! 0.0 initialise QC by flags assuming good obs
!-----------------------------------------------------------------
iv%instid(i)%tb_qc(:,n) = qc_good
! 1.0 reject all channels over land
!-----------------------------------------------------------------
isflg = iv%instid(i)%isflg(n) !model surface type at ob location
if (isflg > 0 .or. iv%instid(i)%landsea_mask(n) == 0) then
iv%instid(i)%tb_qc(:,n) = qc_bad
nrej_land = nrej_land + 1
end if
! 2.0 check sun_glint angle for 6.9, 7.3 and 10.6 GHz
! It should be >= 25 degrees
!-----------------------------------------------------------------
bearaz = 270.0 - iv%instid(i)%satazi(n)
sun_zenith = iv%instid(i)%solzen(n)
sun_azimuth = 90.0 - iv%instid(i)%solazi(n)
sun_glint = acos(coscon * cos( bearaz ) * cos( sun_zenith ) * cos( sun_azimuth ) + &
coscon * sin( bearaz ) * cos( sun_zenith ) * sin( sun_azimuth ) + &
sincon * sin( sun_zenith )) * rad_to_deg
if (sun_glint < 25.0) then
! apply only to 6.9,7.3 and 10.6 GHz for both V & H polarizations
iv%instid(i)%tb_qc(1:6,n) = qc_bad
nrej_glint = nrej_glint + 1
end if
! 3.0 check iuse
!-----------------------------------------------------------------
do k = 1, nchan
if (satinfo(i)%iuse(k) .eq. -1) &
iv%instid(i)%tb_qc(k,n) = qc_bad
end do
! 4.0 check cloud
!-----------------------------------------------------------------
if (.not. crtm_cloud ) then
do k = 1, nchan
! clw check
! channel dependent criteria
if( iv%instid(i)%clw(n) < 0 .or. & !bad obs retrieved clw
iv%instid(i)%clw(n) > clwcutofx(k) .or. & !obs retrieved clw
iv%instid(i)%clwp(n) > clwcutofx(k) ) then !model/guess clwp
iv%instid(i)%tb_qc(k,n) = qc_bad
if (iv%instid(i)%info%proc_domain(1,n)) &
nrej_clw(k) = nrej_clw(k) + 1
end if
end do
end if
! 5.0 check innovation
!-----------------------------------------------------------------
do k = 1, nchan
! absolute departure check
if ( k <= 7 .or. k == 11 .or. k == 12) then
if (abs(iv%instid(i)%tb_inv(k,n)) > 6.0) then
iv%instid(i)%tb_qc(k,n) = qc_bad
if (iv%instid(i)%info%proc_domain(1,n)) &
nrej_omb_abs(k) = nrej_omb_abs(k) + 1
end if
else if ( k == 8 .or. k == 9 ) then
if (abs(iv%instid(i)%tb_inv(k,n)) > 8.0) then
iv%instid(i)%tb_qc(k,n) = qc_bad
if (iv%instid(i)%info%proc_domain(1,n)) &
nrej_omb_abs(k) = nrej_omb_abs(k) + 1
end if
else
if (abs(iv%instid(i)%tb_inv(k,n)) > 10.0) then
iv%instid(i)%tb_qc(k,n) = qc_bad
if (iv%instid(i)%info%proc_domain(1,n)) &
nrej_omb_abs(k) = nrej_omb_abs(k) + 1
end if
end if
! relative departure check
if (use_error_factor_rad) then
iv%instid(i)%tb_error(k,n) = &
satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
else
iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
end if
if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
iv%instid(i)%tb_qc(k,n) = qc_bad
if (iv%instid(i)%info%proc_domain(1,n)) &
nrej_omb_std(k) = nrej_omb_std(k) + 1
end if
! final QC decsion
if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
iv%instid(i)%tb_error(k,n) = 500.0
if (iv%instid(i)%info%proc_domain(1,n)) &
nrej(k) = nrej(k) + 1
else
if (iv%instid(i)%info%proc_domain(1,n)) &
ngood(k) = ngood(k) + 1
end if
end do ! nchan
end do ! end loop pixel
! Do inter-processor communication to gather statistics.
call da_proc_sum_int (num_proc_domain)
call da_proc_sum_ints (nrej_clw)
call da_proc_sum_int (nrej_land)
call da_proc_sum_int (nrej_glint)
call da_proc_sum_ints (nrej_omb_abs)
call da_proc_sum_ints (nrej_omb_std)
call da_proc_sum_ints (nrej)
call da_proc_sum_ints (ngood)
if (rootproc) then
if (num_fgat_time > 1) then
write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
else
write(filename,'(i2.2,a)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)
end if
call da_get_unit(fgat_rad_unit)
open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
if (ios /= 0) then
write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
call da_error(__FILE__,__LINE__,message(1:1))
end if
write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
if(num_proc_domain > 0) write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
if(nrej_land > 0)write(fgat_rad_unit,'(a20,i7)') ' nrej_land = ', nrej_land
if(nrej_glint> 0)write(fgat_rad_unit,'(a20,i7)') ' nrej_glint = ', nrej_glint
write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
write(fgat_rad_unit,'(a20)') ' nrej_clw(:) = '
write(fgat_rad_unit,'(10i7)') nrej_clw(:)
write(fgat_rad_unit,'(a20)') ' nrej(:) = '
write(fgat_rad_unit,'(10i7)') nrej(:)
write(fgat_rad_unit,'(a20)') ' ngood(:) = '
write(fgat_rad_unit,'(10i7)') ngood(:)
close(fgat_rad_unit)
call da_free_unit(fgat_rad_unit)
end if
if (trace_use) call da_trace_exit("da_qc_amsr2")
end subroutine da_qc_amsr2