subroutine prewgt(mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    prewgt
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: setup smoothing and grid transform for background error     
!
! program history log:
!   2000-03-15  wu           
!   2004-02-03  kleist, updated to load background stats according
!               to updated mpi distribution on horizontal slabs
!   2004-03-15  derber, kleist, incorporate variances into this routine
!               stats from single file, additional clean up
!   2004-05-17  kleist, documentation and clean up
!   2004-08-03  treadon - add only to module use, add intent in/out
!   2004-10-26  wu - include factors hzscl in the range of RF table
!   2004-11-02  treadon - add horizontal resolution error check on berror_stats
!   2004-11-16  treadon - add longitude dimension to variance array dssv
!   2004-11-20  derber - modify to make horizontal table more reproducable and  
!               move most of table calculations to berror 
!   2005-01-22  parrish - split out balance variables to subroutine prebal--this
!               to make anisotropic filtering option less confusing
!   2005-02-23  wu - setup background variance for qoption=2
!   2005-03-28  wu - change loop index (mlat+1 to mlat) over varq
!   2005-04-14  treadon - add corq2 to global berror_stats read
!   2005-04-22  treadon - change berror file to 4-byte reals
!   2005-05-27  kleist - add setup call for new patch interpolation
!   2005-08-16  guo - add gmao surface interface
!   2005-09-28  guo - set nrr=nlat to support the GMAO grid
!   2005-09-28  guo - fixed nrr to nlat-2 to avoid the subscript of
!		array rlats being out of the range, and to avoid the
!		singularity of rs at rlats(1)=-pi/2.
!   2005-11-16  wgu - set nolp=nr+1+(ny-nlat)/2 to make sure
!               nmix+nrmxb=nr in routine smoothrf no matter what
!               number nlat is.
!   2005-11-29  derber - unify ozone variance calculation
!   2006-01-10  treadon - replace rdsfull with read_gfssfc_full
!   2006-01-11  kleist - place upper/lower bounds on qoption=2 variance
!   2006-01-31  treadon - invert hzscl
!   2006-02-03  derber - fix up sl array stuff
!   2006-04-12  treadon - remove sigl (not used)
!   2006-04-21  kleist  - add capability to perturb background error parameters
!   2006-07-28  derber  - use r1000 from constants
!   2006-09-18  derber  - change minimum moisture variance
!   2007-05-30  h.liu   - use oz stats from berror file
!   2007-07-03  kleist  - add option for flow-dependent background error variances
!   2008-04-23  safford - rm unused uses and vars
!   2008-07-30  guo     - read stats using m_berror_stats
!   2009-01-12  gayno   - rm use of read_gfssfc_full
!   2010-02-25  zhu     - mv varq to m_berror_stats
!                       - make changes for generalizing control variables,
!                         change interface of berror_read_wgt,use nrf*
!   2010-04-01  treadon - move strip to gridmod
!   2010-04-10  parrish - remove rhgues, no longer needed
!   2010-05-28  todling - obtain variable id's on the fly (add getindex)
!   2010-06-01  todling - bypass calculation when pointer not defined
!                         sort of generalize pert_berr
!                         rename as,tsfc_sdv to as3d,as2d,atsfc_sdv (alloc now)
!   2010-06-03  todling - protect motley in dssvs w/ mvars check
!                       - turn nrf2_loc/nrf3_loc into local variables
!   2010-06-18  todling - add call to write_bkgvars_grid and write_bkgvars2_grid
!   2010-07-07  kokron/todling - fix definition of hwllp to do sfc-only
!   2011-07-03  todling - calculation of bl and bl2 must be done in double-prec
!                         or GSI won'd work when running in single precision; go figure!
!   2012-05-14  wargan - add adjustozvar
!   2012-11-26  parrish - move subroutine blend to module blendmod.f90, and add "use blendmod, only: blend"
!   2012-12-15  zhu     - add two cwoption options
!   2013-10-19  todling - all guess variables in met-guess
!   2013-10-25  todling - reposition ltosi and others to commvars
!   2014-02-01  todling - update interface to berror_read_wgt
!   2014-02-05  mkim/todling - move cw overwrite w/ q to m_berror_stats
!   2014-08-02  zhu     - set up new background error variance and correlation lengths of cw 
!                         for all-sky radiance assimilation
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_single,r_double,r_quad
  use berror, only: dssvs,wtaxs,&
       bw,wtxrs,inaxs,inxrs,nr,ny,nx,mr,ndeg,&
       nf,vs,be,dssv,norh,bl2,bl,init_rftable,hzscl,&
       pert_berr,bkgv_flowdep,slw,slw1,slw2,bkgv_write,nhscrf,&
       adjustozvar
  use m_berror_stats,only : berror_read_wgt
  use mpimod, only: nvar_id,levs_id
  use mpimod, only: mpi_comm_world,ierror,mpi_rtype
  use jfunc, only: varcw,cwoption
  use jfunc, only: varq,qoption
  use control_vectors, only: cvars2d,cvars3d
  use control_vectors, only: cvars => nrf_var
  use control_vectors, only: as2d,as3d,atsfc_sdv
  use control_vectors, only: nrf,nc2d,nc3d,mvars
  use gridmod, only: istart,jstart,lat2,lon2,rlats,nlat,nlon,nsig,&
       nnnn1o,lat1,lon1,itotsub,iglobal,ijn,displs_g,&
       strip
  use general_commvars_mod, only: ltosi,ltosj
  use constants, only: zero,quarter,half,one,two,three,&
       rearth_equator,pi,r1000,r400
  use guess_grids, only: isli2
  use guess_grids, only: ntguessig
  use smooth_polcarf, only: norsp,setup_smooth_polcas
  use mpeu_util, only: getindex
  use blendmod, only: blend
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_metguess_mod, only: gsi_metguess_bundle

  implicit none

! Declare passed variables
  integer(i_kind),intent(in   ) :: mype

! Declare local variables
  integer(i_kind) n,nrr,iii,jjj,nxg,i2,im,jm,j2
  integer(i_kind) i,j,k,ii,nn,nbuf,nmix,nxe,nor,ndx,ndy
  integer(i_kind) nlathh,mm1,nolp,mm,ir,k1
  integer(i_kind) ix,jx,mlat
  integer(i_kind) nf2p,istatus
  integer(i_kind),dimension(0:40):: iblend
  integer(i_kind) nrf3_sf,nrf3_q,nrf3_vp,nrf3_t,nrf3_oz,nrf2_ps,nrf2_sst,nrf3_cw
  integer(i_kind),allocatable,dimension(:) :: nrf3_loc,nrf2_loc

  real(r_kind) wlipi,wlipih,df
  real(r_kind) samp,s2u,df2,pi2
  real(r_quad) y,x,dxx
  real(r_kind),dimension(ndeg):: rate
  real(r_kind),dimension(ndeg,ndeg):: turn
  real(r_kind),dimension(lat2,lon2)::temp
  real(r_kind),dimension(nlat,nlon):: sl,factx
  real(r_kind),dimension(-nf:nf,-nf:nf) :: fact1,fact2
  real(r_kind),dimension(mr:nlat-2):: rs
  real(r_kind),dimension(lat1*lon1)::zsm
  real(r_kind),dimension(itotsub)::work1
  real(r_kind),dimension(ny,nx,3):: scsli
  real(r_kind),dimension(-nf:nf,-nf:nf,3):: scs12
  real(r_single),dimension(nlat,nlon):: corsst
  real(r_kind),dimension(lon2,nsig):: dsv
  real(r_single) hsstmin
  real(r_kind) minhsst
  real(r_kind) my_corz
  real(r_kind),allocatable:: randfct(:)
  real(r_kind),allocatable,dimension(:,:,:,:):: sli,sli1,sli2

  real(r_kind),allocatable,dimension(:,:,:):: vz
  real(r_kind),allocatable,dimension(:,:,:):: hwll
  real(r_kind),allocatable,dimension(:,:)  :: hwllp

  real(r_single),allocatable,dimension(:,:,:):: corz
  real(r_single),allocatable,dimension(:,:,:):: hwllin
  real(r_single),allocatable,dimension(:,:,:):: vscalesin
  real(r_single),allocatable,dimension(:,:)  :: corp
  real(r_single),allocatable,dimension(:,:)  :: hwllinp
  real(r_single),allocatable,dimension(:,:)  :: hsst

  real(r_kind),dimension(lat2,lon2,nsig):: sfvar,vpvar,tvar
  real(r_kind),dimension(lat2,lon2):: psvar
  real(r_kind),dimension(:,:,:),pointer :: ges_oz=>NULL()
! real(r_kind),parameter:: eight_tenths = 0.8_r_kind
! real(r_kind),parameter:: six          = 6.0_r_kind
! real(r_kind),parameter:: r800         = 800.0_r_kind
! real(r_kind),parameter:: r40000       = 40000.0_r_kind
! real(r_kind),parameter:: r25          = one/25.0_r_kind

! Initialize local variables
  pi2=two*pi
  ndy=(nlat-ny)/2
  nxe=nlon/8
  nor=norh*2
  mm1=mype+1
  nlathh=nlat/4
  nf2p=2*nf+1

  if(nc2d>0) then
     allocate(hwllp(0:nlat+1,nc2d))
     allocate(corp(nlat,nc2d))
     allocate(hwllinp(nlat,nc2d))
     hwllp=zero
  endif

  if(nc3d>0)then
     allocate(vz(nsig,0:nlat+1,nc3d))
     allocate(hwll(0:nlat+1,nsig,nc3d))
     allocate(corz(nlat,nsig,nc3d))
     allocate(hwllin(nlat,nsig,nc3d))
     allocate(vscalesin(nsig,nlat,nc3d))
  endif

  allocate(hsst(nlat,nlon))

! Get indexes to required CV variables
  nrf3_oz   = getindex(cvars3d,'oz')
  nrf3_t    = getindex(cvars3d,'t')
  nrf3_sf   = getindex(cvars3d,'sf')
  nrf3_vp   = getindex(cvars3d,'vp')
  nrf3_q    = getindex(cvars3d,'q')
  nrf3_cw   = getindex(cvars3d,'cw')
  nrf2_ps   = getindex(cvars2d,'ps')
  nrf2_sst  = getindex(cvars2d,'sst')
! nrf2_stl  = getindex(cvarsmd,'stl')
! nrf2_sti  = getindex(cvarsmd,'sti')

! Setup blending
  mm=4
  call blend(mm,iblend)

  nolp=nr+1+(ny-nlat)/2
!  nbuf=nolp/4
  nbuf=0
  nmix=nolp-nbuf*2
  dxx=1._r_quad/real((nmix+1),r_quad)
  bl2=0._r_kind
  k=0
  do i=1,nmix
     k=k+1
     x=i*dxx
     y=0._r_quad
     y=real(iblend(mm),r_quad)
     do j=mm-1,0,-1
        y=x*y+real(iblend(j),r_quad)
     enddo
     y=y*x**(mm+1)
     bl2(k)=1._r_quad-y
  enddo
  if(minval(bl2)<zero) then
     write(6,*) 'prewgt: trouble bl2 coeffs negative ', bl2
     call stop2(99)
  endif
  do k=1,nmix    
    bl2(k)=sqrt(bl2(k))
  end do
  
! Modify precision to be consistent with bl2 for higher res grids
! ** NOTE ** bl and bl2 are actually defined to be r_kind
! 
  nmix=(nx-nlon)
  dxx=1._r_quad/real((nmix+1),r_quad)
  ndx=(nx-nlon)
  bl=0._r_kind
  k=ndx-nmix
  do i=1,nmix
     k=k+1
     x=real(i,r_quad)*dxx
     y=0._r_quad
     y=real(iblend(mm),r_quad)
     do j=mm-1,0,-1
        y=x*y+real(iblend(j),r_quad)
     enddo
     y=y*x**real((mm+1),r_quad)
     bl(k)=1._r_quad-y
  enddo

  if(minval(bl)<zero) then
     write(6,*) 'prewgt: trouble bl coeffs negative ', bl
     call stop2(99)
  endif
  do k=1,nmix
     bl(k)=sqrt(bl(k))
  end do

! Setup sea-land mask
  sl=zero
  if(bw /= zero)then
    do j=1,lon1*lat1
       zsm(j)=zero
     end do
     do j=1,lon2
        do i=1,lat2
           temp(i,j)=float(isli2(i,j))
        end do
     end do

     call strip(temp,zsm)

     call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
        work1,ijn,displs_g,mpi_rtype,&
        mpi_comm_world,ierror)

     do k=1,iglobal
        i=ltosi(k) ; j=ltosj(k)
        sl(i,j)=work1(k)
     end do



     do j=1,nlon
        do i=1,nlat
           if(sl(i,j) > one)sl(i,j)=zero
        enddo
     enddo
     call smoothww(nlat,nlon,sl,half,2,1)
     do j=1,nlon
        do i=1,nlat
           sl(i,j)=min(max(sl(i,j),zero),one)
        enddo
     enddo
  end if

! Get background error statistics from a file ("berror_stats").
  call berror_read_wgt(corz,corp,hwllin,hwllinp,vscalesin,corsst,hsst,varq,qoption,varcw,cwoption,mype)
  mlat=nlat

! load the horizontal length scales
  hwll=zero
  do j=1,nc3d
     do k=1,nsig
        do i=1,nlat
           hwll(i,k,j)=hwllin(i,k,j)
        end do
     end do
  end do
  if(nrf3_oz>0) hwll(:,:,nrf3_oz)=hwll(:,:,nrf3_oz)*three   !inflate scale

! surface pressure
  if(nrf2_ps>0) then
     do i=1,nlat
        hwllp(i,nrf2_ps)=hwllinp(i,nrf2_ps)
     end do
  endif


! sea surface temperature, convert from km to m
! also calculate a minimum horizontal length scale for
! sst to be used for land skin temp and ice temp
  hsstmin=1.e10_r_single
  minhsst=1.e10_r_kind
  do j=1,nlon
     do i=1,nlat
        hsst(i,j)=r1000*hsst(i,j)
        hsstmin=min(hsstmin,hsst(i,j))
     end do
  end do
  minhsst=hsstmin


! perturb background error
! Things to perturb: as(1-8), hzscl(1-3) and vs(1)
  if (pert_berr) then
     allocate(randfct(nc2d+nc3d+3+1))

     call get_randoms(nc2d+nc3d+3+1,randfct)
     do i=1,nc3d
        as3d(i)=as3d(i)+as3d(i)*randfct(i)
     end do
     ii=nc3d
     do i=1,nc2d
        ii=ii+1
        as2d(i)=as2d(i)+as2d(i)*randfct(ii)
     end do
     do i=1,nhscrf
        hzscl(i)=hzscl(i)+hzscl(i)*randfct(nc2d+nc3d+i)
     end do
     vs=vs+vs*randfct(nc2d+nc3d+3+1)
     if (mype==0) then
        write(6,*) 'PREWGT: REDEFINE AS = ',as3d,as2d
        write(6,*) 'PREWGT: REDEFINE HZSCL = ',hzscl
        write(6,*) 'PREWGT: REDEFINE VS = ',vs
     end if
     deallocate(randfct)
  end if

! As used in the code, the horizontal length scale
! parameters are used in an inverted form.  Invert
! the parameter values here.
  do i=1,nhscrf
     hzscl(i)=one/hzscl(i)
  end do
! apply scaling (deflate/inflate) to vertical length scales
! note: parameter vs needs to be inverted
  vs=one/vs

! Initialize full array to zero before loading part of array below
  vz=zero

! load vertical length scales
  do j=1,nc3d
     do k=1,nsig
        do i=1,nlat
           vz(k,i,j)=vs*vscalesin(k,i,j)
        end do
     end do
  end do

  call rfdpar1(be,rate,ndeg)
  call rfdpar2(be,rate,turn,samp,ndeg)

! Load background error variances onto subdomains
  if(nrf3_sf>0.and.nrf3_vp>0) then
     do k=1,nsig
        do i=1,lat2
           ix=istart(mm1)+i-2
           ix=max(ix,2)
           ix=min(nlat-1,ix)
           do j=1,lon2
              sfvar(i,j,k)=corz(ix,k,nrf3_sf)
              vpvar(i,j,k)=corz(ix,k,nrf3_vp)
           end do
        end do
     end do
  end if
  if(nrf3_t>0) then
     do k=1,nsig
        do i=1,lat2
           ix=istart(mm1)+i-2
           ix=max(ix,2)
           ix=min(nlat-1,ix)
           do j=1,lon2
              tvar(i,j,k)=corz(ix,k,nrf3_t)
           end do
        end do
     end do
  end if

  if(nrf2_ps>0) then
     do i=1,lat2
        ix=istart(mm1)+i-2
        ix=max(ix,2)
        ix=min(nlat-1,ix)
        do j=1,lon2
           psvar(i,j)=corp(ix,nrf2_ps)
        end do
     end do
  end if

! Reweight the variances based on flow dependence if flag set
  if (bkgv_flowdep) then
      call bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype)
  else
      if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,mype)
  endif

! vertical length scales
!!!$omp parallel do  schedule(dynamic,1) private(i,n,k,j,jx,ix,loc,dsv)
  do n=1,nc3d
     do j=1,lat2         
        jx=istart(mm1)+j-2
        jx=max(jx,2)
        jx=min(nlat-1,jx)
        call smoothzo(vz(1,jx,n),samp,rate,n,j,dsv)

!       load variances onto subdomains
        if (n==nrf3_sf) then
           do k=1,nsig
              do i=1,lon2
                 dssv(j,i,k,n)=dsv(i,k)*sfvar(j,i,k)*as3d(n)   ! streamfunction
              end do
           end do
        else if (n==nrf3_vp) then
           do k=1,nsig
              do i=1,lon2
                 dssv(j,i,k,n)=dsv(i,k)*vpvar(j,i,k)*as3d(n)   ! velocity potential
              end do
           end do
        else if (n==nrf3_t) then
           do k=1,nsig
              do i=1,lon2
                 dssv(j,i,k,n)=dsv(i,k)*tvar(j,i,k)*as3d(n)    ! temperature
              end do
           end do
        else if (n==nrf3_oz.and.adjustozvar) then
           call gsi_bundlegetpointer (gsi_metguess_bundle(ntguessig),'oz',ges_oz,istatus)
           if (istatus==0) then
              do k=1,nsig
                 do i=1,lon2
                    my_corz = max(ges_oz(j,i,k),0.000000002_r_kind)
                    ! Reduce weight in the stratosphere
                    if (my_corz .gt. 0.0000001_r_kind) my_corz = my_corz/4.0 
                    dssv(j,i,k,n)=dsv(i,k)*my_corz*as3d(n)   ! ozone
                 end do
              end do
           endif
        else
           do k=1,nsig
              do i=1,lon2
                 dssv(j,i,k,n)=dsv(i,k)*corz(jx,k,n)*as3d(n)
              end do
           end do
        end if
    enddo
  end do

! Special case of dssv for qoption=2 and cw
  if (qoption==2) call compute_qvar3d

!!!$omp parallel do  schedule(dynamic,1) private(i,n,j,jx,ix,loc)
  do n=1,nc2d
     if (n==nrf2_ps) then
        do j=1,lat2         
           do i=1,lon2
              dssvs(j,i,n)=psvar(j,i)*as2d(n)             ! surface pressure
           end do
        end do
     else if (n==nrf2_sst) then
        do j=1,lat2         
           do i=1,lon2
              if(mvars>=2 .and. isli2(j,i)==1)then
                 dssvs(j,i,nc2d+1)= atsfc_sdv(1)          ! land surface temperature
              else if(mvars>=2 .and. isli2(j,i)==2)then
                 dssvs(j,i,nc2d+2)= atsfc_sdv(2)          ! ice surface temperature
              else
                 jx=istart(mm1)+j-2
                 jx=max(jx,2)
                 jx=min(nlat-1,jx)
                 ix=jstart(mm1)+i-2
                 if (ix==0) ix=nlon
                 ix=max(ix,1)
                 if (ix==nlon+1) ix=1
                 ix=min(nlon,ix)
                 dssvs(j,i,n)=corsst(jx,ix)*as2d(n)        ! sea surface temperature
              end if
           end do
        end do
     end if
  end do

  if (bkgv_write) call write_bkgvars2_grid

! distance of gaussian lat from pole on stereogaphic map
! r=r/(1+z)
  do ir=mr,ubound(rs,1)		! ubound(rs,1) is nlat-2 to skip S.P.
     rs(ir)=cos(rlats(nlat-ir))/(one+sin(rlats(nlat-ir)))
  enddo
  df=tan(pi2/nlon*half)

! set up polcas
  call setwts(wtaxs,wtxrs,inaxs,inxrs,rs,df,nor,nxe,nf,mr,nr)

! set up smooth_polcas if desired (norsp > 0)
  if(norsp>0) call setup_smooth_polcas

! load arrays which are in correct units, to be used to
! define the scales below
  do j=1,nlon
     do i=2,nlat-1
        factx(i,j)=one/(one+(one-sl(i,j))*bw)
     end do
     factx(1,j)=factx(2,j)
     factx(nlat,j)=factx(nlat-1,j)
  end do

  wlipi=nlon/pi2
  wlipih=nlon/pi2*half*samp*samp
  do j=1,nx
     jjj=j-ndx
     if(jjj<1)jjj=nlon+jjj
     if(jjj>nlon)jjj=jjj-nlon
     do i=1,ny
        iii=i+ndy
        scsli(i,j,1)=(rlats(iii+1)-rlats(iii-1))*wlipih*cos(rlats(iii))* &
                     factx(iii,jjj)**2
        scsli(i,j,2)=(rlats(iii     )-rlats(iii-1))*wlipi*factx(iii,jjj)
        scsli(i,j,3)=cos(rlats(iii))*factx(iii,jjj)
     enddo
  enddo

  nxg=nxe+norh
  nrr=ubound(rs,1)	! was nf*3/2
  ndx=(nx-nlon)/2
  call polcasl(factx,fact1,fact2,1,nf,mr,nrr,nor,rs,df,nxe,nxg)
  fact1(0,0)=quarter*(fact1(1,0)+fact1(0,1)+fact1(-1,0)+fact1(0,-1))
  fact2(0,0)=quarter*(fact2(1,0)+fact2(0,1)+fact2(-1,0)+fact2(0,-1))

  df2=df*df
  do j=-nf,nf
     jm=j-1
     j2=j*j
     do i=-nf,nf
        im=i-1
        i2=i*i
        scs12(i,j,1)=(samp/(one+(i2+j2)*df2))**2*fact1(i,j)**2
        scs12(i,j,2)=one/(one+((im*im+i2)*half+j2)*df2)*fact1(i,j)
        scs12(i,j,3)=one/(one+(i2+(j2+jm*jm)*half)*df2)*fact1(i,j)
     enddo
  enddo


! Convert horizontal scales from physical units (m) to grid relative units
! rearth_equator is the equatorial radius from a 1999 IAG report.  The
! horizontal scales are defined at the equator, hence the need for the
! equatorial radius.
  s2u=(two*pi*rearth_equator)/float(nlon)


  allocate(sli(ny,nx,2,nnnn1o),sli1(-nf:nf,-nf:nf,2,nnnn1o), &
                               sli2(-nf:nf,-nf:nf,2,nnnn1o))

  allocate(nrf3_loc(nc3d),nrf2_loc(nc2d))
  do ii=1,nc3d
     nrf3_loc(ii)=getindex(cvars,cvars3d(ii))
  enddo
  do ii=1,nc2d
     nrf2_loc(ii)=getindex(cvars,cvars2d(ii))
  enddo

!!!$omp parallel do  schedule(dynamic,1) private(k,k1,j,ii,iii,jjj,i,n,nn,factx,fact1,fact2)
  do k=1,nnnn1o
     k1=levs_id(k)
     if (k1==0) then
        do j=1,nlon
           do i=2,nlat-1
              factx(i,j)=zero
           end do
        end do
     else 
        n=nvar_id(k)
        nn=-1
        do ii=1,nc3d
           if (nrf3_loc(ii)==n) then
              nn=ii
              do j=1,nlon
                 do i=2,nlat-1
                    factx(i,j)=s2u/hwll(i,k1,nn)
                 end do
              end do
              exit
           end if
        end do

        if (nn==-1) then
           do ii=1,nc2d
              if (nrf2_loc(ii)==n .or. n>nrf) then
                 nn=ii
                 if (n>nrf) nn=n-nc3d
                 if (nn==nrf2_sst) then
                    do j=1,nlon
                       do i=2,nlat-1
                          factx(i,j)=s2u/hsst(i,j)
                       end do
                    end do
                 else if (nn>nc2d) then 
                    do j=1,nlon
                       do i=2,nlat-1
                          factx(i,j)=two*s2u/minhsst
                       end do
                    end do
                 else  
                    do j=1,nlon
                       do i=2,nlat-1
                          factx(i,j)=s2u/hwllp(i,nn)
                       end do
                    end do
                 end if
                 exit
              end if
           end do
        end if
     endif    ! end if over nvar_id
     do j=1,nlon
        factx(1,j)=factx(2,j)
        factx(nlat,j)=factx(nlat-1,j)
     end do

     call polcasl(factx,fact1,fact2,1,nf,mr,nrr,nor,rs,df,nxe,nxg)
     fact1(0,0)=quarter*(fact1(1,0)+fact1(0,1)+fact1(-1,0)+fact1(0,-1))
     fact2(0,0)=quarter*(fact2(1,0)+fact2(0,1)+fact2(-1,0)+fact2(0,-1))
! first sli
     do j=1,nx
        jjj=j-ndx
        if(jjj < 1)jjj=jjj+nlon
        if(jjj > nlon)jjj=jjj-nlon
        do i=1,ny
           iii=i+ndy
           slw((j-1)*ny+i,k)=scsli(i,j,1)*factx(iii,jjj)**2
           sli(i,j,1,k)=scsli(i,j,2)*factx(iii,jjj)        
           sli(i,j,2,k)=scsli(i,j,3)*factx(iii,jjj)        
        enddo
     enddo           
! now load sli1/sli2 
     do j=-nf,nf
        do i=-nf,nf
           slw2((j+nf)*nf2p+nf+1+i,k)=scs12(i,j,1)*fact1(i,j)**2
           slw1((j+nf)*nf2p+nf+1+i,k)=scs12(i,j,1)*fact2(i,j)**2
           sli2(i,j,1,k)=scs12(i,j,2)*fact1(i,j)
           sli1(i,j,1,k)=scs12(i,j,2)*fact2(i,j)
           sli2(i,j,2,k)=scs12(i,j,3)*fact1(i,j)
           sli1(i,j,2,k)=scs12(i,j,3)*fact2(i,j)
        enddo
     enddo
  end do ! end do over nsig1o/loadling of sli arrays

  deallocate(nrf3_loc,nrf2_loc)

! Load tables used in recursive filters
  call init_rftable(mype,rate,nnnn1o,sli,sli1,sli2)

! Clean up: first-in, last-out
  deallocate(sli,sli1,sli2)
  deallocate(hsst)

  if(allocated(vscalesin)) deallocate(vscalesin)
  if(allocated(hwllin))    deallocate(hwllin)
  if(allocated(corz))      deallocate(corz)
  if(allocated(hwll))      deallocate(hwll)
  if(allocated(vz))        deallocate(vz)

  if(allocated(hwllinp)) deallocate(hwllinp)
  if(allocated(corp))    deallocate(corp)
  if(allocated(hwllp))   deallocate(hwllp)

  return
end subroutine prewgt

subroutine get_randoms(count,randnums)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_randoms
!   prgmmr: kleist           org: np22              date: 2006-04-24
!
! abstract: get random numbers for perturbing background error parms
!
! program history log:
!   2006-04-21  kleist
!   2008-04-23  safford - rm unused uses
!
!   input argument list:
!     count    - number or random numbers to generate
!
!   output argument list:
!     randnums - array of scaled random numbers
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use obsmod, only: iadate
  use berror, only: pert_berr_fct
  use constants, only: one, two
  implicit none

  integer(i_kind)              ,intent(in   ) :: count
  real(r_kind),dimension(count),intent(  out) :: randnums

  integer(i_kind),allocatable,dimension(:):: numrnds
  real(r_kind),dimension(count+1):: temps
  real(r_kind):: rseed

  integer(i_kind) i,ksize

  call random_seed(size=ksize)
  allocate(numrnds(ksize))

! set seed as a function of analysis date
  rseed = 1e6_r_kind*iadate(1) + 1e4_r_kind*iadate(2) &
       + 1e2_r_kind*iadate(3) + iadate(4)

  do i=1,ksize
     numrnds(i)=rseed
  end do

  call random_seed(put=numrnds)
  deallocate(numrnds)

! this goes from 0-1, but want -1 to 1
  call random_number(temps)

! Set range to be +/- factor
! and don't use first random number generated based on date
  do i=1,count
     randnums(i) = pert_berr_fct*(one - two*temps(i+1))
  end do

  return
end subroutine get_randoms