subroutine spruv(wdatain,nn_w,ugesin,vgesin, & 1,3
  fact,wtype,nwdta,nsig,ermax,ermin,gross,lat12,&
  nx,nlon,weight,mype,awork,wpres,nlev,ptop,pbot,vwork)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    spruv       save wind infor.        
!   prgmmr: parrish          org: w/nmc22    date: 90-10-13
!
! abstract: save wind information for later use.
!
! program history log:
!   90-10-13  parrish
!   98-04-07  weiyu yang
!   99-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!
!   input argument list:
!     wdatain  - obs info at obs locations in each pe sub-domain
!     nn_w     - number of data in each wind record
!     fact     - reduction to 10,20 m factor
!     wtype    - wind type
!     nwdta    - number of obs
!     nsig     - number of layers on gaussian grid
!     ermax,ermin,gross - parameters for gross error test
!     lat12    - lat1+2
!     nlon     - number of longitudes
!     weight   - location weighting
!     mype     - pe number
!
!   output argument list:
!     j_w      - wind data location indices
!     wgt_w    - weighting information of wind data
!     other_w  - other information of wind data, other_w
!                is multiplied by aerr already.
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
 use type_kinds, only: fp_kind
 use obsmod
 implicit none

 include 'mpi_inc.h'
 include 'ijstart.h'
 integer iout_uv,mype_uv
 integer jlat,jlon,jsig,ict,icts,m1,n1,n2,ier,ilon,ilat,isig
 integer iures,ivres,iii,k,ioff,i,lat12,nn_w,nwdta,nsig,iix
 integer itypx,itype,latlon11_l,nlev,nx,nlon,mype
 real(fp_kind) dx,dy,ds,dx1,dy1,ds1,scale,ratio,obserror,obserrlm,small
 real(fp_kind) residual,val,ressw,ress,val1,val2,ermin,gross,ermax
 real(fp_kind),dimension(nwdta,nn_w):: wdatain
 real(fp_kind),dimension(nwdta):: wtype,fact,ugesin,vgesin
 real(fp_kind),dimension(nwdta):: weight,wpres
 real(fp_kind),dimension(100+7*nsig,6):: awork
 real(fp_kind),dimension(6):: ptop,pbot
 real(fp_kind),dimension(6,101,3):: vwork

 small=1.e-10
 nwdata_dim=max(nwdata_dim,1)
 nspddat_dim=max(nspddat_dim,1)

 call makew

 m1=mype+1; n1=nsig+3; n2=n1+nsig
 ier=1; ilon=2; ilat=3; isig=4; iures=5; ivres=6
 do iii=1,nwdta
   if(wdatain(iii,isig)<1.0) wdatain(iii,isig)=1.0
 end do
 ict = 0
 icts = 0
 do iix=1,2
 do iii=1,nwdta
  itypx=nint(wtype(iii))
  if(itypx/=283.and.iix==1)then
  jlat=wdatain(iii,ilat)
  jlon=wdatain(iii,ilon)
  jsig=wdatain(iii,isig)
  dx=wdatain(iii,ilon)-jlon; dx1=1.0-dx
  dy=wdatain(iii,ilat)-jlat; dy1=1.0-dy
  ds=wdatain(iii,isig)-jsig; ds1=1.0-ds
  jlat=min(max(1,jlat),nx); jlon=min(max(0,jlon),nlon)
  latlon11_l=latlon11_obs
  if(jsig==nsig) latlon11_l=0
  if(jstart(m1)==1.and.jlon==nlon) jlon=0
  if(jstart(m1)+jlon1(m1)==nlon+1.and.jlon==0) jlon=nlon
  jlat=jlat-istart(m1)+2; jlon=jlon-jstart(m1)+2
!-----------------------------------gross error test added here
  obserror=1.0/max(sqrt(wdatain(iii,ier)),small)
  obserrlm=max(ermin,min(ermax,obserror))
  residual=sqrt(wdatain(iii,iures)**2+wdatain(iii,ivres)**2)
  ratio=residual/obserrlm
  if(ratio>gross) then
   awork(4,1)=awork(4,1)+weight(iii); wdatain(iii,ier)=0.0
  end if
  awork(4*nsig+jsig+100,1)=awork(4*nsig+jsig+100,1)+wdatain(iii,ier)*&
   wdatain(iii,iures)**2*weight(iii)
  awork(5*nsig+jsig+100,1)=awork(5*nsig+jsig+100,1)+wdatain(iii,ier)*&
   wdatain(iii,ivres)**2*weight(iii)
  awork(6*nsig+jsig+100,1)=awork(6*nsig+jsig+100,1)+weight(iii)
  wdatain(iii,ier)=sqrt(wdatain(iii,ier))
  if(wdatain(iii,ier) > 1.e-10)then
   ict=ict+1
   j1_w(ict)=jlat+(jlon-1)*lat12+(jsig-1)*latlon11_obs
   j2_w(ict)=j1_w(ict)+1
   j3_w(ict)=j1_w(ict)+lat12
   j4_w(ict)=j3_w(ict)+1
   j5_w(ict)=j1_w(ict)+latlon11_l
   j6_w(ict)=j5_w(ict)+1
   j7_w(ict)=j5_w(ict)+lat12
   j8_w(ict)=j7_w(ict)+1
   w1_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy1*ds1              !wgt000
   w2_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy*ds1               !wgt100
   w3_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy1*ds1               !wgt010
   w4_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy*ds1                !wgt110
   w5_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy1*ds               !wgt001
   w6_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy*ds                !wgt101
   w7_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy1*ds                !wgt011
   w8_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy*ds                 !wgt111
   ures(ict)=wdatain(iii,iures)*wdatain(iii,ier)
   vres(ict)=wdatain(iii,ivres)*wdatain(iii,ier)
   uges(ict)=ugesin(iii)*wdatain(iii,ier)
   vges(ict)=vgesin(iii)*wdatain(iii,ier)
   uverr(ict)=wdatain(iii,ier)
   weight_w(ict)=weight(iii)
  end if
  end if
  if(itypx==283.and.iix==2)then
  jlat=wdatain(iii,ilat)
  jlon=wdatain(iii,ilon)
  dx=wdatain(iii,ilon)-jlon; dx1=1.0-dx
  dy=wdatain(iii,ilat)-jlat; dy1=1.0-dy
  if(jstart(m1)==1.and.jlon==nlon) jlon=0
  if(jstart(m1)+jlon1(m1)==nlon+1.and.jlon==0) jlon=nlon
  jlat=jlat-istart(m1)+2; jlon=jlon-jstart(m1)+2
!-----------------------------------gross error test added here
  obserror=1.0/max(sqrt(wdatain(iii,ier)),small)
  obserrlm=max(ermin,min(ermax,obserror))
  residual=abs(wdatain(iii,iures)-wdatain(iii,ivres))
  ratio=residual/obserrlm
  if(ratio>gross) then
   awork(4,1)=awork(4,1)+weight(iii); wdatain(iii,ier)=0.0
  end if
  awork(5,1)=awork(5,1)+wdatain(iii,ier)*&
   (wdatain(iii,iures)-wdatain(iii,ivres))**2*weight(iii)
  awork(6,1)=awork(6,1)+weight(iii)
  wdatain(iii,ier)=sqrt(wdatain(iii,ier))
  if(wdatain(iii,ier) > 1.e-10)then
   icts=icts+1
   j1_spd(icts)=jlat+(jlon-1)*lat12
   j2_spd(icts)=j1_spd(icts)+1
   j3_spd(icts)=j1_spd(icts)+lat12
   j4_spd(icts)=j3_spd(icts)+1
   w1_spd(icts)=fact(iii)*wdatain(iii,ier)*dx1*dy1              !wgt000
   w2_spd(icts)=fact(iii)*wdatain(iii,ier)*dx1*dy               !wgt100
   w3_spd(icts)=fact(iii)*wdatain(iii,ier)*dx*dy1               !wgt010
   w4_spd(icts)=fact(iii)*wdatain(iii,ier)*dx*dy                !wgt110
   spdres(icts)=wdatain(iii,iures)*wdatain(iii,ier)
   usges(icts)=ugesin(iii)*wdatain(iii,ier)
   vsges(icts)=vgesin(iii)*wdatain(iii,ier)
   spderr(icts)=wdatain(iii,ier)
   weight_spd(icts)=weight(iii)
  end if
  end if
 end do
 end do
 
 nwdata_num=ict
 if(nwdata_num /= nwdata_dim)print *,' wind dim, num ',nwdata_dim,nwdata_num
 nspddat_num=icts
 if(nspddat_num /= nspddat_dim)print *,' spd dim, num ',nspddat_dim,nspddat_num

 ioff=200
 scale=1.0
 do k = 1,nlev
    do i=1,nwdta
       if (wpres(i)>=ptop(k).and.wpres(i)<=pbot(k)) then
          val1=wdatain(i,ier)*wdatain(i,iures)
          val2=wdatain(i,ier)*wdatain(i,ivres)
          if (nint(wdatain(i,9))/=283) then
             ress=scale*sqrt(wdatain(i,5)**2+wdatain(i,6)**2)
             val=.5*(val1*val1+val2*val2)
          else
             ress=scale*(wdatain(i,5)-wdatain(i,6))
             val=(val1-val2)**2
          end if
          val=weight(i)*val
          ressw=ress*ress*weight(i)
          itype=nint(wdatain(i,9))-ioff
          itype=max(1,min(itype,100))
          vwork(k,itype,1) = vwork(k,itype,1)+weight(i)
          vwork(k,itype,2) = vwork(k,itype,2)+ressw
          vwork(k,itype,3) = vwork(k,itype,3)+val
          vwork(k,101,1)   = vwork(k,101,1)+weight(i)
          vwork(k,101,2)   = vwork(k,101,2)+ressw
          vwork(k,101,3)   = vwork(k,101,3)+val

       end if
    end do
 end do


 return
 end