subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& prsl_full) !$$$ subprogram documentation block ! . . . . ! subprogram: read_prepbuf read obs from prepbufr file ! prgmmr: parrish org: np22 date: 1990-10-07 ! ! abstract: This routine reads conventional data found in the prepbufr ! file. Specific observation types read by this routine ! include surface pressure, temperature, winds (components ! and speeds), moisture and total precipitable water. ! ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 1990-10-07 parrish ! 1998-05-15 weiyu yang ! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version ! 2004-02-13 derber, j. - clean up and modify vertical weighting ! 2004-06-16 treadon - update documentation ! 2004-07-23 derber - modify to include conventional sst ! 2004-07-29 treadon - add only to module use, add intent in/out ! 2004-07-30 derber - generalize number of data records per obs type ! 2004-08-26 derber - fix many errors in reading of sst data ! 2004-08-27 kleist - modify pressure calculation ! 2004-10-28 kleist - correct array index bug in hybrid pressure calculation ! 2004-11-16 treadon - deallocate(etabl) prior to exiting routine ! 2005-02-10 treadon - add call destroygrids for obstype = sst ! 2005-05-24 pondeca - add surface analysis option ! 2005-02-24 treadon - remove hardwired setting of metar ps obs error ! 2005-05-27 derber - reduce t, uv, ps error limits ! 2005-05-27 kleist/derber - add option to read in new ob error table ! 2005-07-19 derber - clean up code a bit ! 2005-07-27 derber - add print of monitoring and reject data ! 2005-08-02 derber - modify to use convinfo file ! 2005-09-08 derber - modify to use input group time window ! 2005-10-11 treadon - change convinfo read to free format ! 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-26 treadon - add routine tag to convinfo printout ! 2006-02-03 derber - modify to count read/keep data and new obs control ! 2006-02-03 treadon - use interpolated guess 3d pressure field in errormod ! 2006-02-08 derber - modify to use new convinfo module ! 2006-02-09 treadon - save height for wind observations ! 2006-02-23 kistler - modify to add optional data thinning ! 2006-02-23 kistler - raob instument as subtype and solar elv angle computed ! 2006-02-24 derber - modify to take advantage of convinfo module ! 2006-04-03 derber - modify to properly handle height of surface obs ! 2006-04-05 wu - changes to read in GPS IPW (type 153) ! 2006-05-18 middlecoff/treadon - add huge_i_kind upper limit on nint ! 2006-05-29 treadon - increase nreal to pass more information to setup routines ! 2006-06-08 su - added the option to turn off oiqc ! 2006-06-21 wu - deallocate etabl array ! 2006-07-28 derber - temporarily add subtype for meteosat winds based on sat ID ! 2006-07-31 kleist - change to surface pressure ob error from ln(ps) to ps(cb) ! 2006-10-25 sienkiewicz - add blacklist of raob data ! 2006-12-01 todling - embedded blacklist into a module ! 2007-02-13 parrish - add ability to use obs files with ref time different from analysis time ! 2007-02-20 wu - correct errors in quality mark checks ! 2007-03-01 tremolet - measure time from beginning of assimilation window ! 2007-03-15 su - remove the error table reading part to a subroutine ! 2007-04-24 wu - add TAMDAR (134) to be used as sensible T ! 2007-05-17 kleist - generalize flag for virtual/sensible temperature obs ! 2007-09-28 treadon - truncate/expand obs time to remove extraneous bits ! 2007-10-03 su - Add reading qc mark from satellite wind ! 2007-10-24 Pondeca - add ability to use use_list on mesonet winds ! 2007-11-03 su - modify conventional thinning algorithm ! 2008-03-28 wu - add code to generate optional observation perturbations ! 2008-03-31 li.bi - add ascat ! 2008-05-27 safford - rm unused vars and uses ! 2008-06-02 treadon - check iret from inital readmg and act accordingly ! 2008-09-08 lueken - merged ed's changges into q1fy09 code ! 2008-21-25 todling - adapted Tremolet 2007-03-01 change of time window ! - remove unused vars ! 2009-07-08 pondeca - add ability to convert virtual temperature ! obs into sensible temperature for 2dvar ! 2009-07-08 park,pondeca - add option to use the hilbert curve-based ! cross-validation for 2dvar ! 2009-07-08 pondeca - move handling of "provider use_list" for mesonet winds ! to the new module sfcobsqc ! 2010-03-29 hu - add code to read cloud observation from METAR and NESDIS cloud products ! 2010-05-15 kokron - safety measure: initialize cdata_all to zero ! 2010-08-23 tong - add flg as an input argument of map3grids, so that the subroutine can be used for ! thinning grid with either pressure or height as the vertical coordinate. ! flg=-1 for prepbufr data thinning grid (pressure as the vertical coordinate). ! 2010-09-08 parrish - remove subroutine check_rotate_wind. This was a debug routine introduced when ! the reference wind rotation angle was stored as an angle, beta_ref. This field ! had a discontinuity at the date line (180E), which resulted in erroneous wind ! rotation angles for a small number of winds whose rotation angle was interpolated ! from beta_ref values across the discontinuity. This was fixed by replacing the ! beta_ref field with cos_beta_ref, sin_beta_ref. ! 2010-10-19 wu - add code to limit regional use of MAP winds with P less than 400 mb ! 2010-11-13 su - skip satellite winds from prepbufr ! 2010-11-18 treadon - add check for small POB (if POB=241 & .and. ictype(nc) <260) then cycle else ntmatch=ntmatch+1 ntxall(ntmatch)=nc endif end if if(trim(ioctype(nc)) == trim(obstype) .and. abs(icuse(nc)) <= 1)then if(.not.use_prepb_satwnd .and. trim(ioctype(nc)) == 'uv' .and. ictype(nc) >=241 & .and. ictype(nc) <260) then cycle else ithin=ithin_conv(nc) if(ithin > 0)then ntread=ntread+1 ntx(ntread)=nc end if endif end if end do if(ntmatch == 0)then write(6,*) ' no matching obstype found in obsinfo ',obstype return end if allocate(lmsg(nmsgmax,ntread)) lmsg = .false. maxobs=0 tab=0 nmsg=0 nrep=0 ntb = 0 msg_report: do while (ireadmg(lunin,subset,idate) == 0) if(.not.use_prepb_satwnd .and. trim(subset) == 'SATWND') cycle msg_report ! Time offset if(nmsg == 0) call time_4dvar(idate,toff) nmsg=nmsg+1 if (nmsg>nmsgmax) then write(6,*)'READ_PREPBUFR: messages exceed maximum ',nmsgmax call stop2(50) endif loop_report: do while (ireadsb(lunin) == 0) ntb = ntb+1 nrep(nmsg)=nrep(nmsg)+1 if (ntb>mxtb) then write(6,*)'READ_PREPBUFR: reports exceed maximum ',mxtb call stop2(50) endif ! Extract type information call ufbint(lunin,hdr,4,1,iret,hdstr2) kx=hdr(1) if(twodvar_regional)then ! If running in 2d-var (surface analysis) mode, check to see if observation ! is surface type. If not, read next observation report from bufr file sfctype=(kx>179.and.kx<190).or.(kx>=280.and.kx<=290).or. & (kx>=192.and.kx<=199).or.(kx>=292.and.kx<=299) if (.not.sfctype ) cycle loop_report end if ! temporary specify iobsub until put in bufr file iobsub = 0 if(kx == 280) iobsub=hdr(3) if(use_prepb_satwnd .and. (kx == 243 .or. kx == 253 .or. kx == 254)) iobsub = hdr(2) if(use_prepb_satwnd .and. kx == 245 ) then if(hdr(2) == 259.0_r_kind) iobsub = 15 endif ! For the satellite wind to get quality information and check if it will be used if(use_prepb_satwnd .and. (kx == 243 .or. kx == 253 .or. kx ==254) ) then call ufbint(lunin,satqc,1,1,iret,satqcstr) if(satqc(1) < 85.0_r_double) cycle loop_report ! QI w/o fcst (su's setup !! if(satqc(2) <= 80.0_r_double) cycle loop_report ! QI w/ fcst (old prepdata) endif ! Check for blacklisting of station ID if (blacklst .and. ibcnt > 0) then stnid = transfer(hdr(4),stnid) do i = 1,ibcnt if( kx == blkkx(i) .and. stnid == blkstns(i) ) then write(6,*)'READ_PREPBUFR: blacklist station ',stnid, & 'for obstype ',trim(obstype),' and kx=',kx cycle loop_report endif enddo endif ! Match ob to proper convinfo type ncsave=0 matchloop:do ncx=1,ntmatch nc=ntxall(ncx) if (kx /= ictype(nc))cycle ! Find convtype which match ob type and subtype if(icsubtype(nc) == iobsub) then ncsave=nc exit matchloop else ! Find convtype which match ob type and subtype group (isubtype == ?*) ! where ? specifies the group and icsubtype = ?0) ixsub=icsubtype(nc)/10 iosub=iobsub/10 isubsub=icsubtype(nc)-ixsub*10 if(ixsub == iosub .and. isubsub == 0) then ncsave=nc ! Find convtype which match ob type and subtype is all remaining ! (icsubtype(nc) = 0) else if (ncsave == 0 .and. icsubtype(nc) == 0) then ncsave=nc end if end if end do matchloop ! Save information for next read if(ncsave /= 0) then call ufbint(lunin,levdat,1,255,levs,levstr) maxobs=maxobs+max(1,levs) nx=1 if(ithin_conv(ncsave) > 0)then do ii=2,ntread if(ntx(ii) == ncsave)nx=ii end do end if tab(ntb,1)=ncsave tab(ntb,2)=nx tab(ntb,3)=levs lmsg(nmsg,nx) = .true. end if end do loop_report enddo msg_report if (nmsg==0) goto 900 write(6,*)'READ_PREPBUFR: messages/reports = ',nmsg,'/',ntb,' ntread = ',ntread if(tob) write(6,*)'READ_PREPBUFR: time offset is ',toff,' hours.' !------------------------------------------------------------------------ ! Obtain program code (VTCD) associated with "VIRTMP" step if(tob)call ufbqcd(lunin,'VIRTMP',vtcd) call init_rjlists call init_aircraft_rjlists if (lhilbert) call init_hilbertcurve(maxobs) if (twodvar_regional) call init_ndfdgrid ! loop over convinfo file entries; operate on matches allocate(cdata_all(nreal,maxobs),isort(maxobs)) isort = 0 cdata_all=zero nread=0 ntest=0 nvtest=0 nchanl=0 ilon=2 ilat=3 loop_convinfo: do nx=1, ntread use_all = .true. ithin=0 if(nx > 1) then nc=ntx(nx) ithin=ithin_conv(nc) if (ithin > 0 ) then rmesh=rmesh_conv(nc) pmesh=pmesh_conv(nc) use_all = .false. if(pmesh > zero) then pflag=1 nlevp=r1200/pmesh else pflag=0 nlevp=nsig endif xmesh=rmesh call make3grids(xmesh,nlevp) if (.not.use_all) then allocate(presl_thin(nlevp)) if (pflag==1) then do k=1,nlevp presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth enddo endif endif write(6,*)'READ_PREPBUFR: obstype,ictype(nc),rmesh,pflag,nlevp,pmesh=',& trim(ioctype(nc)),ictype(nc),rmesh,pflag,nlevp,pmesh endif endif call closbf(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) ! Big loop over prepbufr file ntb = 0 nmsg = 0 icntpnt=0 icntpnt2=0 loop_msg: do while (ireadmg(lunin,subset,idate)== 0) if(.not.use_prepb_satwnd .and. trim(subset) =='SATWND') cycle loop_msg nmsg = nmsg+1 if(.not.lmsg(nmsg,nx)) then do i=ntb+1,ntb+nrep(nmsg) icntpnt2=icntpnt2+tab(i,3) end do ntb=ntb+nrep(nmsg) cycle loop_msg ! no useable reports this mesage, skip ahead report count end if loop_readsb: do while(ireadsb(lunin) == 0) ! use msg lookup table to decide which messages to skip ! use report id lookup table to only process matching reports ntb = ntb+1 if(icntpnt < icntpnt2)icntpnt=icntpnt2 icntpnt2=icntpnt2+tab(ntb,3) nc=tab(ntb,1) if(nc <= 0 .or. tab(ntb,2) /= nx) cycle loop_readsb ! Extract type, date, and location information call ufbint(lunin,hdr,8,1,iret,hdstr) if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_readsb if(hdr(2)== r360)hdr(2)=hdr(2)-r360 if(hdr(2) < zero)hdr(2)=hdr(2)+r360 dlon_earth=hdr(2)*deg2rad dlat_earth=hdr(3)*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 loop_readsb ! check to see if outside regional domain else dlat = dlat_earth dlon = dlon_earth call grdcrd(dlat,1,rlats,nlat,1) call grdcrd(dlon,1,rlons,nlon,1) endif !------------------------------------------------------------------------ if(offtime_data) then ! in time correction for observations to account for analysis ! time being different from obs file time. write(date,'( i10)') idate read (date,'(i4,3i2)') iy,im,idd,ihh idate5(1)=iy idate5(2)=im idate5(3)=idd idate5(4)=ihh idate5(5)=0 call w3fs21(idate5,minobs) ! obs ref time in minutes relative to historic date idate5(1)=iadate(1) idate5(2)=iadate(2) idate5(3)=iadate(3) idate5(4)=iadate(4) idate5(5)=0 call w3fs21(idate5,minan) ! analysis ref time in minutes relative to historic date ! Add obs reference time, then subtract analysis time to get obs time relative to analysis time_correction=float(minobs-minan)*r60inv else time_correction=zero end if timeobs=real(real(hdr(4),r_single),r_double) t4dv=timeobs + toff zeps=1.0e-8_r_kind if (t4dv -zeps) t4dv=zero if (t4dv>winlen.and.t4dvwinlen) .and. .not.driftl) cycle loop_readsb ! outside time window else if((real(abs(time)) > real(ctwind(nc)) .or. real(abs(time)) > real(twindin)) & .and. .not. driftl)cycle loop_readsb ! outside time window endif timex=time ! If ASCAT data, determine primary surface type. If not open sea, ! skip this observation. This check must be done before thinning. if (kx==290 .or. kx==289 .or. kx==285) then call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg) if (isflg /= 0) cycle loop_readsb if (tsavg <= 273.0_r_kind) cycle loop_readsb endif sfctype=(kx>179.and.kx<190).or.(kx>=280.and.kx<=290).or. & (kx>=192.and.kx<=199).or.(kx>=292.and.kx<=299) if (sfctype) then call ufbint(lunin,r_prvstg,1,1,iret,prvstr) call ufbint(lunin,r_sprvstg,1,1,iret,sprvstr) else c_prvstg=cspval c_sprvstg=cspval endif ! Extract data information on levels call ufbint(lunin,obsdat,11,255,levs,obstr) call ufbint(lunin,qcmark,8,255,levs,qcstr) call ufbint(lunin,obserr,8,255,levs,oestr) nread=nread+levs if(uvob)then nread=nread+levs else if(sstob)then sstdat=bmiss call ufbint(lunin,sstdat,8,1,levs,sststr) else if(metarcldobs) then metarcld=bmiss metarwth=bmiss metarvis=bmiss call ufbint(lunin,metarcld,2,10,metarcldlevs,metarcldstr) call ufbint(lunin,metarwth,1,10,metarwthlevs,metarwthstr) call ufbint(lunin,metarvis,2,1,iret,metarvisstr) if(levs /= 1 ) then write(6,*) 'READ_PREPBUFR: error in Metar observations, levs sould be 1 !!!' call stop2(110) endif else if(geosctpobs) then geoscld=bmiss call ufbint(lunin,geoscld,4,1,levs,geoscldstr) else if (visob) then metarwth=bmiss call ufbint(lunin,metarwth,1,10,metarwthlevs,metarwthstr) endif ! Check for valid reported pressure (POB). Set POB=bmiss if POB=bmiss) exit ! end of stack end do end do else !peel back events to store sensible temp in case temp is virtual call ufbevn(lunin,tobaux,2,255,20,levs,'TOB TQM') do k=1,levs tvflg(k)=one ! initialize as sensible do j=1,20 if (tpc(k,j)==vtcd) then obsdat(3,k)=tobaux(1,k,j+1) qcmark(3,k)=min(tobaux(2,k,j+1),qcmark_huge) tqm(k)=nint(qcmark(3,k)) end if if (tpc(k,j)>=bmiss) exit ! end of stack end do end do end if end if stnelev=hdr(6) ithin=ithin_conv(nc) ithinp = ithin > 0 .and. pflag /= 0 if(.not. driftl .and. (levs > 1 .or. ithinp))then ! Interpolate guess pressure profile to observation location klon1= int(dlon); klat1= int(dlat) dx = dlon-klon1; dy = dlat-klat1 dx1 = one-dx; dy1 = one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if (klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if (klonp1==nlon+1) klonp1=1 do kk=1,nsig presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) + & w10*prsl_full(klatp1,klon1 ,kk) + & w01*prsl_full(klat1 ,klonp1,kk) + & w11*prsl_full(klatp1,klonp1,kk) end do ! Compute depth of guess pressure layersat observation location if (.not.twodvar_regional .and. levs > 1) then do kk=1,nsig-1 dpres(kk)=presl(kk)-presl(kk+1) end do endif end if LOOP_K_LEVS: do k=1,levs icntpnt=icntpnt+1 ! Extract quality marks if(tob)then qm=tqm(k) else if(uvob) then qm=wqm(k) else if(spdob) then qm=wqm(k) else if(psob) then qm=pqm(k) else if(qob) then if(obsdat(2,k) > r0_01_bmiss)cycle loop_k_levs qm=qqm(k) else if(pwob) then pwq=nint(qcmark(7,k)) qm=pwq else if(sstob) then sstq=100 if (k==1) sstq=nint(min(sstdat(4,k),qcmark_huge)) qm=sstq else if(gustob) then gustqm=0 if (kx==188 .or. kx==288 .or. kx==195 .or. kx==295 ) & call get_gustqm(kx,c_station_id,c_prvstg,c_sprvstg,gustqm) qm=gustqm else if(visob) then visqm=0 ! need to fix this later qm=visqm else if(metarcldobs) then qm=0 else if(geosctpobs) then qm=0 end if ! Check qc marks to see if obs should be processed or skipped if (visob) then if (abs(obsdat(9,k)-bmiss) < r100) then patch_fog=(metarwth(1,1)>= 40.0_r_kind .and. metarwth(1,1)<= 49.0_r_kind) .or. & (metarwth(1,1)>=130.0_r_kind .and. metarwth(1,1)<=135.0_r_kind) .or. & (metarwth(1,1)>=241.0_r_kind .and. metarwth(1,1)<=246.0_r_kind) if (patch_fog) obsdat(9,k)=1000.0_r_kind if (metarwth(1,1)==247.0_r_kind) obsdat(9,k)=75.0_r_kind if (metarwth(1,1)==248.0_r_kind) obsdat(9,k)=45.0_r_kind if (metarwth(1,1)==249.0_r_kind) obsdat(9,k)=15.0_r_kind end if end if if (psob) then cat=nint(min(obsdat(10,k),qcmark_huge)) if ( cat /=0 ) cycle loop_k_levs if ( obsdat(1,k)< r500) qm=100 zqm=nint(qcmark(4,k)) if (zqm>=lim_zqm .and. zqm/=15 .and. zqm/=9) qm=9 endif ! if(convobs .and. pqm(k) >=lim_qm .and. qm/=15 .and. qm/=9 )cycle loop_k_levs ! if(qm >=lim_qm .and. qm /=15 .and. qm /=9)cycle loop_k_levs if(qm > 15 .or. qm < 0) cycle loop_k_levs ! If needed, extract drift information. if(driftl)then if(drfdat(1,k) >= r360)drfdat(1,k)=drfdat(1,k)-r360 if(drfdat(1,k) < zero)drfdat(1,k)=drfdat(1,k)+r360 if(abs(drfdat(2,k)) > r90 .or. drfdat(1,k) > r360 .or. drfdat(1,k) < zero)then drfdat(2,k)=hdr(3) drfdat(1,k)=hdr(2) end if ! Check to ensure header lat and drift lat similar if(abs(drfdat(2,k)-hdr(3)) > r10 .and. & abs(drfdat(1,k)-hdr(2)) > r10)then drfdat(2,k)=hdr(3) drfdat(1,k)=hdr(2) end if ! Check to see if the time is outrageous if so set to header value timeobs = real(real(drfdat(3,k),r_single),r_double) time_drift = timeobs + time_correction if (abs(time_drift-time)>four) time_drift = time ! Check to see if the time is outside range if (l4dvar) then t4dv=toff+time_drift if (t4dvwinlen) then t4dv=toff+timex if (t4dvwinlen) CYCLE LOOP_K_LEVS end if else if(abs(time_drift) > ctwind(nc) .or. abs(time_drift) > twindin)then time_drift=timex if(abs(timex) > ctwind(nc) .or. abs(timex) > twindin) CYCLE LOOP_K_LEVS end if t4dv = toff + time_drift endif dlat_earth = drfdat(2,k) * deg2rad dlon_earth = drfdat(1,k) * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(outside) cycle LOOP_K_LEVS else dlat = dlat_earth dlon = dlon_earth call grdcrd(dlat,1,rlats,nlat,1) call grdcrd(dlon,1,rlons,nlon,1) endif if(levs > 1 .or. ithinp)then ! Interpolate guess pressure profile to observation location klon1= int(dlon); klat1= int(dlat) dx = dlon-klon1; dy = dlat-klat1 dx1 = one-dx; dy1 = one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if (klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if (klonp1==nlon+1) klonp1=1 do kk=1,nsig presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) + & w10*prsl_full(klatp1,klon1 ,kk) + & w01*prsl_full(klat1 ,klonp1,kk) + & w11*prsl_full(klatp1,klonp1,kk) end do ! Compute depth of guess pressure layersat observation location if (.not.twodvar_regional .and. levs > 1) then do kk=1,nsig-1 dpres(kk)=presl(kk)-presl(kk+1) end do endif end if end if ! Special block for data thinning - if requested if (ithin > 0) then ntmp=ndata ! counting moved to map3gridS ! Set data quality index for thinning if (l4dvar) then timedif = zero else timedif=abs(t4dv-toff) endif if(kx == 243 .or. kx == 253 .or. kx ==254) then call ufbint(lunin,satqc,1,1,iret,satqcstr) crit1 = timedif/r6+half + four*(one-satqc(1)/r100)*r3_33 else crit1 = timedif/r6+half endif if (pflag==0) then do kk=1,nsig presl_thin(kk)=presl(kk) end do endif call map3grids(-1,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,& plevs(k),crit1,ithin,ndata,iout,icntpnt,iiout,luse) if (.not. luse) then if(k==levs) then cycle loop_readsb else cycle LOOP_K_LEVS endif endif if(iiout > 0) isort(iiout)=0 if(ndata > ntmp)then nodata=nodata+1 if(uvob)nodata=nodata+1 end if isort(icntpnt)=iout else ndata=ndata+1 nodata=nodata+1 if(uvob)nodata=nodata+1 iout=ndata isort(icntpnt)=iout endif if(ndata > maxobs) then write(6,*)'READ_PREPBUFR: ***WARNING*** ndata > maxobs for ',obstype ndata = maxobs end if ! Set usage variable usage = zero if(icuse(nc) <= 0)usage=100._r_kind if(qm == 15 .or. qm == 12 .or. qm == 9)usage=100._r_kind if(qm >=lim_qm )usage=101._r_kind if(convobs .and. pqm(k) >=lim_qm )usage=102._r_kind if((kx>=192.and.kx<=195) .and. psob )usage=r100 if (gustob .and. abs(obsdat(8,k)-bmiss) < r100) usage=103._r_kind if (visob .and. abs(obsdat(9,k)-bmiss) < r100) usage=103._r_kind if (sfctype) call get_usagerj(kx,obstype,c_station_id,c_prvstg,c_sprvstg, & dlon_earth,dlat_earth,idate,t4dv-toff, & obsdat(5,k),obsdat(6,k),usage) if ((kx>129.and.kx<140).or.(kx>229.and.kx<240) ) then call get_aircraft_usagerj(kx,obstype,c_station_id,usage) endif if(ncnumgrp(nc) > 0 .and. .not.lhilbert )then ! default cross validation on if(mod(ndata+1,ncnumgrp(nc))== ncgroup(nc)-1)usage=ncmiter(nc) end if ! Flag regional MAP wind above 400mb for monitoring if(regional .and. kx==227 .and. obsdat(1,k)<400._r_kind ) usage=r100 ! don't use MESONET psfc obs if 8th character of station id is "x") if( kx==188 .and. psob .and. sidchr(8)=='x' ) usage=r100 ! Get information from surface file necessary for conventional data here call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) if(lhilbert) & call accum_hilbertcurve(usage,c_station_id,c_prvstg,c_sprvstg, & dlat_earth,dlon_earth,dlat,dlon,t4dv,toff,nc,kx,iout) ! Extract pressure level and quality marks dlnpob=log(plevs(k)) ! ln(pressure in cb) ! Set inflate_error logical based on qm flag inflate_error=.false. if (qm==3 .or. qm==7) inflate_error=.true. ! Temperature if(tob) then ppb=obsdat(1,k) call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) toe=obserr(3,k)*errout qtflg=tvflg(k) if (inflate_error) toe=toe*r1_2 if(ppb < r100)toe=toe*r1_2 cdata_all(1,iout)=toe ! temperature error cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=dlnpob ! ln(pressure in cb) if (kx == id_bias_t) then cdata_all(5,iout)=obsdat(3,k)+t0c+conv_bias_t ! temperature ob.+bias else cdata_all(5,iout)=obsdat(3,k)+t0c ! temperature ob. endif cdata_all(6,iout)=rstation_id ! station id cdata_all(7,iout)=t4dv ! time cdata_all(8,iout)=nc ! type cdata_all(9,iout)=qtflg ! qtflg (virtual temperature flag) cdata_all(10,iout)=tqm(k) ! quality mark cdata_all(11,iout)=obserr(3,k) ! original obs error cdata_all(12,iout)=usage ! usage parameter cdata_all(13,iout)=idomsfc ! dominate surface type cdata_all(14,iout)=tsavg ! skin temperature cdata_all(15,iout)=ff10 ! 10 meter wind factor cdata_all(16,iout)=sfcr ! surface roughness cdata_all(17,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(18,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(19,iout)=stnelev ! station elevation (m) cdata_all(20,iout)=obsdat(4,k) ! observation height (m) cdata_all(21,iout)=zz ! terrain height at ob location cdata_all(22,iout)=r_prvstg(1,1) ! provider name cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name cdata_all(24,iout)=obsdat(10,k) ! cat if(perturb_obs)cdata_all(25,iout)=ran01dom()*perturb_fact ! t perturbation if (twodvar_regional) & call adjust_error(cdata_all(17,iout),cdata_all(18,iout),cdata_all(11,iout),cdata_all(1,iout)) ! Winds else if(uvob) then call errormod(pqm,wqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) woe=obserr(5,k)*errout if (inflate_error) woe=woe*r1_2 if(obsdat(1,k) < r50)woe=woe*r1_2 selev=stnelev oelev=obsdat(4,k) if(kx >= 280 .and. kx < 300 )then oelev=r10+selev if (kx == 280 )then it29=nint(hdr(8)) if(it29 == 522 .or. it29 == 523 .or. it29 == 531)then ! oelev=r20+selev oelev=r20 end if end if if (kx == 282) oelev=r20+selev if (kx == 285 .or. kx == 289 .or. kx == 290) then oelev=selev selev=zero endif else if((kx >= 221 .and. kx <= 229) & .and. selev >= oelev) oelev=r10+selev end if ! Rotate winds to rotated coordinate uob=obsdat(5,k) vob=obsdat(6,k) if(regional)then u0=uob v0=vob call rotate_wind_ll2xy(u0,v0,uob,vob,dlon_earth,dlon,dlat) if(diagnostic_reg) then call rotate_wind_xy2ll(uob,vob,u00,v00,dlon_earth,dlon,dlat) nvtest=nvtest+1 disterr=sqrt((u0-u00)**2+(v0-v00)**2) vdisterrmax=max(vdisterrmax,disterr) end if endif cdata_all(1,iout)=woe ! wind error cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=dlnpob ! ln(pressure in cb) cdata_all(5,iout)=oelev ! height of observation cdata_all(6,iout)=uob ! u obs cdata_all(7,iout)=vob ! v obs cdata_all(8,iout)=rstation_id ! station id cdata_all(9,iout)=t4dv ! time cdata_all(10,iout)=nc ! type cdata_all(11,iout)=selev ! station elevation cdata_all(12,iout)=wqm(k) ! quality mark cdata_all(13,iout)=obserr(5,k) ! original obs error cdata_all(14,iout)=usage ! usage parameter cdata_all(15,iout)=idomsfc ! dominate surface type cdata_all(16,iout)=tsavg ! skin temperature cdata_all(17,iout)=ff10 ! 10 meter wind factor cdata_all(18,iout)=sfcr ! surface roughness cdata_all(19,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(20,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(21,iout)=zz ! terrain height at ob location cdata_all(22,iout)=r_prvstg(1,1) ! provider name cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name cdata_all(24,iout)=obsdat(10,k) ! cat if(perturb_obs)then cdata_all(25,iout)=ran01dom()*perturb_fact ! u perturbation cdata_all(26,iout)=ran01dom()*perturb_fact ! v perturbation endif else if(spdob) then woe=obserr(5,k) if (inflate_error) woe=woe*r1_2 elev=r20 oelev=obsdat(4,k) if(kx == 260 .or. kx == 261) elev = oelev ! Nacelle and tower wind speed cdata_all(1,iout)=woe ! wind error cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=dlnpob ! ln(pressure in cb) cdata_all(5,iout)=obsdat(5,k) ! u obs cdata_all(6,iout)=obsdat(6,k) ! v obs cdata_all(7,iout)=rstation_id ! station id cdata_all(8,iout)=t4dv ! time cdata_all(9,iout)=nc ! type cdata_all(10,iout)=elev ! elevation of observation cdata_all(11,iout)=wqm(k) ! quality mark cdata_all(12,iout)=obserr(5,k) ! original obs error cdata_all(13,iout)=usage ! usage parameter cdata_all(14,iout)=idomsfc ! dominate surface type cdata_all(15,iout)=tsavg ! skin temperature cdata_all(16,iout)=ff10 ! 10 meter wind factor cdata_all(17,iout)=sfcr ! surface roughness cdata_all(18,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(19,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(20,iout)=stnelev ! station elevation (m) cdata_all(21,iout)=zz ! terrain height at ob location cdata_all(22,iout)=r_prvstg(1,1) ! provider name cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name ! Surface pressure else if(psob) then poe=obserr(1,k)*one_tenth ! convert from mb to cb if (inflate_error) poe=poe*r1_2 cdata_all(1,iout)=poe ! surface pressure error (cb) cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=exp(dlnpob) ! pressure (in cb) cdata_all(5,iout)=obsdat(4,k) ! surface height cdata_all(6,iout)=obsdat(3,k)+t0c ! surface temperature cdata_all(7,iout)=rstation_id ! station id cdata_all(8,iout)=t4dv ! time cdata_all(9,iout)=nc ! type cdata_all(10,iout)=pqm(k) ! quality mark cdata_all(11,iout)=obserr(1,k)*one_tenth ! original obs error (cb) cdata_all(12,iout)=usage ! usage parameter cdata_all(13,iout)=idomsfc ! dominate surface type cdata_all(14,iout)=tsavg ! skin temperature cdata_all(15,iout)=ff10 ! 10 meter wind factor cdata_all(16,iout)=sfcr ! surface roughness cdata_all(17,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(18,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(19,iout)=stnelev ! station elevation (m) cdata_all(20,iout)=zz ! terrain height at ob location cdata_all(21,iout)=r_prvstg(1,1) ! provider name cdata_all(22,iout)=r_sprvstg(1,1) ! subprovider name if(perturb_obs)cdata_all(23,iout)=ran01dom()*perturb_fact ! ps perturbation if (twodvar_regional) & call adjust_error(cdata_all(17,iout),cdata_all(18,iout),cdata_all(11,iout),cdata_all(1,iout)) ! Specific humidity else if(qob) then qmaxerr=emerr call errormod(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) qoe=obserr(2,k)*one_tenth*errout if (inflate_error) then qmaxerr=emerr*r0_7; qoe=qoe*r1_2 end if qobcon=obsdat(2,k)*convert tdry=r999 if (tqm(k)= 280 .and. kx < 300).or.(kx >= 180 .and. kx < 200))then oelev=r10+selev if ((kx==280).or.(kx==180)) oelev=r20+selev if ((kx==282).or.(kx==182)) oelev=r20+selev if ((kx==285).or.(kx==185)) then oelev=selev selev=zero end if if ((kx==188).or.(kx==288) .or.(kx==195) .or.(kx==295)) then ! gustoe=2.5 gustoe=1.0 windcorr=abs(obsdat(5,k))<1.0 .and. abs(obsdat(6,k))<1.0 .and. obsdat(8,k)>10.0 if (windcorr) gustoe=gustoe*1.5_r_kind if (abs(obsdat(8,k)-sqrt(obsdat(5,k)**2+obsdat(6,k)**2))<1.5) then gustoe=gustoe*1.5_r_kind end if end if end if if (inflate_error) gustoe=gustoe*1.5_r_kind cdata_all(1,iout)=gustoe ! wind gusts error (cb) cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=dlnpob ! ln(pressure in cb) cdata_all(5,iout)=oelev ! observation height cdata_all(6,iout)=obsdat(8,k) ! wind gusts obs cdata_all(7,iout)=rstation_id ! station id cdata_all(8,iout)=t4dv ! time cdata_all(9,iout)=nc ! type cdata_all(10,iout)=gustoe*three ! max error cdata_all(11,iout)=gustqm ! quality mark cdata_all(12,iout)=usage ! usage parameter cdata_all(13,iout)=idomsfc ! dominate surface type cdata_all(14,iout)=tsavg ! skin temperature cdata_all(15,iout)=ff10 ! 10 meter wind factor cdata_all(16,iout)=sfcr ! surface roughness cdata_all(17,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(18,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(19,iout)=selev ! station elevation (m) cdata_all(20,iout)=r_prvstg(1,1) ! provider name cdata_all(21,iout)=r_sprvstg(1,1) ! subprovider name ! Visibility else if(visob) then visoe=4000.0 ! temporarily if ((kx==283).or.(kx==183)) visoe=4500.0 if (inflate_error) visoe=visoe*r1_2 cdata_all(1,iout)=visoe ! visibility error (cb) cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=obsdat(9,k) ! visibility obs cdata_all(5,iout)=rstation_id ! station id cdata_all(6,iout)=t4dv ! time cdata_all(7,iout)=nc ! type cdata_all(8,iout)=visoe*three ! max error cdata_all(9,iout)=visqm ! quality mark cdata_all(10,iout)=usage ! usage parameter cdata_all(11,iout)=idomsfc ! dominate surface type cdata_all(12,iout)=tsavg ! skin temperature cdata_all(13,iout)=ff10 ! 10 meter wind factor cdata_all(14,iout)=sfcr ! surface roughness cdata_all(15,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(16,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(17,iout)=stnelev ! station elevation (m) cdata_all(18,iout)=obsdat(4,k) ! observation height (m) cdata_all(19,iout)=zz ! terrain height at ob location cdata_all(20,iout)=r_prvstg(1,1) ! provider name cdata_all(21,iout)=r_sprvstg(1,1) ! subprovider name ! METAR cloud observation else if(metarcldobs) then cdata_all(1,iout)=rstation_id ! station ID cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=stnelev ! station elevation if(metarvis(1,1) < r0_1_bmiss) then cdata_all(5,iout)=metarvis(1,1) ! visibility (m) else cdata_all(5,iout) = -99999.0_r_kind endif do kk=1, 6 if(metarcld(1,kk) < r0_1_bmiss) then cdata_all(5+kk,iout) =metarcld(1,kk) ! cloud amount else cdata_all(5+kk,iout) = -99999.0_r_kind endif if(metarcld(2,kk) < r0_1_bmiss) then cdata_all(11+kk,iout)=metarcld(2,kk) ! cloud bottom height (m) else cdata_all(11+kk,iout)= -99999.0_r_kind endif enddo do kk=1, 3 if(metarwth(1,kk) < r0_1_bmiss) then cdata_all(17+kk,iout)=metarwth(1,kk) ! weather else cdata_all(17+kk,iout)= -99999.0_r_kind endif enddo cdata_all(21,iout)=timeobs ! time observation cdata_all(22,iout)=usage cdata_all(23,iout)=0.0_r_kind ! reserved for distance between obs and grid ! Calculate dewpoint depression from surface obs, to be used later ! with haze and ceiling logic to exclude dust-caused ceiling obs ! from cloud analysis if(metarvis(2,1) < 1.e10_r_kind) then cdata_all(24,iout)=obsdat(3,1)-metarvis(2,1) ! temperature - dew point else cdata_all(24,iout)=-99999.0_r_kind ! temperature - dew point endif ! NESDIS cloud products else if(geosctpobs) then cdata_all(1,iout)=rstation_id ! station ID cdata_all(2,iout)=dlon ! grid relative longitude cdata_all(3,iout)=dlat ! grid relative latitude cdata_all(4,iout)=geoscld(1,k)/100.0_r_kind ! cloud top pressure (pa) cdata_all(5,iout)=geoscld(2,k) ! cloud cover cdata_all(6,iout)=geoscld(3,k) ! Cloud top temperature (K) cdata_all(7,iout)=timeobs ! time cdata_all(8,iout)=usage end if ! ! End k loop over levs end do LOOP_K_LEVS end do loop_readsb ! ! End of bufr read loop enddo loop_msg ! Close unit to bufr file call closbf(lunin) ! Deallocate arrays used for thinning data if (.not.use_all) then deallocate(presl_thin) call del3grids endif ! Normal exit enddo loop_convinfo! loops over convinfo entry matches deallocate(lmsg) if(lhilbert) & call apply_hilbertcurve(maxobs,cdata_all(10:14,1:maxobs),10,14,& tob,12,uvob,14,spdob,13,psob,12,qob,13,pwob,11,sstob,13, & gustob,12,visob,10) ! Write header record and data to output file for further processing allocate(iloc(ndata)) icount=0 do i=1,maxobs if(isort(i) > 0)then icount=icount+1 iloc(icount)=isort(i) end if end do if(ndata /= icount)then write(6,*) ' PREPBUFR: mix up in read_prepbufr ,ndata,icount ',ndata,icount call stop2(50) end if allocate(cdata_out(nreal,ndata)) do i=1,ndata itx=iloc(i) do k=1,nreal cdata_out(k,i)=cdata_all(k,itx) end do end do deallocate(iloc,isort,cdata_all) ! define a closest METAR cloud observation for each grid point if(metarcldobs .and. ndata > 0) then maxobs=200000 allocate(cdata_all(nreal,maxobs)) call reorg_metar_cloud(cdata_out,nreal,ndata,cdata_all,maxobs,iout) ndata=iout write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((cdata_all(i,j),i=1,nreal),j=1,ndata) deallocate(cdata_all) else write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) cdata_out endif deallocate(cdata_out) call destroy_rjlists call destroy_aircraft_rjlists if (lhilbert) call destroy_hilbertcurve if (twodvar_regional) call destroy_ndfdgrid 900 continue if(diagnostic_reg .and. ntest>0) write(6,*)'READ_PREPBUFR: ',& 'ntest,disterrmax=',ntest,disterrmax if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_PREPBUFR: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax if (ndata == 0) then call closbf(lunin) write(6,*)'READ_PREPBUFR: closbf(',lunin,')' endif close(lunin) close(55) ! End of routine return end subroutine read_prepbufr