subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
prsl_full,nobs,nrec_start)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_prepbuf read obs from prepbufr file
! prgmmr: parrish org: np22 date: 1990-10-07
!
! abstract: This routine reads conventional data found in the prepbufr
! file. Specific observation types read by this routine
! include surface pressure, temperature, winds (components
! and speeds), moisture and total precipitable water.
!
! When running the gsi in regional mode, the code only
! retains those observations that fall within the regional
! domain
!
! program history log:
! 1990-10-07 parrish
! 1998-05-15 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-02-13 derber, j. - clean up and modify vertical weighting
! 2004-06-16 treadon - update documentation
! 2004-07-23 derber - modify to include conventional sst
! 2004-07-29 treadon - add only to module use, add intent in/out
! 2004-07-30 derber - generalize number of data records per obs type
! 2004-08-26 derber - fix many errors in reading of sst data
! 2004-08-27 kleist - modify pressure calculation
! 2004-10-28 kleist - correct array index bug in hybrid pressure calculation
! 2004-11-16 treadon - deallocate(etabl) prior to exiting routine
! 2005-02-10 treadon - add call destroygrids for obstype = sst
! 2005-05-24 pondeca - add surface analysis option
! 2005-02-24 treadon - remove hardwired setting of metar ps obs error
! 2005-05-27 derber - reduce t, uv, ps error limits
! 2005-05-27 kleist/derber - add option to read in new ob error table
! 2005-07-19 derber - clean up code a bit
! 2005-07-27 derber - add print of monitoring and reject data
! 2005-08-02 derber - modify to use convinfo file
! 2005-09-08 derber - modify to use input group time window
! 2005-10-11 treadon - change convinfo read to free format
! 2005-10-17 treadon - add grid and earth relative obs location to output file
! 2005-10-18 treadon - remove array obs_load and call to sumload
! 2005-10-26 treadon - add routine tag to convinfo printout
! 2006-02-03 derber - modify to count read/keep data and new obs control
! 2006-02-03 treadon - use interpolated guess 3d pressure field in errormod
! 2006-02-08 derber - modify to use new convinfo module
! 2006-02-09 treadon - save height for wind observations
! 2006-02-23 kistler - modify to add optional data thinning
! 2006-02-23 kistler - raob instument as subtype and solar elv angle computed
! 2006-02-24 derber - modify to take advantage of convinfo module
! 2006-04-03 derber - modify to properly handle height of surface obs
! 2006-04-05 wu - changes to read in GPS IPW (type 153)
! 2006-05-18 middlecoff/treadon - add huge_i_kind upper limit on nint
! 2006-05-29 treadon - increase nreal to pass more information to setup routines
! 2006-06-08 su - added the option to turn off oiqc
! 2006-06-21 wu - deallocate etabl array
! 2006-07-28 derber - temporarily add subtype for meteosat winds based on sat ID
! 2006-07-31 kleist - change to surface pressure ob error from ln(ps) to ps(cb)
! 2006-10-25 sienkiewicz - add blacklist of raob data
! 2006-12-01 todling - embedded blacklist into a module
! 2007-02-13 parrish - add ability to use obs files with ref time different from analysis time
! 2007-02-20 wu - correct errors in quality mark checks
! 2007-03-01 tremolet - measure time from beginning of assimilation window
! 2007-03-15 su - remove the error table reading part to a subroutine
! 2007-04-24 wu - add TAMDAR (134) to be used as sensible T
! 2007-05-17 kleist - generalize flag for virtual/sensible temperature obs
! 2007-09-28 treadon - truncate/expand obs time to remove extraneous bits
! 2007-10-03 su - Add reading qc mark from satellite wind
! 2007-10-24 Pondeca - add ability to use use_list on mesonet winds
! 2007-11-03 su - modify conventional thinning algorithm
! 2008-03-28 wu - add code to generate optional observation perturbations
! 2008-03-31 li.bi - add ascat
! 2008-05-27 safford - rm unused vars and uses
! 2008-06-02 treadon - check iret from inital readmg and act accordingly
! 2008-09-08 lueken - merged ed's changges into q1fy09 code
! 2008-21-25 todling - adapted Tremolet 2007-03-01 change of time window
! - remove unused vars
! 2009-07-08 pondeca - add ability to convert virtual temperature
! obs into sensible temperature for 2dvar
! 2009-07-08 park,pondeca - add option to use the hilbert curve-based
! cross-validation for 2dvar
! 2009-07-08 pondeca - move handling of "provider use_list" for mesonet winds
! to the new module sfcobsqc
! 2010-03-29 hu - add code to read cloud observation from METAR and NESDIS cloud products
! 2010-05-15 kokron - safety measure: initialize cdata_all to zero
! 2010-08-23 tong - add flg as an input argument of map3grids, so that the subroutine can be used for
! thinning grid with either pressure or height as the vertical coordinate.
! flg=-1 for prepbufr data thinning grid (pressure as the vertical coordinate).
! 2010-09-08 parrish - remove subroutine check_rotate_wind. This was a debug routine introduced when
! the reference wind rotation angle was stored as an angle, beta_ref. This field
! had a discontinuity at the date line (180E), which resulted in erroneous wind
! rotation angles for a small number of winds whose rotation angle was interpolated
! from beta_ref values across the discontinuity. This was fixed by replacing the
! beta_ref field with cos_beta_ref, sin_beta_ref.
! 2010-10-19 wu - add code to limit regional use of MAP winds with P less than 400 mb
! 2010-11-13 su - skip satellite winds from prepbufr
! 2010-11-18 treadon - add check for small POB (if POB<tiny_r_kind then POB=bmiss)
! 2011-02-14 zhu - add gust and visibility
! 2011-07-13 wu - not use mesonet Psfc when 8th character of sid is "x"
! 2011-08-01 lueken - added module use deter_sfc_mod and fixed indentation
! 2011-08-27 todling - add use_prepb_satwnd; cleaned out somebody's left over's
! 2011-11-14 wu - pass CAT to setup routines for raob level enhancement
! 2012-04-03 s.liu - thin new VAD wind
! 2012-11-12 s.liu - identify new VAD wind by vertical resolution
! 2012-08-29 akella - extend nst_gsi option to handle sstobs
! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS)
! 2013-01-26 parrish - WCOSS debug compile error for pflag used before initialized.
! Initialize pflag=0 at beginning of subroutine.
! 2013-02-28 sienkiewicz - put in subset via SAID for kx=290 ASCAT to allow
! separate control of metop-a and metop-b ASCAT if
! needed
! 2013-05-03 sienkiewicz - if ACARS SID == 'ACARS' take ID from ACID instead
! 2013-05-15 zhu - add phase of aircraft flight and vertical velocity for aircraft data
! - match aircraft obs with temperature bias file
! - add new tail number info if there is any
! - add aircraft_t_bc_pof and aircraft_t_bc
! 2013-05-28 wu - add subroutine sonde_ext and call to the subroutine for ext_sonde option
! 2013-06-07 zhu - read aircraft data from prepbufr_profl when aircraft_t_bc=.true.
! 2013-09-08 s.liu - increase nmsgmax to 100000 to read NESDIS cloud product
! 2013-12-08 s.liu - identify VAD wind based on sub type
! 2014-02-28 sienkiewicz - added code for option aircraft_t_bc_ext for external aircraft bias table
! 2014-03-19 pondeca - add 10m wind speed
! 2014-04-10 pondeca - add td2m,mxtm,mitm,pmsl
! 2014-04-15 Su - read errtable and non linear qc b table
! 2014-05-07 pondeca - add significant wave height (howv)
! 2014-06-16 carley/zhu - add tcamt and lcbas
! 2014-06-26 carley - simplify call to apply_hilbertcurve
! 2014-11-20 zhu - added code for aircraft temperature kx=130
! 2014-10-01 Xue - add gsd surface observation uselist
! 2015-02-23 Rancic/Thomas - add thin4d to time window logical
! 2015-03-23 Su -fix array size with maximum message and subset number from fixed number to
! dynamic allocated array
! 2015-07-10 pondeca - add cloud ceiling height (cldch)
! 2015-10-01 guo - consolidate use of ob location (in deg)
! 2016-02-09 Sienkiewicz - explicit KX definition for drifting buoys (formerly ID'd by subtype 562)
! 2016-02-10 s.liu - thin new VAD wind in time level
! 2016-03-15 Su - modified the code so that the program won't stop when no subtype
! is found in non linear qc error tables and b table
! 2016-05-05 pondeca - add 10-m u-wind and v-wind (uwnd10m, vwnd10m)
! 2016-06-01 zhu - use errormod_aircraft
! 2017-06-17 levine - add GLERL program code lookup
! 2017-03-21 Su - add option to thin conventional data in 4 dimension
! 2019-09-27 Su - add hilbert curve application to aircraft winds
! 2018-08-16 akella - explicit KX definition for ships (formerly ID'd by subtype 522/523)
! 2019-02-06 levine - Add lookup of sensor height for mesonet winds
! 2019-06-17 mmorris - Update adjust_goescldobs to reject clear cloud obs over water at night
! 2019-09-27 Su - add hilbert curve application to aircraft winds
! 2019-12-05 mmorris - Update adjust_goescldobs to reject ALL clear cloud obs at night
!
! 2020-05-04 wu - no rotate_wind for fv3_regional
! 2020-09-05 CAPS(C. Tong) - add flag for new vadwind obs to assimilate around the analysis time only
! 2024-10-04 Collard - Saildrone subtype introduced in GSI (not in prepobs) kx=180,280 subtype=02
! 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
! nrec_start - number of subsets without useful information
!
! 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
! 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,deg2rad,fv,t0c,half,&
three,four,rad2deg,tiny_r_kind,huge_r_kind,huge_i_kind,&
r60inv,r10,r100,r2000
use constants,only: rearth,stndrd_atmos_ps,rd,grav
use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,&
tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,&
rlats,rlons,twodvar_regional,fv3_regional
use convinfo, only: nconvtype,ctwind, &
ncmiter,ncgroup,ncnumgrp,icuse,ictype,icsubtype,ioctype, &
ithin_conv,rmesh_conv,pmesh_conv,pmot_conv,ptime_conv, &
use_prepb_satwnd
use convinfo, only: id_drifter,id_ship
use obsmod, only: iadate,oberrflg,perturb_obs,perturb_fact,ran01dom,hilbert_curve
use obsmod, only: blacklst,offtime_data,bmiss,ext_sonde,time_offset, vad_near_analtime
use aircraftinfo, only: aircraft_t_bc,aircraft_t_bc_pof,ntail,taillist,idx_tail,npredt,predt, &
aircraft_t_bc_ext,ntail_update,max_tail,nsort,itail_sort,idx_sort,timelist
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 converr_pw,only: etabl_pw,isuble_pw,maxsub_pw
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,time_4dvar,winlen,thin4d
use convthin, only: make3grids,map3grids,map3grids_m,del3grids,use_all
use convthin_time, only: make3grids_tm,map3grids_tm,map3grids_m_tm,del3grids_tm,use_all_tm
use qcmod, only: errormod,errormod_aircraft,noiqc,newvad,njqc
use qcmod, only: pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres
use qcmod, only: nrand
use nltransf, only: nltransf_forward
use blacklist, only : blacklist_read,blacklist_destroy
use blacklist, only : blkstns,blkkx,ibcnt
use sfcobsqc,only: init_rjlists,get_usagerj,get_gustqm,destroy_rjlists
use sfcobsqc,only: init_gsd_sfcuselist,apply_gsd_sfcuselist,destroy_gsd_sfcuselist
use windht,only: init_windht_lists,readin_windht_list,destroy_windht_lists,find_wind_height
use hilbertcurve,only: init_hilbertcurve, accum_hilbertcurve, &
apply_hilbertcurve,destroy_hilbertcurve
use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error
use jfunc, only: tsensible
use deter_sfc_mod, only: deter_sfc_type,deter_sfc2
use gsi_nstcouplermod, only: nst_gsi,nstinfo
use gsi_nstcouplermod, only: gsi_nstcoupler_deter
use aircraftobsqc, only: init_aircraft_rjlists,get_aircraft_usagerj,&
destroy_aircraft_rjlists
use adjust_cloudobs_mod, only: adjust_convcldobs,adjust_goescldobs
use mpimod, only: npe
use rapidrefresh_cldsurf_mod, only: i_gsdsfc_uselist,i_gsdqc,i_ens_mean
use gsi_io, only: verbose
use phil2, only: denest ! hilbert curve
implicit none
! Declare passed variables
character(len=*) ,intent(in ) :: infile,obstype
character(len=20) ,intent(in ) :: sis
integer(i_kind) ,intent(in ) :: lunout,nrec_start
integer(i_kind) ,intent(inout) :: nread,ndata,nodata
integer(i_kind),dimension(npe) ,intent(inout) :: nobs
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_01 = 0.01_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:: r1_2 = 1.2_r_kind
real(r_kind),parameter:: r1_02 = 1.02_r_kind
real(r_kind),parameter:: r3_33= three + one/three
real(r_kind),parameter:: r6 = 6.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:: 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:: r1200= 1200.0_r_kind
real(r_kind),parameter:: convert= 1.0e-6_r_kind
real(r_kind),parameter:: emerr= 0.2_r_kind
real(r_kind),parameter:: r0_1_bmiss=one_tenth*bmiss
real(r_kind),parameter:: r0_01_bmiss=r0_01*bmiss
character(80),parameter:: cspval= '88888888'
! integer(i_kind),parameter:: mxtb=5000000
! integer(i_kind),parameter:: nmsgmax=100000 ! max message count
! Declare local variables
logical tob,qob,uvob,spdob,sstob,pwob,psob,gustob,visob,tdob,mxtmob,mitmob,pmob,howvob,cldchob
logical metarcldobs,goesctpobs,tcamtob,lcbasob
logical outside,driftl,convobs,inflate_error
logical sfctype
logical luse,ithinp,windcorr
logical patch_fog
logical aircraftset,aircraftobs,aircraftobst,aircrafttype
logical acft_profl_file
logical,allocatable,dimension(:,:):: lmsg ! set true when convinfo entry id found in a message
character(40) drift,hdstr,qcstr,oestr,sststr,satqcstr,levstr,hdstr2
character(40) metarcldstr,goescldstr,metarvisstr,metarwthstr,cldseqstr,cld2seqstr,cldceilhstr
character(40) maxtmintstr,owavestr
character(80) obstr
character(10) date
character(8) subset
character(8) prvstr,sprvstr
character(8) c_prvstg,c_sprvstg
character(8) c_station_id
character(8) cc_station_id
character(1) sidchr(8)
character(8) stnid
character(10) aircraftstr
character(1) cb
character(1) cdummy
logical lhilbert
integer(i_kind) ireadmg,ireadsb,icntpnt,icntpnt2,icount,iiout
integer(i_kind) lunin,i,maxobs,j,idomsfc,it29,nmsgmax,mxtb
integer(i_kind) kk,klon1,klat1,klonp1,klatp1
integer(i_kind) nc,nx,isflg,ntread,itx,ii,ncsave
integer(i_kind) ihh,idd,idate,iret,im,iy,k,levs
integer(i_kind) metarcldlevs,metarwthlevs,cldseqlevs,cld2seqlevs
integer(i_kind) kx,kx0,nreal,nchanl,ilat,ilon,ithin
integer(i_kind) cat,zqm,pwq,sstq,qm,lim_qm,lim_zqm,gustqm,visqm,tdqm,mxtmqm,mitmqm,howvqm,cldchqm
integer(i_kind) lim_tqm,lim_qqm
integer(i_kind) nlevp ! vertical level for thinning
integer(i_kind) ntmp,iout
integer(i_kind) pflag,irec,zflag
integer(i_kind) ntest,nvtest,iosub,ixsub,isubsub,iobsub
integer(i_kind) kl,k1,k2,k1_ps,k1_q,k1_t,k1_uv,k1_pw,k2_q,k2_t,k2_uv,k2_pw,k2_ps
integer(i_kind) itypex,itypey
integer(i_kind) minobs,minan
integer(i_kind) ntb,ntmatch,ncx
integer(i_kind) nmsg ! message index
integer(i_kind) idx ! order index of aircraft temperature bias
integer(i_kind) tcamt_qc,lcbas_qc
integer(i_kind) low_cldamt_qc,mid_cldamt_qc,hig_cldamt_qc
integer(i_kind) iyyyymm
integer(i_kind) jj,start,next,ncount_ps,ncount_q,ncount_uv,ncount_t,ncount_pw
integer(i_kind),dimension(5):: idate5
integer(i_kind),dimension(255):: pqm,qqm,tqm,wqm,pmq
integer(i_kind),dimension(nconvtype)::ntxall
integer(i_kind),dimension(nconvtype+1)::ntx
integer(i_kind),allocatable,dimension(:):: isort,iloc,nrep
integer(i_kind),allocatable,dimension(:,:):: tab
integer(i_kind) ibfms,thisobtype_usage
integer(i_kind) iwmo,ios
integer(i_kind) ntime,itime
integer(i_kind) ierr_ps,ierr_q,ierr_t,ierr_uv,ierr_pw ! the position of error table collum
integer(i_kind) idummy1,idummy2,glret,lindx !glret>0 means GLERL code exists.Others are dummy variables
real(r_kind) time,timex,time_drift,timeobs,toff,t4dv,zeps
real(r_kind) qtflg,tdry,rmesh,ediff,usage,ediff_ps,ediff_q,ediff_t,ediff_uv,ediff_pw
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,gustoe,visoe,qmaxerr
real(r_kind) toe,woe,errout,oelev,dlat,dlon,sstoe,dlat_earth,dlon_earth
real(r_kind) tdoe,mxtmoe,mitmoe,pmoe,howvoe,cldchoe
real(r_kind) dlat_earth_deg,dlon_earth_deg
real(r_kind) selev,elev,stnelev
real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00
real(r_kind) vdisterrmax,u00,v00
real(r_kind) del,terrmin,werrmin,perrmin,qerrmin,pwerrmin,del_ps,del_q,del_t,del_uv,del_pw
real(r_kind) pjbmin,qjbmin,tjbmin,wjbmin
real(r_kind) tsavg,ff10,sfcr,zz
real(r_kind) crit1,timedif,xmesh,pmesh,pmot,ptime ! thinning parameter
real(r_kind) time_correction
real(r_kind) tcamt,lcbas,ceiling
real(r_kind) tcamt_oe,lcbas_oe
real(r_kind) low_cldamt,mid_cldamt,hig_cldamt
real(r_kind),dimension(nsig):: presl
real(r_kind),dimension(nsig-1):: dpres
real(r_kind),dimension(255)::plevs
real(r_kind),dimension(255):: tvflg
real(r_kind),allocatable,dimension(:):: presl_thin
real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out
real(r_kind) :: zob,tref,dtw,dtc,tz_tr
real(r_kind) :: tempvis,visout
real(r_kind) :: tempcldch,cldchout
real(r_kind) :: windsensht
real(r_double) rstation_id,qcmark_huge
real(r_double) vtcd,glcd !virtual temp program code and GLERL program code
real(r_double),dimension(8):: hdr,hdrtsb
real(r_double),dimension(3,255):: hdr3
real(r_double),dimension(8,255):: drfdat,qcmark,obserr,var_jb
real(r_double),dimension(13,255):: obsdat
real(r_double),dimension(8,1):: sstdat
real(r_double),dimension(2,1):: cld2seq
real(r_double),dimension(3,10):: cldseq
real(r_double),dimension(2,10):: metarcld
real(r_double),dimension(1,10):: metarwth
real(r_double),dimension(2,1) :: metarvis
real(r_double),dimension(4,1) :: goescld
real(r_double),dimension(2,255):: maxtmint
real(r_double),dimension(1,255):: owave
real(r_double),dimension(1,255):: cldceilh
real(r_double),dimension(1):: satqc
real(r_double),dimension(1,1):: r_prvstg,r_sprvstg
real(r_double),dimension(1,255):: levdat
real(r_double),dimension(255,20):: tpc
real(r_double),dimension(2,255,20):: tobaux
real(r_double),dimension(2,255):: aircraftwk
! for hilbert curve
integer(i_kind) ndata_hil,nor,ncc,nnrand
integer(i_kind) indexx
real(r_kind) dentrip,dentrip_tmp,vmin,vmax,rmesh_tmp,pmesh_tmp,prest
integer(i_kind) ntime_max,ntime_tmp,itype,ikx
integer(i_kind),dimension(24) :: ntype_arr
integer(i_kind),allocatable,dimension(:,:) :: index_arr
real(r_kind),allocatable,dimension(:,:,:) :: data_hilb
real(r_kind),allocatable,dimension(:) :: rlat_hil,rlon_hil,height,wtob,wght_hilb
! end of block
! 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)
equivalence(rstation_id,sidchr)
! data statements
data hdstr /'SID XOB YOB DHR TYP ELV SAID T29'/
data hdstr2 /'TYP SAID T29 SID'/
data obstr /'POB QOB TOB ZOB UOB VOB PWO MXGS HOVI CAT PRSS TDO PMO' /
data drift /'XDR YDR HRDR '/
data sststr /'MSST DBSS SST1 SSTQM SSTOE '/
data qcstr /'PQM QQM TQM ZQM WQM NUL PWQ PMQ'/
data oestr /'POE QOE TOE NUL WOE NUL PWE '/
! data satqcstr /'RFFL QIFY QIFN EEQF'/
data satqcstr /'QIFN'/
data prvstr /'PRVSTG'/
data sprvstr /'SPRVSTG'/
data levstr /'POB'/
data cld2seqstr /'TOCC HBLCS'/ ! total cloud cover and height above surface of base of lowest cloud seen
data cldseqstr /'VSSO CLAM HOCB'/ ! vertical significance, cloud amount and cloud base height
data metarcldstr /'CLAM HOCB'/ ! cloud amount and cloud base height
data metarwthstr /'PRWE'/ ! present weather
data metarvisstr /'HOVI TDO'/ ! visibility and dew point
data goescldstr /'CDTP TOCC GCDTT CDTP_QM'/ ! NESDIS cloud products: cloud top pressure, total cloud amount,
! cloud top temperature, cloud top temp. qc mark
data aircraftstr /'POAF IALR'/ ! phase of aircraft flight and vertical velocity
data maxtmintstr /'MXTM MITM'/
data owavestr /'HOWV'/
data cldceilhstr /'CEILING'/
data lunin / 13 /
data ithin / -9 /
data rmesh / -99.999_r_kind /
!* test new vad wind
!* for match loction station and time
! character(7*2000) cstn_idtime,cstn_idtime2
! character(7) stn_idtime(2000),stn_idtime2(2000)
! equivalence (stn_idtime(1),cstn_idtime)
! equivalence (stn_idtime2(1),cstn_idtime2)
! integer :: ii1,atmp,btmp,mytimeyy,ibyte
! character(4) stid
! real(8) :: rval
! character(len=8) :: cval
! equivalence (rval,cval)
! character(7) flnm
integer:: icase,klev,ikkk,tkk
real:: diffhgt,diffuu,diffvv
real(r_double),dimension(3,1500):: fcstdat
logical print_verbose
print_verbose=.false.
if(verbose) print_verbose=.true.
! File type
acft_profl_file = index(infile,'_profl')/=0
! Initialize variables
vdisterrmax=zero
zflag=0
nreal=0
satqc=zero
tob = obstype == 't'
uvob = obstype == 'uv' ; if (twodvar_regional) uvob = uvob .or. obstype == 'wspd10m' .or. obstype == 'uwnd10m' .or. obstype == 'vwnd10m'
spdob = obstype == 'spd'
psob = obstype == 'ps'
qob = obstype == 'q'
pwob = obstype == 'pw'
sstob = obstype == 'sst'
gustob = obstype == 'gust'
visob = obstype == 'vis'
tdob = obstype == 'td2m'
mxtmob = obstype == 'mxtm'
mitmob = obstype == 'mitm'
pmob = obstype == 'pmsl'
howvob = obstype == 'howv'
metarcldobs = obstype == 'mta_cld'
goesctpobs = obstype == 'gos_ctp'
tcamtob = obstype == 'tcamt'
lcbasob = obstype == 'lcbas'
cldchob = obstype == 'cldch'
newvad=.false.
convobs = tob .or. uvob .or. spdob .or. qob .or. gustob .or. &
tdob .or. mxtmob .or. mitmob .or. pmob .or. howvob .or. &
tcamtob .or. lcbasob .or. cldchob
aircraftobst=.false.
if(tob)then
nreal=25
else if(uvob) then
nreal=26
else if(spdob) then
nreal=24
else if(psob) then
nreal=20
else if(qob) then
nreal=26
else if(pwob) then
nreal=20
else if(sstob) then
if (nst_gsi > 0) then
nreal=18 + nstinfo
else
nreal=18
end if
else if(gustob) then
nreal=21
else if(visob) then
nreal=18
else if(tdob) then
nreal=25
else if(mxtmob) then
nreal=24
else if(mitmob) then
nreal=24
else if(pmob) then
nreal=24
else if(howvob) then
nreal=23
else if(metarcldobs) then
nreal=27
else if(goesctpobs) then
nreal=8
else if(tcamtob) then
nreal=20
else if(lcbasob) then
nreal=23
else if(cldchob) then
nreal=18
else
write(6,*) ' illegal obs type in READ_PREPBUFR ',obstype
call stop2(94)
end if
! Set qc limits based on noiqc flag
if (noiqc) then
lim_qm=8
if (psob) lim_zqm=7
if (qob.or.tdob) lim_tqm=7
if (tob) lim_qqm=8
else
lim_qm=4
if (psob) lim_zqm=4
if (qob.or.tdob) lim_tqm=4
if (tob) lim_qqm=4
endif
if (tob .and. (aircraft_t_bc_pof .or. aircraft_t_bc .or.&
aircraft_t_bc_ext )) nreal=nreal+3
if(perturb_obs .and. (tob .or. psob .or. qob))nreal=nreal+1
if(perturb_obs .and. uvob )nreal=nreal+2
qcmark_huge = huge_i_kind
lhilbert = twodvar_regional .and. hilbert_curve
if (blacklst) call blacklist_read(obstype)
terrmin=half
werrmin=one
perrmin=0.3_r_kind
qerrmin=0.05_r_kind
pwerrmin=one
tjbmin=zero
qjbmin=zero
wjbmin=zero
pjbmin=zero
!------------------------------------------------------------------------
ntread=1
ntmatch=0
ntx(ntread)=0
ntxall=0
var_jb=zero
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(obstype))then
if(.not.use_prepb_satwnd .and. (trim(ioctype(nc)) == 'uv' .or. trim(ioctype(nc)) == 'wspd10m' .or. &
trim(ioctype(nc)) == 'uwnd10m' .or. trim(ioctype(nc)) == 'vwnd10m') .and. ictype(nc) >=241 &
.and. ictype(nc) <260) then
cycle
else
if (aircraft_t_bc) then
aircrafttype=(ictype(nc) == 130 .or. ictype(nc) == 131 .or. (ictype(nc) >= 133 .and. ictype(nc)<140) .or. &
ictype(nc) == 230 .or. ictype(nc) == 231 .or. (ictype(nc) >= 233 .and. ictype(nc)<240))
if (.not. acft_profl_file .and. aircrafttype) cycle ! skip aircrafttype for prepbufr
if (acft_profl_file .and. (.not. aircrafttype)) cycle ! skip non-aircrafttype for prepbufr_profl
end if
ntmatch=ntmatch+1
ntxall(ntmatch)=nc
endif
end if
if(trim(ioctype(nc)) == trim(obstype) .and. abs(icuse(nc)) <= 1)then
if(.not.use_prepb_satwnd .and. (trim(ioctype(nc)) == 'uv' .or. trim(ioctype(nc)) == 'wspd10m' .or. &
trim(ioctype(nc)) == 'uwnd10m' .or. trim(ioctype(nc)) == 'vwnd10m' ) .and. ictype(nc) >=241 &
.and. ictype(nc) <260) then
cycle
else
if (aircraft_t_bc) then
aircrafttype=(ictype(nc) == 130 .or. ictype(nc) == 131 .or. (ictype(nc) >= 133 .and. ictype(nc)<140).or. &
ictype(nc) == 230 .or. ictype(nc) == 231 .or. (ictype(nc) >= 233 .and. ictype(nc)<240))
if (.not. acft_profl_file .and. aircrafttype) cycle ! skip aircrafttype for prepbufr
if (acft_profl_file .and. (.not. aircrafttype)) cycle ! skip non-aircrafttype for prepbufr_profl
end if
ithin=ithin_conv(nc)
if(ithin > 0 .and. ithin <5)then
ntread=ntread+1
ntx(ntread)=nc
end if
endif
end if
end do
if(ntmatch == 0)then
write(6,*) ' no matching obstype found in obsinfo ',obstype
return
end if
!! get message and subset counts
call getcount_bufr(infile,nmsgmax,mxtb)
allocate(lmsg(nmsgmax,ntread),tab(mxtb,3),nrep(nmsgmax))
lmsg = .false.
maxobs=0
tab=0
nmsg=0
nrep=0
ntb = 0
irec = 0
ncount_ps=0;ncount_q=0;ncount_t=0;ncount_uv=0;ncount_pw=0
! Open, then read date from bufr data
open(lunin,file=trim(infile),form='unformatted')
call openbf(lunin,'IN',lunin)
call datelen(10)
msg_report: do while (ireadmg(lunin,subset,idate) == 0)
irec = irec + 1
if(irec < nrec_start) cycle msg_report
if(.not.use_prepb_satwnd .and. trim(subset) == 'SATWND') cycle msg_report
if (aircraft_t_bc) then
aircraftset = trim(subset)=='AIRCFT' .or. trim(subset)=='AIRCAR'
if (.not. acft_profl_file .and. aircraftset) cycle msg_report
if (acft_profl_file .and. (.not. aircraftset)) cycle msg_report
end if
! Time offset
if(nmsg == 0) call time_4dvar(idate,toff)
nmsg=nmsg+1
if (nmsg>nmsgmax) then
write(6,*)'READ_PREPBUFR: messages exceed maximum ',nmsgmax
call stop2(50)
endif
loop_report: do while (ireadsb(lunin) == 0)
ntb = ntb+1
nrep(nmsg)=nrep(nmsg)+1
if (ntb>mxtb) then
write(6,*)'READ_PREPBUFR: reports exceed maximum ',mxtb
call stop2(50)
endif
! Extract type information
call ufbint(lunin,hdr,4,1,iret,hdstr2)
kx=hdr(1)
if (aircraft_t_bc .and. acft_profl_file) then
kx0=kx
if (.not. uvob) then
if (kx0==330 .or. kx0==430 .or. kx0==530) kx=130
if (kx0==331 .or. kx0==431 .or. kx0==531) kx=131
if (kx0==332 .or. kx0==432 .or. kx0==532) kx=132
if (kx0==333 .or. kx0==433 .or. kx0==533) kx=133
if (kx0==334 .or. kx0==434 .or. kx0==534) kx=134
if (kx0==335 .or. kx0==435 .or. kx0==535) kx=135
else
if (kx0==330 .or. kx0==430 .or. kx0==530) kx=230
if (kx0==331 .or. kx0==431 .or. kx0==531) kx=231
if (kx0==332 .or. kx0==432 .or. kx0==532) kx=232
if (kx0==333 .or. kx0==433 .or. kx0==533) kx=233
if (kx0==334 .or. kx0==434 .or. kx0==534) kx=234
if (kx0==335 .or. kx0==435 .or. kx0==535) kx=235
end if
end if
!* for new vad wind
if(kx==224 .and. .not.newvad) then
call ufbint(lunin,hdrtsb,1,1,iret,'TSB')
if(hdrtsb(1)==2) then
newvad=.true.
else
call ufbint(lunin,obsdat,13,255,levs,obstr)
if(levs>1)then
do k=1, levs-1
diffuu=abs(obsdat(4,k+1)-obsdat(4,k))
if(diffuu==50.0) then
newvad=.true.
exit
end if
end do
end if
end if
if(newvad)write(6,*)'new vad flag::', newvad
end if
!* END new vad wind
! identify drifting buoys - TYP=180/280 T29=562 and last three digits of SID between 500 and 999
! (see https://www.wmo.int/pages/prog/amp/mmop/wmo-number-rules.html) Set kx to 199/299
if (id_drifter .and. (kx==180 .or. kx==280) .and. nint(hdr(3))==562) then
rstation_id=hdr(4)
read(c_station_id,*,iostat=ios) iwmo
if (ios == 0 .and. iwmo > 0) then
if(mod(iwmo,1000) >=500) then
kx = kx + 19
end if
end if
end if
if (id_ship .and. (kx==180) .and. (nint(hdr(3))==522 .or. nint(hdr(3))==523)) then
rstation_id=hdr(4)
kx = kx + 18
end if
if(twodvar_regional)then
! If running in 2d-var (surface analysis) mode, check to see if observation
! is surface type or GOES cloud product(kx=151). If not, read next observation report from bufr file
sfctype=(kx>179.and.kx<190).or.(kx>=280.and.kx<=290).or. &
(kx>=192.and.kx<=199).or.(kx>=292.and.kx<=299) .or. &
(kx==151)
if (.not.sfctype ) cycle loop_report
end if
! temporary specify iobsub until put in bufr file
iobsub = 0
if(kx == 280 .or. kx == 180 ) iobsub=hdr(3)
if(kx == 280 .or. kx ==180) then
if ( hdr(3) >555.0_r_kind .and. hdr(3) <565.0_r_kind ) then
iobsub=00
else
iobsub=01
endif
! Set saildrone to subtype 02
if (nint(hdr(3)) == 560) iobsub = 02
endif
! Su suggested to keep both 289 and 290. But trunk only keep 290
! if(kx == 289 .or. kx == 290) iobsub=hdr(2)
if(kx == 290) iobsub=hdr(2)
if(use_prepb_satwnd .and. (kx >= 240 .and. kx <=260 )) iobsub = hdr(2)
! For the satellite wind to get quality information and check if it will be used
if(use_prepb_satwnd .and. (kx == 243 .or. kx == 253 .or. kx ==254) ) then
call ufbint(lunin,satqc,1,1,iret,satqcstr)
if(satqc(1) < 85.0_r_double) cycle loop_report ! QI w/o fcst (su's setup
!! if(satqc(2) <= 80.0_r_double) cycle loop_report ! QI w/ fcst (old prepdata)
endif
! Check for blacklisting of station ID
if (blacklst .and. ibcnt > 0) then
stnid = transfer(hdr(4),stnid)
do i = 1,ibcnt
if( kx == blkkx(i) .and. stnid == blkstns(i) ) then
write(6,*)'READ_PREPBUFR: blacklist station ',stnid, &
'for obstype ',trim(obstype),' and kx=',kx
cycle loop_report
endif
enddo
endif
! Match ob to proper convinfo type
ncsave=0
matchloop:do ncx=1,ntmatch
nc=ntxall(ncx)
if (kx /= ictype(nc))cycle
! Find convtype which match ob type and subtype
if(icsubtype(nc) == iobsub) then
ncsave=nc
exit matchloop
else
! Find convtype which match ob type and subtype group (isubtype == ?*)
! where ? specifies the group and icsubtype = ?0)
ixsub=icsubtype(nc)/10
iosub=iobsub/10
isubsub=icsubtype(nc)-ixsub*10
if(ixsub == iosub .and. isubsub == 0) then
ncsave=nc
! Find convtype which match ob type and subtype is all remaining
! (icsubtype(nc) = 0)
else if (ncsave == 0 .and. icsubtype(nc) == 0) then
ncsave=nc
end if
end if
end do matchloop
! Save information for next read
if(ncsave /= 0) then
call ufbint(lunin,levdat,1,255,levs,levstr)
maxobs=maxobs+max(1,levs)
nx=1
if(ithin_conv(ncsave) > 0 .and. ithin_conv(ncsave) <5)then
do ii=2,ntread
if(ntx(ii) == ncsave)nx=ii
end do
end if
tab(ntb,1)=ncsave
tab(ntb,2)=nx
tab(ntb,3)=levs
lmsg(nmsg,nx) = .true.
end if
end do loop_report
enddo msg_report
if (nmsg==0) then
call closbf(lunin)
close(lunin)
if(print_verbose)write(6,*)'READ_PREPBUFR: no messages/reports '
return
end if
if(print_verbose)write(6,*)'READ_PREPBUFR: messages/reports = ',nmsg,'/',ntb,' ntread = ',ntread
if(tob .and. print_verbose) write(6,*)'READ_PREPBUFR: time offset is ',toff,' hours.'
!------------------------------------------------------------------------
! Obtain program code (VTCD) associated with "VIRTMP" step
call ufbqcd(lunin,'VIRTMP',vtcd)
!see if file contains GLERL program code (GLCD)
!Obtain code if it exists. Otherwise set to missing (-999)
call status(lunin,lindx,idummy1,idummy2)
call nemtab(lindx,'GLERL',idummy1,cdummy,glret)
if (glret /= 0) then
call ufbqcd(lunin,'GLERL',glcd)
else
!warn that GLERL adjustment is not available.
print*, "WARNING: GLERL program code not in this file."
glcd=-999._r_double
endif
call init_rjlists
call init_aircraft_rjlists
if(i_gsdsfc_uselist==1) call init_gsd_sfcuselist
if (lhilbert) call init_hilbertcurve(maxobs)
if (twodvar_regional) then
call init_ndfdgrid
call init_windht_lists !load wind sensor height provider lists
endif
! loop over convinfo file entries; operate on matches
allocate(cdata_all(nreal,maxobs),isort(maxobs))
isort = 0
cdata_all=zero
nread=0
ntest=0
nvtest=0
nchanl=0
ilon=2
ilat=3
rmesh=zero
pmot=zero
pmesh=zero
ptime=zero
xmesh=zero
pflag=0
loop_convinfo: do nx=1, ntread
use_all_tm = .true.
use_all = .true.
ithin=0
if(nx > 1) then
nc=ntx(nx)
ithin=ithin_conv(nc)
if (ithin > 0 .and. ithin <5) then
rmesh=rmesh_conv(nc)
pmesh=pmesh_conv(nc)
pmot=pmot_conv(nc)
ptime=ptime_conv(nc)
if(pmesh > zero .and. ithin ==1) then
pflag=1
zflag=-1
nlevp=r1200/pmesh
else if(pmesh > zero .and. ithin ==2) then
pflag=1
zflag=1
nlevp=25000.00_r_kind/pmesh
else
zflag=-1
pflag=0
nlevp=nsig
endif
xmesh=rmesh
if( ptime >zero) then
use_all_tm = .false.
ntime=6.0_r_kind/ptime !! 6 hour winddow
call make3grids_tm(xmesh,nlevp,ntime)
allocate(presl_thin(nlevp))
if (zflag==-1 ) then
do k=1,nlevp
presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth
enddo
else if(zflag==1 ) then
do k=1,nlevp
presl_thin(k)=k*pmesh
enddo
endif
else
use_all = .false.
call make3grids(xmesh,nlevp)
allocate(presl_thin(nlevp))
if (zflag==-1 ) then
do k=1,nlevp
presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth
enddo
else if(zflag==1 ) then
do k=1,nlevp
presl_thin(k)=k*pmesh
enddo
endif
endif
if(print_verbose) write(6,*)'READ_PREPBUFR: at line 779: obstype,ictype(nc),rmesh,pflag,nlevp,pmesh,pmot,ptime=',&
trim(ioctype(nc)),ictype(nc),rmesh,pflag,nlevp,pmesh,pmot,ptime,ithin
endif
endif
call closbf(lunin)
close(lunin)
open(lunin,file=infile,form='unformatted')
call openbf(lunin,'IN',lunin)
call datelen(10)
! Big loop over prepbufr file
ntb = 0
nmsg = 0
icntpnt=0
icntpnt2=0
disterrmax=-9999.0_r_kind
irec = 0
loop_msg: do while (ireadmg(lunin,subset,idate)== 0)
irec = irec + 1
if(irec < nrec_start) cycle loop_msg
if(.not.use_prepb_satwnd .and. trim(subset) =='SATWND') cycle loop_msg
if (aircraft_t_bc) then
aircraftset = trim(subset)=='AIRCFT' .or. trim(subset)=='AIRCAR'
if (.not. acft_profl_file .and. aircraftset) cycle loop_msg
if (acft_profl_file .and. (.not. aircraftset)) cycle loop_msg
end if
nmsg = nmsg+1
if(.not.lmsg(nmsg,nx)) then
do i=ntb+1,ntb+nrep(nmsg)
icntpnt2=icntpnt2+tab(i,3)
end do
ntb=ntb+nrep(nmsg)
cycle loop_msg ! no useable reports this mesage, skip ahead report count
end if
loop_readsb: do while(ireadsb(lunin) == 0)
! use msg lookup table to decide which messages to skip
! use report id lookup table to only process matching reports
ntb = ntb+1
if(icntpnt < icntpnt2)icntpnt=icntpnt2
icntpnt2=icntpnt2+tab(ntb,3)
nc=tab(ntb,1)
if(nc <= 0 .or. tab(ntb,2) /= nx) cycle loop_readsb
! Extract type, date, and location information
call ufbint(lunin,hdr,8,1,iret,hdstr)
kx=hdr(5)
if (.not.(aircraft_t_bc .and. acft_profl_file)) then
if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_readsb
if(hdr(2)== r360)hdr(2)=hdr(2)-r360
if(hdr(2) < zero)hdr(2)=hdr(2)+r360
dlon_earth_deg=hdr(2)
dlat_earth_deg=hdr(3)
dlon_earth=hdr(2)*deg2rad
dlat_earth=hdr(3)*deg2rad
!
! identify drifting buoys - TYP=180/280 T29=562 and last three digits of SID between 500 and 999
! (see https://www.wmo.int/pages/prog/amp/mmop/wmo-number-rules.html) Set kx to 199/299
if (id_drifter .and. (kx==180 .or. kx==280) .and. nint(hdr(8))==562 ) then
rstation_id=hdr(1)
read(c_station_id,*,iostat=ios) iwmo
if (ios == 0 .and. iwmo > 0) then
if(mod(iwmo,1000) >=500) then
kx = kx + 19
end if
end if
end if
if (id_ship .and. (kx==180) .and. (nint(hdr(8))==522 .or. nint(hdr(8))==523) ) then
rstation_id=hdr(1)
kx = kx + 18
end if
!
! check VAD subtype. 1--old, 2--new, other--old
if(kx==224) then
call ufbint(lunin,hdrtsb,1,1,iret,'TSB')
if(.not.newvad .and. hdrtsb(1)==2) cycle loop_readsb
if(newvad .and. hdrtsb(1)/=2) cycle loop_readsb
end if
!* thin new VAD in time level
if(kx==224.and.newvad)then
icase=0
if ( vad_near_analtime ) then
if(abs(hdr(4))<0.25_r_kind) icase=1
else
if(abs(hdr(4))>0.17_r_kind.and.abs(hdr(4))<0.32_r_kind) icase=1
if(abs(hdr(4))>0.67_r_kind.and.abs(hdr(4))<0.82_r_kind) icase=1
if(abs(hdr(4))>1.17_r_kind.and.abs(hdr(4))<1.32_r_kind) icase=1
if(abs(hdr(4))>1.67_r_kind.and.abs(hdr(4))<1.82_r_kind) icase=1
if(abs(hdr(4))>2.17_r_kind.and.abs(hdr(4))<2.62_r_kind) icase=1
if(abs(hdr(4))>2.67_r_kind.and.abs(hdr(4))<2.82_r_kind) icase=1
endif
if(icase/=1) cycle
end if
if(regional)then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) ! convert to rotated coordinate
if(diagnostic_reg) then
call txy2ll(dlon,dlat,rlon00,rlat00)
ntest=ntest+1
cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* &
(sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00))
cdist=max(-one,min(cdist,one))
disterr=acos(cdist)*rad2deg
disterrmax=max(disterrmax,disterr)
end if
if(outside) cycle loop_readsb ! check to see if outside regional domain
else
dlat = dlat_earth
dlon = dlon_earth
call grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
else
call ufbint(lunin,hdr3,3,255,levs,'XDR YDR HRDR')
kx0=kx
if (.not. uvob) then
if (kx0==330 .or. kx0==430 .or. kx0==530) kx=130
if (kx0==331 .or. kx0==431 .or. kx0==531) kx=131
if (kx0==332 .or. kx0==432 .or. kx0==532) kx=132
if (kx0==333 .or. kx0==433 .or. kx0==533) kx=133
if (kx0==334 .or. kx0==434 .or. kx0==534) kx=134
if (kx0==335 .or. kx0==435 .or. kx0==535) kx=135
else
if (kx0==330 .or. kx0==430 .or. kx0==530) kx=230
if (kx0==331 .or. kx0==431 .or. kx0==531) kx=231
if (kx0==332 .or. kx0==432 .or. kx0==532) kx=232
if (kx0==333 .or. kx0==433 .or. kx0==533) kx=233
if (kx0==334 .or. kx0==434 .or. kx0==534) kx=234
if (kx0==335 .or. kx0==435 .or. kx0==535) kx=235
end if
endif
!------------------------------------------------------------------------
if(offtime_data) then
! in time correction for observations to account for analysis
! time being different from obs file time.
write(date,'( i10)') idate
read (date,'(i4,3i2)') iy,im,idd,ihh
idate5(1)=iy
idate5(2)=im
idate5(3)=idd
idate5(4)=ihh
idate5(5)=0
call w3fs21(idate5,minobs) ! obs ref time in minutes relative to historic date
idate5(1)=iadate(1)
idate5(2)=iadate(2)
idate5(3)=iadate(3)
idate5(4)=iadate(4)
idate5(5)=0
call w3fs21(idate5,minan) ! analysis ref time in minutes relative to historic date
! Add obs reference time, then subtract analysis time to get obs time relative to analysis
time_correction=float(minobs-minan)*r60inv
else
time_correction=zero
end if
if (.not. (aircraft_t_bc .and. acft_profl_file)) then
timeobs=real(real(hdr(4),r_single),r_double)
t4dv=timeobs + toff
zeps=1.0e-8_r_kind
if (t4dv<zero .and.t4dv> -zeps) t4dv=zero
if (t4dv>winlen.and.t4dv<winlen+zeps) t4dv=winlen
t4dv=t4dv + time_correction
time=timeobs + time_correction
end if
if(use_prepb_satwnd .and. (kx >= 240 .and. kx <= 260)) iobsub = hdr(7)
! Balloon drift information available for these data
driftl=kx==120.or.kx==220.or.kx==221
if (.not. (aircraft_t_bc .and. acft_profl_file)) then
if (l4dvar.or.l4densvar) then
if ((t4dv<zero.OR.t4dv>winlen) .and. .not.driftl) cycle loop_readsb ! outside time window
else
if((real(abs(time)) > real(ctwind(nc)) .or. real(abs(time)) > real(twindin)) &
.and. .not. driftl)cycle loop_readsb ! outside time window
endif
timex=time
end if
! If ASCAT data, determine primary surface type. If not open sea,
! skip this observation. This check must be done before thinning.
if (kx==290 .or. kx==289 .or. kx==285) then
call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg)
if (isflg /= 0) cycle loop_readsb
if (tsavg <= 273.0_r_kind) cycle loop_readsb
endif
sfctype=(kx>179.and.kx<190).or.(kx>=280.and.kx<=290).or. &
(kx>=192.and.kx<=199).or.(kx>=292.and.kx<=299)
if (sfctype) then
call ufbint(lunin,r_prvstg,1,1,iret,prvstr)
call ufbint(lunin,r_sprvstg,1,1,iret,sprvstr)
else if(kx == 120 .and. tob .or. qob .or. psob)then
c_prvstg=cspval
c_sprvstg='PREP'
else if(kx == 220 .and. uvob)then
c_prvstg=cspval
c_sprvstg='PREP'
else
c_prvstg=cspval
c_sprvstg=cspval
endif
! Extract data information on levels
call ufbint(lunin,obsdat,13,255,levs,obstr)
if (twodvar_regional) then
if (mxtmob .or. mitmob) call ufbint(lunin,maxtmint,2,255,levs,maxtmintstr)
if (howvob) call ufbint(lunin,owave,1,255,levs,owavestr)
if (cldchob) call ufbint(lunin,cldceilh,1,255,levs,cldceilhstr)
endif
if(kx==224 .and. newvad) then
call ufbint(lunin,fcstdat,3,255,levs,'UFC VFC TFC ')
end if
call ufbint(lunin,qcmark,8,255,levs,qcstr)
call ufbint(lunin,obserr,8,255,levs,oestr)
call ufbevn(lunin,tpc,1,255,20,levs,'TPC')
! If available, get obs errors from error table
if(oberrflg .and. kx<= 300)then
! Set lower limits for observation errors
terrmin=half
werrmin=one
perrmin=0.3_r_kind
qerrmin=0.05_r_kind
pwerrmin=one
tjbmin=zero
qjbmin=zero
wjbmin=zero
pjbmin=zero
itypey=kx
if( njqc) then
if (psob) then
itypex=itypey
ierr_ps=0
do i =1,maxsub_ps
if(icsubtype(nc)==isuble_ps(itypex,i)) then
ierr_ps=i+1
exit
else if(i== maxsub_ps .and. icsubtype(nc) /= isuble_ps(itypex,i)) then
ncount_ps=ncount_ps+1
do j=1,maxsub_ps
if(isuble_ps(itypex,j) ==0 ) then
ierr_ps=j+1
exit
endif
enddo
if (ncount_ps ==1) then
write(6,*) 'READ_PREPBUFR: WARNING!!psob: cannot find subtyep in the &
error table,itype,iosub=',itypex,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0, error table column=',ierr_ps
endif
endif
enddo
do k=1,levs
ppb=obsdat(1,k)
cat=nint(min(obsdat(10,k),qcmark_huge))
if ( cat /=0 ) cycle
ppb=max(zero,min(ppb,r2000))
if(ppb>=etabl_ps(itypex,1,1)) k1_ps=1
do kl=1,32
if(ppb>=etabl_ps(itypex,kl+1,1).and.ppb<=etabl_ps(itypex,kl,1)) k1_ps=kl
end do
if(ppb<=etabl_ps(itypex,33,1)) k1_ps=5
k2_ps=k1_ps+1
ediff_ps = etabl_ps(itypex,k2_ps,1)-etabl_ps(itypex,k1_ps,1)
if (abs(ediff_ps) > tiny_r_kind) then
del_ps = (ppb-etabl_ps(itypex,k1_ps,1))/ediff_ps
else
del_ps = huge_r_kind
endif
del_ps=max(zero,min(del_ps,one))
if(oberrflg)then
! write(6,*) 'READ_PREPBUFR_PS:',itypex,k1_ps,ierr_ps,k2_ps,ierr_ps
obserr(1,k)=(one-del_ps)*etabl_ps(itypex,k1_ps,ierr_ps)+del_ps*etabl_ps(itypex,k2_ps,ierr_ps)
! Surface pressure error
obserr(1,k)=max(obserr(1,k),perrmin)
endif
! Surface pressure b
var_jb(1,k)=(one-del_ps)*btabl_ps(itypex,k1_ps,ierr_ps)+del_ps*btabl_ps(itypex,k2_ps,ierr_ps)
var_jb(1,k)=max(var_jb(1,k),pjbmin)
if (var_jb(1,k) >=10.0_r_kind) var_jb(1,k)=zero
enddo
endif
if (tob) then
itypex=itypey
ierr_t=0
do i =1,maxsub_t
if( icsubtype(nc) == isuble_t(itypex,i) ) then
ierr_t=i+1
exit
else if( i == maxsub_t .and. icsubtype(nc) /= isuble_t(itypex,i)) then
ncount_t=ncount_t+1
do j=1,maxsub_t
if(isuble_t(itypex,j) ==0 ) then
ierr_t=j+1
exit
endif
enddo
if( ncount_t ==1) then
write(6,*) 'READ_PREPBUFR,WARNING!! tob:cannot find subtyep in the error,&
table,itype,iosub=',itypex,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0,error table column=',ierr_t
endif
endif
enddo
do k=1,levs
ppb=obsdat(1,k)
if(kx==153)ppb=obsdat(11,k)*0.01_r_kind
ppb=max(zero,min(ppb,r2000))
if(ppb>=etabl_t(itypex,1,1)) k1_t=1
do kl=1,32
if(ppb>=etabl_t(itypex,kl+1,1).and.ppb<=etabl_t(itypex,kl,1)) k1_t=kl
end do
if(ppb<=etabl_t(itypex,33,1)) k1_t=5
k2_t=k1_t+1
ediff_t = etabl_t(itypex,k2_t,1)-etabl_t(itypex,k1_t,1)
if (abs(ediff_t) > tiny_r_kind) then
del_t = (ppb-etabl_t(itypex,k1_t,1))/ediff_t
else
del_t = huge_r_kind
endif
del_t=max(zero,min(del_t,one))
! Temperature error
if(oberrflg)then
! write(6,*) 'READ_PREPBUFR_T:',itypex,k1_t,itypey,k2_t,ierr_t,nc,kx,ppb
obserr(3,k)=(one-del_t)*etabl_t(itypex,k1_t,ierr_t)+del_t*etabl_t(itypex,k2_t,ierr_t)
obserr(3,k)=max(obserr(3,k),terrmin)
endif
!Temperature b
var_jb(3,k)=(one-del_t)*btabl_t(itypex,k1_t,ierr_t)+del_t*btabl_t(itypex,k2_t,ierr_t)
var_jb(3,k)=max(var_jb(3,k),tjbmin)
if (var_jb(3,k) >=10.0_r_kind) var_jb(3,k)=zero
enddo
endif
if (qob) then
itypex=itypey
ierr_q=0
do i =1,maxsub_q
if( icsubtype(nc) == isuble_q(itypex,i) ) then
ierr_q=i+1
exit
else if( i == maxsub_q .and. icsubtype(nc) /= isuble_q(itypex,i)) then
ncount_q=ncount_q+1
do j=1,maxsub_q
if(isuble_q(itypex,j) ==0 ) then
ierr_q=j+1
exit
endif
enddo
if(ncount_q ==1 ) then
write(6,*) 'READ_PREPBUFR,WARNING!! qob:cannot find subtyep in the &
error table,itype,iosub=',itypex,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0,error table column=',ierr_q
endif
endif
enddo
do k=1,levs
ppb=obsdat(1,k)
if(kx==153)ppb=obsdat(11,k)*0.01_r_kind
ppb=max(zero,min(ppb,r2000))
if(ppb>=etabl_q(itypex,1,1)) k1_q=1
do kl=1,32
if(ppb>=etabl_q(itypex,kl+1,1).and.ppb<=etabl_q(itypex,kl,1)) k1_q=kl
end do
if(ppb<=etabl_q(itypex,33,1)) k1_q=5
k2_q=k1_q+1
ediff_q = etabl_q(itypex,k2_q,1)-etabl_q(itypex,k1_q,1)
if (abs(ediff_q) > tiny_r_kind) then
del_q = (ppb-etabl_q(itypex,k1_q,1))/ediff_q
else
del_q = huge_r_kind
endif
del_q=max(zero,min(del_q,one))
! Humidity error
if(oberrflg)then
! write(6,*) 'READ_PREPBUFR_Q:',itypex,k1_q,itypey,k2_q,ierr_q,nc,kx,ppb
obserr(2,k)=(one-del_q)*etabl_q(itypex,k1_q,ierr_q)+del_q*etabl_q(itypex,k2_q,ierr_q)
obserr(2,k)=max(obserr(2,k),qerrmin)
endif
!Humidity b
var_jb(2,k)=(one-del_q)*btabl_q(itypex,k1_q,ierr_q)+del_q*btabl_q(itypex,k2_q,ierr_q)
var_jb(2,k)=max(var_jb(2,k),qjbmin)
if (var_jb(2,k) >=10.0_r_kind) var_jb(2,k)=zero
! if(itypey==120 ) then
! write(6,*) 'READ_PREPBUFR:120_q,obserr,var_jb=',obserr(2,k),var_jb(2,k),ppb
! endif
enddo
endif
if (uvob) then
itypex=itypey
ierr_uv=0
do i =1,maxsub_uv
if( icsubtype(nc) == isuble_uv(itypex,i) ) then
ierr_uv=i+1
exit
else if( i == maxsub_uv .and. icsubtype(nc) /= isuble_uv(itypex,i)) then
ncount_uv=ncount_uv+1
do j=1,maxsub_uv
if(isuble_uv(itypex,j) ==0 ) then
ierr_uv=j+1
exit
endif
enddo
if( ncount_uv == 1) then
write(6,*) 'READ_PREPBUFR,WARNING!! uvob:cannot find subtyep in the error,&
table,itype,iosub=',itypex,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0,error table column=',ierr_uv
endif
endif
enddo
do k=1,levs
ppb=obsdat(1,k)
if(kx==153)ppb=obsdat(11,k)*0.01_r_kind
ppb=max(zero,min(ppb,r2000))
if(ppb>=etabl_uv(itypex,1,1)) k1_uv=1
do kl=1,32
if(ppb>=etabl_uv(itypex,kl+1,1).and.ppb<=etabl_uv(itypex,kl,1)) k1_uv=kl
end do
if(ppb<=etabl_uv(itypex,33,1)) k1_uv=5
k2_uv=k1_uv+1
ediff_uv = etabl_uv(itypex,k2_uv,1)-etabl_uv(itypex,k1_uv,1)
if (abs(ediff_uv) > tiny_r_kind) then
del_uv = (ppb-etabl_uv(itypex,k1_uv,1))/ediff_uv
else
del_uv = huge_r_kind
endif
del_uv=max(zero,min(del_uv,one))
! Wind error
! write(6,*) 'READ_PREPBUFR_UV:',itypex,k1_uv,itypey,k2_uv,ierr_uv,nc,kx,ppb
obserr(5,k)=(one-del_uv)*etabl_uv(itypex,k1_uv,ierr_uv)+del_uv*etabl_uv(itypex,k2_uv,ierr_uv)
obserr(5,k)=max(obserr(5,k),werrmin)
!Wind b
var_jb(5,k)=(one-del_uv)*btabl_uv(itypex,k1_uv,ierr_uv)+del_uv*btabl_uv(itypex,k2_uv,ierr_uv)
var_jb(5,k)=max(var_jb(5,k),wjbmin)
if (var_jb(5,k) >=10.0_r_kind) var_jb(5,k)=zero
! if(itypey==220) then
! write(6,*) 'READ_PREPBUFR:220_uv,obserr,var_jb=',obserr(5,k),var_jb(5,k),ppb,k2_uv,del_uv
! endif
enddo
endif
if (pwob) then
itypex=itypey
ierr_pw=0
do i =1,maxsub_pw
if(icsubtype(nc) == isuble_pw(itypex,i) ) then
ierr_pw=i+1
exit
else if( i == maxsub_pw .and. icsubtype(nc) /= isuble_pw(itypex,i)) then
ncount_pw=ncount_pw+1
do j=1,maxsub_pw
if(isuble_pw(itypex,j) ==0 ) then
ierr_pw=j+1
exit
endif
enddo
if(ncount_pw ==1 ) then
write(6,*) 'READ_PREPBUFR,WARNING!! pwob:cannot find subtyep in the error,&
table,itypex,iosub=',itypex,icsubtype(nc)
write(6,*) 'read error table at colomn subtype as 0,error table column=',ierr_pw
endif
endif
enddo
do k=1,levs
ppb=obsdat(1,k)
if(kx==153)ppb=obsdat(11,k)*0.01_r_kind
ppb=max(zero,min(ppb,r2000))
if(ppb>=etabl_pw(itypex,1,1)) k1_pw=1
do kl=1,32
if(ppb>=etabl_pw(itypex,kl+1,1).and.ppb<=etabl_pw(itypex,kl,1)) k1_pw=kl
end do
if(ppb<=etabl_pw(itypex,33,1)) k1_pw=5
k2_pw=k1_pw+1
ediff_pw = etabl_pw(itypex,k2_pw,1)-etabl_pw(itypex,k1_pw,1)
if (abs(ediff_pw) > tiny_r_kind) then
del_pw = (ppb-etabl_pw(itypex,k1_pw,1))/ediff_pw
else
del_pw = huge_r_kind
endif
del_pw=max(zero,min(del_pw,one))
if(oberrflg)then
! Precip water error
! write(6,*) 'READ_PREPBUFR_Pw:',itypex,itypey,ierr_pw,k2_pw,ierr_pw,nc,kx,ppb
obserr(7,k)=(one-del_pw)*etabl_pw(itypex,k1_pw,ierr_pw)+del_pw*etabl_pw(itypex,k2_pw,ierr_pw)
obserr(7,k)=max(obserr(7,k),pwerrmin)
endif
enddo
endif
else
do k=1,levs
itypex=kx
ppb=obsdat(1,k)
if(kx==153)ppb=obsdat(11,k)*0.01_r_kind
ppb=max(zero,min(ppb,r2000))
if(ppb>=etabl(itypex,1,1)) k1=1
do kl=1,32
if(ppb>=etabl(itypex,kl+1,1).and.ppb<=etabl(itypex,kl,1)) k1=kl
end do
if(ppb<=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 ! endif for njqc
endif ! endif for oberrflg
! If data with drift position, get drift information
if(driftl)call ufbint(lunin,drfdat,8,255,iret,drift)
! raob level enhancement on temp and q obs
if(ext_sonde .and. kx==120) call sonde_ext(obsdat,tpc,qcmark,obserr,drfdat,levs,kx,vtcd)
nread=nread+levs
aircraftobs = (kx==130) .or. (kx==131) .or. (kx>=133 .and. kx<140) .or. &
(kx==230) .or. (kx==231) .or. (kx>=233 .and. kx<240)
aircraftobst = .false.
if(uvob)then
nread=nread+levs
else if(tob) then
! aircraft temperature data
! aircraftobst = kx>129.and.kx<140
aircraftobst = (kx==131) .or. (kx>=133 .and. kx<=135) .or. (kx==130) ! for currently known types
aircraftwk = bmiss
if (aircraftobst) then
if (aircraft_t_bc) then
call ufbint(lunin,aircraftwk,2,255,levs,aircraftstr)
if (kx0>=330 .and. kx0<340) aircraftwk(2,:) = zero
else if (aircraft_t_bc_pof) then
call ufbint(lunin,aircraftwk,2,255,levs,aircraftstr)
aircraftwk(2,:) = bmiss
if (kx==130) aircraftwk(1,:) = 3.0_r_kind
else if (aircraft_t_bc_ext) then
call ufbint(lunin,aircraftwk,2,255,levs,aircraftstr)
aircraftwk(2,:) = bmiss
end if
end if
else if(sstob)then
sstdat=bmiss
call ufbint(lunin,sstdat,8,1,levs,sststr)
else if(metarcldobs) then
metarcld=bmiss
metarwth=bmiss
metarvis=bmiss
call ufbint(lunin,metarcld,2,10,metarcldlevs,metarcldstr)
call ufbint(lunin,metarwth,1,10,metarwthlevs,metarwthstr)
call ufbint(lunin,metarvis,2,1,iret,metarvisstr)
if(levs /= 1 ) then
write(6,*) 'READ_PREPBUFR: error in Metar observations, levs sould be 1 !!!'
call stop2(110)
endif
else if(goesctpobs) then
goescld=bmiss
call ufbint(lunin,goescld,4,1,levs,goescldstr)
else if (visob) then
metarwth=bmiss
call ufbint(lunin,metarwth,1,10,metarwthlevs,metarwthstr)
else if(tcamtob .or. lcbasob) then
if (trim(subset) == 'GOESND') then
goescld=bmiss
call ufbint(lunin,goescld,4,1,levs,goescldstr)
if (all(goescld==bmiss)) cycle
else
cldseq=bmiss
metarwth=bmiss
cld2seq =bmiss
call ufbint(lunin,cldseq,3,10,cldseqlevs,cldseqstr)
call ufbrep(lunin,cld2seq,2,1,cld2seqlevs,cld2seqstr)
call ufbint(lunin,metarwth,1,10,metarwthlevs,metarwthstr)
if (all(cldseq==bmiss) .and. all(cld2seq==bmiss) .and. all(metarwth==bmiss)) cycle
endif
endif
! Set station ID
rstation_id=hdr(1)
if ((kx==133 .or. kx==233) .and. c_station_id=='ACARS') then
call ufbint(lunin,hdr,1,1,iret,'ACID')
if(ibfms(hdr(1))==0) rstation_id=hdr(1)
end if
! Check for valid reported pressure (POB). Set POB=bmiss if POB<tiny_r_kind
do k=1,levs
if (obsdat(1,k)<tiny_r_kind) then
write(6,*)'READ_PREPBUFR: ***WARNING*** invalid pressure pob=',&
obsdat(1,k),' at k=',k,' for obstype=',obstype,' kx=',kx,&
' c_station_id=',c_station_id,' reset pob=',bmiss
obsdat(1,k)=bmiss
endif
end do
! Determine tail number for aircraft temperature data
idx = 0
iyyyymm = iadate(1)*100+iadate(2)
if (tob)then
if (aircraftobst .and. (aircraft_t_bc_pof .or. &
aircraft_t_bc .or. aircraft_t_bc_ext)) then
! Determine if the tail number is included in the taillist
! special treatment since kx130 has only flight NO. info, no
! aircraft type info
if (kx==130) then
cc_station_id = 'KX130'
else
cc_station_id = c_station_id
end if
cb = cc_station_id(1:1)
do j=1,nsort
if (cb==itail_sort(j)) then
start = idx_sort(j)
if (j==nsort) then
next = ntail
else
next=idx_sort(j+1)-1
end if
do jj=start,next
if (trim(cc_station_id)==trim(taillist(jj))) then
idx = jj
if (timelist(jj)/=iyyyymm) timelist(jj) = iyyyymm
exit
end if
end do
end if
end do
if (idx==0 .and. ntail_update>ntail) then
do j = ntail+1,ntail_update
if (cc_station_id == trim(taillist(j))) then
idx = j
exit
end if
end do
end if
! Append new tail number at the end of existing tail numbers.
! At 1st analysis, the obs will be used without bias correction,
! patch new tail number;
! At 2nd analysis, bias coefficients will be generated for this new
! tail number.
if (idx == 0) then
ntail_update = ntail_update+1
! print*, cc_station_id, ' ntail_update=',ntail_update,'
! ntail=',ntail
if (ntail_update > max_tail) then
write(6,*)'READ_PREPBUFR: ***ERROR*** tail number exceeds maximum'
write(6,*)'READ_PREPBUFR: stop program execution'
call stop2(341)
end if
idx_tail(ntail_update) = ntail_update
taillist(ntail_update) = cc_station_id
timelist(ntail_update) = iyyyymm
do j = 1,npredt
predt(j,ntail_update) = zero
end do
end if
! Re-set idx if idx>ntail
if (idx>ntail) idx = 0
end if
end if
! Loop over levels
do k=1,levs
do i=1,8
qcmark(i,k) = min(qcmark(i,k),qcmark_huge)
end do
plevs(k)=one_tenth*obsdat(1,k) ! convert mb to cb
if (kx == 290) plevs(k)=101.0_r_kind ! Assume 1010 mb = 101.0 cb
if (goesctpobs) plevs(k)=goescld(1,k)/1000.0_r_kind ! cloud top pressure in 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))
pmq(k)=nint(qcmark(8,k))
end do
! If temperature ob, extract information regarding virtual
! versus sensible temperature
if(tob) then
if (.not. twodvar_regional .or. .not.tsensible) then
do k=1,levs
tvflg(k)=one ! initialize as sensible
do j=1,20
if (tpc(k,j)==vtcd) tvflg(k)=zero ! reset flag if virtual
if (tpc(k,j)>=bmiss) exit ! end of stack
end do
end do
else
!look for GLERL-adjusted ob first in events stack. If not there,
!peel back events to store sensible temp in case temp is virtual
call ufbevn(lunin,tobaux,2,255,20,levs,'TOB TQM')
do k=1,levs
tvflg(k)=one ! initialize as sensible
do j=1,20
if (glret /= 0) then !GLERL adjusted obs possible
if (tpc(k,j)==glcd) then !found GLERL ob - use that and jump out of events stack
obsdat(3,k)=tobaux(1,k,j)
qcmark(3,k)=min(tobaux(2,k,j),qcmark_huge)
tqm(k)=nint(qcmark(3,k))
exit
endif
endif
if (tpc(k,j)==vtcd) then
obsdat(3,k)=tobaux(1,k,j+1)
qcmark(3,k)=min(tobaux(2,k,j+1),qcmark_huge)
tqm(k)=nint(qcmark(3,k))
end if
if (tpc(k,j)>=bmiss) exit ! end of stack
end do
end do
end if
end if
if(i_gsdqc==2) then
! AMV acceptance for all obs (E. James)
if (kx >= 240 .and. kx <= 260) then
do k=1,levs
pqm(k)=2
wqm(k)=2
end do
end if
! END of the AMV acceptance section (E. James)
! USE q from 300-10 mb for aircraft and raobs (E. James)
if(qob .and. (kx==120 .or. kx==131 .or. kx==133 .or. kx==134)) then
do k=1,levs
if( plevs(k)<=30.0_r_kind .and. plevs(k)>=1.0_r_kind ) then
if(qqm(k) == 9) qqm(k)=2
endif
end do
endif
! END use q from 300-10 mb
endif
stnelev=hdr(6)
ithin=ithin_conv(nc)
ithinp = ithin > 0 .and. ithin <5 .and. pflag /= 0
if(.not. (driftl .or. (aircraft_t_bc .and. acft_profl_file)) .and. &
(((tob .or. qob .or. uvob).and. levs > 1) .or. ithinp))then
! Interpolate guess pressure profile to observation location
klon1= int(dlon); klat1= int(dlat)
dx = dlon-klon1; dy = dlat-klat1
dx1 = one-dx; dy1 = one-dy
w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy
klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon)
if (klon1==0) klon1=nlon
klatp1=min(nlat,klat1+1); klonp1=klon1+1
if (klonp1==nlon+1) klonp1=1
do kk=1,nsig
presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) + &
w10*prsl_full(klatp1,klon1 ,kk) + &
w01*prsl_full(klat1 ,klonp1,kk) + &
w11*prsl_full(klatp1,klonp1,kk)
end do
! Compute depth of guess pressure layersat observation location
if (.not.twodvar_regional .and. levs > 1) then
do kk=1,nsig-1
dpres(kk)=presl(kk)-presl(kk+1)
end do
endif
end if
LOOP_K_LEVS: do k=1,levs
if( zflag ==-1) then
ppb=obsdat(1,k)*one_tenth
else if(zflag ==1) then
ppb=obsdat(4,k)
endif
if(kx==224 .and. newvad)then
if(mod(k,6)/=0) cycle LOOP_K_LEVS
end if
icntpnt=icntpnt+1
! Extract quality marks
if(tob)then
qm=tqm(k)
else if(uvob) then
qm=wqm(k)
else if(spdob) then
qm=wqm(k)
else if(psob) then
qm=pqm(k)
else if(qob) then
if(obsdat(2,k) > r0_01_bmiss)cycle loop_k_levs
qm=qqm(k)
else if(pwob) then
pwq=nint(qcmark(7,k))
qm=pwq
else if(sstob) then
sstq=100
if (k==1) sstq=nint(min(sstdat(4,k),qcmark_huge))
qm=sstq
else if(gustob) then
gustqm=0
if (kx==188 .or. kx==288 .or. kx==195 .or. kx==295 ) &
call get_gustqm(kx,c_station_id,c_prvstg,c_sprvstg,gustqm)
qm=gustqm
else if(visob) then
visqm=0 ! need to fix this later
qm=visqm
else if(tdob) then
if(obsdat(12,k) > r0_01_bmiss)cycle loop_k_levs
tdqm=qqm(k)
qm=tdqm
else if(mxtmob) then
mxtmqm=0
qm=mxtmqm
else if(mitmob) then
mitmqm=0
qm=mitmqm
else if(pmob) then
qm=pmq(k)
else if(howvob) then
howvqm=0
qm=howvqm
else if(cldchob) then
cldchqm=0
qm=cldchqm
else if(metarcldobs) then
qm=0
else if(goesctpobs) then
qm=0
else if(tcamtob) then
qm=0
if (kx==151)pqm=0 !Make sure GOESND data are not rejected due to the pressure quality mark
else if(lcbasob) then
qm=0
if (kx==151)pqm=0 !Make sure GOESND data are not rejected due to the pressure quality mark
end if
! Check qc marks to see if obs should be processed or skipped
if (visob) then
if (obsdat(9,k) > r0_1_bmiss) then
patch_fog=(metarwth(1,1)>= 40.0_r_kind .and. metarwth(1,1)<= 49.0_r_kind) .or. &
(metarwth(1,1)>=130.0_r_kind .and. metarwth(1,1)<=135.0_r_kind) .or. &
(metarwth(1,1)>=241.0_r_kind .and. metarwth(1,1)<=246.0_r_kind)
if (patch_fog) obsdat(9,k)=1000.0_r_kind
if (metarwth(1,1)==247.0_r_kind) obsdat(9,k)=75.0_r_kind
if (metarwth(1,1)==248.0_r_kind) obsdat(9,k)=45.0_r_kind
if (metarwth(1,1)==249.0_r_kind) obsdat(9,k)=15.0_r_kind
end if
end if
if (psob) then
cat=nint(min(obsdat(10,k),qcmark_huge))
if ( cat /=0 ) cycle loop_k_levs
if ( obsdat(1,k)< r500) qm=100
zqm=nint(qcmark(4,k))
if (zqm>=lim_zqm .and. zqm/=15 .and. zqm/=9) qm=9
endif
! if(convobs .and. pqm(k) >=lim_qm .and. qm/=15 .and. qm/=9 )cycle loop_k_levs
! if(qm >=lim_qm .and. qm /=15 .and. qm /=9)cycle loop_k_levs
if(qm > 15 .or. qm < 0) cycle loop_k_levs
! extract aircraft profile information
if (aircraft_t_bc .and. acft_profl_file) then
if (nint(obsdat(10,k))==7) cycle LOOP_K_LEVS
if(abs(hdr3(2,k))>r90 .or. abs(hdr3(1,k))>r360) cycle LOOP_K_LEVS
if(hdr3(1,k)== r360)hdr3(1,k)=hdr3(1,k)-r360
if(hdr3(1,k) < zero)hdr3(1,k)=hdr3(1,k)+r360
dlon_earth_deg=hdr3(1,k)
dlat_earth_deg=hdr3(2,k)
dlon_earth=hdr3(1,k)*deg2rad
dlat_earth=hdr3(2,k)*deg2rad
if(regional)then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) ! convert to rotated coordinate
if(outside) cycle loop_readsb ! check to see if outside regional domain
else
dlat = dlat_earth
dlon = dlon_earth
call grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
timeobs=real(real(hdr3(3,k),r_single),r_double)
t4dv=timeobs + toff
zeps=1.0e-8_r_kind
if (t4dv<zero .and.t4dv> -zeps) t4dv=zero
if (t4dv>winlen.and.t4dv<winlen+zeps) t4dv=winlen
t4dv=t4dv + time_correction
time=timeobs + time_correction
if (l4dvar.or.l4densvar) then
if (t4dv<zero.OR.t4dv>winlen) cycle LOOP_K_LEVS
else
if (real(abs(time))>real(ctwind(nc)) .or. real(abs(time))>real(twindin)) cycle LOOP_K_LEVS
endif
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
timeobs = real(real(drfdat(3,k),r_single),r_double)
time_drift = timeobs + time_correction
if (abs(time_drift-time)>four) time_drift = time
! Check to see if the time is outside range
if (l4dvar.or.l4densvar) then
t4dv=toff+time_drift
if (t4dv<zero .or. t4dv>winlen) then
t4dv=toff+timex
if (t4dv<zero .or. t4dv>winlen) CYCLE LOOP_K_LEVS
end if
else
if(abs(time_drift) > ctwind(nc) .or. abs(time_drift) > twindin)then
time_drift=timex
if(abs(timex) > ctwind(nc) .or. abs(timex) > twindin) CYCLE LOOP_K_LEVS
end if
t4dv = toff + time_drift
endif
dlat_earth_deg = drfdat(2,k)
dlon_earth_deg = drfdat(1,k)
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 grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
end if
if((driftl .or. (aircraft_t_bc .and. acft_profl_file)) &
.and. ((tob.or. qob.or. uvob .and. levs > 1) .or. ithinp))then
! Interpolate guess pressure profile to observation location
klon1= int(dlon); klat1= int(dlat)
dx = dlon-klon1; dy = dlat-klat1
dx1 = one-dx; dy1 = one-dy
w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy
klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon)
if (klon1==0) klon1=nlon
klatp1=min(nlat,klat1+1); klonp1=klon1+1
if (klonp1==nlon+1) klonp1=1
do kk=1,nsig
presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) + &
w10*prsl_full(klatp1,klon1 ,kk) + &
w01*prsl_full(klat1 ,klonp1,kk) + &
w11*prsl_full(klatp1,klonp1,kk)
end do
! Compute depth of guess pressure layersat observation location
if (.not.twodvar_regional .and. levs > 1) then
do kk=1,nsig-1
dpres(kk)=presl(kk)-presl(kk+1)
end do
endif
end if
! Missing Values ==> Cycling! In this case for howv only. #ww3
if (howvob .and. owave(1,k) > r0_1_bmiss) cycle LOOP_K_LEVS
! Set usage variable
usage = zero
if(icuse(nc) <= 0)usage=100._r_kind
if(qm == 15 .or. qm == 12 .or. qm == 9)usage=100._r_kind
if(qm >=lim_qm )usage=101._r_kind
if(convobs .and. pqm(k) >=lim_qm )usage=102._r_kind
if((kx>=192.and.kx<=195) .and. psob )usage=r100
if (gustob .and. obsdat(8,k) > r0_1_bmiss) usage=103._r_kind
if (visob .and. obsdat(9,k) > r0_1_bmiss) usage=103._r_kind
if (tdob .and. obsdat(12,k) > r0_1_bmiss) usage=103._r_kind
if (pmob .and. obsdat(13,k) > r0_1_bmiss) usage=103._r_kind
if (mxtmob .and. maxtmint(1,k) > r0_1_bmiss) usage=103._r_kind
if (mitmob .and. maxtmint(2,k) > r0_1_bmiss) usage=103._r_kind
if (howvob .and. owave(1,k) > r0_1_bmiss) usage=103._r_kind
if (cldchob .and. cldceilh(1,k) > r0_1_bmiss) usage=103._r_kind
if (sfctype) then
if (i_gsdsfc_uselist==1 ) then
if (kx==188 .or. kx==195 .or. kx==288.or.kx==295) &
call apply_gsd_sfcuselist(kx,obstype,c_station_id,c_prvstg,c_sprvstg, &
usage)
else
call get_usagerj(kx,obstype,c_station_id,c_prvstg,c_sprvstg, &
dlon_earth,dlat_earth,idate,t4dv-toff, &
obsdat(5,k),obsdat(6,k),usage)
endif
!retrieve wind sensor height
if (twodvar_regional) then
if ( kx==288.or.kx==295 .or. (gustob .and. (kx==188.or.kx==195)) ) then
call find_wind_height(c_prvstg,c_sprvstg,windsensht)
endif
endif
endif
if (sfctype .and. i_gsdqc==2) then ! filter bad 2-m dew point and 0 mesonet wind obs
if (kx==288.or.kx==295) then ! for mesonet wind
if(abs(obsdat(5,k))<0.01_r_kind .and. abs(obsdat(6,k))<0.01_r_kind) usage=115._r_kind
endif
if (qob .and. (kx >=180 .and. kx<=189) .and. obsdat(2,k) < 1.0e10_r_kind) then ! for 2-m dew point
if(obsdat(12,k) < min(-40.0_r_kind,obsdat(3,k)-10.0_r_kind)) usage=116._r_kind ! < min(-40C or T-Td)
if((obsdat(3,k)-obsdat(12,k)) > 70.0_r_kind) usage=117._r_kind ! <70C
if(obsdat(12,k) > 32.2_r_kind) usage=118._r_kind ! > 90F
endif
endif
if ((kx>129.and.kx<140).or.(kx>229.and.kx<240) ) then
call get_aircraft_usagerj(kx,obstype,c_station_id,usage)
endif
if(plevs(k) < 0.0001_r_kind) then
write(*,*) 'warning: obs pressure is too small:',kx,k,plevs(k)
cycle
endif
if(ncnumgrp(nc) > 0 .and. .not.lhilbert )then ! default cross validation on
if(mod(ndata+1,ncnumgrp(nc))== ncgroup(nc)-1)usage=ncmiter(nc)
end if
! Flag regional MAP wind above 400mb for monitoring
if(regional .and. kx==227 .and. obsdat(1,k)<400._r_kind ) usage=r100
! don't use MESONET psfc obs if 8th character of station id is "x")
if( kx==188 .and. psob .and. sidchr(8)=='x' ) usage=r100
! Get information from surface file necessary for conventional data here
! Special block for data thinning - if requested
if (ithin > 0 .and. ithin <5 .and. usage <100.0_r_kind) then
! if (ithin > 0 .and. ithin <5) then
ntmp=ndata ! counting moved to map3gridS
! Set data quality index for thinning
if (thin4d) then
timedif = zero
else
timedif=abs(t4dv-toff)
endif
if(kx == 243 .or. kx == 253 .or. kx ==254) then
call ufbint(lunin,satqc,1,1,iret,satqcstr)
crit1 = timedif/r6+half + four*(one-satqc(1)/r100)*r3_33
else
crit1 = timedif/r6+half
endif
if (pflag==0) then
do kk=1,nsig
presl_thin(kk)=presl(kk)
end do
endif
if (ptime >zero ) then
itime=int((abs(timedif)+three)/ptime)+1
if(itime >ntime) itime=ntime
call map3grids_tm(zflag,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,&
ppb,itime,crit1,ndata,iout,icntpnt,iiout,luse,.false.,.false.)
if (.not. luse) then
if(k==levs) then
cycle loop_readsb
else
cycle LOOP_K_LEVS
endif
endif
if(iiout > 0) isort(iiout)=0
if (ndata > ntmp) then
nodata=nodata+1
if(uvob)nodata=nodata+1
endif
isort(icntpnt)=iout
else
call map3grids(zflag,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,&
ppb,crit1,ndata,iout,icntpnt,iiout,luse,.false.,.false.)
if (.not. luse) then
if(k==levs) then
cycle loop_readsb
else
cycle LOOP_K_LEVS
endif
endif
if(iiout > 0) isort(iiout)=0
if (ndata > ntmp) then
nodata=nodata+1
if(uvob)nodata=nodata+1
endif
isort(icntpnt)=iout
endif
else
ndata=ndata+1
nodata=nodata+1
if(uvob)nodata=nodata+1
iout=ndata
isort(icntpnt)=iout
endif
if(ndata > maxobs) then
write(6,*)'READ_PREPBUFR: ***WARNING*** ndata > maxobs for ',obstype
ndata = maxobs
end if
call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz)
if(lhilbert) &
call accum_hilbertcurve(usage,c_station_id,c_prvstg,c_sprvstg, &
dlat_earth,dlon_earth,dlat,dlon,t4dv,toff,nc,kx,iout)
! Extract pressure level and quality marks
dlnpob=log(plevs(k)) ! ln(pressure in cb)
! Set inflate_error logical based on qm flag
inflate_error=.false.
if (qm==3 .or. qm==7) inflate_error=.true.
! Temperature
if(tob) then
ppb=obsdat(1,k)
if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then
call errormod_aircraft(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3)
else
call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm)
end if
toe=obserr(3,k)*errout
qtflg=tvflg(k)
if (inflate_error) toe=toe*r1_2
if(ppb < r100)toe=toe*r1_2
if (aircraft_t_bc .and. kx==130 .and. ppb>=500.0_r_kind) toe=toe*r10
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)=obsdat(3,k)+t0c ! temperature ob.
cdata_all(6,iout)=rstation_id ! station id
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
cdata_all(9,iout)=qtflg ! qtflg (virtual temperature flag)
cdata_all(10,iout)=tqm(k) ! quality mark
cdata_all(11,iout)=obserr(3,k) ! original obs error
cdata_all(12,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=12 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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)=stnelev ! station elevation (m)
cdata_all(20,iout)=obsdat(4,k) ! observation height (m)
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
cdata_all(25,iout)=var_jb(3,k) ! non linear qc for T
if (aircraft_t_bc_pof .or. aircraft_t_bc .or.aircraft_t_bc_ext) then
cdata_all(26,iout)=aircraftwk(1,k) ! phase of flight
cdata_all(27,iout)=aircraftwk(2,k) ! vertical velocity
cdata_all(28,iout)=idx ! index of temperature bias
end if
if(perturb_obs)cdata_all(nreal,iout)=ran01dom()*perturb_fact ! t perturbation
if (twodvar_regional) &
call adjust_error(cdata_all(17,iout),cdata_all(18,iout),cdata_all(11,iout),cdata_all(1,iout))
! Winds
else if(uvob) then
if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then
call errormod_aircraft(pqm,wqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3)
else
call errormod(pqm,wqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm)
end if
woe=obserr(5,k)*errout
if (inflate_error) woe=woe*r1_2
if(obsdat(1,k) < r50)woe=woe*r1_2
selev=stnelev
oelev=obsdat(4,k)
if(kx >= 280 .and. kx < 300 )then
if (twodvar_regional.and.(kx==288.or.kx==295)) then
oelev=windsensht+selev !windsensht: read in from prepbufr
else
oelev=r10+selev
endif
if (kx == 280 )then
it29=nint(hdr(8))
if(it29 == 522 .or. it29 == 523 .or. it29 == 531)then
! oelev=r20+selev
oelev=r20
end if
end if
if (kx == 282) oelev=r20+selev
if (kx == 285 .or. kx == 289 .or. kx == 290) then
oelev=selev
selev=zero
endif
else
if((kx >= 221 .and. kx <= 229) &
.and. selev >= oelev) oelev=r10+selev
end if
! Rotate winds to rotated coordinate
uob=obsdat(5,k)
vob=obsdat(6,k)
!* thin new VAD wind and generate VAD superob
if(kx==224.and.newvad)then
klev=k+5 !*average over 6 points
! klev=k !* no average
if(klev>levs) cycle loop_readsb
diffuu=obsdat(5,k)-fcstdat(1,k)
diffvv=obsdat(6,k)-fcstdat(2,k)
if(sqrt(diffuu**2+diffvv**2)>10.0_r_kind) cycle loop_k_levs
if(abs(diffvv)>8.0_r_kind) cycle loop_k_levs
!if(abs(diffvv)>5.0.and.oelev<5000.0.and.fcstdat(3,k)>276.3) cycle loop_k_levs
if(oelev>7000.0_r_kind) cycle loop_k_levs
if(abs(diffvv)>5.0_r_kind.and.oelev<5000.0_r_kind) cycle loop_k_levs
! write(6,*)'sliu diffuu,vv::',diffuu, diffvv
uob=0.0
vob=0.0
oelev=0.0
tkk=0
do ikkk=k,klev
diffhgt=obsdat(4,ikkk)-obsdat(4,k)
if(diffhgt<301.0_r_kind)then
uob=uob+obsdat(5,ikkk)
vob=vob+obsdat(6,ikkk)
oelev=oelev+obsdat(4,ikkk)
tkk=tkk+1
end if
end do
uob=uob/tkk
vob=vob/tkk
oelev=oelev/tkk
diffuu=5.0_r_kind;diffvv=5.0_r_kind
diffhgt=0.0_r_kind
do ikkk=k,klev
diffuu=abs(obsdat(5,ikkk)-uob)
if(diffhgt<diffuu)diffhgt=diffuu
diffvv=abs(obsdat(6,ikkk)-vob)
if(diffhgt<diffvv)diffhgt=diffvv
end do
if(diffhgt>5.0_r_kind)cycle LOOP_K_LEVS !* if u-u_avg>5.0, reject
if(tkk<3) cycle LOOP_K_LEVS !* obs numb<3, reject
!* unreasonable observation, will fix this in QC package
if(sqrt(uob**2+vob**2)>60.0_r_kind)cycle LOOP_readsb
end if
if(regional .and. .not. fv3_regional)then
u0=uob
v0=vob
call rotate_wind_ll2xy(u0,v0,uob,vob,dlon_earth,dlon,dlat)
if(diagnostic_reg) then
call rotate_wind_xy2ll(uob,vob,u00,v00,dlon_earth,dlon,dlat)
nvtest=nvtest+1
disterr=sqrt((u0-u00)**2+(v0-v00)**2)
vdisterrmax=max(vdisterrmax,disterr)
end if
endif
cdata_all(1,iout)=woe ! wind error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=dlnpob ! ln(pressure in cb)
cdata_all(5,iout)=oelev ! height of observation
cdata_all(6,iout)=uob ! u obs
cdata_all(7,iout)=vob ! v obs
cdata_all(8,iout)=rstation_id ! station id
cdata_all(9,iout)=t4dv ! time
cdata_all(10,iout)=nc ! type
cdata_all(11,iout)=selev ! station elevation
cdata_all(12,iout)=wqm(k) ! quality mark
cdata_all(13,iout)=obserr(5,k) ! original obs error
cdata_all(14,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=14 ! save INDEX of where usage
! is stored for hilbertcurve cross validation (if requested)
cdata_all(15,iout)=idomsfc ! dominate surface type
cdata_all(16,iout)=tsavg ! skin temperature
cdata_all(17,iout)=ff10 ! 10 meter wind factor
cdata_all(18,iout)=sfcr ! surface roughness
cdata_all(19,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(20,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
cdata_all(25,iout)=var_jb(5,k) ! non linear qc parameter
cdata_all(26,iout)=one ! hilbert curve weight, modified later
if(perturb_obs)then
cdata_all(27,iout)=ran01dom()*perturb_fact ! u perturbation
cdata_all(28,iout)=ran01dom()*perturb_fact ! v perturbation
endif
else if(spdob) then
woe=obserr(5,k)
if (inflate_error) woe=woe*r1_2
elev=r20
if (((kx==295).or.(kx==288)).and.twodvar_regional) then !account for mesonet wind ht
oelev=windsensht+selev
else
oelev=obsdat(4,k)
endif
if(kx == 260 .or. kx == 261) elev = oelev ! Nacelle and tower wind speed
cdata_all(1,iout)=woe ! wind error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=dlnpob ! ln(pressure in cb)
cdata_all(5,iout)=obsdat(5,k) ! u obs
cdata_all(6,iout)=obsdat(6,k) ! v obs
cdata_all(7,iout)=rstation_id ! station id
cdata_all(8,iout)=t4dv ! time
cdata_all(9,iout)=nc ! type
cdata_all(10,iout)=elev ! elevation of observation
cdata_all(11,iout)=wqm(k) ! quality mark
cdata_all(12,iout)=obserr(5,k) ! original obs error
cdata_all(13,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=13 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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 (degrees)
cdata_all(19,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(20,iout)=stnelev ! station elevation (m)
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
! Surface pressure
else if(psob) then
poe=obserr(1,k)*one_tenth ! convert from mb to cb
if (inflate_error) poe=poe*r1_2
cdata_all(1,iout)=poe ! surface pressure error (cb)
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=exp(dlnpob) ! pressure (in cb)
cdata_all(5,iout)=obsdat(4,k) ! surface height
cdata_all(6,iout)=obsdat(3,k)+t0c ! surface temperature
cdata_all(7,iout)=rstation_id ! station id
cdata_all(8,iout)=t4dv ! time
cdata_all(9,iout)=nc ! type
cdata_all(10,iout)=pqm(k) ! quality mark
cdata_all(11,iout)=obserr(1,k)*one_tenth ! original obs error (cb)
cdata_all(12,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=12 ! save INDEX of where usage is stored
! for hilbertcurve cross validation (if requested)
cdata_all(13,iout)=idomsfc ! dominate surface type
cdata_all(14,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(15,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(16,iout)=stnelev ! station elevation (m)
cdata_all(17,iout)=zz ! terrain height at ob location
cdata_all(18,iout)=r_prvstg(1,1) ! provider name
cdata_all(19,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(20,iout)=var_jb(1,k) ! non linear qc b parameter
if(perturb_obs)cdata_all(21,iout)=ran01dom()*perturb_fact ! ps perturbation
if (twodvar_regional) &
call adjust_error(cdata_all(14,iout),cdata_all(15,iout),cdata_all(11,iout),cdata_all(1,iout))
! Specific humidity
else if(qob) then
qmaxerr=emerr
if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then
call errormod_aircraft(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3)
else
call errormod(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm)
end if
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,iout)=qoe ! q 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)=qobcon ! q ob
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
cdata_all(10,iout)=tdry ! dry temperature (obs is tv)
cdata_all(11,iout)=qqm(k) ! quality mark
cdata_all(12,iout)=obserr(2,k)*one_tenth ! original obs error
cdata_all(13,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=13 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(14,iout)=idomsfc ! dominate surface type
cdata_all(15,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(16,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(17,iout)=stnelev ! station elevation (m)
cdata_all(18,iout)=obsdat(4,k) ! observation height (m)
cdata_all(19,iout)=zz ! terrain height at ob location
cdata_all(20,iout)=r_prvstg(1,1) ! provider name
cdata_all(21,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(22,iout)=obsdat(10,k) ! cat
cdata_all(23,iout)=var_jb(2,k) ! non linear qc b parameter
if(perturb_obs)cdata_all(24,iout)=ran01dom()*perturb_fact ! q perturbation
if (twodvar_regional) &
call adjust_error(cdata_all(15,iout),cdata_all(16,iout),cdata_all(12,iout),cdata_all(1,iout))
! Total precipitable water (ssm/i)
else if(pwob) then
pwoe=obserr(7,k)
pwmerr=pwoe*three
cdata_all(1,iout)=pwoe ! pw error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=obsdat(7,k) ! pw obs
cdata_all(5,iout)=rstation_id ! station id
cdata_all(6,iout)=t4dv ! time
cdata_all(7,iout)=nc ! type
cdata_all(8,iout)=pwmerr ! pw max error
cdata_all(9,iout)=pwq ! quality mark
cdata_all(10,iout)=obserr(7,k) ! original obs error
cdata_all(11,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=11 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(12,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(13,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(14,iout)=stnelev ! station elevation (m)
cdata_all(15,iout)=obsdat(1,k) ! observation pressure (hPa)
cdata_all(16,iout)=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
cdata_all(1,iout)=sstoe ! sst error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=sstdat(3,k) ! sst obs
cdata_all(5,iout)=rstation_id ! station id
cdata_all(6,iout)=t4dv ! time
cdata_all(7,iout)=nc ! type
cdata_all(8,iout)=sstoe*three ! pw max error
cdata_all(9,iout)=sstdat(2,k) ! depth of measurement
cdata_all(10,iout)=sstdat(1,k) ! measurement type
cdata_all(11,iout)=sstq ! quality mark
cdata_all(12,iout)=sstdat(5,k) ! original obs error
cdata_all(13,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=13 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(14,iout)=idomsfc ! dominate surface type
cdata_all(15,iout)=tsavg ! skin temperature
cdata_all(16,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(17,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(18,iout)=stnelev ! station elevation (m)
if( nst_gsi > 0) then
zob = sstdat(2,k)
if (zob > 10.0) then
tref = tsavg
dtw = zero
dtc = zero
tz_tr = one
else
call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr)
end if
cdata_all(19,iout) = tref ! foundation temperature
cdata_all(20,iout) = dtw ! dt_warm at zob
cdata_all(21,iout) = dtc ! dt_cool at zob
cdata_all(22,iout) = tz_tr ! d(Tz)/d(Tr)
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
! Wind gusts
else if(gustob) then
! need to find out gustoe
! gustoe=1.8
gustoe=1.0
selev=stnelev
oelev=obsdat(4,k)
if(selev == oelev)oelev=r10+selev
if((kx >= 280 .and. kx < 300).or.(kx >= 180 .and. kx < 200))then
oelev=r10+selev
if ((kx==280).or.(kx==180)) oelev=r20+selev
if ((kx==299).or.(kx==199)) oelev=r20+selev
if ((kx==282).or.(kx==182)) oelev=r20+selev
if (((kx==295).or.(kx==288).or.(kx==195).or.(kx==188)).and.twodvar_regional) then
!account for mesonet wind sensor height
oelev=windsensht+selev
end if
if (kx==198) oelev=r20+selev
if ((kx==285).or.(kx==185)) then
oelev=selev
selev=zero
end if
if ((kx==188).or.(kx==288) .or.(kx==195) .or.(kx==295)) then
! gustoe=2.5
gustoe=1.0
windcorr=abs(obsdat(5,k))<1.0 .and. abs(obsdat(6,k))<1.0 .and. obsdat(8,k)>10.0
if (windcorr) gustoe=gustoe*1.5_r_kind
if (abs(obsdat(8,k)-sqrt(obsdat(5,k)**2+obsdat(6,k)**2))<1.5) then
gustoe=gustoe*1.5_r_kind
end if
end if
end if
if (inflate_error) gustoe=gustoe*1.5_r_kind
cdata_all(1,iout)=gustoe ! wind gusts error (cb)
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=dlnpob ! ln(pressure in cb)
cdata_all(5,iout)=oelev ! observation height
cdata_all(6,iout)=obsdat(8,k) ! wind gusts obs
cdata_all(7,iout)=rstation_id ! station id
cdata_all(8,iout)=t4dv ! time
cdata_all(9,iout)=nc ! type
cdata_all(10,iout)=gustoe*three ! max error
cdata_all(11,iout)=gustqm ! quality mark
cdata_all(12,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=12 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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)=selev ! station elevation (m)
cdata_all(20,iout)=r_prvstg(1,1) ! provider name
cdata_all(21,iout)=r_sprvstg(1,1) ! subprovider name
! Visibility
else if(visob) then
!......................................................................
!NLTRCV: must setup as true
! visoe is in NLTR space, and is read in from the namelist. Is this OK?
!......................................................................
visoe=estvisoe
if ((kx==283).or.(kx==183)) visoe=visoe*r1_02
if (inflate_error) visoe=visoe*r1_02
cdata_all(1,iout)=visoe ! visibility error (cb)
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
!......................................................................
! simple QC check: if an observation vis is negative, assign it as bmiss
! if obs. = zero, reassign it as one_r_kind
! about bmiss:
! #ifdef ibm_sp ! real(r_kind), parameter:: bmiss = 1.0e11_r_kind !#else
! real(r_kind), parameter:: bmiss = 1.0e9_r_kind !#endif
! in setupvis: missing data is checked and assigned not use in muse
! visthres is much smaller than bmiss
! i.e: this holds: (obsdat(9,k)> zero .and. obsdat(9,k)<=vis_thres)
!......................................................................
if(obsdat(9,k) < zero) then
cdata_all(4,iout)=bmiss
elseif(obsdat(9,k)> r0_1_bmiss)then
cdata_all(4,iout)=obsdat(9,k)
elseif(obsdat(9,k)> vis_thres .and. obsdat(9,k)<= r0_1_bmiss )then
obsdat(9,k)=vis_thres
else
obsdat(9,k)=max(obsdat(9,k),one)
endif
if(obsdat(9,k)> zero .and. obsdat(9,k)<=vis_thres)then
tempvis=obsdat(9,k)
call nltransf_forward(tempvis,visout,pvis,scale_cv)
cdata_all(4,iout) = visout
endif
cdata_all(5,iout)=rstation_id ! station id
cdata_all(6,iout)=t4dv ! time
cdata_all(7,iout)=nc ! type
cdata_all(8,iout)=visoe*three ! max error
cdata_all(9,iout)=visqm ! quality mark
cdata_all(10,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=10 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(11,iout)=idomsfc ! dominate surface type
cdata_all(12,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(13,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(14,iout)=stnelev ! station elevation (m)
cdata_all(15,iout)=obsdat(4,k) ! observation height (m)
cdata_all(16,iout)=zz ! terrain height at ob location
cdata_all(17,iout)=r_prvstg(1,1) ! provider name
cdata_all(18,iout)=r_sprvstg(1,1) ! subprovider name
! 2m-Dewpoint
else if(tdob) then
tdoe=obserr(3,k)*r1_2
qobcon=obsdat(2,k)*convert
tdry=r999
if (tqm(k)<lim_tqm) tdry=(obsdat(3,k)+t0c)/(one+fv*qobcon)
cdata_all(1,iout)=tdoe ! td error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=dlnpob ! ln(pressure in cb)
cdata_all(5,iout)=obsdat(12,k)+t0c ! td ob
cdata_all(6,iout)=rstation_id ! station id
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
cdata_all(9,iout)=tdoe*three ! td max error
cdata_all(10,iout)=tdry ! dry temperature (obs is tv)
cdata_all(11,iout)=tdqm ! quality mark
cdata_all(12,iout)=tdoe ! original obs error
cdata_all(13,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=13 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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 (degrees)
cdata_all(19,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(20,iout)=stnelev ! station elevation (m)
cdata_all(21,iout)=obsdat(4,k) ! observation height (m)
cdata_all(22,iout)=zz ! terrain height at ob location
cdata_all(23,iout)=r_prvstg(1,1) ! provider name
cdata_all(24,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(25,iout)=obsdat(10,k) ! cat
! Maximum temperature
else if(mxtmob) then
mxtmoe=obserr(3,k)
qtflg=one
cdata_all(1,iout)=mxtmoe ! maximum 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)=maxtmint(1,k) ! maximum temperature ob.
cdata_all(6,iout)=rstation_id ! station id
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
cdata_all(9,iout)=qtflg ! qtflg (virtual temperature flag)
cdata_all(10,iout)=mxtmqm ! quality mark
cdata_all(11,iout)=obserr(3,k) ! original obs error
cdata_all(12,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=12 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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)=stnelev ! station elevation (m)
cdata_all(20,iout)=obsdat(4,k) ! observation height (m)
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
! Minimum temperature
else if(mitmob) then
mitmoe=obserr(3,k)
qtflg=one
cdata_all(1,iout)=mitmoe ! minimum 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)=maxtmint(2,k) ! minimum temperature ob.
cdata_all(6,iout)=rstation_id ! station id
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
cdata_all(9,iout)=qtflg ! qtflg (virtual temperature flag)
cdata_all(10,iout)=mitmqm ! quality mark
cdata_all(11,iout)=obserr(3,k) ! original obs error
cdata_all(12,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=12 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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)=stnelev ! station elevation (m)
cdata_all(20,iout)=obsdat(4,k) ! observation height (m)
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
! Pressure at mean sea level
else if(pmob) then
pmoe=obserr(1,k)*one_tenth ! convert from mb to cb
if (inflate_error) pmoe=pmoe*r1_2
cdata_all(1,iout)=pmoe ! pressure at mean sea level error (cb)
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=exp(dlnpob) ! pressure (in cb)
cdata_all(5,iout)=one_tenth*obsdat(13,k) ! pressure at mean sea level (in cb)
cdata_all(6,iout)=obsdat(4,k) ! surface height
cdata_all(7,iout)=obsdat(3,k)+t0c ! surface temperature
cdata_all(8,iout)=rstation_id ! station id
cdata_all(9,iout)=t4dv ! time
cdata_all(10,iout)=nc ! type
cdata_all(11,iout)=pmq(k) ! quality mark
cdata_all(12,iout)=obserr(1,k)*one_tenth ! original obs error (cb)
cdata_all(13,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=13 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
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 (degrees)
cdata_all(19,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(20,iout)=stnelev ! station elevation (m)
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
! Significant wave height
else if(howvob) then
zz=0_r_kind ! #ww3
howvoe=0.3_r_kind ! use temporarily
cdata_all(1,iout)=howvoe ! significant wave height error (m)
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=exp(dlnpob) ! pressure (in cb)
cdata_all(5,iout)=owave(1,k) ! significant wave height (in m)
cdata_all(6,iout)=rstation_id ! station id
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
cdata_all(9,iout)=howvqm ! quality mark
cdata_all(10,iout)=howvoe ! original obs error (m)
cdata_all(11,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=11 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(12,iout)=idomsfc ! dominate surface type
cdata_all(13,iout)=tsavg ! skin temperature
cdata_all(14,iout)=ff10 ! 10 meter wind factor
cdata_all(15,iout)=sfcr ! surface roughness
cdata_all(16,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(17,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(18,iout)=stnelev ! station elevation (m)
cdata_all(19,iout)=obsdat(4,k) ! observation height (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)=obsdat(10,k) ! cat
! METAR cloud observation
else if(metarcldobs) then
cdata_all(1,iout)=rstation_id ! station ID
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=stnelev ! station elevation
if(metarvis(1,1) < r0_1_bmiss) then
cdata_all(5,iout)=metarvis(1,1) ! visibility (m)
else
cdata_all(5,iout) = -99999.0_r_kind
endif
do kk=1, 6
if(metarcld(1,kk) < r0_1_bmiss) then
cdata_all(5+kk,iout) =metarcld(1,kk) ! cloud amount
else
cdata_all(5+kk,iout) = -99999.0_r_kind
endif
if(metarcld(2,kk) < r0_1_bmiss) then
cdata_all(11+kk,iout)=metarcld(2,kk) ! cloud bottom height (m)
else
cdata_all(11+kk,iout)= -99999.0_r_kind
endif
enddo
do kk=1, 3
if(metarwth(1,kk) < r0_1_bmiss) then
cdata_all(17+kk,iout)=metarwth(1,kk) ! weather
else
cdata_all(17+kk,iout)= -99999.0_r_kind
endif
enddo
cdata_all(21,iout)=timeobs ! time observation
cdata_all(22,iout)=usage
if (lhilbert) thisobtype_usage=22 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(23,iout)=0.0_r_kind ! reserved for distance between obs and grid
! Calculate dewpoint depression from surface obs, to be used later
! with haze and ceiling logic to exclude dust-caused ceiling obs
! from cloud analysis
if(metarvis(2,1) < 1.e10_r_kind) then
cdata_all(24,iout)=obsdat(3,1)-metarvis(2,1) ! temperature - dew point
else
cdata_all(24,iout)=-99999.0_r_kind ! temperature - dew point
endif
cdata_all(25,iout)=nc ! type
cdata_all(26,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(27,iout)=dlat_earth_deg ! earth relative latitude (degrees)
! NESDIS cloud products
else if(goesctpobs) then
cdata_all(1,iout)=rstation_id ! station ID
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=goescld(1,k)/100.0_r_kind ! cloud top pressure (pa)
cdata_all(5,iout)=goescld(2,k) ! cloud cover
cdata_all(6,iout)=goescld(3,k) ! Cloud top temperature (K)
cdata_all(7,iout)=timeobs ! time
cdata_all(8,iout)=usage
if (lhilbert) thisobtype_usage=8 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
! Total cloud amount
else if(tcamtob) then
if (k==1) then
! adjust quality mark/usage parameter
if (trim(subset) == 'GOESND') then
call adjust_goescldobs(goescld(2,1),timeobs,dlat_earth,dlon_earth, &
low_cldamt,low_cldamt_qc,mid_cldamt,mid_cldamt_qc, &
hig_cldamt,hig_cldamt_qc,tcamt,tcamt_qc)
else
call adjust_convcldobs(cld2seq,cld2seqlevs,cldseq,cldseqlevs,metarwth,metarwthlevs, &
low_cldamt,low_cldamt_qc,mid_cldamt,mid_cldamt_qc, &
hig_cldamt,hig_cldamt_qc,tcamt,lcbas,tcamt_qc,lcbas_qc,ceiling,stnelev)
end if
if(tcamt_qc==15 .or. tcamt_qc==12 .or. tcamt_qc==9) usage=100._r_kind
tcamt_oe=20.0_r_kind
if(tcamt_qc==1) tcamt_oe=tcamt_oe*1.25_r_kind
if(tcamt_qc==2) tcamt_oe=tcamt_oe*1.50_r_kind
if(tcamt_qc==3) tcamt_oe=tcamt_oe*1.75_r_kind
cdata_all( 1,iout)=tcamt_oe ! obs error
cdata_all( 2,iout)=dlon ! grid relative longitude
cdata_all( 3,iout)=dlat ! grid relative latitude
cdata_all( 4,iout)=tcamt ! total cloud amount (%)
cdata_all( 5,iout)=rstation_id ! station ID
cdata_all( 6,iout)=t4dv ! time
cdata_all( 7,iout)=nc ! type
cdata_all( 8,iout)=tcamt_qc ! quality mark
cdata_all( 9,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=9 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(10,iout)=idomsfc ! dominate surface type
cdata_all(11,iout)=tsavg ! skin temperature
cdata_all(12,iout)=ff10 ! 10 meter wind factor
cdata_all(13,iout)=sfcr ! surface roughness
cdata_all(14,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(15,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(16,iout)=stnelev ! station elevation (m)
cdata_all(17,iout)=obsdat(4,k) ! observation height (m)
cdata_all(18,iout)=zz ! terrain height at ob location
cdata_all(19,iout)=r_prvstg(1,1) ! provider name
cdata_all(20,iout)=r_sprvstg(1,1) ! subprovider name
end if
! Base height of the lowest cloud seen
else if(lcbasob) then
if (k==1) then
! adjust quality mark/usage parameter
call adjust_convcldobs(cld2seq,cld2seqlevs,cldseq,cldseqlevs,metarwth,metarwthlevs, &
low_cldamt,low_cldamt_qc,mid_cldamt,mid_cldamt_qc, &
hig_cldamt,hig_cldamt_qc,tcamt,lcbas,tcamt_qc,lcbas_qc,ceiling,stnelev)
if(lcbas_qc==15 .or. lcbas_qc==12 .or. lcbas_qc==9) usage=100._r_kind
lcbas_oe=4500.0_r_kind
if(lcbas_qc==3) lcbas_oe=lcbas_oe*1.25_r_kind
if(lcbas_qc==4) lcbas_oe=lcbas_oe*1.5_r_kind
cdata_all( 1,iout)=lcbas_oe ! obs error
cdata_all( 2,iout)=dlon ! grid relative longitude
cdata_all( 3,iout)=dlat ! grid relative latitude
cdata_all( 4,iout)=lcbas ! base height of lowest cloud (m)
cdata_all( 5,iout)=rstation_id ! station ID
cdata_all( 6,iout)=t4dv ! time
cdata_all( 7,iout)=nc ! type
cdata_all( 8,iout)=lcbas_qc ! quality mark
cdata_all( 9,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=9 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(10,iout)=idomsfc ! dominate surface type
cdata_all(11,iout)=tsavg ! skin temperature
cdata_all(12,iout)=ff10 ! 10 meter wind factor
cdata_all(13,iout)=sfcr ! surface roughness
cdata_all(14,iout)=dlon_earth_deg ! earth relative longitude (degrees)
cdata_all(15,iout)=dlat_earth_deg ! earth relative latitude (degrees)
cdata_all(16,iout)=stnelev ! station elevation (m)
cdata_all(17,iout)=obsdat(4,k) ! observation height (m)
cdata_all(18,iout)=zz ! terrain height at ob location
cdata_all(19,iout)=ceiling ! cloud ceiling obs
if (trim(subset) == 'GOESND') then
! cdata_all(20,iout)=goescld(1,k)/100.0_r_kind ! cloud top pressure (pa)
cdata_all(20,iout)=goescld(1,k) ! cloud top pressure
cdata_all(21,iout)=goescld(3,k) ! Cloud top temperature (K)
else
cdata_all(20,iout)=bmiss
cdata_all(21,iout)=bmiss
end if
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
end if
! Cloud ceiling height
else if(cldchob) then
!......................................................................
!NLTRCV: must setup as true
! cldchoe is in NLTR space, and is read in via the namelist. Is this OK?
!......................................................................
cldchoe=estcldchoe
if (inflate_error) cldchoe=cldchoe*r1_02
cdata_all(1,iout)=cldchoe ! cloud ceiling height error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
!......................................................................
! NLTRCV:
! simple QC check and designate bad observation.
! if obs. cldch < zero, assign it bmiss
! if obs/first cldch guess is zero, assigne it as one
! about bmiss:
! #ifdef ibm_sp ! real(r_kind), parameter:: bmiss = 1.0e11_r_kind !#else
! real(r_kind), parameter:: bmiss = 1.0e9_r_kind !#endif
! in setupcldch: missing data is checked and assigned as not-use
! cldchthres is much smaller than bmiss
! i.e: this holds: (obsdat(x,k)> zero .and. obsdat(x,k)<=cldch_thres)
!......................................................................
if(cldceilh(1,k) < zero) then
cdata_all(4,iout)=bmiss
elseif (cldceilh(1,k)> r0_1_bmiss) then
cdata_all(4,iout)=cldceilh(1,k)
elseif (cldceilh(1,k)>=cldch_thres .and. cldceilh(1,k)<= r0_1_bmiss) then
cldceilh(1,k)=cldch_thres
else
cldceilh(1,k)=max(cldceilh(1,k),one) !consider cldceilh(1,k)=zero a valid data
endif
if (cldceilh(1,k)> zero .and. cldceilh(1,k)<=cldch_thres)then
tempcldch=cldceilh(1,k)
call nltransf_forward(tempcldch,cldchout,pcldch,scale_cv)
cdata_all(4,iout) = cldchout
endif ! ceiling height obs
cdata_all(5,iout)=rstation_id ! station id
cdata_all(6,iout)=t4dv ! time
cdata_all(7,iout)=nc ! type
cdata_all(8,iout)=cldchoe*three ! max error
cdata_all(9,iout)=cldchqm ! quality mark
cdata_all(10,iout)=usage ! usage parameter
if (lhilbert) thisobtype_usage=10 ! save INDEX of where usage is stored for hilbertcurve cross validation (if requested)
cdata_all(11,iout)=idomsfc ! dominate surface type
cdata_all(12,iout)=dlon_earth*rad2deg ! earth relative longitude (degrees)
cdata_all(13,iout)=dlat_earth*rad2deg ! earth relative latitude (degrees)
cdata_all(14,iout)=stnelev ! station elevation (m)
cdata_all(15,iout)=obsdat(4,k) ! observation height (m)
cdata_all(16,iout)=zz ! terrain height at ob location
cdata_all(17,iout)=r_prvstg(1,1) ! provider name
cdata_all(18,iout)=r_sprvstg(1,1) ! subprovider name
end if
!
! End k loop over levs
end do LOOP_K_LEVS
end do loop_readsb
!
! End of bufr read loop
enddo loop_msg
! Deallocate arrays used for thinning data
if (.not.use_all) then
deallocate(presl_thin)
call del3grids
endif
if (.not.use_all_tm) then
deallocate(presl_thin)
call del3grids_tm
endif
! Normal exit
enddo loop_convinfo! loops over convinfo entry matches
deallocate(lmsg,tab,nrep)
! Close unit to bufr file
call closbf(lunin)
close(lunin)
! Apply hilbert curve for cross validation if requested
if(lhilbert) &
call apply_hilbertcurve(maxobs,obstype,cdata_all(thisobtype_usage,1:maxobs))
! Write header record and data to output file for further processing
allocate(iloc(ndata))
icount=0
do i=1,maxobs
if(isort(i) > 0)then
icount=icount+1
iloc(icount)=isort(i)
end if
end do
if(ndata /= icount)then
write(6,*) ' PREPBUFR: mix up in read_prepbufr ,ndata,icount ',ndata,icount
call stop2(50)
end if
allocate(cdata_out(nreal,ndata))
do i=1,ndata
itx=iloc(i)
do k=1,nreal
cdata_out(k,i)=cdata_all(k,itx)
end do
end do
deallocate(iloc,isort,cdata_all)
! the following is gettin the types which will be applied hilbert curve to
! estimate the density
if(obstype == 'uv') then
vmin=-10.00_r_kind
vmax=18000.00_r_kind
nor=0
ithin=0
allocate(wght_hilb(ndata))
wght_hilb=one
pmesh=zero
rmesh=zero
dentrip=0
pmesh_tmp=zero
rmesh_tmp=zero
dentrip_tmp=0
ntime_max=0
ntime_tmp=0
ntime=0
nnrand=nrand
do ncc=1,nconvtype
if( trim(ioctype(ncc)) == 'uv') then
itype=ictype(ncc)
if( itype ==230 .or. itype ==231 .or. itype ==233 ) then
if(itype ==230 .and. ithin_conv(ncc) >5) then
if(ptime_conv(ncc) >zero) then
ntime_max= int(6.0_r_kind/ptime_conv(ncc))
else
ntime_max=1
endif
pmesh=pmesh_conv(ncc)
rmesh=rmesh_conv(ncc)
dentrip=ithin_conv(ncc)/10
else if(itype == 231 .and. ithin_conv(ncc) >5) then
if(ptime_conv(ncc) >zero) then
ntime_tmp=int(6.0_r_kind/ptime_conv(ncc))
else
ntime_tmp=1
endif
pmesh_tmp=pmesh_conv(ncc)
rmesh_tmp=rmesh_conv(ncc)
dentrip_tmp=ithin_conv(ncc)/10
if(pmesh >zero .and. rmesh >zero) then
if(ntime_tmp /= ntime_max .or. pmesh_tmp /= pmesh .or. &
rmesh_tmp /= rmesh .or. dentrip_tmp /= dentrip) then
write(6,*) 'READ_PREPBUFR:WARING convinfo file settings are not right,use first one itype=',itype
write(6,*) 'READ_PREPBUFR: ntime_max,pmesh,rmesh,ntime_tmp,pmesh_tmp,rmesh_tmp,rmesh_tmp=',&
ntime_max,pmesh,rmesh,ntime_tmp,pmesh_tmp,rmesh_tmp,dentrip_tmp
endif
else if(pmesh_tmp >zero .and. rmesh_tmp >zero) then
pmesh=pmesh_tmp
rmesh=rmesh_tmp
ntime_max=ntime_tmp
endif
else if(itype == 233 .and. ithin_conv(ncc) >5) then
if(ptime_conv(ncc) >zero) then
ntime_tmp=int(6.0_r_kind/ptime_conv(ncc))
else
ntime_tmp=1
endif
pmesh_tmp=pmesh_conv(ncc)
rmesh_tmp=rmesh_conv(ncc)
dentrip_tmp=ithin_conv(ncc)/10
if(pmesh >zero .and. rmesh >zero) then
if(ntime_tmp /= ntime_max .or. pmesh_tmp /= pmesh .or. &
rmesh_tmp /= rmesh .or. dentrip_tmp /= dentrip) then
write(6,*) 'READ_PREPBUFR:WARING convinfo file settings are not right,use first one itype=',itype
write(6,*) 'READ_PREPBUFR: ntime_max,pmesh,rmesh,ntime_tmp,pmesh_tmp,rmesh_tmp,rmesh_tmp=',&
ntime_max,pmesh,rmesh,ntime_tmp,pmesh_tmp,rmesh_tmp,dentrip_tmp
endif
else if(pmesh_tmp >zero .and. rmesh_tmp >zero) then
pmesh=pmesh_tmp
rmesh=rmesh_tmp
ntime_max=ntime_tmp
endif
endif
endif
endif
enddo
write(6,*),'READ_PREPBUFR:dentrip,pmesh,rmesh,ndata=',dentrip,pmesh,rmesh,ntime_max,ndata
if(dentrip >= one .and. pmesh >zero .and. rmesh >zero) then
allocate(data_hilb(3,ndata,6),index_arr(ndata,ntime_max))
ndata_hil=0
ntype_arr=0
ntime=1
index_arr=0
do k=1,ndata
ikx=nint(cdata_out(10,k))
if (ikx>0) then
itype=ictype(ikx)
else
itype=0
endif
if( itype ==230 .or. itype ==231 .or. itype ==233) then
prest=r10*exp(cdata_out(4,k))
if (prest <100.0_r_kind) cycle
if(ithin_conv(ikx) >=5) then
if(ptime_conv(ikx) >zero) then
ntime=int(((cdata_out(9,k)-time_offset)+three)/ptime_conv(ikx))+1
endif
if(ntime >ntime_max) ntime=ntime_max
if(ntime <0) ntime=1
ntype_arr(ntime)=ntype_arr(ntime)+1
ndata_hil=ntype_arr(ntime)
data_hilb(1,ndata_hil,ntime)=cdata_out(20,k)
data_hilb(2,ndata_hil,ntime)=cdata_out(19,k)
prest=prest*100.0_r_kind
if(prest >stndrd_atmos_ps) then
prest=zero
else
prest=rd*265.00_r_kind*log(stndrd_atmos_ps/prest)/grav
endif
data_hilb(3,ndata_hil,ntime)=prest
index_arr(ndata_hil,ntime)=k
if(data_hilb(1,ndata_hil,ntime) >90.0_r_kind .or. &
data_hilb(1,ndata_hil,ntime) <-90.0_r_kind .or. &
data_hilb(2,ndata_hil,ntime) <zero .or. &
data_hilb(2,ndata_hil,ntime) >360.0_r_kind .or. &
data_hilb(3,ndata_hil,ntime) <vmin .or. &
data_hilb(3,ndata_hil,ntime) >vmax ) then
write(6,*),'READ_PREPBUFR :something is wrong,lat,lon,prest=',&
data_hilb(1,ndata_hil,ntime),&
data_hilb(2,ndata_hil,ntime),&
cdata_out(4,k),data_hilb(3,ndata_hil,ntime)
endif
endif
endif
enddo
rmesh=rmesh*1000.0_r_kind
do kk=1,ntime_max
ndata_hil=ntype_arr(kk)
if(ndata_hil >=2) then
allocate(rlat_hil(ndata_hil),rlon_hil(ndata_hil),height(ndata_hil),wtob(ndata_hil))
rlat_hil(1:ndata_hil)=data_hilb(1,1:ndata_hil,kk)
rlon_hil(1:ndata_hil)=data_hilb(2,1:ndata_hil,kk)
height(1:ndata_hil)=data_hilb(3,1:ndata_hil,kk)
call denest(ndata_hil,nnrand,nor,rearth,dentrip,rmesh,pmesh,&
vmin,vmax,rlat_hil,rlon_hil,height,wtob)
do i=1,ndata_hil
indexx=index_arr(i,kk)
wght_hilb(indexx)=wtob(i)
enddo
ndata_hil=0
deallocate(rlat_hil,rlon_hil,height,wtob)
endif
enddo
deallocate(data_hilb,index_arr)
endif
do i=1,ndata
cdata_out(26,i)=wght_hilb(i)
enddo
deallocate(wght_hilb)
endif
! end of hilbert curve
! define a closest METAR cloud observation for each grid point
if(metarcldobs .and. ndata > 0) then
if(i_ens_mean /= 1) then
maxobs=2000000
allocate(cdata_all(nreal,maxobs))
call reorg_metar_cloud(cdata_out,nreal,ndata,cdata_all,maxobs,iout)
ndata=iout
deallocate(cdata_out)
allocate(cdata_out(nreal,ndata))
do i=1,nreal
do j=1,ndata
cdata_out(i,j)=cdata_all(i,j)
end do
end do
deallocate(cdata_all)
endif
endif
call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs)
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata
write(lunout) cdata_out
deallocate(cdata_out)
call destroy_rjlists
call destroy_aircraft_rjlists
if(i_gsdsfc_uselist==1) call destroy_gsd_sfcuselist
if (lhilbert) call destroy_hilbertcurve
if (twodvar_regional) then
call destroy_ndfdgrid
call destroy_windht_lists
endif
if(diagnostic_reg .and. ntest>0) write(6,*)'READ_PREPBUFR: ',&
'ntest,disterrmax=',ntest,disterrmax
if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_PREPBUFR: ',&
'nvtest,vdisterrmax=',ntest,vdisterrmax
if(print_verbose)write(6,*)'READ_PREPBUFR: closbf(',lunin,')'
! End of routine
return
end subroutine read_prepbufr
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!
! !ROUTINE: sonde_ext -level enhancemnt for raob
!
! !INTERFACE:
!
subroutine sonde_ext(obsdat,tpc,qcmark,obserr,drfdat,levsio,kx,vtcd)
!$$$ subprogram documentation block
! . . . .
! subprogram: sonde_ext level enhancemnt for raob
! prgmmr: wu org: np22 date: 2013-05-17
!
! abstract: This routine adds bogus raob so that at least one report
! at each model layer, by interpolate between a significant
! report and the neighboring obs
!
! program history log:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
! !USES:
use kinds, only: r_kind,r_single,r_double,i_kind
use constants, only: zero,one,one_tenth
use guess_grids, only: ges_psfcavg,ges_prslavg
use gridmod, only: nsig
use obsmod, only: bmiss
implicit none
! !INPUT PARAMETERS:
integer(i_kind) , intent(in ) ::kx
real(r_double) , intent(in ) ::vtcd
! !INPUT/OUTPUT PARAMETERS:
integer(i_kind) , intent(inout) ::levsio
real(r_double),dimension(13,255), intent(inout) :: obsdat
real(r_double),dimension(8,255), intent(inout) :: drfdat,qcmark,obserr
real(r_double),dimension(255,20), intent(inout) :: tpc
real(r_kind) wim,wi
real(r_kind),dimension(nsig) :: prsltmp,dpmdl
integer(i_kind) i,j,k,levs
integer(i_kind) ku,kl,ll,im
real rsig(nsig)
integer(i_kind),dimension(255):: pqm,qqm,tqm,wqm,cat,zqm
real(r_kind),dimension(255):: dpres,tvflg,dpobs
!!! find averaged sigma levels !!!!!!!!
levs=levsio
ll=levsio
do k=1,nsig
rsig(k)=ges_prslavg(k)/ges_psfcavg
enddo
do k=1,levs
cat(k)=nint(obsdat(10,k))
enddo
!!! find model levels at obs location in log(cb) !!!!!!!!
do k=1,nsig
dpmdl(k)=obsdat(1,1)*rsig(k)
prsltmp(k)=log(dpmdl(k)*one_tenth)
enddo
!!! find obs levels in log(cb) !!!!!!!!
do k=1,levs
dpres(k)=log(obsdat(1,k)*one_tenth)
dpobs(k)=dpres(k)
enddo
if(kx==120)then
pqm(1)=nint(min(qcmark(1,1),10000.0))
qqm(1)=nint(min(qcmark(2,1),10000.0))
tqm(1)=nint(min(qcmark(3,1),10000.0))
zqm(1)=nint(min(qcmark(4,1),10000.0))
call grdcrd(dpres,levs,prsltmp(1),nsig,-1)
do k=1,levs
tvflg(k)=one ! initialize as sensible
do j=1,20
if (tpc(k,j)==vtcd) tvflg(k)=zero ! reset flag if virtual
if (tpc(k,j)>=bmiss) exit ! end of stack
end do
end do
do i=2,levs
im=i-1
pqm(i)=nint(min(qcmark(1,i),10000.0))
qqm(i)=nint(min(qcmark(2,i),10000.0))
tqm(i)=nint(min(qcmark(3,i),10000.0))
zqm(i)=nint(min(qcmark(4,i),10000.0))
if ( (cat(i)==2 .or. cat(im)==2 .or. cat(i)==5 .or. cat(im)==5) .and. &
pqm(i)<4 .and. pqm(im)<4 )then
ku=dpres(i)-1
ku=min(nsig,ku)
kl=dpres(im)+2
kl=max(2,kl)
do k = kl,ku
ll=ll+1
if(ll>255)then
write(6,*)'error in SONDE_EXT levs > 255'
return
endif
obsdat(1,ll)=dpmdl(k)
qcmark(1,ll) =max (qcmark(1,i),qcmark(1,im)) !PQM
qcmark(2,ll) = bmiss
qcmark(3,ll) = bmiss
qcmark(4,ll) = bmiss
qcmark(5,ll) = bmiss
qcmark(7,ll) = bmiss
do j=1,20
tpc(ll,j)=tpc(i,j)
end do
wim=(prsltmp(k)-dpobs(i))/(dpobs(im)-dpobs(i))
wi=(dpobs(im)-prsltmp(k))/(dpobs(im)-dpobs(i))
!!! find tob, only bogus if both good obs and of the same type (sensible/virtual)
if( tqm(i)<4 .and. tqm(im)<4 .and. tvflg(i)==tvflg(im) ) then
obsdat(3,ll)=obsdat(3,im)*wim + obsdat(3,i)*wi
drfdat(1,ll) = drfdat(1,im)*wim + drfdat(1,i)*wi
drfdat(2,ll) = drfdat(2,im)*wim + drfdat(2,i)*wi
drfdat(3,ll) = drfdat(3,im)*wim + drfdat(3,i)*wi
qcmark(3,ll) =max (qcmark(3,i),qcmark(3,im)) !TQM
obserr(3,ll) =max (obserr(3,i),obserr(3,im)) ! TOE
endif
!!! find qob
if( qqm(i)<4 .and. qqm(im)<4 ) then
obsdat(2,ll)=obsdat(2,im)*wim + obsdat(2,i)*wi
drfdat(1,ll) = drfdat(1,im)*wim + drfdat(1,i)*wi
drfdat(2,ll) = drfdat(2,im)*wim + drfdat(2,i)*wi
drfdat(3,ll) = drfdat(3,im)*wim + drfdat(3,i)*wi
qcmark(2,ll) =max (qcmark(2,i),qcmark(2,im)) !QQM
obserr(2,ll) =max (obserr(2,i),obserr(2,im)) ! QOE
endif
!!! define zob
if( zqm(i)<4 .and. zqm(im)<4 ) then
obsdat(4,ll)=obsdat(4,im)*wim + obsdat(4,i)*wi
else
obsdat(4,ll)=max(obsdat(4,im),obsdat(4,i))
endif
qcmark(4,ll) =max (qcmark(4,i),qcmark(4,im)) !ZQM
enddo !kl,ku
endif !cat
enddo !levs
!!!!!!!!! w (not used) !!!!!!!!!!!!!!!!!!!!!!!!!!!
elseif(kx==220)then
pqm(1)=nint(min(qcmark(1,1),10000.0))
wqm(1)=nint(min(qcmark(5,1),10000.0))
call grdcrd(dpres,levs,prsltmp(1),nsig,-1)
do i=2,levs
im=i-1
wqm(i)=nint(min(qcmark(5,i),10000.0))
zqm(i)=nint(min(qcmark(4,i),10000.0))
pqm(i)=nint(min(qcmark(1,i),10000.0))
if( wqm(i)<4 .and. wqm(im)<4 .and. pqm(i)<4 .and. pqm(im)<4 .and.&
(cat(i)==2 .or. cat(im)==2 .or. cat(i)==5 .or. cat(im)==5) )then
ku=dpres(i)-1
ku=min(nsig,ku)
kl=dpres(im)+2
kl=max(2,kl)
do k = kl,ku
ll=ll+1
if(ll>255)then
write(6,*)'error in SONDE_EXT levs > 255'
return
endif
obsdat(1,ll)=dpmdl(k)
qcmark(1,ll) =max (qcmark(1,i),qcmark(1,im)) !PQM
qcmark(2,ll) = bmiss
qcmark(3,ll) = bmiss
qcmark(4,ll) = bmiss
qcmark(5,ll) = bmiss
qcmark(7,ll) = bmiss
wim=(prsltmp(k)-dpobs(i))/(dpobs(im)-dpobs(i))
wi=(dpobs(im)-prsltmp(k))/(dpobs(im)-dpobs(i))
!!! find wob (wint)
obsdat(5,ll)=obsdat(5,im)*wim + obsdat(5,i)*wi
obsdat(6,ll)=obsdat(6,im)*wim + obsdat(6,i)*wi
drfdat(1,ll) = drfdat(1,im)*wim + drfdat(1,i)*wi
drfdat(2,ll) = drfdat(2,im)*wim + drfdat(2,i)*wi
drfdat(3,ll) = drfdat(3,im)*wim + drfdat(3,i)*wi
qcmark(5,ll) =max (qcmark(5,i),qcmark(5,im)) !WQM
obserr(5,ll) =max (obserr(5,i),obserr(5,im)) ! WOE
qcmark(1,ll) =max (qcmark(1,i),qcmark(1,im))
!!! find zob
if( zqm(i)<4 .and. zqm(im)<4 ) then
obsdat(4,ll)=obsdat(4,im)*wim + obsdat(4,i)*wi
else
obsdat(4,ll)=max(obsdat(4,im),obsdat(4,i))
endif
enddo !kl,ku
endif !cat
enddo !levs
endif ! 120,220
!11 change the number of levels and output
levsio=ll
! End of routine
return
end subroutine sonde_ext