subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& rmesh,jsatid,gstime,infile,lunout,obstype,& nread,ndata,nodata,twind,sis, & mype_root,mype_sub,npe_sub,mpi_comm_sub, nobs, & nrec_start,nrec_start_ears,nrec_start_db,dval_use,radmod) !$$$ subprogram documentation block ! . . . . ! subprogram: read_bufrtovs read bufr tovs 1b data ! prgmmr: treadon org: np23 date: 2003-09-13 ! ! abstract: This routine reads BUFR format TOVS 1b radiance ! (brightness temperature) files. Optionally, the data ! are thinned to a specified 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: ! 2003-09-13 treadon ! 2004-05-28 kleist - update subroutine call ! 2004-06-16 treadon - update documentation ! 2004-07-23 derber - make changes to eliminate obs. earlier in thinning ! 2004-07-29 treadon - add only to module use, add intent in/out ! 2004-10-15 derber - various changes to "quality" prediction used ! in data selection algorithm ! 2005-01-26 derber - land/sea determination and weighting for data selection ! 2005-02-10 treadon - correct spelling in runtime print message; specify ! _r_kind precision for real constants ! 2005-07-06 derber - add mhs and hirs/4 from NOAA-18, clean up code and ! modify data selection criteria ! 2005-09-08 derber - modify to use input group time window ! 2005-09-28 derber - modify to produce consistent surface info ! 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-10-24 derber - use more precise MHS scan properties (step & start) ! 2005-10-26 treadon - clean up formatting, add clarifying comments, correct ! ndata,ndata1 printout ! 2005-11-22 derber - include mean in bias correction ! 2005-11-29 parrish - modify getsfc to work for different regional options ! 2006-02-01 parrish - remove getsfc (different version called now in read_obs) ! 2006-02-03 derber - modify for new obs control and obs count ! 2006-02-11 liu - add ssu ! 2006-03-07 derber - correct error in nodata count ! 2006-04-27 derber - clean up code ! 2006-05-02 treadon - move close lnbufr to after 900 continue ! 2006-05-19 eliu - add logic to reset relative weight when all channels not used ! 2006-07-28 derber - add solar and satellite azimuth angles remove isflg from output ! 2007-01-18 derber - add kidsat for metop ! 2007-03-01 tremolet - measure time from beginning of assimilation window ! 2007-08-31 h.liu - add kidsat for n05, n06, n07, n08, n09, n10, n11, n12,tiros-n ! 2008-05-27 safford - rm unused vars ! 2008-09-08 lueken - merged ed's changes into q1fy09 code ! 2008-10-14 derber - properly use EARS data and use MPI_IO ! 2009-01-02 todling - remove unused vars ! 2009-01-09 gayno - add option to calculate surface fields based on ! size/shape of field of view. ! 2009-04-17 todling - zero out azimuth angle when unavailable (otherwise can't use old files) ! 2009-04-18 woollen - improve mpi_io interface with bufrlib routines ! 2009-04-21 derber - add ithin to call to makegrids ! 2009-12-20 gayno - modify for updated version of FOV surface code which ! calculates relative antenna power for some instruments. ! 2010-02-25 collard - changes to call to crtm_init for CRTM v2.0 ! 2010-06-29 zhu - add newpc4pred option ! 2010-07-12 zhu - include global offset for amsua in bias correction for adp_anglebc option ! 2010-09-02 zhu - add use_edges option ! 2010-10-12 zhu - use radstep and radstart from radinfo ! 2011-04-07 todling - newpc4pred now in radinfo ! 2011-04-08 li - (1) use nst_gsi, nstinfo, and add NSST vars ! (2) get zob, tz_tr (call skindepth and cal_tztr) ! (3) interpolate NSST Variables to Obs. location (call deter_nst) ! (4) add more elements (nstinfo) in data array ! 2011-05-05 todling - merge in Min-Jeong update for cloudy-radiance work ! 2011-05-20 mccarty - updated to read ATMS data ! 2011-07-04 todling - fixes to run either single or double precision ! 2011-08-01 lueken - removed deter_sfc subroutines, placed in new module deter_sfc_mod ! 2011-09-13 gayno - improve error handling for FOV-based sfc calculation ! (isfcalc=1) ! 2011-12-13 collard Replace find_edges code to speed up execution. ! 2011-12-14 collard Remove ATMS ! 2012-03-05 akella - nst now controlled via coupler ! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) ! 2013-12-20 zhu - change icw4crtm>0 to icw4crtm>10 (bug fix) ! 2014-01-31 mkim - added iql4crtm for all-sky mw radiance data assimilation ! 2014-04-27 eliu/zhu - add thinning options for AMSU-A under allsky condition ! 2015-02-23 Rancic/Thomas - add thin4d to time window logical ! 2015-08-20 zhu - use centralized radmod for all-sky and aerosol usages in radiance assimilation ! 2016-04-28 jung - added logic for RARS and direct broadcast from NESDIS/UW ! 2016-10-20 collard - fix to allow monitoring and limited assimilation of spectra when key ! channels are missing. ! 2018-04-19 eliu - allow data selection for precipitation-affected data ! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90. ! 2020-04021 s.sieron - change converting brightness temperatures to antenna temperautres to the ! other way around. added support for multiple SpcCoeff files and ACCoeffs ! ! input argument list: ! mype - mpi task id ! val_tovs - weighting factor applied to super obs ! ithin - flag to thin data ! isfcalc - method to calculate surface fields within FOV ! when one, calculate accounting for size/shape of FOV. ! otherwise, use bilinear interpolation. ! rmesh - thinning mesh size (km) ! jsatid - satellite to read ! gstime - analysis time in minutes from reference date ! infile - unit from which to read BUFR data ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! twind - input group time window(hours) ! sis - sensor/instrument/satellite indicator ! 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 ! dval_use - logical for using dval ! nrec_start - first subset with useful information ! nrec_start_ears - first ears subset with useful information ! nrec_start_db - first db subset with useful information ! ! output argument list: ! nread - number of BUFR TOVS 1b observations read ! ndata - number of BUFR TOVS 1b profiles retained for further processing ! nodata - number of BUFR TOVS 1b observations retained for further processing ! 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 satthin, only: super_val,itxmax,makegrids,destroygrids,checkob, & finalcheck,map2tgrid,score_crit use satthin, only: radthin_time_info,tdiff2crit use obsmod, only: time_window_max, ta2tb use radinfo, only: iuse_rad,newchn,cbias,predx,nusis,jpch_rad,air_rad,ang_rad, & use_edges,radedge1, radedge2, radstart,radstep,newpc4pred use radinfo, only: crtm_coeffs_path,adp_anglebc use gridmod, only: diagnostic_reg,regional,nlat,nlon,tll2xy,txy2ll,rlats,rlons use constants, only: deg2rad,zero,one,two,three,five,rad2deg,r60inv,r1000,h300,r100 use crtm_module, only: success, & crtm_kind => fp, & MAX_SENSOR_ZENITH_ANGLE use crtm_spccoeff, only: sc,crtm_spccoeff_load,crtm_spccoeff_destroy use ACCoeff_Define, only: ACCoeff_type use calc_fov_crosstrk, only : instrument_init, fov_cleanup, fov_check use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen use antcorr_application, only: remove_antcorr, apply_antcorr use mpeu_util, only: getindex use deter_sfc_mod, only: deter_sfc_fov,deter_sfc use gsi_nstcouplermod, only: nst_gsi,nstinfo use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter use mpimod, only: npe use radiance_mod, only: rad_obs_type use gsi_io, only: verbose implicit none ! Declare passed variables character(len=*),intent(in ) :: infile,obstype,jsatid character(len=20),intent(in ) :: sis integer(i_kind) ,intent(in ) :: mype,lunout,ithin integer(i_kind) ,intent(in ) :: nrec_start,nrec_start_ears,nrec_start_db integer(i_kind) ,intent(inout) :: isfcalc integer(i_kind) ,intent(inout) :: nread integer(i_kind),dimension(npe) ,intent(inout) :: nobs integer(i_kind) ,intent( out) :: ndata,nodata real(r_kind) ,intent(in ) :: rmesh,gstime,twind real(r_kind) ,intent(inout) :: val_tovs 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 logical, intent(in ) :: dval_use type(rad_obs_type),intent(in ) :: radmod ! Declare local parameters character(8),parameter:: fov_flag="crosstrk" ! change from 13 to 14 for sacv integer(i_kind),parameter:: n1bhdr=14 integer(i_kind),parameter:: n2bhdr=4 real(r_kind),parameter:: r360=360.0_r_kind real(r_kind),parameter:: tbmin=50.0_r_kind real(r_kind),parameter:: tbmax=550.0_r_kind ! Declare local variables logical hirs,msu,amsua,amsub,mhs,hirs4,hirs3,hirs2,ssu logical outside,iuse,assim,valid character(len=40):: infile2 character(len=8) :: subset character(len=80):: hdr1b,hdr2b integer(i_kind) ireadsb,ireadmg,irec,next,nrec_startx integer(i_kind) i,j,k,ifov,ntest,llll integer(i_kind) sacv integer(i_kind) iret,idate,nchanl,n,idomsfc(1) integer(i_kind) ich1,ich2,ich8,ich15,ich16,ich17 integer(i_kind) kidsat,instrument,maxinfo integer(i_kind) nmind,itx,nreal,nele,itt,ninstruments integer(i_kind) iskip,ichan2,ichan1,ichan15 integer(i_kind) lnbufr,ksatid,ichan8,isflg,ichan3,ich3,ich4,ich6 integer(i_kind) ilat,ilon,ifovmod integer(i_kind),dimension(5):: idate5 integer(i_kind) instr,ichan integer(i_kind) error_status,irecx,ierr integer(i_kind) radedge_min, radedge_max integer(i_kind),allocatable,dimension(:)::nrec character(len=20),dimension(1):: sensorlist real(r_kind) cosza,sfcr real(r_kind) ch1,ch2,ch3,ch8,d0,d1,d2,ch15,qval real(r_kind) ch1flg real(r_kind) expansion real(r_kind),dimension(0:3):: sfcpct real(r_kind),dimension(0:3):: ts real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10 real(r_kind) :: zob,tref,dtw,dtc,tz_tr real(r_kind) pred, pred_water, pred_not_water real(r_kind) rsat,dlat,panglr,dlon,rato,sstime,tdiff,t4dv real(r_kind) dlon_earth_deg,dlat_earth_deg,sat_aziang real(r_kind) dlon_earth,dlat_earth,r01 real(r_kind) crit1,step,start,ch8flg,dist1 real(r_kind) terrain,lza,df2,tt,lzaest real(r_kind),dimension(0:4):: rlndsea real(r_kind),allocatable,dimension(:,:):: data_all real(crtm_kind),allocatable,dimension(:):: data1b4 real(r_double),allocatable,dimension(:):: data1b8,data1b8x real(r_double),dimension(n1bhdr):: bfr1bhdr real(r_double),dimension(n2bhdr):: bfr2bhdr real(r_kind) disterr,disterrmax,cdist,dlon00,dlat00 real(r_kind) :: ptime,timeinflat,crit0 integer(i_kind) :: ithin_time,n_tbin,it_mesh logical :: critical_channels_missing,quiet logical :: print_verbose logical :: spc_coeff_found integer(i_kind) :: spc_coeff_versions character(len=80) :: spc_filename type(ACCoeff_type),dimension(3) :: accoeff_sets !************************************************************************** ! Initialize variables print_verbose=.false. if(verbose) print_verbose=.true. maxinfo=31 lnbufr = 15 disterrmax=zero ntest=0 ndata = 0 nodata = 0 nread = 0 ilon=3 ilat=4 if(nst_gsi>0) then call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype endif 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 ! Make thinning grids call makegrids(rmesh,ithin,n_tbin=n_tbin) ! Set various variables depending on type of data to be read hirs2 = obstype == 'hirs2' hirs3 = obstype == 'hirs3' hirs4 = obstype == 'hirs4' hirs = hirs2 .or. hirs3 .or. hirs4 msu = obstype == 'msu' amsua = obstype == 'amsua' amsub = obstype == 'amsub' mhs = obstype == 'mhs' ssu = obstype == 'ssu' ! instrument specific variables d1 = 0.754_r_kind d2 = -2.265_r_kind r01 = 0.01_r_kind ich1 = 1 !1 ich2 = 2 !2 ich3 = 3 !3 ich4 = 4 !4 ich6 = 6 !6 ich8 = 8 !8 ich15 = 15 !15 ich16 = 16 !16 ich17 = 17 !17 ! Set array index for surface-sensing channels ichan1 = newchn(sis,ich1) ichan2 = newchn(sis,ich2) ichan3 = newchn(sis,ich3) if (hirs) then ichan8 = newchn(sis,ich8) endif if (amsua) ichan15 = newchn(sis,ich15) if(jsatid == 'n05')kidsat=705 if(jsatid == 'n06')kidsat=706 if(jsatid == 'n07')kidsat=707 if(jsatid == 'tirosn')kidsat=708 if(jsatid == 'n08')kidsat=200 if(jsatid == 'n09')kidsat=201 if(jsatid == 'n10')kidsat=202 if(jsatid == 'n11')kidsat=203 if(jsatid == 'n12')kidsat=204 if(jsatid == 'n14')kidsat=205 if(jsatid == 'n15')kidsat=206 if(jsatid == 'n16')kidsat=207 if(jsatid == 'n17')kidsat=208 if(jsatid == 'n18')kidsat=209 if(jsatid == 'n19')kidsat=223 if(jsatid == 'metop-a')kidsat=4 if(jsatid == 'metop-b')kidsat=3 if(jsatid == 'metop-c')kidsat=5 if(jsatid == 'npp')kidsat=224 radedge_min = 0 radedge_max = 1000 do i=1,jpch_rad if (trim(nusis(i))==trim(sis)) then step = radstep(i) start = radstart(i) if (radedge1(i)/=-1 .and. radedge2(i)/=-1) then radedge_min=radedge1(i) radedge_max=radedge2(i) end if exit endif end do rato=1.1363987_r_kind if ( hirs ) then nchanl=19 ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 15._r_kind rlndsea(2) = 10._r_kind rlndsea(3) = 15._r_kind rlndsea(4) = 30._r_kind if (isfcalc == 1) then expansion=one ! use one for ir sensors ichan=-999 ! not used for hirs if (hirs2) then if(kidsat==203.or.kidsat==205)then instr=5 ! hirs-2i on noaa 11 and 14 elseif(kidsat==708.or.kidsat==706.or.kidsat==707.or. & kidsat==200.or.kidsat==201.or.kidsat==202.or. & kidsat==204)then instr=4 ! hirs-2 on tiros-n,noaa6-10,noaa12 else ! sensor/sat mismatch write(6,*) 'READ_BUFRTOVS: *** WARNING: HIRS2 SENSOR/SAT MISMATCH' instr=5 ! set to something to prevent abort endif elseif (hirs4) then instr=8 elseif (hirs3) then if(kidsat==206) then instr=6 ! noaa 15 elseif(kidsat==207.or.kidsat==208)then instr=7 ! noaa 16/17 else write(6,*) 'READ_BUFRTOVS: *** WARNING: HIRS3 SENSOR/SAT MISMATCH' instr=7 endif endif endif ! isfcalc == 1 else if ( msu ) then nchanl=4 if (isfcalc==1) then instr=10 ichan=-999 expansion=2.9_r_kind endif ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 20._r_kind rlndsea(2) = 15._r_kind rlndsea(3) = 20._r_kind rlndsea(4) = 100._r_kind else if ( amsua ) then nchanl=15 if (isfcalc==1) then instr=11 ichan=15 ! pick a surface sens. channel expansion=2.9_r_kind ! use almost three for microwave sensors. endif ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 15._r_kind rlndsea(2) = 10._r_kind rlndsea(3) = 15._r_kind rlndsea(4) = 100._r_kind else if ( amsub ) then nchanl=5 if (isfcalc==1) then instr=12 ichan=-999 expansion=2.9_r_kind endif ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 15._r_kind rlndsea(2) = 20._r_kind rlndsea(3) = 15._r_kind rlndsea(4) = 100._r_kind else if ( mhs ) then nchanl=5 if (isfcalc==1) then instr=13 ichan=1 ! all channels give similar result. expansion=2.9_r_kind endif ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 15._r_kind rlndsea(2) = 20._r_kind rlndsea(3) = 15._r_kind rlndsea(4) = 100._r_kind else if ( ssu ) then nchanl=3 if (isfcalc==1) then instr=9 ichan=-999 ! not used for this sensor expansion=one endif ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 15._r_kind rlndsea(2) = 10._r_kind rlndsea(3) = 15._r_kind rlndsea(4) = 30._r_kind end if ! If all channels of a given sensor are set to monitor or not ! assimilate mode (iuse_rad<1), reset relative weight to zero. ! We do not want such observations affecting the relative ! weighting between observations within a given thinning group. assim=.false. search: do i=1,jpch_rad if ((nusis(i)==sis) .and. (iuse_rad(i)>0)) then assim=.true. exit search endif end do search if (.not.assim) val_tovs=zero ! Initialize variables for use by FOV-based surface code. if (isfcalc == 1) then call instrument_init(instr,jsatid,expansion,valid) if (.not. valid) then if (assim) then write(6,*)'READ_BUFRTOVS: ***ERROR*** IN SETUP OF FOV-SFC CODE. STOP' call stop2(71) else call fov_cleanup isfcalc = 0 write(6,*)'READ_BUFRTOVS: ***ERROR*** IN SETUP OF FOV-SFC CODE' endif endif endif if (isfcalc==1) then if (amsub.or.mhs)then rlndsea(4) = max(rlndsea(0),rlndsea(1),rlndsea(2),rlndsea(3)) else rlndsea=0 endif endif ! Allocate arrays to hold all data for given satellite if(dval_use) maxinfo=maxinfo+2 nreal = maxinfo + nstinfo nele = nreal + nchanl hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS SACV' hdr2b ='SAZA SOZA BEARAZ SOLAZI' allocate(data_all(nele,itxmax),data1b8(nchanl),data1b4(nchanl),nrec(itxmax)) next=0 irec=0 ! Big loop over standard data feed and possible ears/db data ! llll=1 is normal feed, llll=2 EARS/RARS data, llll=3 DB/UW data) ears_db_loop: do llll= 1, 3 if(llll == 1)then nrec_startx=nrec_start infile2=trim(infile) ! Set bufr subset names based on type of data to read elseif(llll == 2) then nrec_startx=nrec_start_ears infile2=trim(infile)//'ears' ! Set bufr subset names based on type of data to read if(amsua .and. kidsat >= 200 .and. kidsat <= 207) cycle ears_db_loop elseif(llll == 3) then nrec_startx=nrec_start_db infile2=trim(infile)//'_db' ! Set bufr subset names based on type of data to read if(amsua .and. kidsat >= 200 .and. kidsat <= 207) cycle ears_db_loop end if ! Reopen unit to satellite bufr file open(lnbufr,file=trim(infile2),form='unformatted',status = 'old',iostat=ierr) if(ierr /= 0) cycle ears_db_loop call openbf(lnbufr,'IN',lnbufr) ! support multiple spc coefficient files for any given sensor if(amsua .or. amsub .or. mhs)then quiet=.not.verbose allocate(data1b8x(nchanl)) spc_coeff_versions = 0 spc_coeff_found = .true. do while (spc_coeff_found) if (spc_coeff_versions == 0) then sensorlist(1)=sis else i = spc_coeff_versions+1 write(sensorlist(1),'(a,a,i1)') trim(sis),'_v',i end if if( crtm_coeffs_path /= "" ) then if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_BUFRTOVS: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"' end if spc_filename = trim(crtm_coeffs_path) // trim(sensorlist(1)) // '.SpcCoeff.bin' INQUIRE(FILE=trim(spc_filename), EXIST=spc_coeff_found) if (.NOT. spc_coeff_found) then if (spc_coeff_versions == 0) then write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& ' TERMINATE PROGRAM EXECUTION' call stop2(71) else write(6,*)'READ_BUFRTOVS: ', spc_coeff_versions, ' versions of SpcCoeff found for ', trim(sis) end if else spc_coeff_versions = spc_coeff_versions+1 if( crtm_coeffs_path /= "" ) then error_status = crtm_spccoeff_load(sensorlist,& File_Path = crtm_coeffs_path, quiet=quiet ) else error_status = crtm_spccoeff_load(sensorlist,quiet=quiet) endif if (error_status /= success) then write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& ' despite file ',trim(spc_filename),' existing, TERMINATE PROGRAM EXECUTION' call stop2(71) endif ninstruments = size(sc) instrument_loop: do n=1,ninstruments if(sis == sc(n)%sensor_id)then instrument=n exit instrument_loop end if end do instrument_loop if(instrument == 0)then write(6,*)' failure to find instrument in read_bufrtovs ',sis end if accoeff_sets(spc_coeff_versions) = sc(instrument)%ac ! deallocate crtm info error_status = crtm_spccoeff_destroy() if (error_status /= success) & write(6,*)'OBSERVER: ***ERROR*** crtm_spccoeff_destroy error_status=',error_status end if end do end if ! Loop to read bufr file irecx=0 read_subset: do while(ireadmg(lnbufr,subset,idate)>=0) irecx=irecx+1 if(irecx < nrec_startx) cycle read_subset irec=irec+1 next=next+1 if(next == npe_sub)next=0 if(next/=mype_sub)cycle read_subset read_loop: do while (ireadsb(lnbufr)==0) ! Read header record. (llll=1 is normal feed, 2=EARS data, 3=DB data) call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b) ! Extract satellite id. If not the one we want, read next record ksatid=nint(bfr1bhdr(1)) if(ksatid /= kidsat) cycle read_subset rsat=bfr1bhdr(1) ! If msu, drop obs from first (1) and last (11) scan positions ifov = nint(bfr1bhdr(2)) if (use_edges) then if (msu .and. (ifov==1 .or. ifov==11)) cycle read_loop else if ((ifov < radedge_min .OR. ifov > radedge_max )) then cycle read_loop end if ! Check that ifov is not out of range of cbias dimension if (ifov < 1 .OR. ifov > 90) cycle read_loop ! Extract observation location and other required information if(abs(bfr1bhdr(11)) <= 90._r_kind .and. abs(bfr1bhdr(12)) <= r360)then dlat_earth = bfr1bhdr(11) dlon_earth = bfr1bhdr(12) elseif(abs(bfr1bhdr(9)) <= 90._r_kind .and. abs(bfr1bhdr(10)) <= r360)then dlat_earth = bfr1bhdr(9) dlon_earth = bfr1bhdr(10) else cycle read_loop end if if(dlon_earth=r360) dlon_earth = dlon_earth-r360 dlat_earth_deg = dlat_earth dlon_earth_deg = dlon_earth dlat_earth = dlat_earth*deg2rad dlon_earth = dlon_earth*deg2rad ! Regional case if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(diagnostic_reg) then call txy2ll(dlon,dlat,dlon00,dlat00) ntest=ntest+1 cdist=sin(dlat_earth)*sin(dlat00)+cos(dlat_earth)*cos(dlat00)* & (sin(dlon_earth)*sin(dlon00)+cos(dlon_earth)*cos(dlon00)) cdist=max(-one,min(cdist,one)) disterr=acos(cdist)*rad2deg disterrmax=max(disterrmax,disterr) end if ! Check to see if in domain if(outside) cycle read_loop ! Global case 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) = bfr1bhdr(3) !year idate5(2) = bfr1bhdr(4) !month idate5(3) = bfr1bhdr(5) !day idate5(4) = bfr1bhdr(6) !hour idate5(5) = bfr1bhdr(7) !minute call w3fs21(idate5,nmind) t4dv= (real((nmind-iwinbgn),r_kind) + bfr1bhdr(8)*r60inv)*r60inv ! add in seconds sstime= real(nmind,r_kind) + bfr1bhdr(8)*r60inv ! add in seconds tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle read_loop else if(abs(tdiff) > twind) cycle read_loop endif nread=nread+nchanl terrain = 50._r_kind if(llll == 1)terrain = 0.01_r_kind*abs(bfr1bhdr(13)) crit0 = 0.01_r_kind + terrain if (llll > 1 ) crit0 = crit0 + r100 * float(llll) timeinflat=two 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_loop ! Extract satellite antenna corrections version number if (llll > 1) then sacv = nint(bfr1bhdr(14)) if (sacv > spc_coeff_versions) then write(6,*) 'READ_BUFRTOVS WARNING sacv greater than spc_coeff_versions' end if else ! normal feed doesn't have antenna correction, so set sacv to 0 sacv = 0 end if call ufbint(lnbufr,bfr2bhdr,n2bhdr,1,iret,hdr2b) ! Set common predictor parameters ifovmod=ifov ! Account for assymetry due to satellite build error if(hirs .and. ((jsatid == 'n16') .or. (jsatid == 'n17'))) & ifovmod=ifovmod+1 panglr=(start+float(ifovmod-1)*step)*deg2rad lzaest = asin(rato*sin(panglr)) if( msu .or. hirs2 .or. ssu)then lza = lzaest else lza = bfr2bhdr(1)*deg2rad ! local zenith angle if((amsua .and. ifovmod <= 15) .or. & (amsub .and. ifovmod <= 45) .or. & (mhs .and. ifovmod <= 45) .or. & (hirs .and. ifovmod <= 28) ) lza=-lza end if ! Check for errors in satellite zenith angles if(abs(lzaest-lza)*rad2deg > one) then write(6,*)' READ_BUFRTOVS WARNING uncertainty in lza ', & lza*rad2deg,lzaest*rad2deg,sis,ifovmod,start,step cycle read_loop end if if(abs(lza)*rad2deg > MAX_SENSOR_ZENITH_ANGLE) then write(6,*)'READ_BUFRTOVS WARNING lza error ',bfr2bhdr(1),panglr cycle read_loop end if sat_aziang=bfr2bhdr(3) if (abs(sat_aziang) > r360) then sat_aziang=zero !_RT write(6,*) 'READ_BUFRTOVS: bad azimuth angle ',sat_aziang !_RT cycle read_loop endif ! Read data record. Increment data counter ! TMBR is actually the antenna temperature for most microwave sounders. ! Changed from accepting TMBR and converting TMBRST to ! antenna temperature (remove_antcorr), to converting TMBR to ! brightness temperature (add_antcorr) and accepting TMBRST ! if (ta2tb) then if (llll == 1) then call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') if ( (amsua .or. amsub .or. mhs) .and. & .not.(jsatid == 'n15' .or. jsatid == 'n16') )then ! convert antenna temperature to brightness temperature, ! unless the satellite is n15 or n16, because tranamsua ! does this conversion because the coefficient files exist ! for it to use data1b8x=data1b8 data1b4=data1b8 !call apply_antcorr(accoeff_sets(spc_coeff_versions),ifov,data1b4) call apply_antcorr(accoeff_sets(1),ifov,data1b4) data1b8=data1b4 do j=1,nchanl if(data1b8x(j) > r1000) data1b8(j) = 1000000._r_kind end do end if else ! EARS / DB call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') !if ( amsua .or. amsub .or. mhs .AND. sacv .ne. spc_coeff_versions)then if ( amsua .or. amsub .or. mhs .AND. spc_coeff_versions /= 1 .AND. sacv /= 1)then ! convert brightness temperature to antenna temperature using ! the satellite antenna correction version (sacv) used by the ! data originator, ! then convert back to brightness temperature using the version ! of parameters used by the CRTM data1b8x=data1b8 data1b4=data1b8 call remove_antcorr(accoeff_sets(sacv),ifov,data1b4) !call apply_antcorr(accoeff_sets(spc_coeff_versions),ifov,data1b4) call apply_antcorr(accoeff_sets(1),ifov,data1b4) data1b8=data1b4 do j=1,nchanl if(data1b8x(j) > r1000) data1b8(j) = 1000000._r_kind end do end if end if else if (llll == 1) then call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') else ! EARS / DB call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') if ( amsua .or. amsub .or. mhs )then data1b8x=data1b8 data1b4=data1b8 call remove_antcorr(accoeff_sets(1),ifov,data1b4) data1b8=data1b4 do j=1,nchanl if(data1b8x(j) > r1000)data1b8(j) = 1000000._r_kind end do end if end if endif ! Transfer observed brightness temperature to work array. If any ! temperature exceeds limits, reset observation to "bad" value iskip=0 critical_channels_missing = .false. do j=1,nchanl if (data1b8(j) < tbmin .or. data1b8(j) > tbmax) then iskip = iskip + 1 ! Flag profiles where key channels are bad if(( msu .and. j == ich1) .or. & (amsua .and. (j == ich1 .or. j == ich2 .or. j == ich3 .or. & j == ich4 .or. j == ich6 .or. j == ich15 )) .or.& (hirs .and. (j == ich8 )) .or. & (amsub .and. j == ich1) .or. & (mhs .and. (j == ich1 .or. j == ich2)) ) critical_channels_missing = .true. endif end do if (iskip >= nchanl) cycle read_loop ! Map obs to thinning grid crit1 = crit1 + 10._r_kind*float(iskip) call checkob(dist1,crit1,itx,iuse) if(.not. iuse)cycle read_loop ! Determine surface properties based on ! sst and land/sea/ice mask ! ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! FOV-based surface code requires fov number. if out-of-range, then ! skip this ob. if (isfcalc == 1) then call fov_check(ifov,instr,ichan,valid) if (.not. valid) cycle read_loop ! When isfcalc is one, calculate surface fields based on size/shape of fov. ! Otherwise, use bilinear method. call deter_sfc_fov(fov_flag,ifov,instr,ichan,sat_aziang,dlat_earth_deg,& dlon_earth_deg,expansion,t4dv,isflg,idomsfc(1), & sfcpct,vfr,sty,vty,stp,sm,ff10,sfcr,zz,sn,ts,tsavg) else call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg, & idomsfc(1),sfcpct,ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) endif crit1 = crit1 + rlndsea(isflg) call checkob(dist1,crit1,itx,iuse) if(.not. iuse)cycle read_loop if (critical_channels_missing) then pred=1.0e8_r_kind else ! Set data quality predictor if (msu) then if (newpc4pred) then ch1 = data1b8(ich1)-ang_rad(ichan1)*cbias(ifov,ichan1)- & predx(1,ichan1)*air_rad(ichan1) else ch1 = data1b8(ich1)-ang_rad(ichan1)*cbias(ifov,ichan1)- & r01*predx(1,ichan1)*air_rad(ichan1) end if ch1flg = tsavg-ch1 if(isflg == 0)then pred = 100._r_kind-min(ch1flg,100.0_r_kind) else pred = abs(ch1flg) end if else if (hirs) then if (newpc4pred) then ch8 = data1b8(ich8) -ang_rad(ichan8)*cbias(ifov,ichan8)- & predx(1,ichan8)*air_rad(ichan8) else ch8 = data1b8(ich8) -ang_rad(ichan8)*cbias(ifov,ichan8)- & r01*predx(1,ichan8)*air_rad(ichan8) end if ch8flg = tsavg-ch8 pred = 10.0_r_kind*max(zero,ch8flg) else if (amsua) then ! Remove angle dependent pattern (not mean) if (adp_anglebc .and. newpc4pred) then ch1 = data1b8(ich1)-ang_rad(ichan1)*cbias(ifov,ichan1) ch2 = data1b8(ich2)-ang_rad(ichan2)*cbias(ifov,ichan2) ch15= data1b8(ich15)-ang_rad(ichan15)*cbias(ifov,ichan15) else ch1 = data1b8(ich1)-ang_rad(ichan1)*cbias(ifov,ichan1)+ & air_rad(ichan1)*cbias(15,ichan1) ch2 = data1b8(ich2)-ang_rad(ichan2)*cbias(ifov,ichan2)+ & air_rad(ichan2)*cbias(15,ichan2) ch15= data1b8(ich15)-ang_rad(ichan15)*cbias(ifov,ichan15) end if if (isflg == 0 .and. ch1<285.0_r_kind .and. ch2<285.0_r_kind) then cosza = cos(lza) d0 = 8.24_r_kind - 2.622_r_kind*cosza + 1.846_r_kind*cosza*cosza qval=cosza*(d0+d1*log(285.0_r_kind-ch1)+d2*log(285.0_r_kind-ch2)) if (radmod%lcloud_fwd) then ! no preference in selecting clouds/precipitation ! qval=zero ! favor non-precipitating clouds qval=-113.2_r_kind+(2.41_r_kind-0.0049_r_kind*ch1)*ch1 + & 0.454_r_kind*ch2-ch15 if (qval>=9.0_r_kind) then qval=1000.0_r_kind*qval else qval=zero end if if (radmod%lprecip) qval=zero ! favor thinner clouds ! cosza = cos(lza) ! d0= 8.24_r_kind - 2.622_r_kind*cosza + 1.846_r_kind*cosza*cosza ! qval=cosza*(d0+d1*log(285.0_r_kind-ch1)+d2*log(285.0_r_kind-ch2)) ! if (qval>0.2_r_kind) then ! qval=1000.0_r_kind*qval ! else ! qval=zero ! end if end if pred = max(zero,qval)*100.0_r_kind else if (adp_anglebc .and. newpc4pred) then ch3 = data1b8(ich3)-ang_rad(ichan3)*cbias(ifov,ichan3) ch15 = data1b8(ich15)-ang_rad(ichan15)*cbias(ifov,ichan15) else ch3 = data1b8(ich3)-ang_rad(ichan3)*cbias(ifov,ichan3)+ & air_rad(ichan3)*cbias(15,ichan3) ch15 = data1b8(ich15)-ang_rad(ichan15)*cbias(ifov,ichan15)+ & air_rad(ichan15)*cbias(15,ichan15) end if pred = abs(ch1-ch15) if(ch1-ch15 >= three) then df2 = 5.10_r_kind +0.78_r_kind*ch1-0.96_r_kind*ch3 tt = 168._r_kind-0.49_r_kind*ch15 if(ch1 > 261._r_kind .or. ch1 >= tt .or. & (ch15 <= 273._r_kind .and. df2 >= 0.6_r_kind))then pred = 100._r_kind end if end if endif ! sval=-113.2_r_kind+(2.41_r_kind-0.0049_r_kind*ch1)*ch1 + & ! 0.454_r_kind*ch2-ch15 else if (amsub .or. mhs) then if (newpc4pred) then ch1 = data1b8(ich1)-ang_rad(ichan1)*cbias(ifov,ichan1)- & predx(1,ichan1)*air_rad(ichan1) ch2 = data1b8(ich2)-ang_rad(ichan2)*cbias(ifov,ichan2)- & predx(1,ichan2)*air_rad(ichan2) else ch1 = data1b8(ich1)-ang_rad(ichan1)*cbias(ifov,ichan1)- & r01*predx(1,ichan1)*air_rad(ichan1) ch2 = data1b8(ich2)-ang_rad(ichan2)*cbias(ifov,ichan2)- & r01*predx(1,ichan2)*air_rad(ichan2) end if pred_water = zero if(sfcpct(0) > zero)then cosza = cos(lza) if(ch2 < h300)then pred_water = (0.13_r_kind*(ch1+33.58_r_kind*log(h300-ch2)- & 341.17_r_kind))*five else pred_water = 100._r_kind end if end if pred_not_water = 42.72_r_kind + 0.85_r_kind*ch1-ch2 pred = (sfcpct(0)*pred_water) + ((one-sfcpct(0))*pred_not_water) pred = max(zero,pred) endif end if ! Compute "score" for observation. All scores>=0.0. Lowest score is "best" crit1 = crit1+pred call finalcheck(dist1,crit1,itx,iuse) if(.not. iuse)cycle read_loop ! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr if(nst_gsi>0) then tref = ts(0) dtw = zero dtc = zero tz_tr = one if(sfcpct(0)>zero) then call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr) endif endif ! Load selected observation into data array data_all(1 ,itx)= rsat ! satellite ID data_all(2 ,itx)= t4dv ! time data_all(3 ,itx)= dlon ! grid relative longitude data_all(4 ,itx)= dlat ! grid relative latitude data_all(5 ,itx)= lza ! local zenith angle data_all(6 ,itx)= bfr2bhdr(3) ! local azimuth angle data_all(7 ,itx)= panglr ! look angle data_all(8 ,itx)= ifov ! scan position data_all(9 ,itx)= bfr2bhdr(2) ! solar zenith angle data_all(10,itx)= bfr2bhdr(4) ! solar azimuth angle data_all(11,itx) = sfcpct(0) ! sea percentage of data_all(12,itx) = sfcpct(1) ! land percentage data_all(13,itx) = sfcpct(2) ! sea ice percentage data_all(14,itx) = sfcpct(3) ! snow percentage data_all(15,itx)= ts(0) ! ocean skin temperature data_all(16,itx)= ts(1) ! land skin temperature data_all(17,itx)= ts(2) ! ice skin temperature data_all(18,itx)= ts(3) ! snow skin temperature data_all(19,itx)= tsavg ! average skin temperature data_all(20,itx)= vty ! vegetation type data_all(21,itx)= vfr ! vegetation fraction data_all(22,itx)= sty ! soil type data_all(23,itx)= stp ! soil temperature data_all(24,itx)= sm ! soil moisture data_all(25,itx)= sn ! snow depth data_all(26,itx)= zz ! surface height data_all(27,itx)= idomsfc(1) + 0.001_r_kind ! dominate surface type data_all(28,itx)= sfcr ! surface roughness data_all(29,itx)= ff10 ! ten meter wind factor data_all(30,itx) = dlon_earth_deg ! earth relative longitude (deg) data_all(31,itx) = dlat_earth_deg ! earth relative latitude (deg) if(dval_use) then data_all(32,itx)= val_tovs data_all(33,itx)= itt end if if(nst_gsi>0) then data_all(maxinfo+1,itx) = tref ! foundation temperature data_all(maxinfo+2,itx) = dtw ! dt_warm at zob data_all(maxinfo+3,itx) = dtc ! dt_cool at zob data_all(maxinfo+4,itx) = tz_tr ! d(Tz)/d(Tr) endif do i=1,nchanl data_all(i+nreal,itx)=data1b8(i) end do nrec(itx)=irec ! End of bufr read loops enddo read_loop enddo read_subset call closbf(lnbufr) close(lnbufr) if (allocated(data1b8x)) deallocate(data1b8x) end do ears_db_loop deallocate(data1b8,data1b4) call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) ! if(mype_sub==mype_root)then do n=1,ndata do i=1,nchanl if(data_all(i+nreal,n) > tbmin .and. & data_all(i+nreal,n) < tbmax)nodata=nodata+1 end do end do if(dval_use .and. assim)then do n=1,ndata itt=nint(data_all(33,n)) super_val(itt)=super_val(itt)+val_tovs end do end if ! Write final set of "best" observations to output file call count_obs(ndata,nele,ilat,ilon,data_all,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata) end if ! Deallocate local arrays deallocate(data_all,nrec) ! Deallocate satthin arrays call destroygrids ! Deallocate FOV surface code arrays and nullify pointers. if (isfcalc == 1) call fov_cleanup if(diagnostic_reg.and.ntest>0) write(6,*)'READ_BUFRTOVS: ',& 'mype,ntest,disterrmax=',mype,ntest,disterrmax ! End of routine return end subroutine read_bufrtovs