subroutine read_rapidscat(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_rapidscat read scatterometer winds ! prgmmr: Ling Liu date: 2015-04-03 ! ! abstract: This routine reads RapidScat scatterometer winds from BUFR dump. ! It also has options to thin the data by using conventional ! thinning programs ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 2015-04-03 Ling Liu ! 2015-09-17 Thomas - add l4densvar and thin4d to data selection procedure ! 2016-03-11 j. guo - Fixed {dlat,dlon}_earth_deg in the obs data stream ! 2020-05-04 wu - no rotate_wind for fv3_regional ! ! input argument list: ! ithin - flag to thin data ! rmesh - thinning mesh size (km) ! gstime - analysis time in minutes from reference date ! infile - unit from which to read BUFR data ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! twind - input group time window (hours) ! sis - satellite/instrument/sensor indicator ! ! output argument list: ! nread - number of satellite winds read ! ndata - number of satellite winds retained for further processing ! nodata - number of satellite winds retained for further processing ! nobs - array of observations on each subdomain for each processor ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_double,i_kind,r_single use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,& tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,& rlats,rlons,fv3_regional use qcmod, only: errormod,noiqc use convthin, only: make3grids,map3grids,del3grids,use_all use constants, only: deg2rad,zero,rad2deg,one_tenth,& tiny_r_kind,huge_r_kind,r60inv,one_tenth,& one,two,three,four,five,half,quarter,r60inv,r10,r100,r2000 ! use converr,only: etabl use obsmod, only: ran01dom,bmiss use convinfo, only: nconvtype, & icuse,ictype,icsubtype,ioctype, & ithin_conv,rmesh_conv,pmesh_conv use gsi_4dvar, only: l4dvar,iwinbgn,winlen,time_4dvar,l4densvar,thin4d use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 use mpimod, only: npe implicit none ! Declare passed variables character(len=*) ,intent(in ) :: infile,obstype character(len=*) ,intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata,nodata real(r_kind) ,intent(in ) :: twind real(r_kind),dimension(nlat,nlon,nsig),intent(in ) :: prsl_full integer(i_kind),dimension(npe) ,intent(inout) :: nobs ! Declare local parameters real(r_kind),parameter:: r1_2= 1.2_r_kind real(r_kind),parameter:: r3_33= 3.33_r_kind real(r_kind),parameter:: r6= 6.0_r_kind real(r_kind),parameter:: r50= 50.0_r_kind real(r_kind),parameter:: r54= 54.0_r_kind real(r_kind),parameter:: r55= 55.0_r_kind real(r_kind),parameter:: r56= 56.0_r_kind real(r_kind),parameter:: r70= 70.0_r_kind real(r_kind),parameter:: r85= 85.0_r_kind real(r_kind),parameter:: r90= 90.0_r_kind real(r_kind),parameter:: r105= 105.0_r_kind real(r_kind),parameter:: r110= 110.0_r_kind real(r_kind),parameter:: r125=125.0_r_kind real(r_kind),parameter:: r200=200.0_r_kind real(r_kind),parameter:: r250=250.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind real(r_kind),parameter:: r700=700.0_r_kind real(r_kind),parameter:: r199=199.0_r_kind real(r_kind),parameter:: r299=299.0_r_kind real(r_kind),parameter:: r421=421.0_r_kind real(r_kind),parameter:: r296=296.0_r_kind real(r_kind),parameter:: r799=799.0_r_kind real(r_kind),parameter:: r1200= 1200.0_r_kind real(r_kind),parameter:: r10000= 10000.0_r_kind ! Declare local variables logical outside logical luse,ithinp logical,allocatable,dimension(:,:):: lmsg ! set true when convinfo entry id found in a message character(70) obstr,hdrtr,wndstr character(8) subset character(8) c_prvstg,c_sprvstg integer(i_kind) ireadmg,ireadsb,iuse,mxtb,nmsgmax integer(i_kind) i,maxobs,idomsfc,nsattype integer(i_kind) nc,nx,isflg,itx,nchanl integer(i_kind) ntb,ntmatch,ncx,ncsave,ntread integer(i_kind) kk,klon1,klat1,klonp1,klatp1 integer(i_kind) nmind,lunin,idate,ilat,ilon,iret,k integer(i_kind) nreal,ithin,iout,ntmp,icount,iiout,ii integer(i_kind) itype,iosub,ixsub,isubsub,iobsub integer(i_kind) lim_qm integer(i_kind) nlevp ! vertical level for thinning integer(i_kind) pflag integer(i_kind) ntest,nvtest integer(i_kind) kl,k1,k2 integer(i_kind) nmsg ! message index integer(i_kind) qc1 integer(i_kind),dimension(nconvtype) :: ntxall integer(i_kind),dimension(nconvtype+1) :: ntx integer(i_kind),dimension(5):: idate5 integer(i_kind),allocatable,dimension(:):: isort,iloc,nrep integer(i_kind),allocatable,dimension(:,:)::tab integer(i_kind) ietabl,itypex,lcount,iflag,m real(r_single),allocatable,dimension(:,:,:) :: etabl real(r_kind) toff,t4dv real(r_kind) rmesh,ediff,usage,tdiff real(r_kind) u0,v0,uob,vob,dx,dy,dx1,dy1,w00,w10,w01,w11 real(r_kind) dlnpob,ppb,ppb2,qifn,qify,ee real(r_kind) woe,dlat,dlon,dlat_earth,dlon_earth,oelev real(r_kind) dlat_earth_deg,dlon_earth_deg real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00 real(r_kind) vdisterrmax,u00,v00,uob1,vob1 real(r_kind) del,werrmin,obserr,ppb1 real(r_kind) tsavg,ff10,sfcr,sstime,gstime,zz real(r_kind) crit1,timedif,xmesh,pmesh real(r_kind),dimension(nsig):: presl real(r_kind) uob_1,uob_2,uob_3,uob_4,vob_1,vob_2,vob_3,vob_4 real(r_kind) lkcs_1,lkcs_2,lkcs_3,lkcs_4 real(r_double),dimension(9):: hdrdat real(r_double),dimension(2):: satqc real(r_double),dimension(5):: obsdat real(r_double),dimension(3,4):: wnddat real(r_double),dimension(1,1):: r_prvstg,r_sprvstg real(r_kind),allocatable,dimension(:):: presl_thin real(r_kind),allocatable,dimension(:,:):: cdata_all real(r_kind),allocatable,dimension(:,:):: cdata_out ! equivalence to handle character names equivalence(r_prvstg(1,1),c_prvstg) equivalence(r_sprvstg(1,1),c_sprvstg) !******** Modify below from the bufrtable: data hdrtr /'SAID YEAR MNTH DAYS HOUR MINU SECO WS10 WD10'/ data obstr/'CLAT CLON WS10 WD10 SWVQ'/ data wndstr/'WS10 WD10 SWVQ'/ data ithin / -9 / data lunin / 21 / data rmesh / -99.999_r_kind / !************************************************************************** ! Return when SATWND are coming from prepbufr file ! if(use_prepb_satwnd) return ! Read observation error table ! itype 291 has been modified in the error table allocate(etabl(300,33,6)) etabl=1.e9_r_kind ietabl=19 open(ietabl,file='errtable',form='formatted') rewind ietabl etabl=1.e9_r_kind lcount=0 loopd : do read(ietabl,100,IOSTAT=iflag) itypex if( iflag /= 0 ) exit loopd lcount=lcount+1 do k=1,33 read(ietabl,110)(etabl(itypex,k,m),m=1,6) end do end do loopd 100 format(1x,i3) 110 format(1x,6e12.5) if(lcount<=0 ) then write(6,*)'READ_RAPIDSCAT: obs error table not available to 3dvar. the program will stop' call stop2(49) else write(6,*)'READ_RAPIDSCAT: observation errors provided by local file errtable' endif close(ietabl) ! Set lower limits for observation errors ! ** keep this way for now ! nreal keep the dimension of cdata_all werrmin=one nsattype=0 nreal=23 if (noiqc) then lim_qm=8 else lim_qm=4 endif ! ** read convtype from convinfo file ! ** only read in rapidsat 296 for now ** ntread=1 ntmatch=0 ntx(ntread)=0 ntxall=0 do nc=1,nconvtype if(trim(ioctype(nc)) == 'uv' .and. ictype(nc) ==296) then ntmatch=ntmatch+1 ntxall(ntmatch)=nc ithin=ithin_conv(nc) if(ithin > 0)then ntread=ntread+1 ntx(ntread)=nc end if end if end do if(ntmatch == 0)then write(6,*) ' READ_RAPIDSCAT: no matching obstype found in convinfo ',obstype return end if !! go through the satedump to find out how many subset to process !** Open and read data from bufr data file open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) !! 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 msg_report: do while (ireadmg(lunin,subset,idate) == 0) ! Time offset if(nmsg == 0) call time_4dvar(idate,toff) nmsg=nmsg+1 if (nmsg>nmsgmax) then write(6,*)'READ_RAPIDSCAT: messages exceed maximum ',nmsgmax call stop2(49) endif loop_report: do while (ireadsb(lunin) == 0) ntb = ntb+1 maxobs=maxobs+1 nrep(nmsg)=nrep(nmsg)+1 if (ntb>mxtb) then write(6,*)'READ_RAPIDSCAT: reports exceed maximum ',mxtb call stop2(49) endif !** Extract sat ID information from BUFR and assign type ! This part will extract information from the bufrtable ! bufrtable need to be modified including the rapidscat data entry ! iobsub is the prepbufr subtype, still part of the convinfo file !********* iobsub=0 for rapidscat* call ufbint(lunin,hdrdat,9,1,iret,hdrtr) ! determine the satellite wind type ! 296: rapidscat data iobsub=0 itype=-1 if(trim(subset) == 'NC012255') then if( hdrdat(1) == r296 ) then ! rapidscat data itype=296 endif endif ! Match ob to proper convinfo type ncsave=0 matchloop:do ncx=1,ntmatch nc=ntxall(ncx) if (itype /= ictype(nc)) cycle matchloop ! 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 maxobs=maxobs+1 nx=1 if(ithin_conv(ncsave) > 0)then do ii=2,ntread if(ntx(ii) == ncsave)nx=ii end do end if tab(ntb,1)=ncsave tab(ntb,2)=nx tab(ntb,3)=1 lmsg(nmsg,nx) = .true. end if enddo loop_report enddo msg_report ! 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 !! read satellite winds one type a time ! same as in the read_prepbufr.f90 file loop_convinfo: do nx=1,ntread use_all = .true. ithin=0 if(nx >1) then nc=ntx(nx) ithin=ithin_conv(nc) if (ithin > 0 ) then rmesh=rmesh_conv(nc) pmesh=pmesh_conv(nc) use_all = .false. if(pmesh > zero) then pflag=1 nlevp=r1200/pmesh else pflag=0 nlevp=nsig endif xmesh=rmesh call make3grids(xmesh,nlevp) if (.not.use_all) then allocate(presl_thin(nlevp)) if (pflag==1) then do k=1,nlevp presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth enddo endif endif write(6,*)'READ_RAPIDSCAT: ictype(nc),rmesh,pflag,nlevp,pmesh,nc ',& ioctype(nc),ictype(nc),rmesh,pflag,nlevp,pmesh,nc endif endif call closbf(lunin) close(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) ! Big loop over BUFR file ntb = 0 nmsg = 0 loop_msg: do while(ireadmg(lunin,subset,idate) == 0) nmsg = nmsg+1 if(.not.lmsg(nmsg,nx)) then 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 nc=tab(ntb,1) if(nc <= 0 .or. tab(ntb,2) /= nx) cycle loop_readsb hdrdat=bmiss obsdat=bmiss wnddat=bmiss satqc=bmiss iobsub=0 itype=-1 uob=bmiss vob=bmiss ppb=bmiss ppb1=bmiss ppb2=bmiss uob1=bmiss vob1=bmiss ee=r110 qifn=r110 qify=r110 uob_1=bmiss vob_1=bmiss uob_2=bmiss vob_2=bmiss uob_3=bmiss vob_3=bmiss uob_4=bmiss vob_4=bmiss lkcs_1=bmiss lkcs_2=bmiss lkcs_3=bmiss lkcs_4=bmiss ! Extract type, date, and location information call ufbint(lunin,hdrdat,9,1,iret,hdrtr) call ufbint(lunin,obsdat,5,1,iret,obstr) call ufbrep(lunin,wnddat,3,4,iret,wndstr) !** potential can reject bad cell, etc, place holder for now ! reject the data with bad quality mark from SDM !** cycle loop means skip to the next record if(hdrdat(1) /= r296) cycle loop_readsb ! Compare relative obs time with window. If obs ! falls outside of window, don't use this obs idate5(1) = hdrdat(2) !year idate5(2) = hdrdat(3) ! month idate5(3) = hdrdat(4) ! day idate5(4) = hdrdat(5) ! hours idate5(5) = hdrdat(6) ! minutes call w3fs21(idate5,nmind) t4dv = real((nmind-iwinbgn),r_kind)*r60inv sstime = real(nmind,r_kind) tdiff=(sstime-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle loop_readsb else if (abs(tdiff)>twind) cycle loop_readsb endif ! determine the satellite wind type as in prepbufr ! 296: rapidscat winds iosub=0 if(abs(obsdat(1)) >r90 ) cycle loop_readsb if(obsdat(2) r360) cycle loop_readsb if(obsdat(5) >1) cycle loop_readsb if(trim(subset) == 'NC012255') then ! rapidscat wind if( hdrdat(1) == r296) then itype=296 endif endif nread=nread+2 dlon_earth_deg = obsdat(2) dlat_earth_deg = obsdat(1) dlon_earth=obsdat(2)*deg2rad dlat_earth=obsdat(1)*deg2rad ! If regional, map obs lat,lon to rotated grid. if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if(diagnostic_reg) then call txy2ll(dlon,dlat,rlon00,rlat00) ntest=ntest+1 cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* & (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00)) cdist=max(-one,min(cdist,one)) disterr=acos(cdist)*rad2deg disterrmax=max(disterrmax,disterr) end if if(outside) cycle loop_readsb ! check to see if outside regional else dlon=dlon_earth dlat=dlat_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! If rapidscat data, determine primary surface type. If not open sea, ! skip this observation. This check must be done before thinning. ! isflg - surface flag 0:sea 1:land 2:sea ice 3:snow 4:mixed if (itype==296) 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 !! convert from wind direction and speed to u,v component uob=obsdat(3)*sin(obsdat(4)*deg2rad) vob=obsdat(3)*cos(obsdat(4)*deg2rad) qc1=obsdat(5) !! Get observation error from PREPBUFR observation error table ! only need read the 4th column for type 291 from the right ppb=max(zero,min(ppb,r2000)) if(ppb>=etabl(itype,1,1)) k1=1 do kl=1,32 if(ppb>=etabl(itype,kl+1,1).and.ppb<=etabl(itype,kl,1)) k1=kl end do if(ppb<=etabl(itype,33,1)) k1=5 k2=k1+1 ediff = etabl(itype,k2,1)-etabl(itype,k1,1) if (abs(ediff) > tiny_r_kind) then del = (ppb-etabl(itype,k1,1))/ediff else del = huge_r_kind endif del=max(zero,min(del,one)) obserr=(one-del)*etabl(itype,k1,4)+del*etabl(itype,k2,4) obserr=max(obserr,werrmin) ! Set usage variable usage = 0 iuse=icuse(nc) if(iuse <= 0)usage=r100 ! Get information from surface file necessary for conventional data here ! This is different from the previous sfc_type call call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) !! process the thining procedure ithin=ithin_conv(nc) ithinp = ithin > 0 .and. pflag /= 0 if(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 end if ! The following two lines are new in SATWND: dlnpob=log(one_tenth*ppb) ! ln(pressure in cb) ppb=one_tenth*ppb ! from mb to cb ! Special block for data thinning - if requested if (ithin > 0 .and. iuse >=0) then ntmp=ndata ! counting moved to map3gridS ! Set data quality index for thinning if (thin4d) then timedif = zero else timedif=abs(t4dv-toff) endif crit1 = timedif/r6+half if (pflag==0) then do kk=1,nsig presl_thin(kk)=presl(kk) end do endif call map3grids(-1,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,& ppb,crit1,ndata,iout,ntb,iiout,luse,.false.,.false.) if (.not. luse) cycle loop_readsb if(iiout > 0) isort(iiout)=0 if (ndata > ntmp) then nodata=nodata+2 endif isort(ntb)=iout else ndata=ndata+1 nodata=nodata+2 iout=ndata isort(ntb)=iout endif woe=obserr oelev=r10 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 ! Output result to array, need restruct from here. 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 ! index of height cdata_all(6,iout)=uob ! u obs cdata_all(7,iout)=vob ! v obs cdata_all(8,iout)=ndata ! station id cdata_all(9,iout)=t4dv ! time cdata_all(10,iout)=nc ! index of type in convinfo file cdata_all(11,iout)=0 ! index of station elevation cdata_all(12,iout)=qc1 ! index of quality mark cdata_all(13,iout)=obserr ! original obs error cdata_all(14,iout)=usage ! usage parameter cdata_all(15,iout)=idomsfc ! dominate surface type cdata_all(16,iout)=9999999 ! 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 enddo loop_readsb enddo loop_msg ! Deallocate arrays used for thinning data if (.not.use_all) then deallocate(presl_thin) call del3grids endif ! Normal exit enddo loop_convinfo! loops over convinfo entry matches deallocate(lmsg,tab,nrep) call closbf(lunin) ! 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,*) ' READ_RAPIDSCAT: mix up in read_satwnd ,ndata,icount ',ndata,icount call stop2(49) 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) deallocate(etabl) call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) cdata_out deallocate(cdata_out) 900 continue if(diagnostic_reg .and. ntest>0) write(6,*)'READ_RAPIDSCAT: ',& 'ntest,disterrmax=',ntest,disterrmax if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_RAPIDSCAT: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax write(6,*)'READ_RAPIDSCAT: closbf(',lunin,')' write(6,*) 'READ_RAPIDSCAT,nread,ndata,nreal,nodata=',nread,ndata,nreal,nodata close(lunin) ! End of routine return end subroutine read_rapidscat