         subroutine get_streamlined_stats_cv_mpi(cgrid, & 
                        ista,iend,jsta,jend,ribuffer_km,mype,npe)
!
!  to compile:
!    xlf90 -c get_streamlined_stats_cv.f -L/nwprod/lib -lw3_8
!
!  abstract: 
!    read in from original unformatted gsi diagnostic files 
!    and write out formatted streamlined files. Valid for
!    gsi release 200609.
  
         use mpi
         use cressanl_common, only: nobsmax,nmax

         implicit none

!Declare passed variables
         integer(4),intent(in):: mype,npe
         integer(4),intent(in):: ista,iend,jsta,jend
         real(4),intent(in):: ribuffer_km
         character(60),intent(in):: cgrid

!Declare local variables
         real(8) alat1,elon1,dx,elonv,alatan
         real(4) dg2rad,elonv4,alatan4,xn

         integer(4),parameter::lun=10

         integer(4) idate,nchar,nreal,i,ii,mypegsi
         character(8),allocatable,dimension(:):: cdiagbuf
         character(8),allocatable,dimension(:):: cprvstg
         character(8),allocatable,dimension(:):: csprvstg
         real(4),allocatable,dimension(:,:)::rdiagbuf
         real(4) rlat,rlon,oberr,ob,ob_model,ddiff,ddifforig, & 
                 dudiff,dvdiff,rmuse,xx,yy
         real(4) uob,uob_model,vob,vob_model,rfactor,qtflg
         real(4) uobe,uobe_model,vobe,vobe_model
         real(4) wobe,wobe_model,wdire,wdire_model
         real(4) t0,t0_model,td0,td0_model
         real(4) p0,p0_model
         real(4) shgt0  !station height
         real(4) hgt0  !observation elevation 
         real(4) hgt   !model terrain at ob location
         real(4) slm
         real(4) dtime
         real(4) rtvts
         real(4) bb
         real(4) pgesorig
         character(3) otype
         character(8) obstype
         real(8) rlat8,rlon8,xx8,yy8
         character(8) cstation,cstationmax
         character(8) cprovider,csubprovider
         integer(4) itype
         integer(4) imax
         logical fexist, r15887_diagfile_fmt
         logical muse
         logical lboundtd
         logical lambconform
         logical polarstereo

         integer(4),dimension(npe):: nts,nqs, & 
                    nps,nws,nspds,ngusts,nvis,npblhs,ndists
         integer(4),dimension(npe):: nus,nvs

         integer(4),dimension(npe):: nts2,nqs2, & 
                    nps2,nws2,nspds2,ngusts2,nvis2,npblhs2
         integer(4) n,ntsmax,nqsmax,npsmax,nwsmax, & 
                    nspdsmax,ngustsmax,nvismax,npblhsmax

         integer(4) n_t,n_p,n_q,n_w,n_spd,n_gust,n_vis,n_pblh
         integer(4) nmax0

         real(4) xsta,xend,ysta,yend
         real(4) xsta2,xend2

         integer(4) ibuffer

         character(2) clun
         logical insubdom
         logical insubdom2

         integer(4) ierror

         integer(4)   nx,ny
         real(4) da
         character(3) gesoranl
         real(4),allocatable,dimension(:,:)::pgrid1

         namelist/diagfile_fmt/r15887_diagfile_fmt
!*******************************************************************************
         data r15887_diagfile_fmt/.true./

         inquire(file='diagfile_fmt_input',exist=fexist)
         if (fexist) then
            open (55,file='diagfile_fmt_input',form='formatted')
            read(55,*) r15887_diagfile_fmt
            close(55)
         endif
         if (mype==0) print*,'in get_streamlined_stats_cv: r15887_diagfile_fmt=',r15887_diagfile_fmt

         call proj_info(cgrid,dx,alat1,elon1,elonv,alatan)

         if (mype==0) then
            print*,'in get_streamlined_stats_cv: cgrid=',trim(cgrid)
            print*,'in get_streamlined_stats_cv: alat1=',alat1
            print*,'in get_streamlined_stats_cv: elon1=',elon1
            print*,'in get_streamlined_stats_cv: dx=',dx
            print*,'in get_streamlined_stats_cv: elonv=',elonv
            print*,'in get_streamlined_stats_cv: alatan=',alatan
         endif

         call domain_dims(cgrid,nx,ny,da)

         if (mype==0) then
            print*,'in get_streamlined_stats_cv: cgrid=',trim(cgrid)
            print*,'in get_streamlined_stats_cv: nx,ny,da=',nx,ny,da
         endif

         lambconform=trim(cgrid)=='conus'.or.trim(cgrid)=='cohres'.or.trim(cgrid)=='hrrr'.or.trim(cgrid)=='cohresext'
         polarstereo=trim(cgrid)=='alaska'.or.trim(cgrid)=='akhres'.or.trim(cgrid)=='juneau'

         dg2rad=atan(1.)/45.
         elonv4=elonv
         alatan4=alatan

         xn=1.
         if (lambconform) xn=sin(alatan4*dg2rad)

         if (mype==0) print*,'in get_streamlined_stats_cv: lambconform,polarstereo,xn=', & 
                                                           lambconform,polarstereo,xn

         print*,'in get_streamlined_stats_cv: mype,npe,ista,iend,jsta,jend=', & 
                                              mype,npe,ista,iend,jsta,jend

         if (mype==0) print*,'in get_streamlined_stats_cv: ribuffer_km=',ribuffer_km

         xsta=float(ista)
         xend=float(iend)
         ysta=float(jsta)
         yend=float(jend)

!====================================================================================
!==> 
!====================================================================================
         allocate(pgrid1(nx,ny))

         gesoranl='ges'
         call getguess_or_anal('psfc',gesoranl,pgrid1,nx,ny)

         if (mype==0) print*,'in get_streamlined_stats_cv: pgrid1,min,max=',minval(pgrid1),maxval(pgrid1)

         ibuffer=1000.*ribuffer_km/da

         if (mype==0) print*,'in get_streamlined_stats_cv: number of ibuffer points in x: ', ibuffer

         xsta2=float(max(1,ista-ibuffer))
         xend2=float(min(nx,iend+ibuffer))

         print*,'in get_streamlined_stats_cv: mype,xsta,xend,xsta2,xend2=', &
                                              mype,xsta,xend,xsta2,xend2

!====================================================================================
         write(clun,'(i2.2)') mype

         open (10,file='stats_ges_cv_'//clun,form='formatted')          !output file

         open (7,file='diag_conv_ges.dat',form='unformatted')           !input file

         read(7) idate

         nts(:)=0 ; nqs(:)=0 ; nps(:)=0 ; nws(:)=0 ; nspds(:)=0
         ngusts(:)=0 ; nvis(:)=0 ; npblhs(:)=0 ; ndists(:)=0

         n_t=0 ; n_p=0 ; n_q=0 ; n_w=0 ; n_spd=0
         n_gust=0 ; n_vis=0 ; n_pblh=0  

         nobsmax=0

         bb=0. !dummy variable
100      continue
         read(7,end=200)otype,nchar,nreal,ii,mypegsi

         print*,'in get_streamlined_stats_cv, idate=',idate
         print*,'in get_streamlined_stats_cv, otype,nchar,nreal,ii,mypegsi=', & 
                    otype,nchar,nreal,ii,mypegsi

         if (ii==0) then
            read(7,end=200)
            if (r15887_diagfile_fmt) read(7,end=200)
            goto 100
         endif
  
         allocate(cdiagbuf(ii))
         allocate(cprvstg(ii))
         allocate(csprvstg(ii))
         allocate(rdiagbuf(nreal,ii))

         if (r15887_diagfile_fmt) then
             read(7) cdiagbuf,rdiagbuf
             read(7) cprvstg,csprvstg
           else
             read(7) cdiagbuf,rdiagbuf,cprvstg,csprvstg
         endif

         do i=1,ii
           if (cprvstg(i)(1:4)=='B7Hv')   cprvstg(i)(5:8)='   ' !this provider name comes with weird
           if (csprvstg(i)(1:4)=='B7Hv') csprvstg(i)(5:8)='   ' !trailing characters
         enddo


        if (otype(1:3)=='  t') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(20,i)
           hgt=rdiagbuf(21,i)
           rtvts=qtflg
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
              nobsmax=nobsmax+1
              n_t=n_t+1
              if (insubdom) nts(mype+1)=nts(mype+1)+1
              obstype='t-ob'
              write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
              write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                           ob,ob_model,qtflg,rtvts,shgt0,hgt0,hgt,slm,rmuse,insubdom
            endif
          enddo
        end if


        if (otype(1:3)=='  q') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           p0=rdiagbuf(6,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             n_q=n_q+1
             if (insubdom) nqs(mype+1)=nqs(mype+1)+1
             obstype='q-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                          ob,ob_model,qtflg,p0,shgt0,hgt0,hgt,slm,rmuse,insubdom

             call bilinear_2d0v2(pgrid1,1,nx,1,ny,p0_model,yy,xx)!Note the reverse order "yy,xx"

             t0=-9999. ; t0_model=-9999.
             lboundtd=.false.
             call get_stndewpt_r4(100.*p0,ob,t0,td0,lboundtd)
             call get_stndewpt_r4(100.*p0_model,ob_model,t0_model,td0_model,lboundtd)

             nobsmax=nobsmax+1
             obstype='td-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, &
                           td0,td0_model,qtflg,p0,shgt0,hgt0,hgt,slm,rmuse,insubdom

           endif
          enddo
        end if

        if (otype(2:3)=='ps') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i) !ges used in analysis (coverted to hPa)
           ddifforig=rdiagbuf(19,i)
           pgesorig=rdiagbuf(17,i)-rdiagbuf(19,i) !ges w/o adjustment to guess surface pressure (hPa)
           slm=rdiagbuf(20,i)
           hgt=rdiagbuf(21,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             n_p=n_p+1
             if (insubdom) nps(mype+1)=nps(mype+1)+1
             obstype='p-ob'
             oberr=(1./oberr)*100. !convert to Pa
             ob=ob*100. !convert to Pa
             ob_model=ob_model*100.! convert to Pa
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,oberr, & 
                          ob,ob_model,qtflg,pgesorig,shgt0,hgt0,hgt,slm,rmuse,insubdom
           endif
          enddo
        end if


        if (otype(2:3)=='uv') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           uob=rdiagbuf(17,i)
           dudiff=rdiagbuf(18,i)
           uob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           vob=rdiagbuf(20,i)
           dvdiff=rdiagbuf(21,i)
           vob_model=rdiagbuf(20,i)-rdiagbuf(21,i)

           uobe=uob
           vobe=vob
           uobe_model=uob_model
           vobe_model=vob_model
          
!          no need for wind rotation anymore since wind-components in 
!          diagnostic file are now earth-relative !MPondeca, 23Jan2012
!
!          if (lambconform .or. polarstereo) then
!               call get_earth_winds(uobe,vobe,uob,vob,xn,elonv4,rlon,dg2rad)
!               call get_earth_winds(uobe_model,vobe_model,uob_model,vob_model,xn,elonv4,rlon,dg2rad)
!          endif

           call speeddir(uobe,vobe,wdire,wobe)
           call speeddir(uobe_model,vobe_model,wdire_model,wobe_model)

           slm=rdiagbuf(24,i)
           hgt=rdiagbuf(25,i)
           rfactor=rdiagbuf(23,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             if (insubdom) nws(mype+1)=nws(mype+1)+1
             obstype='u-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*)  itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                           uob,uob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom

             nobsmax=nobsmax+1
             obstype='v-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                           vob,vob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom

             nobsmax=nobsmax+1
             n_w=n_w+1
             obstype='w-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                           wobe,wobe_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom

             nobsmax=nobsmax+1
             obstype='w2-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                           wobe,wobe_model,uob,vob,shgt0,hgt0,hgt,slm,rmuse,insubdom

             nobsmax=nobsmax+1
             obstype='wd-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                           wdire,wdire_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom

             if (lambconform .or. polarstereo) then
                  uobe=uob
                  vobe=vob
                  uobe_model=uob_model
                  vobe_model=vob_model

                  call get_grid_winds(uobe,vobe,uob,vob,xn,elonv4,rlon,dg2rad)
                  call get_grid_winds(uobe_model,vobe_model,uob_model,vob_model,xn,elonv4,rlon,dg2rad)

                  nobsmax=nobsmax+1
                  obstype='ugrel-ob'
                  write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
                  write(lun,*)  itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                                uob,uob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom

                  nobsmax=nobsmax+1
                  obstype='vgrel-ob'
                  write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
                  write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                                vob,vob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom
             endif

           endif
          enddo
        end if

        if (otype(1:3)=='spd') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           rfactor=rdiagbuf(20,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             n_spd=n_spd+1
             if (insubdom) nspds(mype+1)=nspds(mype+1)+1
             obstype='spd-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                          ob,ob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom
           endif
          enddo
        end if

        if (otype(1:3)=='gst') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           rfactor=rdiagbuf(20,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             n_gust=n_gust+1
             if (insubdom) ngusts(mype+1)=ngusts(mype+1)+1
             obstype='gust-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                          ob,ob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom
           endif
          enddo
        end if

        if (otype(1:3)=='vis') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           rfactor=rdiagbuf(20,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             n_vis=n_vis+1
             if (insubdom) nvis(mype+1)=nvis(mype+1)+1
             obstype='vis-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                          ob,ob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom
           endif
          enddo
        end if

        if (otype(1:3)=='pbl') then
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider='xxxxxxxx'         !cprvstg(i)
           csubprovider='xxxxxxxx'      !=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           rfactor=rdiagbuf(20,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           slm=-9999.                 !rdiagbuf(21,i)
           hgt=-9999.                 !rdiagbuf(22,i)
           call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8

           insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
           insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend

           if (insubdom2 .and. oberr.gt.1.e-05) then
             nobsmax=nobsmax+1
             n_pblh=n_pblh+1
             if (insubdom) npblhs(mype+1)=npblhs(mype+1)+1
             obstype='pblh-ob'
             write(lun,'(4a8)') cstation,cprovider,csubprovider,obstype
             write(lun,*) itype,rlat,rlon,dtime,xx,yy,1./oberr, & 
                          ob,ob_model,rfactor,bb,shgt0,hgt0,hgt,slm,rmuse,insubdom
           endif
          enddo
        end if
 
        deallocate(cdiagbuf,rdiagbuf)
        deallocate(cprvstg,csprvstg)
        goto 100
200     continue
        close(lun)
        close(7)

        print*,'in get_streamlined_stats_cv: mype,n_t,n_p,n_q,n_w,n_spd,n_gust,n_vis,n_pblh=',& 
                                             mype,n_t,n_p,n_q,n_w,n_spd,n_gust,n_vis,n_pblh

        nmax0=max(n_t,n_p,n_q,n_w,n_spd,n_gust,n_vis,n_pblh)

        call mpi_allreduce(nmax0,nmax,npe,mpi_integer,mpi_max,mpi_comm_world,ierror)
        if (mype==0) print*,'in get_streamlined_stats_cv: nmax=',nmax

        print*,'in get_streamlined_stats_cv: mype,nobsmax=',mype,nobsmax
        print*
 
        nts2(:)=nts(:)
        nqs2(:)=nqs(:)
        nps2(:)=nps(:)
        nws2(:)=nws(:)
        nspds2(:)=nspds(:)
        ngusts2(:)=ngusts(:)
        nvis2(:)=nvis(:)
        npblhs2(:)=npblhs(:)

        call mpi_allreduce(nts2,    nts,    npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(nqs2,    nqs,    npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(nps2,    nps,    npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(nws2,    nws,    npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(nspds2,  nspds,  npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(ngusts2, ngusts, npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(nvis2,   nvis,   npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)
        call mpi_allreduce(npblhs2, npblhs, npe, mpi_integer, mpi_sum, mpi_comm_world, ierror)

        ntsmax=maxval(nts)
        nqsmax=maxval(nqs)
        npsmax=maxval(nps)
        nwsmax=maxval(nws)
        nspdsmax=maxval(nspds)
        ngustsmax=maxval(ngusts)
        nvismax=maxval(nvis)
        npblhsmax=maxval(npblhs)
 
        nus(:)=nws(:)
        nvs(:)=nws(:)

        if (mype==0) then
           print*,'----------in get_streamlined_stats_cv:-------------'
           print*,'nts(:),ntsmax,all = '       ,(nts(n),n=1,npe),ntsmax,sum(nts(:))
           print*,'nqs(:),nqsmax,all = '       ,(nqs(n),n=1,npe),nqsmax,sum(nqs(:))
           print*,'nps(:),npsmax,all = '       ,(nps(n),n=1,npe),npsmax,sum(nps(:))
           print*,'nws(:),nwsmax,all = '       ,(nws(n),n=1,npe),nwsmax,sum(nws(:))
           print*,'nspds(:),nspdsmax,all ='    ,(nspds(n),n=1,npe),nspdsmax,sum(nspds(:))
           print*,'ngusts(:),ngustsmax,all = ' ,(ngusts(n),n=1,npe),ngustsmax,sum(ngusts(:))
           print*,'nvis(:),nvismax,all = '     ,(nvis(n),n=1,npe),nvismax,sum(nvis(:))
           print*,'npblhs(:),npblhsmax,all = ' ,(npblhs(n),n=1,npe),npblhsmax,sum(npblhs(:))
!          print*,'nmax=',nmax
        endif

        deallocate(pgrid1)

end subroutine get_streamlined_stats_cv_mpi
!!***************************************************************
!      subroutine get_earth_winds(un,vn,ut,vt,ROTCON_P,LON_XX_P, &
!                                 olon,dg2rad)
!!
!!     Abstract: Rotates wind relative to grid on Lambert Conformal
!!               or polar streo projection to get winds with respect
!!               to the true north.
!!     Adapted from Stan Benjamin's code snippet
!
!      implicit none
!
!!Declare passed variables
!      real(4),intent(in)::  ut,vt,ROTCON_P,LON_XX_P,olon,dg2rad
!      real(4),intent(out):: un,vn
!
!!Declare local variables
!      real(4)  angle2,sinx2,cosx2
!
!!**  ROTCON_P          R  WIND ROTATION CONSTANT, = 1 FOR POLAR STEREO
!!**                         AND SIN(LAT_TAN_P) FOR LAMBERT CONFORMAL
!!**  LAT_TAN_P         R  LATITUDE AT LAMBERT CONFORMAL PROJECTION
!!**                         IS TRUE (DEG)
!!**  LON_XX_P          R  MERIDIAN ALIGNED WITH CARTESIAN X-AXIS(DEG)
!
!           angle2 = rotcon_p*(olon-lon_xx_p)*dg2rad
!           sinx2 = sin(angle2)
!           cosx2 = cos(angle2)
!           un = cosx2*ut+sinx2*vt
!           vn =-sinx2*ut+cosx2*vt
!
!           return
!           end
!!***************************************************************
!***************************************************************
!note:

!***************************************************************
!        if (otype(3:3)=='t') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
! 
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!   
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = prest              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(iobshgt,i)    ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)

!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = data(iqt,i)        ! setup qc or event mark (currently qtflg only)
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one
!       endif
!
!       err_input = data(ier2,i)
!       err_adjst = data(ier,i)
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input=one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst=one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final=one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (K**-1)
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (K**-1)
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (K**-1)
!
!       rdiagbuf(17,ii) = data(itob,i)       ! temperature observation (K)
!       rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis (K)
!       rdiagbuf(19,ii) = tob-tges           ! obs-ges w/o bias correction (K) (future slot)
!***************************************************************

 
!***************************************************************
!       if (otype(3:3)=='q') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!   
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = presq              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(iobshgt,i)    ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
! 
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark 
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one                    
!       endif
! 
!       err_input = data(ier2,i)*qsges            ! convert rh to q
!       err_adjst = data(ier,i)*qsges             ! convert rh to q
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse observation error
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error
!
!       rdiagbuf(17,ii) = data(iqob,i)       ! observation
!       rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis
!       rdiagbuf(19,ii) = qob-qges           ! obs-ges w/o bias correction (future slot)
!
!       rdiagbuf(20,ii) = qsges              ! guess saturation specific humidity
!       
!***************************************************************


!***************************************************************
!       if (otype(2:3)=='ps') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = data(ipres,i)*r10  ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = dhgt               ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
!
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one                    
!       endif
!
!       pob      = pob*r10
!       pges     = pges*r10
!       pgesorig = pgesorig*r10
!
!       err_input = data(ier2,i)*r10   ! r10 converts cb to mb
!       err_adjst = data(ier,i)*r10
!       if (ratio_errors*error/r10>tiny_r_kind) then
!          err_final = r10/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (hPa**-1)
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (hPa**-1)
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (hPa**-1)
!
!       rdiagbuf(17,ii) = pob                ! surface pressure observation (hPa)
!       rdiagbuf(18,ii) = pob-pges           ! obs-ges used in analysis (coverted to hPa)
!       rdiagbuf(19,ii) = pob-pgesorig       ! obs-ges w/o adjustment to guess surface pressure (hPa)
!***************************************************************


!***************************************************************
!       if (otype(2:3)=='uv') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!   
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(ielev,i)      ! station elevation (meters)
!       rdiagbuf(6,ii)  = presw              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(ihgt,i)       ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
!
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one
!       endif
!
!       err_input = data(ier2,i)
!       err_adjst = data(ier,i)
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (m/s)**-1
!
!       rdiagbuf(17,ii) = data(iuob,i)       ! u wind component observation (m/s)
!       rdiagbuf(18,ii) = dudiff             ! u obs-ges used in analysis (m/s)
!       rdiagbuf(19,ii) = uob-ugesin         ! u obs-ges w/o bias correction (m/s) (future slot)
!
!       rdiagbuf(20,ii) = data(ivob,i)       ! v wind component observation (m/s)
!       rdiagbuf(21,ii) = dvdiff             ! v obs-ges used in analysis (m/s)
!       rdiagbuf(22,ii) = vob-vgesin         ! v obs-ges w/o bias correction (m/s) (future slot)
!
!       rdiagbuf(23,ii) = factw              ! 10m wind reduction factor
!***************************************************************


!***************************************************************
!       if (otype(1:3)=='spd') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!    
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = presw              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(ihgt,i)       ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
!
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one
!       endif
!
!       spdob0    = sqrt(data(iuob,i)*data(iuob,i)+data(ivob,i)*data(ivob,i))
!       err_input = data(ier2,i)
!       err_adjst = data(ier,i)
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (m/s)**-1
!
!       rdiagbuf(17,ii) = spdob              ! wind speed observation (m/s)
!       rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis (m/s)
!       rdiagbuf(19,ii) = spdob0-spdges      ! obs-ges w/o bias correction (m/s) (future slot)
!
!       rdiagbuf(20,ii) = factw              ! 10m wind reduction factor
!***************************************************************
