subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis,&
prsl_full,nobs)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_fl_hdob read obs from hdob bufr file
! prgmmr: eliu org: np22 date: 2013-02-05
!
! abstract: This routine reads high-density flight-level observations and surface data
! from Stepped Frequency Microwave Radiometer (SFMR). The observation
! types read by this routine include temperature, dew-point temperature,
! wind direction, wind speed, surface wind speed, and total rain rate
!
! When running the gsi in regional mode, the code only
! retains those observations that fall within the regional
! domain
! program history log:
! 2013-02-05 eliu - initial coding
! 2015-02-23 Rancic/Thomas - add thin4d to time window logical
! 2015-02-26 su - add njqc as an option to choose new non linear qc
! 2016-03-15 Su - modified the code so that the program won't stop when no subtype is found in non
! linear qc error table and b table
! 2015-10-01 guo - calc ob location once in deg
!
! 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
! gstime - analysis time in minutes from reference date
! twind - input group time window (hours)
! sis - satellite/instrument/sensor indicator
! 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
! nobs - array of observations on each subdomain for each processor
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_single,r_kind,r_double,i_kind
use constants, only: zero,one_tenth,one,two,ten,deg2rad,t0c,half,&
three,four,rad2deg,tiny_r_kind,huge_r_kind,r0_01,&
r60inv,r10,r100,r2000,hvap,eps,omeps,rv,grav
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, &
icuse,ictype,icsubtype,ioctype, &
ithin_conv,rmesh_conv,pmesh_conv
use obsmod, only: perturb_obs,perturb_fact,ran01dom
use obsmod, only: bmiss
use converr,only: etabl
use converr_ps,only: etabl_ps,isuble_ps,maxsub_ps
use converr_q,only: etabl_q,isuble_q,maxsub_q
use converr_t,only: etabl_t,isuble_t,maxsub_t
use converr_uv,only: etabl_uv,isuble_uv,maxsub_uv
use convb_ps,only: btabl_ps
use convb_q,only: btabl_q
use convb_t,only: btabl_t
use convb_uv,only: btabl_uv
use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,time_4dvar,winlen,thin4d
use qcmod, only: errormod,njqc
use convthin, only: make3grids,map3grids,del3grids,use_all
use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error
use deter_sfc_mod, only: deter_sfc_type,deter_sfc2
use mpimod, only: npe
implicit none
! Declare passed variables
character(len=*), intent(in ) :: infile,obstype
character(len=20),intent(in ) :: sis
integer(i_kind) , intent(in ) :: lunout
integer(i_kind) , dimension(npe), intent(inout) :: nobs
integer(i_kind) , intent(inout) :: nread,ndata,nodata
real(r_kind) , intent(in ) :: twind
real(r_kind) , intent(in ) :: gstime
real(r_kind) , intent(in ) :: prsl_full(nlat,nlon,nsig)
! Declare local variables
! Logical variables
logical :: outside
logical :: inflate_error
logical :: ltob,lqob,luvob,lspdob,lpsob
logical :: luse
! Character variables
character(40) :: timestr,locstr,tmpstr,mststr,wndstr,sfmrstr
character(40) :: psfstr,prsstr,g10str,qcmstr
character(40) :: obs_region
character( 8) :: subset
character( 8) :: c_prvstg,c_sprvstg
character( 8) :: c_station_id
character( 6) :: bulstr1,bulstr2
character( 6) :: obsbul(2,1)
! Integer variables
integer(i_kind), parameter :: mxib = 31
integer(i_kind), parameter :: ietabl= 19
integer(i_kind) :: i,k,kl,k1,k2,j
integer(i_kind) :: lunin
integer(i_kind) :: ireadmg,ireadsb
integer(i_kind) :: idate
integer(i_kind) :: ilat,ilon
integer(i_kind) :: nlv
integer(i_kind) :: nreal,nchanl
integer(i_kind) :: idomsfc,isflg
integer(i_kind) :: ithin,iout
integer(i_kind) :: nc,ncsave
integer(i_kind) :: ntmatch,ntb
integer(i_kind) :: nmsg
integer(i_kind) :: maxobs
integer(i_kind) :: itype,itypey,iecol
integer(i_kind) :: ierr_ps,ierr_q,ierr_t,ierr_uv,ncount_ps,ncount_q,ncount_t,ncount_uv
integer(i_kind) :: qcm,lim_qm
integer(i_kind) :: p_qm,g_qm,t_qm,q_qm,uv_qm,wspd_qm,ps_qm
integer(i_kind) :: ntest,nvtest
! integer(i_kind) :: m,itypex,lcount,iflag
integer(i_kind) :: nlevp ! vertical level for thinning
integer(i_kind) :: pflag
integer(i_kind) :: ntmp,iiout,igood
integer(i_kind) :: kk,klon1,klat1,klonp1,klatp1
integer(i_kind) :: iuse
integer(i_kind) :: nmind
integer(i_kind) :: nib
integer(i_kind) :: ibit(mxib)
integer(i_kind) :: idate5(5)
integer(i_kind), allocatable,dimension(:) :: isort
! Real variables
real(r_kind), parameter :: r0_001 = 0.001_r_kind
real(r_kind), parameter :: r1_2 = 1.2_r_kind
real(r_kind), parameter :: r0_7 = 0.7_r_kind
real(r_kind), parameter :: r6 = 6.0_r_kind
real(r_kind), parameter :: r50 = 50.0_r_kind
real(r_kind), parameter :: r1200 = 1200.0_r_kind
real(r_kind), parameter :: emerr = 0.2_r_kind ! RH
real(r_kind), parameter :: missing = 1.0e+11_r_kind
real(r_kind) :: toff,t4dv
real(r_kind) :: rmesh
real(r_kind) :: usage
real(r_kind) :: woe,toe,qoe,psoe,obserr,var_jb
real(r_kind) :: dlat,dlon,dlat_earth,dlon_earth
real(r_kind) :: dlat_earth_deg,dlon_earth_deg
real(r_kind) :: cdist,disterr,disterrmax,rlon00,rlat00
real(r_kind) :: vdisterrmax,u00,v00,u0,v0
real(r_kind) :: dx,dy,dx1,dy1,w00,w10,w01,w11
real(r_kind) :: wdir,wspd
real(r_kind) :: tob,uob,vob,qob,spdob,rrob
real(r_kind) :: rhob,tdob
real(r_kind) :: pob_mb,pob_cb,pob_pa,gob
real(r_kind) :: psob_mb,psob_cb,psob_pa
real(r_kind) :: qmaxerr
real(r_kind) :: dlnpsob,dlnpob,ppb
real(r_kind) :: crit1,timedif,xmesh,pmesh
real(r_kind) :: sstime,tdiff
real(r_kind) :: tsavg,ff10,sfcr,zz
real(r_kind) :: es,qsat,rhob_calc,tdob_calc,tdry
real(r_kind) :: dummy
real(r_kind) :: del,ediff,errmin,jbmin
real(r_kind) :: tvflg
real(r_kind) :: presl(nsig)
real(r_kind) :: obstime(6,1)
real(r_kind) :: obsloc(2,1)
real(r_kind) :: obstmp(2,1)
real(r_kind) :: obswnd(4,1)
real(r_kind) :: obsfmr(2,1)
real(r_kind) :: obsmst(3,1)
real(r_kind) :: obsprs(1,1)
real(r_kind) :: obspsf(1,1)
real(r_kind) :: obsg10(1,1)
real(r_kind) :: obsqcm(2,1)
real(r_double) :: rstation_id
real(r_double) :: r_prvstg(1,1),r_sprvstg(1,1)
real(r_kind), allocatable,dimension(:,:) :: cdata_all,cdata_out
real(r_kind), allocatable,dimension(:) :: presl_thin
! Equivalence to handle character names
equivalence(r_prvstg(1,1),c_prvstg)
equivalence(r_sprvstg(1,1),c_sprvstg)
equivalence(rstation_id,c_station_id)
! Data
data bulstr1 / 'BUHD' /
data bulstr2 / 'BORG' /
data timestr / 'YEAR MNTH DAYS HOUR MINU SECO' /
data sfmrstr / 'PKSWSP TRRT' /
data locstr / 'CLAT CLON' /
data tmpstr / 'QMAT TMDB' /
data mststr / 'QMDD TMDP REHU' /
data wndstr / 'QMWN WDIR WSPD PKWDSP' /
data prsstr / 'PRLC' /
data psfstr / '' / ! nor in the bufr yet
data g10str / 'GP10' /
data qcmstr / 'QHDOP QHDOM'/
data lunin / 13 /
data ithin / -9 /
data rmesh / -99.999_r_kind /
!------------------------------------------------------------------------------------------------
write(6,*)'READ_FL_HDOB: begin to read flight-level high density data ...'
! Initialize parameters
! Set common variables
ltob = obstype == 't'
luvob = obstype == 'uv'
lspdob = obstype == 'spd'
lpsob = obstype == 'ps'
lqob = obstype == 'q'
nreal = 0
iecol = 0
ierr_ps = 0
ierr_q = 0
ierr_t = 0
ierr_uv = 0
var_jb=zero
jbmin=zero
lim_qm = 4
iecol=0
if (ltob) then
nreal = 25
iecol = 2
errmin = half ! set lower bound of ob error for T or Tv
else if (luvob) then
nreal = 26
iecol = 4
errmin = one ! set lower bound of ob error for u,v winds
else if (lspdob) then
nreal = 23
iecol = 4
errmin = one
else if (lqob) then ! set lower bound of ob err for surface wind speed
nreal = 26
iecol = 3
errmin = half ! set lower bound of ob error for moisture (RH)
else if (lpsob) then
nreal = 23
iecol = 5
errmin = one_tenth ! set lower bound of ob error for moisture (RH)
else
write(6,*) ' illegal obs type in read_fl_hdob '
call stop2(94)
end if
if (perturb_obs .and. luvob) nreal = nreal+2
if (perturb_obs .and. (ltob .or. lpsob .or. lqob)) nreal = nreal+1
inflate_error = .false.
! Read in entire observation error table
! 33 pressure levels & 6 variables
! 1: pressure levels [mb]
! 2: temperature error [K]
! 3: relative humidity error*10
! 4: wind speed error [m/s]
! 5: surface pressure error [mb]
! 6: total precipitable water error [?]
! open(ietabl,file='errtable',form='formatted')
! rewind ietabl
! lcount = 0
! etabl = 1.e9_r_kind
! loopd : do
! read(ietabl,100,IOSTAT=iflag) itypex
! if( iflag /= 0 ) exit loopd
!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 loopd
! if (lcount <= 0) then
! write(6,*)'READ_FL_HDOB: obs error table not available'
! call stop2(49)
! else
! write(6,*)'READ_FL_HDOB: obs errors provided by local file errtable'
! endif
!
! Check if the obs type specified in the convinfo is in the fl hdob bufr file
! If found, get the index (nc) from the convinfo for the specified type
ntmatch = 0
ncsave = 0
do nc = 1, nconvtype
if (trim(ioctype(nc)) == trim(obstype))then
if (trim(ioctype(nc)) == 'uv' .and. ictype(nc) == 236 .or. &
trim(ioctype(nc)) == 'spd' .and. ictype(nc) == 292 .or. &
trim(ioctype(nc)) == 't' .and. ictype(nc) == 136 .or. &
trim(ioctype(nc)) == 'q' .and. ictype(nc) == 136 .or. &
trim(ioctype(nc)) == 'ps' .and. ictype(nc) == 136 ) then
ntmatch = ntmatch+1
ncsave = nc
ithin = ithin_conv(nc) ! 0: no thinning 1: thinning
itype = ictype(nc)
end if
end if
enddo
if(ntmatch == 0)then ! Return if not specified in convinfo
write(6,*) ' READ_FL_HDOB: No matching obstype found in obsinfo ',obstype
return
else
nc = ncsave
write(6,*) ' READ_FL_HDOB: Processing FL HDOB data : ', ntmatch, nc, ioctype(nc), ictype(nc), itype
end if
ncount_ps=0;ncount_q=0;ncount_t=0;ncount_uv=0
! Setup thinning parameters
use_all = .true.
ithin = ithin_conv(nc)
if (ithin > 0) then
rmesh = rmesh_conv(nc) ! horizontal mesh size
pmesh = pmesh_conv(nc) ! vertical mesh size
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_FL_HDOB: ictype(nc),rmesh,pflag,nlevp,pmesh,nc ',&
ioctype(nc),ictype(nc),rmesh,pflag,nlevp,pmesh,nc
endif
!------------------------------------------------------------------------------------------------
! Go through the bufr file to find out how mant subsets to process
nmsg = 0
maxobs = 0
call closbf(lunin)
open(lunin,file=trim(infile),form='unformatted')
call openbf(lunin,'IN',lunin)
call datelen(10)
loop_msg1: do while(ireadmg(lunin,subset,idate) >= 0)
if(nmsg == 0) call time_4dvar(idate,toff) ! time offset (hour)
nmsg = nmsg+1
loop_readsb1: do while(ireadsb(lunin) == 0)
maxobs = maxobs+1
end do loop_readsb1
end do loop_msg1
call closbf(lunin)
write(6,*) 'READ_FL_HDOB: total number of data found in the bufr file ',maxobs,obstype
write(6,*) 'READ_FL_HDOB: time offset is ',toff,' hours'
!---------------------------------------------------------------------------------------------------
! Allocate array to hold data
allocate(cdata_all(nreal,maxobs))
allocate(isort(maxobs))
! Initialize
cdata_all = zero
isort = 0
nread = 0
nchanl = 0
ntest = 0
nvtest = 0
ilon = 2
ilat = 3
! Open bufr file again for reading
call closbf(lunin)
open(lunin,file=trim(infile),form='unformatted')
call openbf(lunin,'IN',lunin)
call datelen(10)
ntb = 0
igood = 0
! Loop through BUFR file
loop_msg2: do while(ireadmg(lunin,subset,idate) >= 0)
loop_readsb2: do while(ireadsb(lunin) == 0)
ntb = ntb+1
! Extract observation bulletin info
! URNT15 --- data in Atlantic region
! URPN15 --- data in Eastern and Central Pacific region
! URPA15 --- data in West Pacitic region
! KNHC --- Air Force product
! KWBC --- NOAA product
! KBIX --- Air Force backup product
call readlc(lunin,obsbul(1,1),bulstr1)
call readlc(lunin,obsbul(2,1),bulstr2)
obs_region = 'Unknown'
if (obsbul(1,1) == 'URNT15') obs_region = 'Atlantic'
if (obsbul(1,1) == 'URPN15') obs_region = 'East and Central Pacific'
if (obsbul(1,1) == 'URPA15') obs_region = 'West Pacific'
c_station_id = 'FL_HDOB'
c_prvstg = obsbul(2,1)
c_sprvstg = obsbul(1,1)
! QC mark 9: will be monitored but not assimilated
! QC mark 4: reject - will not be monitored nor assimilated
! QC mark 3: suspect
! QC mark 2: neutral or not checked
! QC mark 1: good
! QC mark 0: keep - will be always assimilated
qcm = 0
p_qm = 0
g_qm = 0
t_qm = 0
q_qm = 0
uv_qm = 0
wspd_qm = 0
ps_qm = 0
! Read observation time
call ufbint(lunin,obstime,6,1,nlv,timestr)
idate5(1) = obstime(1,1) ! year
idate5(2) = obstime(2,1) ! month
idate5(3) = obstime(3,1) ! day
idate5(4) = obstime(4,1) ! hour
idate5(5) = obstime(5,1) ! minute
call w3fs21(idate5,nmind)
t4dv = real((nmind-iwinbgn),r_kind)*r60inv
sstime = real(nmind,r_kind)
tdiff = (sstime-gstime)*r60inv
if (l4dvar.or.l4densvar) then
if (t4dv < zero .OR. t4dv > winlen) cycle loop_readsb2
else
if (abs(tdiff)>twind) cycle loop_readsb2
endif
nread = nread+2
! Read QC control flag for HDOB positional data
! QHDOP: 0 all parameters of nominal accuracy
! 1 lat/lon questionable
! 2 geopotential altitude or static pressure questionable
! 3 both lat/lon abd GA/PS questionable
call ufbint(lunin,obsqcm,2,1,nlv,qcmstr)
call upftbv(lunin,"QHDOP",obsqcm(1,1),mxib,ibit,nib)
if (nib > 0) then
ibit(1:nib) = ibit(1:nib)-1
if (any(ibit(1:nib) > 0)) then
write(6,*) 'READ_FL_HDOB: bad positional data ... toss away'
cycle loop_readsb2
endif
else ! will keep for further QC check
write(6,*) 'READ_FL_HDOB: missing QC info '
cycle loop_readsb2
endif
! Read QC mark for HDOB meteorological status
! QHDOP: 0 all parameters of nominal accuracy
! 1 T or Td questionable
! 2 flight-level winds questionable
! 3 SFMR parameters questionable
! 4 T/Td and FL winds questionable
! 5 T/Td and SFMR questionable
! 6 FL winds and SFMR questionable
! 9 T/Td, FL winds, and SFMR questionable
call upftbv(lunin,"QHDOM",obsqcm(2,1),mxib,ibit,nib)
if (nib > 0) then
ibit(1:nib) = ibit(1:nib)-1
if (any(ibit(1:nib) == 1)) then ! for T/Td
t_qm = 0
q_qm = 4
endif
if (any(ibit(1:nib) == 2)) then ! for uv
uv_qm = 4
endif
if (any(ibit(1:nib) == 3)) then ! for SFMR data
wspd_qm = 4
endif
if (any(ibit(1:nib) == 4)) then ! for T/Td and uv
t_qm = 4
q_qm = 4
uv_qm = 4
endif
if (any(ibit(1:nib) == 5)) then ! for T/Td and SFMR data
t_qm = 4
q_qm = 4
wspd_qm = 4
endif
if (any(ibit(1:nib) == 6)) then ! for uv and SFMR data
uv_qm = 4
wspd_qm = 4
endif
if (any(ibit(1:nib) == 9)) then ! for T/Td, uv, and SFMR data
t_qm = 4
q_qm = 4
uv_qm = 4
wspd_qm = 4
endif
else
write(6,*) 'READ_FL_HDOB: missing QC info'
endif
usage = zero ! will be considered for assimilation
! subject to further QC in setupt subroutine
iuse = icuse(nc) ! assimilation flag
if (iuse <=0) usage = r100 ! will be monitored but not assimilated
! Read observation location (lat/lon degree)
call ufbint(lunin,obsloc,2,1,nlv,locstr)
if (obsloc(1,1) == missing .or. abs(obsloc(1,1)) > 90.0_r_kind .or. &
obsloc(2,1) == missing .or. abs(obsloc(2,1)) > 360.0_r_kind) then
write(6,*) 'READ_FL_HDOB: bad lat/lon values: ', obsloc(1,1),obsloc(2,1)
cycle loop_readsb2
endif
dlon_earth_deg = obsloc(2,1)
dlat_earth_deg = obsloc(1,1)
dlon_earth = obsloc(2,1)*deg2rad ! degree to radian
dlat_earth = obsloc(1,1)*deg2rad ! degree to radian
! Convert obs lat/lon to rotated coordinate and check
! if the obs is outside of the regional domain
if (regional) then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
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_readsb2
else
dlon = dlon_earth
dlat = dlat_earth
call grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
! Read flight-level pressure [Pa] and convert to [cb]
call ufbint(lunin,obsprs,1,1,nlv,prsstr)
if (obsprs(1,1) >= missing .or. &
obsprs(1,1) .gt. 110000.0_r_kind .or. obsprs(1,1) .lt. 5000.0_r_kind) then
write(6,*) 'READ_FL_HDOB: bad flight-level pressure [pa] values: ', obsprs(1,1)
cycle loop_readsb2
endif
pob_pa = obsprs(1,1) ! [Pa]
pob_mb = obsprs(1,1)*r0_01 ! convert {Pa] to [mb]
pob_cb = obsprs(1,1)*r0_001 ! convert [Pa] to [cb]
dlnpob = log(pob_cb) ! [cb]
! Read flight-level geopotential height [(m/s)**2] and convert to height [m]
call ufbint(lunin,obsg10,1,1,nlv,g10str)
if (obsg10(1,1) == missing) then
write(6,*) 'READ_FL_HDOB: bad flight-level geopotential height [(m/s)**2] values: ', obsg10(1,1)
cycle loop_readsb2
endif
gob = obsg10(1,1)/grav ! convert to height [m]
! Get observation error from error table
ppb = max(zero,min(pob_mb,r2000))
if(.not. njqc) then
if(ppb >= etabl(itype,1,1)) k1 = 1
do kl = 1,32
if(ppb >= etabl(itype,kl+1,1) .and. ppb <= etabl(itype,kl,1)) k1 = kl
end do
if(ppb <= etabl(itype,33,1)) k1 = 5
k2 = k1+1
ediff = etabl(itype,k2,1)-etabl(itype,k1,1)
if (abs(ediff) > tiny_r_kind) then
del = (ppb-etabl(itype,k1,1))/ediff
else
del = huge_r_kind
endif
del = max(zero,min(del,one))
obserr = (one-del)*etabl(itype,k1,iecol)+del*etabl(itype,k2,iecol)
obserr = max(obserr,errmin)
endif
! Read extrapolated surface pressure [pa] and convert to [cb]
if (lpsob) then
call ufbint(lunin,obspsf,1,1,nlv,psfstr)
if (obspsf(1,1) >= missing .or. &
obspsf(1,1) .gt. 110000.0_r_kind .or. obspsf(1,1) .lt. 5000.0_r_kind) then
write(6,*) 'READ_FL_HDOB: bad surface pressure [pa] values: ', obspsf(1,1)
cycle loop_readsb2
endif
psob_pa = obspsf(1,1) ! [Pa]
psob_mb = obspsf(1,1)*r0_01 ! convert {Pa] to [mb]
psob_cb = obspsf(1,1)*r0_001 ! convert [Pa] to [cb]
dlnpsob = log(psob_cb) ! [cb]
! Get observation error from error table
if (njqc) then
ppb = max(zero,min(pob_mb,r2000))
itypey=itype
ierr_ps=0
do i =1,maxsub_ps
if( icsubtype(nc) == isuble_ps(itypey,i) ) then
ierr_ps=i+1
exit
else if( i == maxsub_ps .and. icsubtype(nc) /= isuble_ps(itypey,i)) then
ncount_ps=ncount_ps+1
do j=1,maxsub_ps
if(isuble_ps(itypey,j) ==0 ) then
ierr_ps=j+1
exit
endif
enddo
if (ncount_ps ==1) then
write(6,*) 'READ_FL_HDOB,WARNING!!psob: cannot find subtyep in the error,&
table,itype,iosub=',itypey,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0, error table column= ',ierr_ps
endif
endif
enddo
if(ppb >= etabl_ps(itypey,1,1)) k1 = 1
do kl = 1,32
if(ppb >= etabl_ps(itypey,kl+1,1) .and. ppb <= etabl_ps(itypey,kl,1)) k1 = kl
end do
if(ppb <= etabl_ps(itypey,33,1)) k1 = 5
k2 = k1+1
ediff = etabl_ps(itypey,k2,1)-etabl_ps(itypey,k1,1)
if (abs(ediff) > tiny_r_kind) then
del = (ppb-etabl_ps(itypey,k1,1))/ediff
else
del = huge_r_kind
endif
del = max(zero,min(del,one))
obserr = (one-del)*etabl_ps(itypey,k1,ierr_ps)+del*etabl_ps(itypey,k2,ierr_ps)
var_jb=(one-del)*btabl_ps(itypey,k1,ierr_ps)+del*btabl_ps(itypey,k2,ierr_ps)
obserr = max(obserr,errmin)
var_jb=max(var_jb,jbmin)
endif
endif
! Convert raw temperature data to T or Tv
! Read temperature [K],dew pointer temperature [K] and related QC marks
! If both T and Td are available and in good condition, then calculate
! virtual temperature and pass it to setupt. Otherwise, The dry airs
! temperature will be used
if (ltob) then
qcm = t_qm
if (qcm > 3) usage = r100
call ufbint(lunin,obstmp,2,1,nlv,tmpstr)
call ufbint(lunin,obsmst,3,1,nlv,mststr)
tob = obstmp(2,1) ! airs temperature [K]
tdob = obsmst(2,1) ! dew point temperature [K]
rhob = obsmst(3,1) ! relative humidity (%)
if (tob >= missing .or. tob <= 170.0_r_kind .or. tob >= 320.0_r_kind) then
cycle loop_readsb2
endif
if (tdob >= missing) then
tvflg = one ! tob is sensible temperature
qob = bmiss
rhob = bmiss
tob = obstmp(2,1)
else
tvflg = 0 ! tob is virtual temperature temperature
! Calculate specific humidity from tob and td
if (rhob >= missing) then
! Calculate RH [%] since rhob is missing
rhob_calc = exp((one-tob/tdob)*(hvap/rv)/tob) ! e.g. rh=0.98
call fpvsx_ad(tob,es,dummy,dummy,.false.)
qsat = eps*es/(pob_cb-omeps*es)
rhob = rhob_calc ! calculate RH (%) since rhob is missing
qob = rhob*qsat
else
call fpvsx_ad(tob,es,dummy,dummy,.false.)
qsat = eps*es/(pob_cb-omeps*es)
tdob_calc = tob*(one-tob*log(rhob/100)) ! for comparison
qob = rhob*qsat
endif
tob = tob*(1.0_r_kind+0.61_r_kind*qob) ! conver t to tv
endif
! Get observation error from error table
if (njqc) then
ppb = max(zero,min(pob_mb,r2000))
itypey=itype
ierr_t=0
do i =1,maxsub_t
if( icsubtype(nc) == isuble_t(itypey,i) ) then
ierr_t=i+1
exit
else if( i == maxsub_t .and. icsubtype(nc) /= isuble_t(itypey,i)) then
ncount_t=ncount_t+1
do j=1,maxsub_t
if(isuble_t(itypey,j) ==0 ) then
ierr_t=j+1
exit
endif
enddo
if(ncount_t ==1) then
write(6,*) 'READ_FL_HDOB,WARNING!! tob:cannot find subtyep in the error table,&
itype,iosub=',itype,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0,error table column=',ierr_t
endif
endif
enddo
if(ppb >= etabl_t(itypey,1,1)) k1 = 1
do kl = 1,32
if(ppb >= etabl_t(itypey,kl+1,1) .and. ppb <= etabl_t(itypey,kl,1)) k1 = kl
end do
if(ppb <= etabl_t(itypey,33,1)) k1 = 5
k2 = k1+1
ediff = etabl_t(itypey,k2,1)-etabl_t(itypey,k1,1)
if (abs(ediff) > tiny_r_kind) then
del = (ppb-etabl_t(itypey,k1,1))/ediff
else
del = huge_r_kind
endif
del = max(zero,min(del,one))
obserr = (one-del)*etabl_t(itypey,k1,ierr_t)+del*etabl_t(itypey,k2,ierr_t)
var_jb = (one-del)*btabl_t(itypey,k1,ierr_t)+del*btabl_t(itypey,k2,ierr_t)
obserr = max(obserr,errmin)
var_jb=max(var_jb,jbmin)
endif
endif
! Convert raw moisture data from dew point temperature to specific humidity
! Read temperature,dew point temperature,relative humidity,
! related QC mark and then calculate specific humidity
if (lqob) then
qcm = q_qm
if (qcm > 3) usage = r100
call ufbint(lunin,obstmp,2,1,nlv,tmpstr)
call ufbint(lunin,obsmst,3,1,nlv,mststr)
tob = obstmp(2,1) ! dry airs temperature [K]
tdob = obsmst(2,1) ! dew point temperature [K]
rhob = obsmst(3,1) ! relative humidity (%)
tdry = tob
if (tob >= missing .or. tdob >= missing .or. &
tob <= 170.0_r_kind .or. tob >= 320.0_r_kind .or. &
tdob <= 170.0_r_kind .or. tdob >= 320.0_r_kind) then
cycle loop_readsb2
endif
! Calculate specific humidity from relative humidity if abailable
if (rhob >= missing) then
rhob_calc = exp((one-tob/tdob)*(hvap/rv)/tob) ! e.g. rh=0.98
call fpvsx_ad(tob,es,dummy,dummy,.false.)
qsat = eps*es/(pob_cb-omeps*es)
rhob = rhob_calc ! calculate RH (%) since rhob is missing
qob = rhob*qsat
else
call fpvsx_ad(tob,es,dummy,dummy,.false.)
qsat = eps*es/(pob_cb-omeps*es)
tdob_calc = tob*(one-tob*log(rhob/100)) ! for comparison
qob = rhob*qsat
endif
! write(4000,1004) nread,pob_mb,tob,tdob,qob,qsat,rhob,q_qm,usage
!1004 format(i6,6(1x,e20.12),1x,i5,1x,f5.0)
! Get observation error from error table
if (njqc) then
ppb = max(zero,min(pob_mb,r2000))
itypey=itype
ierr_q=0
do i =1,maxsub_q
if( icsubtype(nc) == isuble_q(itypey,i) ) then
ierr_q=i+1
exit
else if( i == maxsub_q .and. icsubtype(nc) /= isuble_q(itypey,i)) then
ncount_q=ncount_q+1
do j=1,maxsub_q
if(isuble_q(itypey,j) ==0 ) then
ierr_q=j+1
exit
endif
enddo
if( ncount_q ==1 ) then
write(6,*) 'READ_FL_HDOB,WARNING!! qob:cannot find subtyep in the error table,&
itype,iosub=',itype,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0,error table column=',ierr_q
endif
endif
enddo
if(ppb >= etabl_q(itypey,1,1)) k1 = 1
do kl = 1,32
if(ppb >= etabl_q(itypey,kl+1,1) .and. ppb <= etabl_q(itypey,kl,1)) k1 = kl
end do
if(ppb <= etabl_q(itypey,33,1)) k1 = 5
k2 = k1+1
ediff = etabl_q(itypey,k2,1)-etabl_q(itypey,k1,1)
if (abs(ediff) > tiny_r_kind) then
del = (ppb-etabl_q(itypey,k1,1))/ediff
else
del = huge_r_kind
endif
del = max(zero,min(del,one))
obserr = (one-del)*etabl_q(itypey,k1,ierr_q)+del*etabl_q(itypey,k2,ierr_q)
var_jb = (one-del)*btabl_q(itypey,k1,ierr_q)+del*btabl_q(itypey,k2,ierr_q)
obserr = max(obserr,errmin)
var_jb=max(var_jb,jbmin)
endif
endif
! Convert raw wind data from wind direction/speed to u & v winds
! Read wind direction [degree true], speed [m/s] and related QC mark
! Convert wind direction and spped to u and v wind components
if (luvob) then
qcm = uv_qm
if (qcm > 3) usage=r100
call ufbint(lunin,obswnd,4,1,nlv,wndstr)
wdir = obswnd(2,1) ! degree true
wspd = obswnd(3,1) ! m/s
if (wdir >= missing .or. wspd >= missing) cycle loop_readsb2
uob = -wspd*sin(wdir*deg2rad) ! u-wind component
vob = -wspd*cos(wdir*deg2rad) ! v-wind component
endif
! Read surface wind speed [m/s] and total rain rate [mm/hr] from SFMR
if (lspdob) then
qcm = wspd_qm
if (qcm > 3) usage=r100
call ufbint(lunin,obsfmr,2,1,nlv,sfmrstr)
! print*, 'PKSWSP = ', obstype,obsfmr(1,1)
! print*, 'TRRP = ', obstype,obsfmr(2,1)
spdob = obsfmr(1,1) ! surface wind speed
rrob = obsfmr(2,1) ! rain rate
if (spdob >= missing .or. rrob >=missing) cycle loop_readsb2
endif
if( lspdob .or. luvob) then
! Get observation error from error table
if (njqc) then
ppb = max(zero,min(pob_mb,r2000))
itypey=itype
ierr_uv=0
do i =1,maxsub_uv
if( icsubtype(nc) == isuble_uv(itypey,i) ) then
ierr_uv=i+1
exit
else if( i == maxsub_uv .and. icsubtype(nc) /= isuble_uv(itypey,i)) then
ncount_uv=ncount_uv+1
do j=1,maxsub_uv
if(isuble_uv(itypey,j) ==0 ) then
ierr_uv=j+1
exit
endif
enddo
if(ncount_uv ==1) then
write(6,*) 'READ_FL_HDOB,WARNING!! uvob:cannot find subtyep in the error table,&
itype,iosub=',itype,icsubtype(nc)
write(6,*) 'read error table at colomn subtype 0,error table column=',ierr_uv
endif
endif
enddo
if(ppb >= etabl_uv(itypey,1,1)) k1 = 1
do kl = 1,32
if(ppb >= etabl_uv(itypey,kl+1,1) .and. ppb <= etabl_uv(itypey,kl,1)) k1 = kl
end do
if(ppb <= etabl_uv(itypey,33,1)) k1 = 5
k2 = k1+1
ediff = etabl_uv(itypey,k2,1)-etabl_uv(itypey,k1,1)
if (abs(ediff) > tiny_r_kind) then
del = (ppb-etabl_uv(itypey,k1,1))/ediff
else
del = huge_r_kind
endif
del = max(zero,min(del,one))
obserr = (one-del)*etabl_uv(itypey,k1,ierr_uv)+del*etabl_uv(itypey,k2,ierr_uv)
var_jb = (one-del)*btabl_uv(itypey,k1,ierr_uv)+del*btabl_uv(itypey,k2,ierr_uv)
obserr=max(obserr,errmin)
var_jb=max(var_jb,jbmin)
endif
endif
! Obtain information necessary for conventional data assimilation
! Detect surface type (isfg) and skin temperature (tsavg)
! isflg - surface flag
! 0 sea
! 1 land
! 2 sea ice
! 3 snow
! 4 mixed
if ( .not. twodvar_regional) then
call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg)
endif
! Get information from surface file necessary for conventional data
call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz)
! Process data passed quality control
igood = igood+1
! Process data thinning procedure on good data
if (ithin > 0) then
if (pflag == 0) 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
endif ! pflag
ntmp = ndata ! counting moved into map3grids
! Set data quality index for thinning
if (thin4d) then
timedif = zero
else
timedif = abs(t4dv-toff)
endif
crit1 = timedif/r6+half
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,&
pob_cb,crit1,ndata,iout,igood,iiout,luse,.false.,.false.)
if (.not. luse) cycle loop_readsb2
if(iiout > 0) isort(iiout) = 0
if (ndata > ntmp) then
nodata = nodata+2
if (luvob) &
nodata = nodata+2
endif
isort(igood) = iout
else
ndata = ndata+1
nodata = nodata+2
if (luvob) &
nodata = nodata+2
iout = ndata
isort(igood) = iout
endif ! ithin
!-------------------------------------------------------------------------------------------------
! Write data into output arrays
if (var_jb >=10.0_r_kind) var_jb=zero
if (qcm == 3) inflate_error = .true.
if (lpsob) then
qcm = ps_qm
psoe = obserr*one_tenth ! convert from mb to cb
if (inflate_error) psoe = psoe*r1_2
if (qcm > lim_qm ) psoe = psoe*1.0e6_r_kind
cdata_all( 1,iout)=psoe ! 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(dlnpsob) ! pressure (in cb)
cdata_all( 5,iout)=zz ! surface height ! use model terrian elevation from model surface file
cdata_all( 6,iout)=bmiss ! surface temperature ! this is not provided
cdata_all( 7,iout)=rstation_id ! station id
cdata_all( 8,iout)=t4dv ! time
cdata_all( 9,iout)=nc ! type
cdata_all(10,iout)=qcm ! quality mark
cdata_all(11,iout)=obserr*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_deg ! earth relative longitude (degree)
cdata_all(18,iout)=dlat_earth_deg ! earth relative latitude (degree)
cdata_all(19,iout)=gob ! 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
cdata_all(23,iout)=var_jb ! non linear qc b
if(perturb_obs)cdata_all(24,iout)=ran01dom()*perturb_fact ! ps perturbation
endif
! Winds --- u, v components
if (luvob) then
woe = obserr
if (pob_mb < r50) woe = woe*r1_2
if (inflate_error) woe = woe*r1_2
if (qcm > lim_qm ) woe = woe*1.0e6_r_kind
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)=gob ! height of observation (m)
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 ! index of type in convinfo
cdata_all(11,iout)=gob ! station elevation (m)
cdata_all(12,iout)=qcm ! quality mark
cdata_all(13,iout)=obserr ! original obs error
cdata_all(14,iout)=usage ! usage parameter
cdata_all(15,iout)=idomsfc ! dominate surface
cdata_all(16,iout)=tsavg ! skin temperature
cdata_all(17,iout)=ff10 ! 10 meter wind
cdata_all(18,iout)=sfcr ! surface roughness
cdata_all(19,iout)=dlon_earth_deg ! earth relative longitude (degree)
cdata_all(20,iout)=dlat_earth_deg ! earth relative latitude (degree)
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)=qcm ! cat
cdata_all(25,iout)=var_jb ! non linear qc
cdata_all(26,iout)=one
if(perturb_obs)then
cdata_all(26,iout)=ran01dom()*perturb_fact ! u perturbation
cdata_all(27,iout)=ran01dom()*perturb_fact ! v perturbation
endif
write(3000,1003) nread,pob_mb,uob,vob,qcm,usage
1003 format(i12,3e25.18,f5.0,f5.0)
endif
! Temperature
if(ltob) then
toe = obserr
if (pob_mb < r100) toe = toe*r1_2
if (inflate_error) toe = toe*r1_2
if (qcm > lim_qm ) toe = toe*1.0e6_r_kind
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)
cdata_all( 5,iout)=tob ! temperature ob.
cdata_all( 6,iout)=rstation_id ! station id
cdata_all( 7,iout)=t4dv ! time
cdata_all( 8,iout)=nc ! ob type
cdata_all( 9,iout)=tvflg ! qtflg (0=virtual temperature 1=sensible temperature)
cdata_all(10,iout)=qcm ! quality mark
cdata_all(11,iout)=obserr ! 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_deg ! earth relative longitude (degrees)
cdata_all(18,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(19,iout)=gob ! station elevation (m)
cdata_all(20,iout)=gob ! 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)=qcm ! cat
cdata_all(25,iout)=var_jb ! non linear qc
if(perturb_obs) &
cdata_all(26,iout)=ran01dom()*perturb_fact ! t perturbation
write(1000,1001) nread,tdiff,tvflg,pob_mb,qob,rhob,obstmp(2,1),tob,qcm,usage
1001 format(i12,f8.3,f5.0,5e25.18,f5.0,f5.0)
endif
! Specific humidity
if(lqob) then
qoe = obserr*one_tenth ! RH (e.g. 0.98)
qmaxerr = emerr
if (inflate_error) then
qmaxerr = emerr*r0_7
qoe = qoe*r1_2
end if
if (qcm > lim_qm ) qoe = qoe*1.0e6_r_kind
cdata_all( 1,iout)=qoe ! q error (RH e.g. 0.98)
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)=qob ! q ob (specific humidity)
cdata_all( 6,iout)=rstation_id ! station id
cdata_all( 7,iout)=t4dv ! time
cdata_all( 8,iout)=nc ! type
cdata_all( 9,iout)=qmaxerr ! q max error (RH e.g. 0.2)
cdata_all(10,iout)=tdry ! dry temperature (obsis tv)
cdata_all(11,iout)=qcm ! quality mark
cdata_all(12,iout)=obserr*one_tenth ! original obs error (RH e.g. 0.98)
cdata_all(13,iout)=usage ! usage parameter
cdata_all(14,iout)=idomsfc ! dominate surface type
cdata_all(15,iout)=dlon_earth_deg ! earth relative longitude (degree)
cdata_all(16,iout)=dlat_earth_deg ! earth relative latitude (degree)
cdata_all(17,iout)=gob ! station elevation (m)
cdata_all(18,iout)=gob ! 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
cdata_all(22,iout)=qcm ! cat
cdata_all(26,iout)=var_jb ! non linear qc
if(perturb_obs) &
cdata_all(27,iout)=ran01dom()*perturb_fact ! q perturbation
! write(2000,1002)nread,tdiff,tvflg,pob_mb,tob,qob,rhob,qoe,obserr*one_tenth,qcm,usage
1002 format(i12,f8.3,f5.0,6e20.12,i5,1x,f5.0)
endif
! Winds --- surface wind speed
if (lspdob) then
woe = obserr
if (inflate_error) woe = woe*r1_2
if (qcm > lim_qm ) woe = woe*1.0e6_r_kind
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)=dlnpsob ! ln(surface pressure in cb)
cdata_all( 5,iout)=spdob*sqrt(two)*half ! u obs
cdata_all( 6,iout)=spdob*sqrt(two)*half ! 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)=r10 ! elevation of observation ! 10-m wind
cdata_all(11,iout)=qcm ! quality mark
cdata_all(12,iout)=obserr ! 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_deg ! earth relative longitude (degree)
cdata_all(19,iout)=dlat_earth_deg ! earth relative latitude (degree)
cdata_all(20,iout)=gob ! 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
endif
end do loop_readsb2
end do loop_msg2
! 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
! Write header record and data to output file for further processing
allocate(cdata_out(nreal,ndata))
do i=1,ndata
do k=1,nreal
cdata_out(k,i)=cdata_all(k,i)
end do
end do
deallocate(cdata_all)
! deallocate(etabl)
call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs)
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
write(lunout) cdata_out
deallocate(cdata_out)
900 continue
if(diagnostic_reg .and. ntest>0) write(6,*)'READ_FL_HDOB: ',&
'ntest, disterrmax=', ntest,disterrmax
if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_FL_HDOB: ',&
'nvtest,vdisterrmax=',ntest,vdisterrmax
if (ndata == 0) then
call closbf(lunin)
write(6,*)'READ_FL_HDOB: no data to process'
endif
write(6,*)'READ_FL_HDOB: nreal=',nreal
write(6,*)'READ_FL_HDOB: ntb,nread,ndata,nodata=',ntb,nread,ndata,nodata
close(lunin)
! End of routine
return
end subroutine read_fl_hdob