subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& 2,20
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)
!
! 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_double,i_kind
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
character(40) drift,hdstr,qcstr,oestr,sststr
character(80) obstr
character(10) date
character(8) subset
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_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_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
real(r_kind) del,terrmin,werrmin,perrmin,qerrmin,pwerrmin
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 lunin / 10 /
data ithin / -9 /
data rmesh / -99.999 /
!**************************************************************************
! 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
! 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))
! 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
! 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(icuse(ikx) < 0)usage=100.
if(qm == 15 .or. qm == 9)usage=100.
if(ncnumgrp(ikx) > 0 )then ! cross validation on
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
! Temperature
if(tob) then
ppb=obsdat(1,k)
qtflg=zero
if(qqm(k)>=lim_qqm.and.ppb>r300) 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)
! 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)<lim_tqm) tdry=(obsdat(3,k)+t0c)/(one+fv*qobcon)
cdata_all(1,ndata)=qoe ! q 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)=qobcon ! q 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)=qmaxerr ! q max error
cdata_all(10,ndata)=tdry ! dry temperature (obs is tv)
cdata_all(11,ndata)=qqm(k) ! quality mark
cdata_all(12,ndata)=obserr(2,k)*one_tenth ! 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)
cdata_all(17,ndata)=obsdat(4,k) ! observation height (m)
! Total precipitable water (ssm/i)
else if(pwob) then
pwoe=obserr(7,k)
pwmerr=pwoe*three
cdata_all(1,ndata)=pwoe ! pw error
cdata_all(2,ndata)=dlon ! grid relative longitude
cdata_all(3,ndata)=dlat ! grid relative latitude
cdata_all(4,ndata)=obsdat(7,k) ! pw obs
cdata_all(5,ndata)=rstation_id ! station id
cdata_all(6,ndata)=time ! time
cdata_all(7,ndata)=ikx ! type
cdata_all(8,ndata)=pwmerr ! pw max error
cdata_all(9,ndata)=pwq ! quality mark
cdata_all(10,ndata)=obserr(7,k) ! original obs error
cdata_all(11,ndata)=usage ! usage parameter
cdata_all(12,ndata)=dlon_earth*rad2deg ! earth relative longitude (degrees)
cdata_all(13,ndata)=dlat_earth*rad2deg ! earth relative latitude (degrees)
cdata_all(14,ndata)=stnelev ! station elevation (m)
cdata_all(15,ndata)=obsdat(1,k) ! observation pressure (hPa)
cdata_all(16,ndata)=obsdat(4,k) ! observation height (m)
! Conventional sst observations
else if(sstob) then
! Locate the observation on the analysis grid. Get land/sea/ice
! mask at nearest analysis grid points.
sstoe=r0_75
! sstoe=sstdat(5,k)
cdata_all(1,ndata)=sstoe ! sst error
cdata_all(2,ndata)=dlon ! grid relative longitude
cdata_all(3,ndata)=dlat ! grid relative latitude
cdata_all(4,ndata)=sstdat(3,k) ! sst obs
cdata_all(5,ndata)=rstation_id ! station id
cdata_all(6,ndata)=time ! time
cdata_all(7,ndata)=ikx ! type
cdata_all(8,ndata)=sstoe*three ! pw max error
cdata_all(9,ndata)=sstdat(2,k) ! depth of measurement
cdata_all(10,ndata)=sstdat(1,k) ! measurement type
cdata_all(11,ndata)=sstq ! quality mark
cdata_all(12,ndata)=sstdat(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)
! if(kx == 120 .or. kx == 282)then
! write(6,*)'READ_PREPBUFR: kx,rstation_id,sstq=',&
! kx,rstation_id,sstq
! do i=1,10
! write(6,*)'READ_PREPBUFR: i,cdata_all=',i,cdata_all(i,ndata)
! end do
! end if
! Measurement types
! 0 Ship intake
! 1 Bucket
! 2 Hull contact sensor
! 3 Reversing Thermometer
! 4 STD/CTD sensor
! 5 Mechanical BT
! 6 Expendable BT
! 7 Digital BT
! 8 Thermistor chain
! 9 Infra-red scanner
! 10 Micro-wave scanner
! 11-14 Reserved
end if
!
! End k loop over levs
end do LOOP_K_LEVS
!
! End of bufr read loop
go to 10
! Normal exit
1000 continue
! Write header record and data to output file for further processing
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
write(lunout) ((cdata_all(k,i),k=1,nreal),i=1,ndata)
deallocate(cdata_all)
! Close unit to bufr file
1020 continue
if (oberrflg) deallocate(etabl)
call closbf(lunin)
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
! End of routine
return
end subroutine read_prepbufr