subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & obstype,twind,sis,ithin,rmesh,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_ozone read ozone data ! prgmmr: yang org: np23 date: 1998-05-15 ! ! abstract: This routine reads SBUV/2 ozone observations. Both layer ! and total column values are read in. The routine has ! the ability to read both IEEE and BUFR format SBUV/2 ! ozone data files. OMI and GOME data is optionally thinned ! to a specific resolution using simple quality control checks. ! ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 1998-05-15 yang, weiyu ! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version ! 2004-06-16 treadon - update documentation ! 2004-07-29 treadon - add only to module use, add intent in/out ! 2004-09-17 todling - fixed intent of jsatid ! 2004-12-02 todling - compilation in OSF1 forces big_endian for bufr files; ! need to force little_endian for ieee files ! 2004-12-22 kokron - change cpp tokens to add support for ifort compiler ! efc does not have a convert option so it should use ! the other 'open' ! 2005-03-14 treadon - define numeric constants to r_kind precision ! 2005-05-12 wu - add OMI total ozone ! 2005-06-27 guo - bug fix: hour read from header was incorrect ! 2005-09-08 derber - modify to use input group time window ! 2005-09-19 treadon - add check on NOAA-17 sbuv data (toss bad data) ! 2005-10-17 treadon - add grid and earth relative obs location to output file ! 2005-10-18 treadon - remove array obs_load and call to sumload ! 2005-12-23 treadon - bound longitude to be less than 360.0 ! 2006-01-26 treadon - remove ieee sbuv option ! 2006-02-03 derber - modify for new obs control and obs count ! 2007-03-01 tremolet - measure time from beginning of assimilation window ! 2007-07-10 zhou - modify to read version 8 SBUV/2 BUFR data(keep ! option to read version 6 data), also add ! total ozone and ozone profile quality control. ! 2007-09-11 h.liu - add kidsat for nimbus-7, n09, n11, n14 ! 2007-10-16 zhou - organize ozone flag control for all satellites ! 2008-04-16 h.liu - thin OMI and read in GOME data ! 2008-05-27 safford - rm unused vars and uses ! 2008-05-30 treadon - accept version8 poq=7 obs for further processing ! 2008-06-01 treadon - adjust logic to correctly handle zero length BUFR files ! 2008-06-03 treadon - add use_poq7 flag ! 2008-09-08 lueken - merged ed's changes into q1fy09 code ! 2009-01-20 sienkiewicz - merge in changes for MLS ozone ! 2009-04-21 derber - add ithin to call to makegrids ! 2009-3-05 h.liu - read in OMI bufr, QC GOME2 and OMI ! 2009-7-02 h.liu - toss the OMI data with AFBO=3 (c-pair correction) and clean up codes ! 2010-05-26 treadon - add timedif=zero for l4dvar (used in thinning) ! 2010-06-02 sienkiewicz - care for closing bufr other than for o3lev ! 2011-07-04 todling - fixes to run either single or double precision ! 2011-08-01 lueken - replaced F90 with f90 (no machine logic) ! 2012-10-12 h.liu - read in MLS v2 Near Real Time (NRT) and v2.2 standard bufr data ! 2013-01-17 h.liu - read in MLS v3 Near Real Time (NRT) ! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) ! 2014-02-03 guo - removed unused "o3lev" handling, which can (and should) be ! implemented again in module m_extOzone, if ever needed. ! 2015-02-23 Rancic/Thomas - add thin4d to time window logical ! 2015-10-01 guo - consolidate use of ob location (in deg ! 2017-12-05 wargan - implement OMPS nadir capability ! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90. ! 2018-06-26 todling - total column nadir OMPS of Wargan and Liu handling redundant NM/TC8 names ! 2018-08-13 H. Liu - add capability to use OMPS nadir profiler and nadir mapper data ! ! input argument list: ! obstype - observation type to process ! jsatid - satellite id to read ! infile - unit from which to read ozone 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) ! ! output argument list: ! nread - number of sbuv/omi ozone observations read ! ndata - number of sbuv/omi ozone profiles retained for further processing ! nodata - number of sbuv/omi ozone observations retained for further processing ! nobs - array of observations on each subdomain for each processor ! ! remarks: ! NCEP stopped producing IEEE format sbuv ozone files in April 2004. ! Hence, the IEEE portion of this routine no future application. It ! is retained in the GSI package for use with retrospective runs. The ! IEEE portion of this routine may be removed from the GSI at a later date. ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_double,i_kind use satthin, only: makegrids,map2tgrid,destroygrids, & finalcheck,itxmax use satthin, only: radthin_time_info,tdiff2crit use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons use constants, only: deg2rad,zero,one_tenth,r60inv,two use obsmod, only: nloz_v6,nloz_v8, ompslp_mult_fact use obsmod, only: time_window_max use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen use radinfo, only: dec2bin use qcmod, only: use_poq7 use ozinfo, only: jpch_oz,nusis_oz,iuse_oz use mpimod, only: npe implicit none ! Declare passed variables character(len=*),intent(in ) :: obstype,infile,jsatid character(len=20),intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout,ithin integer(i_kind) ,intent(inout) :: nread integer(i_kind),dimension(npe) ,intent(inout) :: nobs integer(i_kind) ,intent(inout) :: ndata,nodata real(r_kind) ,intent(in ) :: gstime,twind,rmesh ! Declare local parameters real(r_kind),parameter:: r6 = 6.0_r_kind real(r_kind),parameter:: r76 = 76.0_r_kind real(r_kind),parameter:: r84 = 84.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind real(r_kind),parameter:: rmiss = -9999.9_r_kind real(r_kind),parameter:: badoz = 10000.0_r_kind real(r_kind),parameter:: montoz = 100.0_r_kind !monitored ozone ! Declare local variables logical outside,version6,version8,iuse character(2) version character(8) subset,subset6,subset8,subset8_ompsnp character(49) ozstr,ozostr character(63) lozstr character(51) ozgstr character(27) ozgstr2 character(42) ozostr2 character(64) mlstr character(14) mlstrl integer(i_kind) maxobs,nozdat,nloz integer(i_kind) idate,jdate,ksatid,kk,iy,iret,im,ihh,idd,lunin integer(i_kind) nmind,i,j integer(i_kind) nmrecs,k,ilat,ilon,nreal,nchanl ! integer(i_kind) ithin,kidsat integer(i_kind) kidsat integer(i_kind) idate5(5) integer(i_kind) JULIAN,IDAYYR,IDAYWK integer(i_kind) ikx integer(i_kind) decimal,binary(14),binary_mls(18) integer(i_kind) itx,itt,ipoq7 real(r_kind) tdiff,sstime,dlon,dlat,t4dv,crit1,dist1 real(r_kind) slons0,slats0,rsat,solzen,solzenp,dlat_earth,dlon_earth real(r_kind) dlat_earth_deg,dlon_earth_deg real(r_kind),allocatable,dimension(:):: poz ! maximum number of observations set to real(r_kind),allocatable,dimension(:,:):: ozout real(r_double) toq,poq real(r_double),dimension(nloz_v6):: ozone_v6 real(r_double),dimension(29,nloz_v8):: ozone_v8 real(r_double),dimension(10):: hdroz real(r_double),dimension(10):: hdrozg real(r_double),dimension(5):: hdrozg2 real(r_double),dimension(10):: hdrozo real(r_double),dimension(8) :: hdrozo2 real(r_double),dimension(13):: hdrmls real(r_double),allocatable,dimension(:,:) :: hdrmlsl real(r_kind),allocatable,dimension(:):: mlspres,mlsoz,mlsozpc,usage1 integer(i_kind),allocatable,dimension(:):: ipos real(r_double) totoz,hdrmls13 integer(i_kind) :: k0 logical :: first,read_success real(r_double),allocatable,dimension(:,:):: olpdtsq,lpsdvals real(r_double),allocatable,dimension(:):: press,omr,omrstd real(r_double) said, lat, lon, year, month, day, hour, minu real(r_double) soza ! MLS data version: mlsv=22 is version 2.2 standard data; ! mlsv=20 is v2 near-real-time data ! mlsv=30 is v3 near-real-time data integer(i_kind) :: mlsv data lozstr & / 'OSP12 OSP11 OSP10 OSP9 OSP8 OSP7 OSP6 OSP5 OSP4 OSP3 OSP2 OSP1 ' / data ozgstr & / 'SAID CLAT CLON YEAR DOYR HOUR MINU SECO SOZA SOLAZI' / data ozgstr2 & / 'CLDMNT SNOC ACIDX STKO FOVN' / data ozostr & / 'SAID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SOZA' / ! since 2009020412, the omi bufr contains fovn data ozostr2 & / 'CLDMNT ACIDX STKO VZAN TOQC TOQF FOVN AFBO' / data mlstr & / 'SAID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SOZA CONV MLST PCCF' / data mlstrl & / 'PRLC OZMX OZMP' / data lunin / 10 / data subset6 / 'NC008010' / data subset8 / 'NC008011' / data subset8_ompsnp / 'NC008017'/ real(r_kind) :: ptime,timeinflat,crit0 integer(i_kind) :: ithin_time,n_tbin,it_mesh !************************************************************************** ! Set constants. Initialize variables rsat=999._r_kind maxobs=1e6 ilon=3 ilat=4 ipoq7=0 if (use_poq7) ipoq7=7 ! Separately process sbuv or omi ozone if (obstype == 'sbuv2' .or. obstype == 'ompsnp') then nreal=9 open(lunin,file=trim(infile),form='unformatted') nmrecs=0 call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) version6 = .false. version8 = .false. if (subset == subset6) then version6 = .true. nloz = nloz_v6 version = 'v6' elseif (subset == subset8 .or. subset == subset8_ompsnp) then ! OMPS-NP processed with V8 algorithm version8 = .true. nloz = nloz_v8 version = 'v8' else write(6,*)'READ_OZONE: *** WARNING: unknown layer ozone version type, subset=',subset write(6,*)' infile=',trim(infile), ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid write(6,*)' SKIP PROCESSING OF THIS OZONE LAYER FILE' call closbf(lunin) close(lunin) return endif ! Set dependent variables and allocate arrays nchanl=nloz+1 nozdat=nreal+nchanl allocate (ozout(nozdat,maxobs)) allocate ( poz(nloz+1)) ! Set BUFR string based on sbuv version if (version6) then ozstr='SAID CLAT CLON YEAR MNTH DAYS HOUR MINU OSZA OPSZ' else if (version8) then ozstr='SAID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SOZA' endif read_loop1: do call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) if (iret/=0) exit read_loop1 cycle read_loop1 endif ! extract header information ! BUFR code values for satellite identifiers are listed in ! Dennis Keyser's website, ! http://www.emc.ncep.noaa.gov/mmb/papers/keyser/Satellite_Historical.txt call ufbint(lunin,hdroz,10,1,iret,ozstr) rsat = hdroz(1); ksatid=rsat if(jsatid == 'nim07') kidsat = 767 if(jsatid == 'n09') kidsat = 201 if(jsatid == 'n11') kidsat = 203 if(jsatid == 'n14') kidsat = 205 if(jsatid == 'n16') kidsat = 207 if(jsatid == 'n17') kidsat = 208 if(jsatid == 'n18') kidsat = 209 if(jsatid == 'n19') kidsat = 223 if(jsatid == 'npp') kidsat = 224 if(jsatid == 'n20') kidsat = 225 if(jsatid == 'n21') kidsat = 226 if (ksatid /= kidsat) cycle read_loop1 nmrecs=nmrecs+nloz+1 ! Convert observation location to radians slats0= hdroz(2) slons0= hdroz(3) if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle read_loop1 if(slons0< zero) slons0=slons0+r360 if(slons0==r360) slons0=zero dlat_earth_deg = slats0 dlon_earth_deg = slons0 dlat_earth = slats0 * deg2rad dlon_earth = slons0 * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(outside) cycle read_loop1 else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! Special check for NOAA-17 version 6 ! Before July 2007 NOAA-17 SBUV/2 has a stray light problem which produces ! erroneous ozone profile retrievals for a limited portion ! of its measurements. The contaminated signals only occur ! in the Southern Hemisphere and only for Solar Zenith ! Angles (SZA) greater than 76 Degrees. if (version6) then solzen = hdroz(9) ! solar zenith angle solzenp= hdroz(10) ! profile solar zenith angle if (ksatid==208 .and. dlat_earth r76) cycle read_loop1 else if(version8)then solzen = hdroz(10) ! solar zenith angle endif ! Convert observation time to relative time idate5(1) = hdroz(4) !year idate5(2) = hdroz(5) !month idate5(3) = hdroz(6) !day idate5(4) = hdroz(7) !hour idate5(5) = hdroz(8) !minute call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv sstime=real(nmind,r_kind) tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if(t4dvwinlen) cycle read_loop1 else if(abs(tdiff) > twind) cycle read_loop1 end if ! Extract layer ozone values and compute profile total ozone if (version8) then call ufbseq(lunin,ozone_v8,29,21,iret,'OZOPQLSQ') totoz=zero do k=1,nloz kk=nloz-k+1 poz(k) = ozone_v8(6,kk) totoz=totoz+ozone_v8(6,k) end do poz(nloz+1) = totoz endif if (version6) then call ufbint(lunin,ozone_v6,nloz,1,iret,lozstr) do k=1,nloz kk=nloz-k+1 poz(k) = ozone_v6(kk) end do ! extract total ozone call ufbint(lunin,totoz,1,1,iret,'OTSP') poz(nloz+1) = totoz endif ! Extract and apply version 8 total and profile ozone quaility information ! Toss observations for which the total ozone error code is neither 0 nor 2 ! Toss observations for which the profile ozone error code is neither 0 nor 1 ! NOTES: ! 1) Profile ozone error code 0 identifies good data; 1 identifies good ! data with a solar zenith angle > 84 degrees; 7 identifies profile ! for which stray light correction applied ! 2) Total ozone error code 0 indentifies good data; 2 identifies good ! data with a solar zenith angle > 84 degrees. ! 3) We do not use the version 6 error flags. Thus, initialize toq and ! poq to 0 (use the data) toq=0._r_double poq=0._r_double if (version8) then call ufbint(lunin,toq,1,1,iret,'SBUVTOQ') call ufbint(lunin,poq,1,1,iret,'SBUVPOQ') if (toq/=0 .and. toq/=2) cycle read_loop1 if (poq/=0 .and. poq/=1 .and. poq/=ipoq7) cycle read_loop1 endif ! Check ozone layer values. If any layer value is bad, toss entire profile do k=1,nloz if (poz(k)>badoz) cycle read_loop1 end do ! Write ozone record to output file ndata=min(ndata+1,maxobs) nodata=nodata+nloz+1 ozout(1,ndata)=rsat ozout(2,ndata)=t4dv ozout(3,ndata)=dlon ! grid relative longitude ozout(4,ndata)=dlat ! grid relative latitude ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) ozout(7,ndata)=toq ! total ozone error flag ozout(8,ndata)=poq ! profile ozone error flag ozout(9,ndata)=solzen ! solar zenith angle do k=1,nloz+1 ozout(k+9,ndata)=poz(k) end do ! Loop back to read next profile end do read_loop1 call closbf(lunin) close(lunin) ! End of bufr ozone block ! Process GOME-2 data else if ( obstype == 'gome') then call radthin_time_info(obstype, jsatid, sis, ptime, ithin_time) if( ptime > 0.0_r_kind) then n_tbin=nint(2*time_window_max/ptime) else n_tbin=1 endif open(lunin,file=trim(infile),form='unformatted') nmrecs=0 call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) if (subset == 'NC008012') then write(6,*)'READ_OZONE: GOME-2 data type, subset=',subset else write(6,*)'READ_OZONE: *** WARNING: unknown ozone data type, subset=',subset write(6,*)' infile=',trim(infile), ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid call closbf(lunin) close(lunin) return endif ! Make thinning grids call makegrids(rmesh,ithin,n_tbin=n_tbin) ! Set dependent variables and allocate arrays nreal=14 nloz=0 nchanl=1 nozdat=nreal+nchanl allocate (ozout(nozdat,itxmax)) do k=1,itxmax do i=1,nozdat ozout(i,k)=rmiss end do end do iy=0 idd=0 ihh=0 obsloop: do call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) if (iret/=0) exit obsloop cycle obsloop endif ! extract header information call ufbint(lunin,hdrozg,10,1,iret,ozgstr) call ufbint(lunin,hdrozg2,5,1,iret,ozgstr2) rsat = hdrozg(1); ksatid=rsat if(jsatid == 'metop-a')kidsat = 4 if(jsatid == 'metop-b')kidsat = 3 if(jsatid == 'metop-c')kidsat = 5 if (ksatid /= kidsat) cycle obsloop ! NESDIS does not put a flag for high SZA gome-2 data (SZA > 84 degree) if ( hdrozg(9) > r84 ) cycle obsloop nmrecs=nmrecs+nloz+1 ! Convert observation location to radians slats0= hdrozg(2) slons0= hdrozg(3) if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle obsloop if(slons0< zero) slons0=slons0+r360 if(slons0==r360) slons0=zero dlat_earth_deg = slats0 dlon_earth_deg = slons0 dlat_earth = slats0 * deg2rad dlon_earth = slons0 * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(outside) cycle obsloop else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! Convert observation time to relative time idate5(1) = hdrozg(4) !year IDAYYR = hdrozg(5) ! Day of year JULIAN = -31739 + 1461 * (idate5(1) + 4799) /4 & -3 * ((idate5(1) + 4899) / 100) / 4 + IDAYYR call w3fs26(JULIAN,idate5(1),idate5(2),idate5(3),IDAYWK,IDAYYR) ! idate5(2) month ! idate5(3) day idate5(4) = hdrozg(6) !hour idate5(5) = hdrozg(7) !minute call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv sstime=real(nmind,r_kind) tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if(t4dvwinlen) cycle obsloop else if(abs(tdiff) > twind) cycle obsloop end if ! extract total ozone call ufbint(lunin,totoz,1,1,iret,'OZON') if (totoz > badoz ) cycle obsloop ! only accept flag 0 (good) data toq=0._r_double call ufbint(lunin,toq,1,1,iret,'GOMEEF') if (toq/=0) cycle obsloop ! only accept scan positions from 2 to 25 if( hdrozg2(5) < two .or. hdrozg2(5) > 25.0_r_kind ) cycle obsloop ! thin GOME data ! GOME data has bias when the satellite looks to the east. Consider QC out this data. crit0 = 0.01_r_kind timeinflat=r6 call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh) call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh) if(.not. iuse) cycle obsloop call finalcheck(dist1,crit1,itx,iuse) if(.not. iuse) cycle obsloop ndata=ndata+1 nodata=ndata ozout(1,itx)=rsat ozout(2,itx)=t4dv ozout(3,itx)=dlon ! grid relative longitude ozout(4,itx)=dlat ! grid relative latitude ozout(5,itx)=dlon_earth_deg ! earth relative longitude (degrees) ozout(6,itx)=dlat_earth_deg ! earth relative latitude (degrees) ozout(7,itx)=toq ! total ozone error flag ozout(8,itx)=hdrozg(9) ! solar zenith angle ozout(9,itx)=hdrozg(10) ! solar azimuth angle ozout(10,itx)=hdrozg2(1) ! CLOUD AMOUNT IN SEGMENT ozout(11,itx)=hdrozg2(2) ! SNOW COVER ozout(12,itx)=hdrozg2(3) ! AEROSOL CONTAMINATION INDEX ozout(13,itx)=hdrozg2(4) ! ASCENDING/DESCENDING ORBIT QUALIFIER ozout(14,itx)=hdrozg2(5) ! scan position (fovn) ozout(15,itx)=totoz end do obsloop call closbf(lunin) close(lunin) ! End of GOME bufr block ! Process OMI/OMPS data without efficiency factors else if ( obstype == 'omi' .or. obstype == 'ompsnm' .or. obstype == 'ompstc8') then call radthin_time_info(obstype, jsatid, sis, ptime, ithin_time) if( ptime > 0.0_r_kind) then n_tbin=nint(2*time_window_max/ptime) else n_tbin=1 endif nmrecs=0 open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) select case(subset) case('NC008013') write(6,*)'READ_OZONE: OMI data type, subset=',subset case('NC008018') write(6,*)'READ_OZONE: OMPS Nadir Mapper data type, subset=',subset case default write(6,*)'READ_OZONE: *** WARNING: unknown ozone data type, subset=',subset write(6,*)' infile=',trim(infile), ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid call closbf(lunin) close(lunin) return end select ! Make thinning grids call makegrids(rmesh,ithin,n_tbin=n_tbin) ! Set dependent variables and allocate arrays nreal=14 nloz=0 nchanl=1 nozdat=nreal+nchanl allocate (ozout(nozdat,itxmax)) do k=1,itxmax do i=1,nozdat ozout(i,k)=rmiss end do end do iy=0 im=0 idd=0 ihh=0 read_loop2: do call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) if (iret/=0) exit read_loop2 cycle read_loop2 endif ! extract header information call ufbint(lunin,hdrozo,10,1,iret,ozostr) call ufbint(lunin,hdrozo2,8,1,iret,ozostr2) rsat = hdrozo(1); ksatid=rsat if(jsatid == 'aura')kidsat = 785 if(jsatid == 'npp') kidsat = 224 if(jsatid == 'n20') kidsat = 225 if(jsatid == 'n21') kidsat = 226 if (ksatid /= kidsat) cycle read_loop2 nmrecs=nmrecs+nloz+1 ! Convert observation location to radians slats0= hdrozo(2) slons0= hdrozo(3) if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle read_loop2 if(slons0< zero) slons0=slons0+r360 if(slons0==r360) slons0=zero dlat_earth_deg = slats0 dlon_earth_deg = slons0 dlat_earth = slats0 * deg2rad dlon_earth = slons0 * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(outside) cycle read_loop2 else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! convert observation time to relative time idate5(1) = hdrozo(4) !year idate5(2) = hdrozo(5) !month idate5(3) = hdrozo(6) !day idate5(4) = hdrozo(7) !hour idate5(5) = hdrozo(8) !minute call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv sstime=real(nmind,r_kind) tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle read_loop2 else if(abs(tdiff) > twind) cycle read_loop2 end if ! extract total ozone call ufbint(lunin,totoz,1,1,iret,'OZON') if (totoz > badoz ) cycle read_loop2 ! QC for omi_aura if (obstype == 'omi') then ! Bit 10 in TOQF represents row anomaly. decimal=int(hdrozo2(6)) call dec2bin(decimal,binary,14) if (binary(10) == 1 ) cycle read_loop2 ! remove the bad scan position data: fovn beyond 25 if (hdrozo2(7) >=25.0_r_double) cycle read_loop2 end if ! only accept flag 0 1, flag 2 is high SZA data which is not used for now toq=hdrozo2(5) if (toq/=0 .and. toq/=1) cycle read_loop2 ! remove the data in which the C-pair algorithm ((331 and 360 nm) is used. if (hdrozo2(8) == 3_r_double .or. hdrozo2(8) == 13_r_double) cycle read_loop2 ! thin OMI/OMPS-NM(or TC8) data crit0 = 0.01_r_kind timeinflat=r6 call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh) call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh) if(.not. iuse)cycle read_loop2 call finalcheck(dist1,crit1,itx,iuse) if(.not. iuse)cycle read_loop2 ndata=ndata+1 nodata=ndata ozout(1,itx)=rsat ozout(2,itx)=t4dv ozout(3,itx)=dlon ! grid relative longitude ozout(4,itx)=dlat ! grid relative latitude ozout(5,itx)=dlon_earth_deg ! earth relative longitude (degrees) ozout(6,itx)=dlat_earth_deg ! earth relative latitude (degrees) ozout(7,itx)=hdrozo2(5) ! total ozone quality code ozout(8,itx)=hdrozo(10) ! solar zenith angle if (obstype == 'omi') then ozout(9,itx)=binary(10) ! row anomaly flag end if ozout(10,itx)=hdrozo2(1) ! cloud amount ozout(11,itx)=hdrozo2(4) ! vzan ozout(12,itx)=hdrozo2(2) ! aerosol index ozout(13,itx)=hdrozo2(3) ! ascending/descending ozout(14,itx)=hdrozo2(7) ! scan position ozout(15,itx)=totoz ! End of loop over observations end do read_loop2 call closbf(lunin) close(lunin) ! End of OMI/OMPS-NM(or TC8) block ! Process MLS bufr data else if ( index(obstype,'mls')/=0 ) then nmrecs=0 open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) if (subset == 'NC008015') then write(6,*)'READ_OZONE: MLS data type, subset=',subset else write(6,*)'READ_OZONE: *** WARNING: unknown ozone data type, subset=',subset write(6,*)' infile=',trim(infile), ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid call closbf(lunin) close(lunin) return endif read_success=.false. if(iret==0) then call readsb(lunin,iret) if (iret/=0) then do call readmg(lunin,subset,jdate,iret) if (iret/=0) exit call readsb(lunin,iret) if(iret == 0)then read_success=.true. exit end if end do endif endif if(read_success)then ! Get # of vertical pressure levels nloz and MLS NRT data version which depends on nloz allocate(hdrmlsl(3,100)) call ufbrep(lunin,hdrmlsl,3,100,iret,mlstrl) nloz=iret ! for NRT data, mlsv=20 or 30 depending on the nloz mlsv=-999 if(nloz==37) then if(index(sis,'mls22')/=0 ) then !mls v2.2 data mlsv=22 else if(index(sis,'mls20')/=0 ) then !mls v2 nrt data mlsv=20 end if else if (nloz==55) then !mls v3 nrt data if (index(sis,'mls30')/=0 ) then mlsv=30 endif else write(6,*) 'invalid vertical level number: ', nloz write(6,*) '******STOP*******: error reading MLS vertical levels in read_ozone.f90' call stop2(338) end if deallocate(hdrmlsl) write(6,*) 'READ_OZONE: MLS data version=',mlsv write(6,*) 'READ_OZONE: MLS vertical level number=',nloz if (mlsv<0) then write(6,*) 'inconsistent MLS versions. bufr nloz=',nloz,' obsinput sis= ',trim(sis) write(6,*) '******STOP*******: error bufr and specified MLS versions' call stop2(338) end if ! Allocate arrays allocate(hdrmlsl(3,nloz)) allocate (mlspres(nloz)) allocate (mlsoz(nloz)) allocate (mlsozpc(nloz)) allocate(ipos(nloz)) allocate (usage1(nloz)) ! Set dependent variables and allocate arrays nreal=12 nchanl=1 nozdat=nreal+nchanl allocate (ozout(nozdat,maxobs)) do k=1,maxobs do i=1,nozdat ozout(i,k)=rmiss end do end do ikx=0 k0=0 ipos=999 first=.false. do k=1,jpch_oz if( (.not. first) .and. index(nusis_oz(k),sis)/=0 ) then k0=k first=.true. end if if(first .and. index(nusis_oz(k),sis)/=0 ) then ikx=ikx+1 ipos(ikx)=k0+ikx-1 end if end do ! Reopen unit to bufr file call closbf(lunin) close(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) end if read_loop4: do if(.not. read_success) exit call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) if (iret/=0) exit read_loop4 cycle read_loop4 endif do k=1,nloz if (iuse_oz(ipos(k)) < 0) then usage1(k) = montoz else usage1(k) = zero endif end do ! extract header information call ufbint(lunin,hdrmls,13,1,iret,mlstr) rsat = hdrmls(1); ksatid=rsat if(jsatid == 'aura')kidsat = 785 if (ksatid /= kidsat) cycle read_loop4 nmrecs=nmrecs+nloz ! Convert observation location to radians slats0= hdrmls(2) slons0= hdrmls(3) if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle read_loop4 if(slons0< zero) slons0=slons0+r360 if(slons0==r360) slons0=zero dlat_earth_deg = slats0 dlon_earth_deg = slons0 dlat_earth = slats0 * deg2rad dlon_earth = slons0 * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(outside) cycle read_loop4 else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! convert observation time to relative time idate5(1) = hdrmls(4) !year idate5(2) = hdrmls(5) !month idate5(3) = hdrmls(6) !day idate5(4) = hdrmls(7) !hour idate5(5) = hdrmls(8) !minute call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv sstime=real(nmind,r_kind) tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle read_loop4 else if(abs(tdiff) > twind) cycle read_loop4 end if ! v2.2 data screening, only accept: ! Pressure range(PRLC): 215-0.02mb (lev5-27) ! Precision(OZMP): positive OZMP; ! Status flag(MLST): only use even number ! Quality(PCCF): use >1.2 for data at 215-100mb & low latitude, ! use >0.4 for data elsewhere ! Convergence(CONV): use <1.8 ! v2 NRT data screening, only accept: ! Pressure range(PRLC): 68-0.2mb (lev8-23) ! Precision(OZMP): positive OZMP; ! Status flag(MLST): only use even number ! Quality(PCCF): do NOT use <1.2 or >3.0 ! v3 NRT data screening, only accept: ! Pressure range(PRLC): 261-0.1mb (lev8-43) ! Precision(OZMP): positive OZMP; ! Status flag(MLST): only use even number ! Quality(PCCF): only use if >0.4 ! Convergence(CONV): use <1.2 ! status: Bit 1 in MLST represents data should not be used ! Note: in BUFR bits are defined from left to right as: 123456789... ! whereas in HDF5 (and the nasa document) bits are defined from right to left as: ...876543210 decimal=int(hdrmls(12)) call dec2bin(decimal,binary_mls,18) if (binary_mls(1) == 1 ) cycle read_loop4 ! v2.2 data, remove data when convergence>1.8 ! v3 NRT data,remove data when convergence>1.2 if(mlsv==22) then if(hdrmls(11) >= 1.8_r_kind) cycle read_loop4 else if(mlsv==30) then if(hdrmls(11) >= 1.2_r_kind) cycle read_loop4 end if ! extract pressure, ozone mixing ratio and precision call ufbrep(lunin,hdrmlsl,3,nloz,iret,mlstrl) do k=1,nloz mlspres(k)=log(hdrmlsl(1,k)*0.001_r_kind) ! mls pressure in Pa, coverted to log(cb) mlsoz(k)=hdrmlsl(2,k) ! ozone mixing ratio in ppmv mlsozpc(k)=hdrmlsl(3,k) ! ozone mixing ratio precision in ppmv ! there is possibility that mlsoz in bufr is 0 or negative or larger than 100 which are not reasonable values. if(mlsoz(k)<1.0e-8_r_kind .or. mlsoz(k)>100.0_r_kind ) then usage1(k)=badoz ! for v2.2 data, if this unreasonable value happens between 215mb (lev5) and 0.02mb (lev27), throw the whole profile ! for v2 NRT data, if this unreasonable value happens between 68mb (lev8) and 0.2mb (lev23), throw the whole profile ! for v3 NRT data, if this unreasonable value happens between 261mb (lev8) and 0.1mb (lev43), throw the whole profile if(mlsv==22 .and. (k<=27 .and. k>=5)) cycle read_loop4 if(mlsv==20 .and. (k<=23 .and. k>=8)) cycle read_loop4 if(mlsv==30 .and. (k<=43 .and. k>=8)) cycle read_loop4 end if end do do k=1,nloz ! pressure range if(mlsv==22) then if(hdrmlsl(1,k)>21700._r_kind .or. hdrmlsl(1,k)<1._r_kind) usage1(k)=badoz else if(mlsv==20) then if(hdrmlsl(1,k)>6900._r_kind .or. hdrmlsl(1,k)<10._r_kind) usage1(k)=badoz else if(mlsv==30) then if(hdrmlsl(1,k)>26500._r_kind .or. hdrmlsl(1,k)<10._r_kind) usage1(k)=badoz end if ! only positive precision accepted if(hdrmlsl(3,k)<=0._r_kind) usage1(k)=badoz end do ! status screening hdrmls13=hdrmls(13)*0.1_r_kind if(mlsv==22) then if (abs(slats0)<30._r_kind) then do k=1,nloz if(hdrmlsl(1,k)>10100._r_kind .and. hdrmlsl(1,k)<21700._r_kind) then if(hdrmls13 <= 1.2_r_kind) usage1(k)=badoz else if(hdrmls13 <= 0.4_r_kind) usage1(k)=badoz endif end do else if(hdrmls13 <= 0.4_r_kind) then do k=1,nloz usage1(k)=badoz end do end if end if else if(mlsv==20) then if(hdrmls13 <= 1.2_r_kind .or. hdrmls13 >= 3.0_r_kind) then do k=1,nloz usage1(k)=badoz end do end if else if(mlsv==30) then if(hdrmls13 <= 0.4_r_kind) then do k=1,nloz usage1(k)=badoz end do end if end if do k=1,nloz ndata=min(ndata+1,maxobs) nodata=ndata ! if(ndata >= nloz) cycle read_loop4 ozout(1,ndata)=rsat ozout(2,ndata)=t4dv ozout(3,ndata)=dlon ! grid relative longitude ozout(4,ndata)=dlat ! grid relative latitude ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) ozout(7,ndata)=hdrmls(10) ! solar zenith angle ozout(8,ndata)=usage1(k) ! ozout(9,ndata)=mlspres(k) ! mls pressure in log(cb) ozout(10,ndata)=mlsozpc(k) ! ozone mixing ratio precision in ppmv ozout(11,ndata)=float(ipos(k)) ! pointer of obs level index in ozinfo.txt ozout(12,ndata)=nloz ! # of mls vertical levels ozout(nreal+1,ndata)=mlsoz(k) ! ozone mixing ratio in ppmv end do end do read_loop4 call closbf(lunin) close(lunin) ! End of MLS bufr loop !Process OMPS LP data elseif(index(obstype,'ompslp') /= 0 )then nloz = 81 nreal=15 nchanl=1 nozdat=nreal+nchanl read_success=.false. open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) if (iret == 0 .and. subset == 'NC008019') then read_success=.true. else write(6,*)'READ_OZONE: *** WARNING: unknown ozone data type, subset=',subset write(6,*)' infile=',trim(infile), ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid call closbf(lunin) close(lunin) return endif ! Allocate arrays allocate(olpdtsq(12,81)) allocate(lpsdvals(3,81)) allocate(press(nloz)) allocate(omr(nloz)) allocate(omrstd(nloz)) allocate(usage1(nloz)) allocate(ipos(nloz)) allocate(ozout(nozdat,maxobs)) do k=1,maxobs do i=1,nozdat ozout(i,k)=rmiss enddo enddo ikx=0 k0=0 ipos=999 first=.false. do k=1,jpch_oz if( (.not. first) .and. index(nusis_oz(k),sis)/=0 ) then k0=k first=.true. end if if(first .and. index(nusis_oz(k),sis)/=0 ) then ikx=ikx+1 ipos(ikx)=k0+ikx-1 endif enddo read_loop5: do if(.not. read_success) exit call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) if (iret/=0) exit read_loop5 cycle read_loop5 endif do k=1,nloz if (iuse_oz(ipos(k)) < 0) then usage1(k) = montoz else usage1(k) = zero endif end do call ufbint(lunin,said,1,1,iret,"SAID") !Convert observation location to radians call ufbint(lunin,lat,1,1,iret,"CLATH") call ufbint(lunin,lon,1,1,iret,"CLONH") if(abs(lat)>90._r_kind .or. abs(lon)>r360) cycle read_loop5 if(lon< zero) lon=lon+r360 if(lon==r360) lon=zero dlat_earth_deg = lat dlon_earth_deg = lon dlat_earth = lat * deg2rad dlon_earth = lon * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(outside) cycle read_loop5 else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif !Convert observation time to relative time call ufbint(lunin,year,1,1,iret,"YEAR") call ufbint(lunin,month,1,1,iret,"MNTH") call ufbint(lunin,day,1,1,iret,"DAYS") call ufbint(lunin,hour,1,1,iret,"HOUR") call ufbint(lunin,minu,1,1,iret,"MINU") idate5(1) = year idate5(2) = month idate5(3) = day idate5(4) = hour idate5(5) = minu call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv sstime=real(nmind,r_kind) tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle read_loop5 else if(abs(tdiff) > twind) cycle read_loop5 endif !Read solar zenith angle call ufbint(lunin,soza,1,1,iret,"SOZA") !Read Pressure and Ozone Mixing Ratio call ufbseq(lunin, olpdtsq, 12, 81, iret, "OLPDTSQ") !Read Ozone Mixing Ratio Standard Deviation call ufbseq(lunin, lpsdvals,3,81,iret,"LPSDVALS") j = 0 do k = 1, nloz press(k) = olpdtsq(2,k)*0.001_r_double ! centibars omr(k) = olpdtsq(11,k) ! ppmv omrstd(k) = lpsdvals(3,k) !omr std j = j + 1 if(omr(k) < 0._r_double .or. omr(k) > 100._r_double) then usage1(k) = badoz j = j - 1 endif enddo do k=1,nloz if(omr(k) > 0._r_double .and. omr(k) < 100._r_double)then ndata=min(ndata+1,maxobs) nodata=ndata nmrecs= ndata ozout(1,ndata)=said ozout(2,ndata)=t4dv ozout(3,ndata)=dlon ! grid relative longitude ozout(4,ndata)=dlat ! grid relative latitude ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) ozout(7,ndata)=soza ! solar zenith angle ozout(8,ndata)=usage1(k) ! ozout(9,ndata)=log(press(k)) ! ompslp pressure in log(cb) ozout(10,ndata)=omrstd(k)*ompslp_mult_fact ! ozone mixing ratio precision in ppmv ozout(11,ndata)=float(ipos(k)) ! pointer of obs level index in ! ozinfo.txt ozout(12,ndata)=j !nloz ! # of ompslp vertical levels ozout(13,ndata)=omr(k) ! ozone mixing ratio in ppmv ozout(14,ndata)=olpdtsq(3,k) ! log10 numberdensity of air ozout(15,ndata)=olpdtsq(6,k) !log10 number density of ozone from UV ozout(16,ndata)=olpdtsq(8,k) !log10 number density of ozone from VIS endif enddo enddo read_loop5 call closbf(lunin) close(lunin) ! end of OMPS LP bufr loop endif if(nmrecs > 0)then ! If gome, omps-nm/tc8 or omi data, compress ozout array to thinned data if (obstype=='omi' .or. obstype=='gome' .or. obstype=='ompsnm' .or. obstype == 'ompstc8') then kk=0 do k=1,itxmax if (ozout(1,k)>zero) then kk=kk+1 do i=1,nozdat ozout(i,kk)=ozout(i,k) end do endif end do ndata=kk nodata=ndata endif ! Write header record and data to output file for further processing call count_obs(ndata,nozdat,ilat,ilon,ozout,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((ozout(k,i),k=1,nozdat),i=1,ndata) nread=nmrecs end if ! Deallocate local arrays if(allocated(ozout))deallocate(ozout) if (obstype == 'sbuv2' .or. obstype == 'ompsnp') deallocate(poz) if (index(obstype,'mls')/=0) then if(allocated(hdrmlsl))deallocate(hdrmlsl) if(allocated(mlspres))deallocate(mlspres) if(allocated(mlsoz))deallocate(mlsoz) if(allocated(mlsozpc))deallocate(mlsozpc) if(allocated(ipos))deallocate(ipos) if(allocated(usage1))deallocate(usage1) end if if(index(obstype,'ompslp')/=0) then if(allocated(olpdtsq))deallocate(olpdtsq) if(allocated(lpsdvals))deallocate(lpsdvals) if(allocated(press))deallocate(press) if(allocated(omr))deallocate(omr) if(allocated(omrstd))deallocate(omrstd) if(allocated(ipos))deallocate(ipos) if(allocated(usage1))deallocate(usage1) endif ! Deallocate satthin arrays if (obstype == 'omi' .or. obstype == 'gome' .or. obstype=='ompsnm' .or. obstype == 'ompstc8' )call destroygrids return end subroutine read_ozone