subroutine read_aerosol(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
obstype,twind,sis,ithin,rmesh, &
mype,mype_root,mype_sub,npe_sub,mpi_comm_sub)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_aerosol read aerosol data
! prgmmr: hchuang org: np23 date: 2009-01-26
!
! abstract: This routine reads MODIS aerosol total column AOD observations.
! ONLY total column values are read in. The routine has
! the ability to read both IEEE and BUFR format MODIS
! aerosol data files.
!
! When running the gsi in regional mode, the code only
! retains those observations that fall within the regional
! domain
!
! program history log:
! 2009-04-08 Huang - modified from read_ozone to read in MODIS AEROSOL data
! 2010-10-20 hclin - modified for total aod in channels
! 2011-01-05 hclin - added three more BUFR records (STYP DBCF QAOD)
!
! input argument list:
! obstype - observation type to process
! jsatid - satellite id to read
! infile - unit from which to read aerosol data
! gstime - analysis time in minutes from reference date
! lunout - unit to which to write data for further processing
! obstype - observation type to process
! twind - input group time window (hours)
! sis - satellite/instrument/sensor indicator
! ithin - flag to thin data
! rmesh - thinning mesh size (km)
! mype - mpi task id
! mype_root - "root" task for sub-communicator
! mype_sub - mpi task id within sub-communicator
! npe_sub - number of data read tasks
! mpi_comm_sub - sub-communicator for data read
!
! output argument list:
! nread - number of modis aerosol observations read
! ndata - number of modis aerosol profiles retained for further processing
! nodata - number of modis aerosol observations retained for further processing
!
! remarks:
!
! attributes:
! language: f90
! machine: IBM AIX Cirrus
!
!$$$
use kinds, only: r_kind, r_double, i_kind
use gridmod, only: nlat, nlon, regional, tll2xy, rlats, rlons
use aod_mod, only: aod_qa_limit, luse_deepblue
use constants, only: deg2rad, zero, two, three, four, rad2deg, r60inv
use obsmod, only: iadate, rmiss_single
use gsi_4dvar, only: l4dvar,iwinbgn,winlen
use satthin, only: itxmax,makegrids,destroygrids,checkob, &
finalcheck,map2tgrid,score_crit
implicit none
!
! Declare local parameters
real(r_kind), parameter :: r6 = 6.0_r_kind
real(r_kind), parameter :: r360 = 360.0_r_kind
!
! Declare passed variables
!
character(10), intent(in) :: obstype, infile, jsatid
character(20), intent(in) :: sis
integer(i_kind), intent(in) :: lunout, ithin
integer(i_kind), intent(inout) :: nread
integer(i_kind), intent(inout) :: ndata, nodata
integer(i_kind) ,intent(in) :: mype
integer(i_kind) ,intent(in) :: mype_root
integer(i_kind) ,intent(in) :: mype_sub
integer(i_kind) ,intent(in) :: npe_sub
integer(i_kind) ,intent(in) :: mpi_comm_sub
real(r_kind), intent(in) :: gstime, twind, rmesh
!
! Declare local variables
!
logical :: outside, iuse
character (len= 8) :: subset
character (len=10) :: date
integer(i_kind) :: naerodat, next, ireadmg, ireadsb
integer(i_kind) :: idate, jdate, ksatid, kk, iy, iret, im, ihh, idd
integer(i_kind) :: lunin = 10
integer(i_kind) :: nmind, i, n
integer(i_kind) :: imin, isec
integer(i_kind) :: k, ilat, ilon, nreal, nchanl
integer(i_kind) :: kidsat
integer(i_kind) :: JULIAN, IDAYYR, IDAYWK
integer(i_kind), dimension(5) :: idate5
!
!| NC008041 | SAID AEROSOL CLONH CLATH YYMMDD HHMMSS SOZA SOLAZI |
!| NC008041 | SCATTA OPTD AEROTP |
!
!| YYMMDD | YEAR MNTH DAYS |
!| | |
!| HHMMSS | HOUR MINU SECO |
!
! SAID Satellite identifier code table (eg, 783 == 'TERRA')
! AEROSOL Aerosol Optical Depth (AOD) source code table (eg, 5 == 'AATSR' )
! YEAR Year
! MNTH Month
! DAYS Day
! HOUR Hour
! MINU Minute
! SECO Second
! CLATH Latitude (high accuracy) degree (5 decimal precision)
! CLONH Longitude (high accuracy) degree (5 decimal precision)
! SOLAZI Solar azimuth degree (2 decimal precision)
! SOZA Solar zenith angle degree (2 decimal precision)
! OPTD Optical depth numeric
! SCATTA Scattering angle degree (2 decimal precsion)
! AEROTP Aerosol type land code table (eg, 1 == 'DUST', 2 == 'SULFATE')
!
! 0-15-195 - AEROTP (Aerosol land type)
!
! CODE DESCRIPTION
! ==== ===========
! 0 Mixed
! 1 Dust
! 2 Sulfate
! 3 Smoke
! 4 Heavy absorbing smoke
! 5-14 Reserved
! 15 Missing value
!
character (len= 4) :: aerostr = 'OPTD'
character (len=53) :: aerogstr = &
'SAID CLATH CLONH YEAR MNTH DAYS HOUR MINU SOZA SOLAZI'
character (len=14) :: flagstr = 'STYP DBCF QAOD'
integer(i_kind) :: itx, itt
real(r_kind) :: tdiff, sstime, slons, slats, dlon, dlat, t4dv, toq, poq, timedif, crit1, dist1
real(r_kind) :: slons0, slats0, rsat, toto3, solzen, azimuth, dlat_earth, dlon_earth
real(r_kind) :: styp, dbcf, qaod
real(r_kind),dimension(0:4):: rlndsea
real(r_kind), allocatable, dimension(:,:) :: aeroout
real(r_kind), allocatable, dimension(:) :: dataaod
real(r_double), dimension( 10) :: hdraerog
real(r_double) :: aod_550
real(r_double), dimension(3) :: aod_flags
!**************************************************************************
! Set constants. Initialize variables
rsat=999._r_kind
! output position of LON and LAT
ilon=3
ilat=4
nread = 0
ndata = 0
nodata = 0
! Set rlndsea for types we would prefer selecting
rlndsea(0) = zero ! styp 0: water
rlndsea(1) = 15._r_kind ! styp 1: coast
rlndsea(2) = 20._r_kind ! styp 2: desert
rlndsea(3) = 10._r_kind ! styp 3: land
rlndsea(4) = 25._r_kind ! styp 4: deep blue
! Make thinning grids
call makegrids(rmesh,ithin)
if ( obstype == 'modis_aod' ) then
!
open(lunin,file=infile,form='unformatted')
call openbf(lunin,'IN',lunin)
call datelen(10)
call readmg(lunin,subset,idate,iret)
if ( iret == 0 ) then
!
if (subset == 'NC008041') then
write(6,*)'READ_AEROSOL: MODIS data type, subset = ',subset
! Set dependent variables and allocate arrays
nreal=11 !9
nchanl=19
naerodat=nreal+nchanl
allocate (aeroout(naerodat,itxmax))
allocate (dataaod(nchanl))
iy = 0
im = 0
idd= 0
ihh= 0
write(date,'( i10)') idate
read (date,'(i4,3i2)') iy,im,idd,ihh
write(6,'(''READ_AEROSOL: aerosol bufr file '',a,'' date is '',i4,4i2.2,a)')trim(infile),iy,im,idd,ihh
read_modis: do
call readsb(lunin,iret)
if (iret/=0) then
call readmg(lunin,subset,jdate,iret)
if (iret/=0) exit read_modis
cycle read_modis
endif
! extract header information
call ufbint(lunin,hdraerog,10,1,iret,aerogstr)
rsat = hdraerog(1); ksatid=rsat
if ( jsatid == 'terra' ) kidsat = 783
if ( jsatid == 'aqua' ) kidsat = 784
if ( ksatid /= kidsat ) cycle read_modis
! Convert observation location to radians
slats0= hdraerog(2)
slons0= hdraerog(3)
if(slons0< zero) slons0=slons0+r360
if(slons0>=r360) slons0=slons0-r360
dlat_earth = slats0 * deg2rad
dlon_earth = slons0 * deg2rad
if(regional)then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
if(outside) cycle read_modis
else
dlat = dlat_earth
dlon = dlon_earth
call grdcrd(dlat,1,rlats,nlat,1)
call grdcrd(dlon,1,rlons,nlon,1)
endif
solzen = hdraerog(9)
azimuth = hdraerog(10)
! Convert observation time to relative time
idate5(1) = hdraerog(4) !year
idate5(2) = hdraerog(5) !month
idate5(3) = hdraerog(6) !day
idate5(4) = hdraerog(7) !hour
idate5(5) = hdraerog(8) !minute
! extract total column aod 1 value 'OPTD' as defined in aerostr
call ufbint(lunin,aod_550,1,1,iret,aerostr)
call w3fs21(idate5,nmind)
t4dv=real((nmind-iwinbgn),r_kind)*r60inv
if (l4dvar) then
if(t4dv<zero .OR. t4dv>winlen) cycle read_modis
else
sstime=real(nmind,r_kind)
tdiff=(sstime-gstime)*r60inv
if ( abs(tdiff) > twind ) cycle read_modis
end if
nread = nread + 1 !nread = nread + nchanl
if (l4dvar) then
timedif = zero
else
timedif = two*abs(tdiff) ! range: 0 to 6
endif
crit1 = 0.01_r_kind + timedif
if ( aod_550 > 1.0e+10_r_double ) cycle read_modis
! extract STYP, DBCF, and QAOD
call ufbint(lunin,aod_flags,3,1,iret,flagstr)
styp = rmiss_single
dbcf = rmiss_single
qaod = rmiss_single
if ( aod_flags(1) < 1.0e+10_r_double ) styp = aod_flags(1)
if ( aod_flags(2) < 1.0e+10_r_double ) dbcf = aod_flags(2)
if ( aod_flags(3) < 1.0e+10_r_double ) qaod = aod_flags(3)
if ( .not. luse_deepblue .and. nint(styp)==4 ) cycle read_modis
if ( qaod < aod_qa_limit ) cycle read_modis
! Map obs to thinning grid
call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis)
if ( .not. iuse ) cycle read_modis
if ( (styp > rmiss_single) .and. (styp >= zero .and. styp <= four) ) then
crit1 = crit1 + rlndsea(nint(styp))
end if
!if ( (dbcf > rmiss_single) .and. (dbcf >= zero .and. dbcf <= three) ) then
! crit1 = crit1 + 10.0_r_kind*(four-dbcf)
!end if
if ( (qaod > rmiss_single) .and. (qaod >= aod_qa_limit .and. qaod <= three) ) then
crit1 = crit1 + 10.0_r_kind*(four-qaod)
end if
call checkob(dist1,crit1,itx,iuse)
if ( .not. iuse ) cycle read_modis
! Compute "score" for observation. All scores>=0.0. Lowest score is "best"
call finalcheck(dist1,crit1,itx,iuse)
if ( .not. iuse ) cycle read_modis
dataaod = rmiss_single
dataaod(4) = aod_550
aeroout( 1,itx) = rsat
aeroout( 2,itx) = tdiff
aeroout( 3,itx) = dlon ! grid relative longitude
aeroout( 4,itx) = dlat ! grid relative latitude
aeroout( 5,itx) = dlon_earth*rad2deg ! earth relative longitude (degrees)
aeroout( 6,itx) = dlat_earth*rad2deg ! earth relative latitude (degrees)
aeroout( 7,itx) = qaod ! total column AOD error flag
aeroout( 8,itx) = solzen ! solar zenith angle
aeroout( 9,itx) = azimuth ! solar azimuth angle
aeroout(10,itx) = styp ! surface type
aeroout(11,itx) = dbcf ! deep blue confidence flag
do i = 1, nchanl
aeroout(i+nreal,itx) = dataaod(i)
end do
end do read_modis
call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
naerodat,itxmax,nread,ndata,aeroout,score_crit)
if ( mype_sub == mype_root ) then
do n = 1, ndata
do i = 1, nchanl
if ( aeroout(i+nreal,n) > rmiss_single ) nodata = nodata + 1
end do
end do
! Write final set of "best" observations to output file
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
write(lunout) ((aeroout(k,n),k=1,naerodat),n=1,ndata)
end if
! Deallocate local arrays
deallocate(aeroout)
deallocate(dataaod)
! End of MODIS bufr block
else ! subset /= NC008041
write(6,*)'READ_AEROSOL: *** WARNING: unknown aerosol data type, subset=',subset
write(6,*)' infile=',infile, ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid
write(6,*)' SKIP PROCESSING OF THIS MODIS FILE'
endif
else ! read subset iret /= 0
write(6,*)'READ_AEROSOL: *** WARNING: read subset error, obstype=',obstype,', iret=',iret
end if
call closbf(lunin)
close(lunin)
else ! obstype /= 'modis'
write(6,*)'READ_AEROSOL: *** WARNING: unknown aerosol input type, obstype=',obstype
endif
! Deallocate satthin arrays
call destroygrids
end subroutine read_aerosol