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-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-09-29 park - add routine to use the HIlber curve for Cross validation ! 2006-10-02 pondeca,park - add routine to use uselist and rejectlist of Mesonet from Prepbufr ! ! input argument list: ! infile - unit from which to read BUFR data ! obstype - observation type to process ! lunout - unit to which to write data for further processing ! prsl_full- 3d pressure on full domain grid ! ! output argument list: ! nread - number of type "obstype" observations read ! nodata - number of individual "obstype" observations read ! ndata - number of type "obstype" observations retained for further processing ! twindin - input group time window (hours) ! sis - satellite/instrument/sensor indicator ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_single,r_double,i_kind !MPondeca, added r_single use constants, only: zero,one_tenth,one,deg2rad,fv,t0c,half,& three,four,rad2deg,tiny_r_kind,huge_r_kind,huge_i_kind use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,& tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,& rlats,rlons,twodvar_regional use convinfo, only: nconvtype,ctwind,cgross,cermax,cermin,cvar_b,cvar_pg, & ncmiter,ncgroup,ncnumgrp,icuse,ictype,icsubtype,ioctype use obsmod, only: iadate,oberrflg use qcmod, only: errormod,noiqc implicit none ! Declare passed variables character(10),intent(in):: infile,obstype character(20),intent(in):: sis integer(i_kind),intent(in):: lunout integer(i_kind),intent(inout):: nread,ndata,nodata real(r_kind),intent(in):: twindin real(r_kind),dimension(nlat,nlon,nsig),intent(in):: prsl_full ! Declare local parameters real(r_kind),parameter:: r0_001=0.001_r_kind real(r_kind),parameter:: r0_5 = 0.5_r_kind real(r_kind),parameter:: r0_75 = 0.75_r_kind real(r_kind),parameter:: r0_7 = 0.7_r_kind real(r_kind),parameter:: r0_8 = 0.8_r_kind real(r_kind),parameter:: r1_2 = 1.2_r_kind real(r_kind),parameter:: r1_5 = 1.5_r_kind real(r_kind),parameter:: r10 = 10.0_r_kind real(r_kind),parameter:: r20 = 20.0_r_kind real(r_kind),parameter:: r50 = 50.0_r_kind real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r100 = 100.0_r_kind real(r_kind),parameter:: r300 = 300.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind real(r_kind),parameter:: r500 = 500.0_r_kind real(r_kind),parameter:: r999 = 999.0_r_kind real(r_kind),parameter:: r1100= 1100.0_r_kind real(r_kind),parameter:: r2000= 2000.0_r_kind real(r_kind),parameter:: convert= 1.0e-6_r_kind real(r_kind),parameter:: emerr= 0.2_r_kind ! Declare local variables logical tob,qob,uvob,spdob,sstob,pwob,psob logical outside,driftl,convobs,inflate_error logical lgood !MPondeca logical listexist !MPondeca logical listexist_mass_rj !MPondeca logical listexist_wind_rj !MPondeca character(40) drift,hdstr,qcstr,oestr,sststr character(80) obstr character(10) date character(8) subset !park start character(8) prvstr,sprvstr character(8) c_prvstg,c_sprvstg character(8) c_station_id character(16) cprovider(200) character(8) cprovider_rj(10000) character(80) cprovider_mass_rj(10000) !MPondeca character(80) cprovider_wind_rj(10000) !MPondeca character(8) ch8 !MPondeca character(60) cgrid !MPondeca character(80) cstring integer(i_kind) meso_unit logical hilbert_curve !park end integer(i_kind) lunin,i,maxobs,sstmeas,j integer(i_kind) kk,klon1,klat1,klonp1,klatp1 integer(i_kind) jdate,ihh,idd,idate,iret,im,iy,k,levs integer(i_kind) kx,nreal,nchanl,ilat,ilon,ithin integer(i_kind) cat,zqm,pwq,sstq,qm,lim_qm,lim_zqm integer(i_kind) lim_tqm,lim_qqm integer(i_kind),dimension(255):: pqm,qqm,tqm,wqm real(r_kind) time,timex,time_drift real(r_kind) qtflg,tdry,rmesh,ediff,usage real(r_kind) u0,v0,uob,vob,dx,dy,dx1,dy1,w00,w10,w01,w11 real(r_kind) qoe,qobcon,pwoe,pwmerr,dlnpob,ppb,poe,qmaxerr,wtype real(r_kind) toe,woe,errout,oelev,dlat,dlon,sstoe,dlat_earth,dlon_earth real(r_kind) selev,elev,stnelev real(r_kind),dimension(nsig):: presl real(r_kind),dimension(nsig-1):: dpres real(r_kind),allocatable,dimension(:,:):: cdata_all real(r_kind) usage_rj !MPondeca !park start real(r_kind),allocatable,dimension(:):: hil_dlon real(r_kind),allocatable,dimension(:):: hil_dlat integer(i_kind),allocatable,dimension(:):: hil_ikx integer(i_kind),allocatable,dimension(:):: hil_kx integer(i_kind),allocatable,dimension(:):: hil_i integer(i_kind),allocatable,dimension(:):: test_set integer(i_kind) ncross !park end real(r_double) rstation_id real(r_double),dimension(7):: hdr real(r_double),dimension(8,255):: drfdat,qcmark,obserr real(r_double),dimension(9,255):: obsdat real(r_double),dimension(8,1):: sstdat real(r_double),dimension(1,1):: r_prvstg,r_sprvstg !MPondeca real(r_kind),dimension(255)::plevs real(r_kind),dimension(1000):: twind,gross,ermax,ermin,var_b,var_pg integer(i_kind),dimension(1000):: itype,iuse,numgrp,ngroup,nmiter real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00 real(r_kind) vdisterr,vdisterrmax,u00,v00 integer(i_kind) ntest,nvtest,iosub,ixsub,isubsub,iosubsub,iobsub real(4),allocatable::etabl(:,:,:) integer(i_kind) ietabl,lcount,itypex,ierrcode,numbcast,kl,k1,k2 integer(i_kind) l,m,ikx integer(i_kind) nprov! MPondeca integer(i_kind) mblack! MPondeca integer(i_kind) nblack! MPondeca integer(i_kind) nx,ny! MPondeca real(r_kind) del,terrmin,werrmin,perrmin,qerrmin,pwerrmin real(r_single),allocatable::slmask(:,:) !MPondeca real(r_single) ds !MPondeca ! include '../../fix.rtma/param.incl' !MPondeca include 'param.incl' !MPondeca namelist/gridname/cgrid !MPondeca equivalence(r_prvstg(1,1),c_prvstg) !MPondeca equivalence(r_sprvstg(1,1),c_sprvstg) !MPondeca equivalence(rstation_id,c_station_id) !Park data hdstr /'SID XOB YOB DHR TYP ELV SAID '/ data obstr /'POB QOB TOB ZOB UOB VOB PWO CAT PRSS' / data drift /'XDR YDR HRDR '/ data sststr /'MSST DBSS SST1 SSTQM SSTOE '/ data qcstr /'PQM QQM TQM ZQM WQM NUL PWQ '/ data oestr /'POE QOE TOE NUL WOE NUL PWE '/ data prvstr /'PRVSTG'/ !MPondeca data sprvstr /'SPRVSTG'/ !MPondeca data lunin / 10 / data ithin / -9 / data rmesh / -99.999 / data meso_unit / 20 / !park data cgrid/'conus'/ !MPondeca !************************************************************************** ! Initialize variables disterrmax=zero vdisterrmax=zero ntest=0 nvtest=0 maxobs=2e6 nread=0 nchanl=0 ilon=2 ilat=3 nreal=0 tob = obstype == 't' uvob = obstype == 'uv' spdob = obstype == 'spd' psob = obstype == 'ps' qob = obstype == 'q' pwob = obstype == 'pw' sstob = obstype == 'sst' convobs = tob .or. uvob .or. spdob .or. qob if(tob)then nreal=16 else if(uvob) then nreal=16 else if(spdob) then nreal=16 else if(psob) then nreal=15 else if(qob) then nreal=17 else if(pwob) then nreal=16 else if(sstob) then nreal=16 else write(6,*) ' illegal obs type in READ_PREPBUFR ' call stop2(94) end if if(oberrflg)then allocate(etabl(300,33,6)) ietabl=19 open(ietabl,file='errtable',form='formatted') rewind ietabl etabl=1.e9_r_kind lcount=0 do l=1,300 read(ietabl,100,end=120,err=120)itypex 100 format(1x,i3) lcount=lcount+1 do k=1,33 read(ietabl,110)(etabl(itypex,k,m),m=1,6) 110 format(1x,6e12.5) end do end do 120 continue if(lcount.le.0) then write(6,*)'READ_PREPBUFR: ***WARNING*** obs error table not available to 3dvar.' oberrflg=.false. deallocate(etabl) else write(6,*)'READ_PREPBUFR: using observation errors from user provided table' endif close(ietabl) ! Set lower limits for observation errors terrmin=r0_5 werrmin=one perrmin=r0_5 qerrmin=one_tenth pwerrmin=one endif !park start hilbert_curve=.true. !hardwired for now ! Read mesonet provider names from the uselist or rejectlist if it exists. inquire(file='mesonetuselist',exist=listexist) if(listexist) then open (meso_unit,file='mesonetuselist',form='formatted') nprov=0 do m=1,3 read(meso_unit,*,end=135) cstring enddo 130 continue nprov=nprov+1 read(meso_unit,*,end=135) cprovider(nprov) goto 130 135 continue nprov=nprov-1 print*,'mesonetuselist: nprov=',nprov endif close(meso_unit) !MPondeca inquire(file='mass_rejectlist',exist=listexist_mass_rj) if(listexist_mass_rj) then open (meso_unit,file='mass_rejectlist',form='formatted') mblack=0 do m=1,3 read(meso_unit,*,end=155) cstring enddo 150 continue mblack=mblack+1 read(meso_unit,*,end=155) cprovider_mass_rj(mblack) goto 150 155 continue mblack=mblack-1 print*,'mass rejectlist mblack=',mblack endif close(meso_unit) !MPondeca inquire(file='wind_rejectlist',exist=listexist_wind_rj) if(listexist_wind_rj) then open (meso_unit,file='wind_rejectlist',form='formatted') nblack=0 do m=1,3 read(meso_unit,*,end=157) cstring enddo 152 continue nblack=nblack+1 read(meso_unit,*,end=157) cprovider_wind_rj(nblack) goto 152 157 continue nblack=nblack-1 print*,'wind rejectlist nblack=',nblack endif close(meso_unit) !MPondeca !park end open (55,file='gridname_input',form='formatted') !MPondeca read(55,gridname) !MPondeca close(55) !MPondeca print*,'in read_prepbufr: cgrid is =',trim(cgrid) !MPondeca call domain_dims(cgrid,nx,ny,ds) !MPondeca allocate(slmask(nx,ny)) open (89,file='rtma_slmask.dat',form='unformatted') !MPondeca read(89) slmask !MPondeca close(89) !MPondeca ! Open, then read date from bufr data open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) call readmg(lunin,subset,idate,iret) if(iret/=0) goto 1020 ! On temperature task, write out date information. If date in prepbufr ! file does not agree with analysis date, print message and stop program ! execution. if(tob) then write(date,'( i10)') idate read (date,'(i4,3i2)') iy,im,idd,ihh write(6,*)'READ_PREPBUFR: bufr file date is ',iy,im,idd,ihh if(iy/=iadate(1).or.im/=iadate(2).or.idd/=iadate(3).or.& ihh/=iadate(4)) then write(6,*)'***READ_PREPBUFR ERROR*** incompatable analysis ',& 'and observation date/time' write(6,*)' year anal/obs ',iadate(1),iy write(6,*)' month anal/obs ',iadate(2),im write(6,*)' day anal/obs ',iadate(3),idd write(6,*)' hour anal/obs ',iadate(4),ihh call stop2(94) end if end if allocate(cdata_all(nreal,maxobs)) !park start !if hilbert curve for cross validation is on, allocate some variables if(hilbert_curve) then allocate(hil_dlon(maxobs)) allocate(hil_dlat(maxobs)) allocate(hil_ikx(maxobs)) allocate(hil_kx(maxobs)) allocate(hil_i(maxobs)) ncross=0 endif !park end ! Big loop over buffer file 10 call readsb(lunin,iret) if(iret/=0) then call readmg(lunin,subset,jdate,iret) if(iret/=0) go to 1000 go to 10 end if ! Extract type, date, and location information call ufbint(lunin,hdr,7,1,iret,hdstr) rstation_id=hdr(1) 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=acosd(cdist) disterrmax=max(disterrmax,disterr) end if if(outside) go to 10 ! 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 time=hdr(4) kx=hdr(5) stnelev=hdr(6) iobsub = 0 ! temporary until put in bufr file if(kx == 243 .or. kx == 253 .or. kx == 254) then iobsub = hdr(7) end if ikx=0 do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. kx == ictype(i) & .and. abs(icuse(i))== 1)then if(icsubtype(i) == iobsub)ikx=i end if end do if(ikx == 0)then do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. kx == ictype(i) & .and. abs(icuse(i))== 1)then ixsub=icsubtype(i)/10 iosub=iobsub/10 isubsub=icsubtype(i)-ixsub*10 if(ixsub == iosub .and. isubsub == 0)ikx=i end if end do end if if(ikx == 0)then do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. kx == ictype(i) & .and. abs(icuse(i))== 1)then if(icsubtype(i) == 0)ikx=i end if end do end if if(ikx == 0) go to 10 ! not ob type used ! 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 if ( twodvar_regional .and. & (kx<180 .or. kx>289 .or. (kx>189 .and. kx<280)) ) go to 10 !park start !-------------------------------------------------------------------------- ! For the mesonet winds, check to see if mesonet provider is in the uselist ! or obs station is in the rejectlist. If not, read next observation report from bufr file. usage_rj=0. if (kx.lt.190 .and. listexist_mass_rj) then do m=1,mblack ch8(1:8)=cprovider_mass_rj(m)(1:8) if (trim(c_station_id) == trim(ch8)) usage_rj=5000. enddo endif if (kx.gt.190 .and. listexist_wind_rj) then do m=1,nblack ch8(1:8)=cprovider_wind_rj(m)(1:8) if (trim(c_station_id) == trim(ch8)) usage_rj=5000. enddo endif if (kx.eq.288 .and. listexist) then call ufbint(lunin,r_prvstg,1,1,iret,prvstr) call ufbint(lunin,r_sprvstg,1,1,iret,sprvstr) if (.not.listexist_wind_rj) usage_rj=5000. do m=1,nprov if (trim(c_prvstg//c_sprvstg) == trim(cprovider(m))) usage_rj=0. enddo endif !park end ! Balloon drift information available for these data driftl=kx==120.or.kx==220.or.kx==221 if((real(abs(time)) > real(ctwind(ikx)) .or. real(abs(time)) > real(twindin)) & .and. .not. driftl)go to 10 ! outside time window timex=time ! Extract data information on levels call ufbint(lunin,obsdat,9,255,levs,obstr) call ufbint(lunin,qcmark,8,255,levs,qcstr) call ufbint(lunin,obserr,8,255,levs,oestr) nread=nread+levs if(uvob)nread=nread+levs sstdat=1.e11 if(sstob)call ufbint(lunin,sstdat,8,1,levs,sststr) ! If available, get obs errors from error table if(oberrflg)then do k=1,levs itypex=kx ppb=obsdat(1,k) if(kx==153)ppb=obsdat(9,k)*0.01 ppb=max(zero,min(ppb,r2000)) if(ppb.ge.etabl(itypex,1,1)) k1=1 do kl=1,32 if(ppb.ge.etabl(itypex,kl+1,1).and.ppb.le.etabl(itypex,kl,1)) k1=kl end do if(ppb.le.etabl(itypex,33,1)) k1=5 k2=k1+1 ediff = etabl(itypex,k2,1)-etabl(itypex,k1,1) if (abs(ediff) > tiny_r_kind) then del = (ppb-etabl(itypex,k1,1))/ediff else del = huge_r_kind endif del=max(zero,min(del,one)) obserr(3,k)=(one-del)*etabl(itypex,k1,2)+del*etabl(itypex,k2,2) obserr(2,k)=(one-del)*etabl(itypex,k1,3)+del*etabl(itypex,k2,3) obserr(5,k)=(one-del)*etabl(itypex,k1,4)+del*etabl(itypex,k2,4) obserr(1,k)=(one-del)*etabl(itypex,k1,5)+del*etabl(itypex,k2,5) obserr(7,k)=(one-del)*etabl(itypex,k1,6)+del*etabl(itypex,k2,6) obserr(3,k)=max(obserr(3,k),terrmin) obserr(2,k)=max(obserr(2,k),qerrmin) obserr(5,k)=max(obserr(5,k),werrmin) obserr(1,k)=max(obserr(1,k),perrmin) obserr(7,k)=max(obserr(7,k),pwerrmin) enddo endif ! If data with drift position, get drift information if(driftl)call ufbint(lunin,drfdat,8,255,iret,drift) ! Loop over levels do k=1,levs do i=1,8 qcmark(i,k) = min(qcmark(i,k),huge_i_kind) end do plevs(k)=one_tenth*obsdat(1,k) ! convert mb to cb pqm(k)=nint(qcmark(1,k)) qqm(k)=nint(qcmark(2,k)) tqm(k)=nint(qcmark(3,k)) wqm(k)=nint(qcmark(5,k)) end do LOOP_K_LEVS: do k=1,levs ! 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) > 1.0e9)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),huge_i_kind)) qm=sstq end if ! Set inflate_error logical and qc limits based on noiqc flag inflate_error=.false. if (noiqc) then lim_qm=8 if (qm==3 .or. qm==7) inflate_error=.true. if (psob) lim_zqm=7 if (qob) lim_tqm=7 if (tob) lim_qqm=8 else lim_qm=4 if (qm==3) inflate_error=.true. if (psob) lim_zqm=4 if (qob) lim_tqm=4 if (tob) lim_qqm=4 endif ! Check qc marks to see if obs should be processed or skipped if (psob) then zqm=nint(qcmark(4,k)) cat=nint(min(obsdat(8,k),huge_i_kind)) if (zqm > lim_zqm .or. cat /=0 .or. obsdat(1,k) < r500)cycle loop_k_levs endif if(convobs .and. pqm(k) >=lim_qm)cycle loop_k_levs if(qm >=lim_qm .and. qm /=15 .and. qm /=9)cycle loop_k_levs ! Set usage variable usage = 0. if (usage_rj > 0.) usage=5000. if(icuse(ikx) < 0)usage=100. if(qm == 15 .or. qm == 9)usage=100. lgood=.true. if (usage==100. .or. usage==5000.) lgood=.false. if(ncnumgrp(ikx) > 0 .and.lgood .and. .not.hilbert_curve)then ! default cross validation on(park) if(mod(ndata+1,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if ! 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 time_drift = drfdat(3,k) if (abs(time_drift)>four) time_drift = time ! Check to see if the time is outside range if(abs(time_drift) > ctwind(ikx) .or. abs(time_drift) > twindin)then time_drift=timex if(abs(timex) > ctwind(ikx) .or. abs(timex) > twindin) CYCLE LOOP_K_LEVS end if time=time_drift 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 end if ! 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) then do kk=1,nsig-1 dpres(kk)=presl(kk)-presl(kk+1) end do endif ! Extract pressure level and quality marks dlnpob=log(plevs(k)) ! ln(pressure in cb) ndata=ndata+1 nodata=nodata+1 if(uvob)nodata=nodata+1 if(ndata > maxobs) then write(6,*)'READ_PREPBUFR: ***WARNING*** ndata > maxobs for ',obstype ndata = maxobs end if !park start !---------------------------------------------------------------- ! prepare for "hilbert-curve" type cross validation !---------------------------------------------------------------- if(hilbert_curve .and. ncnumgrp(ikx) > 0 .and.lgood) then ncross=ncross+1 hil_dlat(ncross)=dlat/nlat hil_dlon(ncross)=dlon/nlon hil_ikx(ncross)=ikx hil_kx(ncross)=kx hil_i(ncross)=ndata ! write(6,*) 'CHECK0:',i,ndata,ncross,hil_ikx(ncross),hil_kx(ncross),ncnumgrp(hil_ikx(ncross)),ncgroup(hil_ikx(ncross)) endif !park end ! Temperature if(tob) then ppb=obsdat(1,k) qtflg=zero if(qqm(k)>=lim_qqm.and.ppb>r300) qtflg=one if(kx==130 .or. kx == 131 .or. kx == 133) qtflg=one call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig) toe=obserr(3,k)*errout if (inflate_error) toe=toe*r1_2 if(ppb < r100)toe=toe*r1_2 cdata_all(1,ndata)=toe ! temperature error cdata_all(2,ndata)=dlon ! grid relative longitude cdata_all(3,ndata)=dlat ! grid relative latitude cdata_all(4,ndata)=dlnpob ! ln(pressure in cb) cdata_all(5,ndata)=obsdat(3,k)+t0c ! temperature ob. cdata_all(6,ndata)=rstation_id ! station id cdata_all(7,ndata)=time ! time cdata_all(8,ndata)=ikx ! type cdata_all(9,ndata)=qtflg ! qtflg cdata_all(10,ndata)=tqm(k) ! quality mark cdata_all(11,ndata)=obserr(3,k) ! original obs error cdata_all(12,ndata)=usage ! usage parameter cdata_all(13,ndata)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(14,ndata)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(15,ndata)=stnelev ! station elevation (m) cdata_all(16,ndata)=obsdat(4,k) ! observation height (m) ! PRINT*,'k,ikx,obserr(3,k),dlon_earth=',k,ikx,obserr(3,k),dlon_earth*rad2deg call adjust_error(cdata_all(13,ndata),cdata_all(14,ndata), & cdata_all(11,ndata),cdata_all(1,ndata),ikx,cgrid,slmask,nx,ny) !MPondeca ! Winds else if(uvob) then call errormod(pqm,wqm,levs,plevs,errout,k,presl,dpres,nsig) 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(selev == oelev)oelev=r10+selev if(kx >= 280 .and. kx < 290 )then if (kx == 280) oelev=r20+selev if (kx == 281) oelev=r10+selev if (kx == 282) oelev=r20+selev if (kx == 284) oelev=r10+selev if (kx == 285) oelev=selev if (kx == 285) selev=zero if (kx == 286) oelev=r10+selev if (kx == 287) oelev=r10+selev if (kx == 288) oelev=r10+selev if (kx == 289) 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,dlat_earth,dlon,dlat) if(diagnostic_reg) then call rotate_wind_xy2ll(uob,vob,u00,v00,dlon_earth,dlat_earth,dlon,dlat) nvtest=nvtest+1 disterr=sqrt((u0-u00)**2+(v0-v00)**2) vdisterrmax=max(vdisterrmax,disterr) end if endif cdata_all(1,ndata)=woe ! wind error cdata_all(2,ndata)=dlon ! grid relative longitude cdata_all(3,ndata)=dlat ! grid relative latitude cdata_all(4,ndata)=dlnpob ! ln(pressure in cb) cdata_all(5,ndata)=oelev ! height of observation cdata_all(6,ndata)=uob ! u obs cdata_all(7,ndata)=vob ! v obs cdata_all(8,ndata)=rstation_id ! station id cdata_all(9,ndata)=time ! time cdata_all(10,ndata)=ikx ! type cdata_all(11,ndata)=selev ! station elevation cdata_all(12,ndata)=wqm(k) ! quality mark cdata_all(13,ndata)=obserr(5,k) ! original obs error cdata_all(14,ndata)=usage ! usage parameter cdata_all(15,ndata)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(16,ndata)=dlat_earth*rad2deg ! earth relative latitude (degrees) else if(spdob) then woe=obserr(5,k) if (inflate_error) woe=woe*r1_2 elev=r20 cdata_all(1,ndata)=woe ! wind error cdata_all(2,ndata)=dlon ! grid relative longitude cdata_all(3,ndata)=dlat ! grid relative latitude cdata_all(4,ndata)=dlnpob ! ln(pressure in cb) cdata_all(5,ndata)=obsdat(5,k) ! u obs cdata_all(6,ndata)=obsdat(6,k) ! v obs cdata_all(7,ndata)=rstation_id ! station id cdata_all(8,ndata)=time ! time cdata_all(9,ndata)=ikx ! type cdata_all(10,ndata)=elev ! elevation of observation cdata_all(11,ndata)=wqm(k) ! quality mark cdata_all(12,ndata)=obserr(5,k) ! original obs error cdata_all(13,ndata)=usage ! usage parameter cdata_all(14,ndata)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(15,ndata)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(16,ndata)=stnelev ! station elevation (m) ! 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,ndata)=poe ! surface pressure error (cb) cdata_all(2,ndata)=dlon ! grid relative longitude cdata_all(3,ndata)=dlat ! grid relative latitude cdata_all(4,ndata)=exp(dlnpob) ! pressure (in cb) cdata_all(5,ndata)=obsdat(4,k) ! surface height cdata_all(6,ndata)=obsdat(3,k)+t0c ! surface temperature cdata_all(7,ndata)=rstation_id ! station id cdata_all(8,ndata)=time ! time cdata_all(9,ndata)=ikx ! type cdata_all(10,ndata)=pqm(k) ! quality mark cdata_all(11,ndata)=obserr(1,k)*one_tenth ! original obs error (cb) cdata_all(12,ndata)=usage ! usage parameter cdata_all(13,ndata)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata_all(14,ndata)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata_all(15,ndata)=stnelev ! station elevation (m) ! Specific humidity else if(qob) then qmaxerr=emerr call errormod(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig) 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)0) write(6,*)'READ_PREPBUFR: ',& 'ntest,disterrmax=',ntest,disterrmax if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_PREPBUFR: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax deallocate(slmask) !MPondeca ! End of routine return end subroutine read_prepbufr !************************************************************ subroutine adjust_error(alon,alat,oberr,oberr2,iobtype,cgrid,slmask,nx,ny) use kinds, only: r_kind,r_single,r_double,i_kind use constants, only: zero implicit none character(60) cgrid real(r_kind), parameter::rfact=5.0 integer(i_kind), parameter::idelta=3 integer(i_kind), parameter::jdelta=3 integer(i_kind), intent(in)::nx,ny,iobtype real (r_single), intent(in)::slmask(nx,ny) real (r_kind), intent(in)::alon,alat real (r_kind), intent(inout)::oberr,oberr2 real(8) da8,alat18,elon18,elonv8,alatan8 !keep "8" rather than "r_kind" so real(8) rlat8,rlon8,xx8,yy8 !to avoid surprises from subroutine !outside GSI integer(i_kind) i,j,istart,jstart,is,ie,js,je real (r_single) rsign1,rsign2 logical lcase1 !exception for Islands off of Southern California rlon8=alon rlat8=alat if (rlon8.gt.180._8) rlon8=rlon8-360._8 lcase1=(rlon8.ge.-122._8 .and. rlon8.le.-117._8 & .and. rlat8.ge.32._8 .and. rlat8.le.35._8) if (trim(cgrid)=='conus') then if(lcase1) return endif call proj_info(cgrid,da8,alat18,elon18,elonv8,alatan8) if (rlon8 .lt. zero) rlon8=rlon8+360._r_kind call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) istart=floor(xx8) jstart=floor(yy8) print*,'da8,alat18,elon18,elonv8,alatan8,rlat8,rlon8,cgrid,xx8,yy8,istart,jstart=',& da8,alat18,elon18,elonv8,alatan8,rlat8,rlon8,cgrid,xx8,yy8,istart,jstart ! print*,'alon,rlon8,alat,rlat8,istart,jstart=',alon,rlon8,alat,rlat8,istart,jstart ! print*,'slmask,min,max=',minval(slmask),maxval(slmask) ! print*,'before, oberr,oberr2=',oberr,oberr2 is=max(1,(istart-idelta)) ie=min((istart+idelta),nx) js=max(1,(jstart-jdelta)) je=min((jstart+jdelta),ny) if (slmask(is,js) .le. 0.5_r_single) rsign1=-1. if (slmask(is,js) .gt. 0.5_r_single) rsign1=+1. do j=js,je do i=is,ie if (slmask(i,j) .le. 0.5_r_single) rsign2=-1. if (slmask(i,j) .gt. 0.5_r_single) rsign2=+1. if (rsign1*rsign2 .lt.0._r_single) then oberr=oberr*rfact oberr2=oberr2*rfact goto 100 endif enddo enddo 100 continue ! print*,'after, oberr,oberr2=',oberr,oberr2 return end subroutine adjust_error !**************************************************************** subroutine domain_dims(cgrid,nx,ny,ds) implicit none character(60) cgrid ! include '../../fix.rtma/param.incl' include 'param.incl' integer(4) nx,ny real(4) ds if (trim(cgrid) == 'conus') then nx=nx_conus ny=ny_conus ds=da_conus elseif (trim(cgrid) == 'alaska') then nx=nx_alaska ny=ny_alaska ds=da_alaska elseif (trim(cgrid) == 'hawaii') then nx=nx_hawaii ny=ny_hawaii ds=da_hawaii elseif (trim(cgrid) == 'puerto_rico') then nx=nx_prico ny=ny_prico ds=da_prico elseif (trim(cgrid) == 'guam') then nx=nx_guam ny=ny_guam ds=da_guam else print*,'in domain_dims: unknown grid ',cgrid,'...aborting' call abort endif return end !**************************************************************** subroutine proj_info(cgrid,da8,alat18,elon18,elonv8,alatan8) implicit none character(60) cgrid ! include '../../fix.rtma/param.incl' include 'param.incl' real(8) da8,alat18,elon18,elonv8,alatan8 if (trim(cgrid) == 'conus') then alat18=alat1_conus elon18=elon1_conus da8=da_conus elonv8=elonv_conus alatan8=alatan_conus elseif (trim(cgrid) == 'alaska') then alat18=alat1_alaska elon18=elon1_alaska da8=da_alaska elonv8=elonv_alaska alatan8=alatan_alaska elseif (trim(cgrid) == 'hawaii') then alat18=alat1_hawaii elon18=elon1_hawaii da8=da_hawaii elonv8=9999._8 alatan8=9999._8 elseif (trim(cgrid) == 'puerto_rico') then alat18=alat1_prico elon18=elon1_prico da8=da_prico elonv8=9999._8 alatan8=9999._8 elseif (trim(cgrid) == 'guam') then alat18=alat1_guam elon18=elon1_guam da8=da_guam elonv8=9999._8 alatan8=9999._8 else print*,'in proj_info: unknown grid ',cgrid,'...aborting' call abort endif return end !**************************************************************** subroutine latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) implicit none character(60) cgrid real(8) da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,xx8,yy8 if (trim(cgrid)=='conus') then call w3fb11(rlat8,rlon8,alat18,elon18,da8,elonv8,alatan8,xx8,yy8) endif if (trim(cgrid)=='alaska') then call w3fb06(rlat8,rlon8,alat18,elon18,da8,elonv8,xx8,yy8) endif ! Add hawaii, puerto_rico, guam options here return end !****************************************************************