subroutine read_anowbufr(nread,ndata,nodata,gstime,& infile,obstype,lunout,twindin,sis,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_anowbufr read pm2_5 obs from AIRNow prepbufr file (based on other bufr readers) ! prgmmr: pagowski date: 2010-09-13 ! ! abstract: This routine reads conventional pm2_5 data from AIRNow ! prepbufr can be done similarly for ozone ! ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 2010-09-13 pagowski adopted prepbufr reader code for ! AIRNow bufr for pm2_5 ! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) ! 2013-11-01 pagowski - make code compatible with ncep/mhu airnow bufr ! 2015-02-23 Rancic/Thomas - add l4densvar to time window logical ! 2015-10-01 guo - consolidate use of ob location (in deg) ! ! input argument list: ! infile - unit from which to read BUFR data ! obstype - observation type to process ! lunout - unit to which to write data for further processing ! twindin - input group time window (hours) ! ! output argument list: ! nread - number of type "obstype" observations read ! ndata - number of type "obstype" observations retained for further processing ! nodata - number of individual "obstype" observations retained for !further processing ! sis - satellite/instrument/sensor indicator ! nobs - array of observations on each subdomain for each processor ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_double,i_kind use constants, only: zero,one,half,deg2rad,& rad2deg,r60inv use gridmod, only: diagnostic_reg,regional,nlon,nlat,& tll2xy,txy2ll,rlats,rlons,region_dx use convinfo, only: nconvtype,ctwind, & icuse,ioctype,ictype,cermin,cermax use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen use chemmod, only : obs2model_anowbufr_pm,& iconc,ierror,ilat,ilon,itime,iid,ielev,isite,iikx,ilate,ilone,& elev_missing,site_scale,tunable_error,& code_pm25_ncbufr,code_pm25_anowbufr,& code_pm10_ncbufr,code_pm10_anowbufr use mpimod, only: npe implicit none ! declare passed variables character(len=*),intent(in ) :: infile,obstype integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata,nodata integer(i_kind),dimension(npe) ,intent(inout) :: nobs real(r_kind) ,intent(in ) :: gstime,twindin character(len=*),intent(in ) :: sis ! declare local parameters integer(i_kind), parameter :: maxobs=1e5 integer(i_kind), parameter :: nsid=1,nxob=2,& nyob=3,ndhr=4,ntyp=5,ncopopm=6 !see headr input format below integer(i_kind), parameter :: nfields=6 !output format parameters integer(i_kind), parameter:: nchanl=0,nreal=ilone real(r_kind),parameter :: r360 = 360.0_r_kind real(r_kind),parameter :: percent=1.e-2_r_kind real(r_kind), parameter :: anow_missing=1.0e11_r_kind,& conc_missing = anow_missing-1 ! declare local variables logical outside integer(i_kind) lunin,i integer(i_kind) idate,iret,k integer(i_kind) kx,ikx integer(i_kind) nmind integer(i_kind) :: iy,im,id,ih,imin,site_char,site_id ! (site_char=1,2,3,4: unknown,urban,suburban,rural) integer(i_kind) :: ireadmg,ireadsb integer(i_kind), dimension(5) :: idate5 integer(i_kind), dimension(8) :: idate8 real(r_kind), dimension(5) :: rinc character(len=8) :: subset character(len=80) :: headr real(r_double), dimension(nfields) :: indata real(r_kind) :: tdiff,obstime,t4dv real(r_kind) :: dlat,dlon,error_1,error_2,obserror,dlat_earth,dlon_earth real(r_kind) :: dlat_earth_deg,dlon_earth_deg real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00 integer(i_kind) ntest,ios real(r_kind) :: conc,site_elev real(r_kind), dimension(nreal,maxobs):: cdata_all character(len=8):: sid character(len=10) :: cdate logical :: ncbufr,anowbufr equivalence (sid,indata(nsid)) data lunin / 10 / ncbufr=.false. anowbufr=.false. site_char=1 ! set unknown site character site_elev=elev_missing ! set unknown site elevation !************************************************************************** ! initialize variables disterrmax=zero ntest=0 nread=0 ndata = 0 nodata = 0 ! open, then read date from bufr data open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) ! reading each report from bufr do while (ireadmg(lunin,subset,idate) == 0) if (trim(obstype)=='pm2_5') then if ( (subset == 'NC008031') .or. (subset == 'NC008032' ) ) then headr='PTID CLONH CLATH TPHR TYPO COPOPM' ncbufr=.true. write(6,*)'READ_PM2_5: AIRNOW data type, subset=',subset else if (subset == 'ANOWPM') then headr='SID XOB YOB DHR TYP COPOPM' anowbufr=.true. write(6,*)'READ_PM2_5: AIRNOW data type, subset=',subset else cycle endif else if (trim(obstype)=='pm10') then if (subset == 'NC008033') then headr='PTID CLONH CLATH TPHR TYPO COPOPM' ncbufr=.true. write(6,*)'READ_PM10: AIRNOW data type, subset=',subset else cycle endif else !keep option for ozone cycle ! subset='AIRNOW' for ozone endif write(cdate,'(i10)') idate read (cdate,'(i4,3i2)') iy,im,id,ih imin=0 do while (ireadsb(lunin) == 0) call ufbint(lunin,indata,nfields,1,iret,headr) if (anowbufr) then kx=indata(ntyp) read(sid,'(Z8)')site_id else if (ncbufr) then kx=indata(ntyp) if (kx/=code_pm25_ncbufr .or. kx/=code_pm10_ncbufr) then cycle else if (trim(obstype)=='pm2_5') then kx=code_pm25_anowbufr else kx=code_pm10_anowbufr endif endif read(sid,'(Z8)',iostat=ios)site_id if (ios/=0) site_id=nint(indata(nsid)) endif nread = nread + 1 conc=indata(ncopopm) if ( iret > 0 .and. (conc < conc_missing ) .and. & (conc >= zero)) then if(indata(nxob) >= r360) indata(nxob) = indata(nxob) - r360 if(indata(nxob) < zero) indata(nxob) = indata(nxob) + r360 dlon_earth_deg=indata(nxob) dlat_earth_deg=indata(nyob) dlon_earth=indata(nxob)*deg2rad dlat_earth=indata(nyob)*deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) ! convert to rotated coordinate if(diagnostic_reg) then call txy2ll(dlon,dlat,rlon00,rlat00) ntest=ntest+1 cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* & (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00)) cdist=max(-one,min(cdist,one)) disterr=acos(cdist)*rad2deg disterrmax=max(disterrmax,disterr) end if if(outside) cycle ! check to see if outside regional domain else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! extract date information. if time outside window, skip this obs idate5(1) = iy !year idate5(2) = im !month idate5(3) = id !day idate5(4) = ih !hour idate5(5) = imin !minute idate8(1:3)=idate5(1:3) idate8(4)=0 idate8(5:6)=idate5(4:5) idate8(7:8)=0 rinc(1)=zero ! half below accounts for the fact that in bufr file measurement ! at 5:30 (ie. average from 5:00 to 6:00) is assigned time 6:00 rinc(2)=indata(ndhr)-half rinc(3:5)=zero call w3movdat(rinc,idate8,idate8) idate5(1:3)=idate8(1:3) idate5(4:5)=idate8(5:6) ikx = 0 do i = 1, nconvtype if (obstype == ioctype(i) .and. kx == ictype(i) .and. & abs(icuse(i))== 1) ikx=i end do if(ikx == 0) cycle ! not ob type used call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv obstime=real(nmind,r_kind) tdiff=(obstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dv < zero .or. t4dv > winlen) cycle else if(abs(tdiff) > twindin .or. & abs(tdiff) > ctwind(ikx)) cycle ! outside time window endif !at a later stage search for site character and elevation !based on sid when reference table available !for now assign default sitecharacter and unknown elevation !calculate pm2_5 obs error !obs error for pm2_5/pm10 is calculated as !obserror=sqrt(error_1^2+error_2^2) !measurement error: !error_1=cermax+cermin*conc !(1.5ug/m3+0.0075*concentration) !cermin, cermax in convinfo.txt are nonstandard !cermax is a constant part, cermin is a percentage part !representativeness error: !error_2=tunable_error*error_1*sqrt(dx/site_scale) !similar for ozone conc=conc*obs2model_anowbufr_pm() error_1=cermax(ikx)+cermin(ikx)*percent*conc error_2=tunable_error*error_1*& sqrt(region_dx(nint(dlat),nint(dlon))/& site_scale(site_char)) obserror=sqrt(error_1**2+error_2**2) ndata=ndata+1 nodata=nodata+1 cdata_all(iconc,ndata) = conc ! pm2_5 obs cdata_all(ierror,ndata) = obserror ! pm2_5 obs error cdata_all(ilat,ndata) = dlat ! grid relative latitude cdata_all(ilon,ndata) = dlon ! grid relative longitude cdata_all(itime,ndata) = tdiff ! time of obs cdata_all(iid,ndata) = site_id ! site id cdata_all(ielev,ndata) = site_elev ! elevation cdata_all(isite,ndata) = site_char ! site character cdata_all(iikx,ndata) = ikx ! ordered number in convinfo table cdata_all(ilate,ndata) = dlat_earth_deg ! earth relative latitude (degrees) cdata_all(ilone,ndata) = dlon_earth_deg ! earth relative longitude (degrees) end if enddo enddo 1000 continue call closbf(lunin) close(lunin) if(diagnostic_reg .and. & ntest > 0) write(6,*)'read_airnow_bufr: ',& 'ntest,disterrmax=',ntest,disterrmax if (nodata == 0) then write(6,*)'did not find pm2_5 in airnow_bufr ' write(6,*)'check input airnow_bufr file' endif ! write header record and data to output file for further processing call count_obs(ndata,nreal,ilat,ilon,cdata_all,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((cdata_all(k,i),k=1,nreal),i=1,ndata) return end subroutine read_anowbufr