      subroutine solvedrvl ( mv_begin,mv_end,NumMinVols,
     &                       im,jnp,mlev,lev,dlon,nq,nr,nv,region,
     &                       nreports,nhgp,np,nw,water,crscov,obserr,
     &                       hcorfe, radei, pi180, rhobar, omega2,plon,
     &                       plat, pstd, wlat, wlon, wstd, wdrg, udel,
     &                       vdel,pdel,pkx,wkx,testing_loc1,npd_flag)

c.....................................................................c
!                                                                     :
!  Program name:  SOLVEDRV                                            :
!                                                                     :
!  Programmer: M. Seablom                                             :
!  Date:       March, 1993                                            :
!  Computer:   Cray YMP C90                                           :
!  Language:   Cray cft 5.0 fortran                                   :
!                                                                     :
! Program SOLVEDRV is the driver for the matrix solver. Looping is    :
! performed over all grid volumes. The Choleski algorithm is called   :
! for each matrix, while the looping is spread over the number of     :
! available CPUs.                                                     :
!                                                                     :
c.....................................................................c
!     NASA/DAO OPTIMUM INTERPOLATION ANALYSIS  -  VERSION 1.4         :
c.....................................................................c

      include 'params.h'
      parameter (maxdat=mxreps)

      integer lev,mv_begin,mv_end
      integer  augment
      real radei, pi180, rhobar, omega, crscov

      logical first, final, onestage, qobs, offlev, solved

      real  
     &      clam(144),bc(700),bf(700), 
     &      slam(144), weights(150,600,3),
     &      delta(150), symat(12000), force(150,600,3), 
     &      force2(270000), npd_mark, aug(3), symat_temp(12000)
!hpf$ distribute clam(block)
!hpf$ distribute slam(block) 
!hpf$ distribute bc(block)
!hpf$ distribute bf(block)
!hpf$ distribute weights(block,block,block)
!hpf$ distribute delta(block)
!hpf$ distribute symat(block)
!hpf$ distribute force(block,block,block)
!hpf$ distribute force2(block)
!hpf$ distribute symat_temp(block)

      logical waterl
      real HCORFE(1000,3)
      real OBSERR(200,4,19)
      real, dimension(NumMinVols)::  region
!hpf$ distribute region(block)

      data aug /1.01, 1.1, 10./
 
      integer  im, jnp, mlev
 
      real,dimension(NumMinVols,700):: plat,plon,pstd,
     &     pdel,pkx,wlat,wlon,wstd,wkx,
     &     udel,vdel,wdrg

      real,dimension(700):: platl,plonl,pstdl,
     &     pdell,pkxl,wlatl,wlonl,wstdl,wkxl,
     &     udell,vdell,wdrgl
!hpf$ distribute     platl(block)
!hpf$ distribute plonl(block)
!hpf$ distribute pstdl(block)
!hpf$ distribute     pdell(block)
!hpf$ distribute pkxl(block)
!hpf$ distribute wlatl(block)
!hpf$ distribute     wlonl(block)
!hpf$ distribute wstdl(block)
!hpf$ distribute wkxl(block)
!hpf$ distribute     udell(block)
!hpf$ distribute vdell(block)
!hpf$ distribute wdrgl(block)
!hpf$ distribute     udel(block,block)
!hpf$ distribute     vdel(block,block)
!hpf$ distribute     plat(block,block)
!hpf$ distribute plon(block,block)
!hpf$ distribute     pstd(block,block)
!hpf$ distribute     pdel(block,block)
!hpf$ distribute     pkx(block,block)
!hpf$ distribute wlat(block,block)
!hpf$ distribute     wlon(block,block)
!hpf$ distribute wstd(block,block)
!hpf$ distribute     wkx(block,block)
!hpf$ distribute     wdrg(block,block)

      integer,dimension(NumMinVols)::np,nw,nhgp,nreports
      integer,dimension(NumMinVols)::  testing_loc1
      integer,dimension(NumMinVols,50):: npd_flag
!hpf$ distribute npd_flag(block,block)
!hpf$ distribute testing_loc1(block)
!hpf$ distribute np(block)
!hpf$ distribute nw(block)
!hpf$ distribute nhgp(block)

      logical WATER(NumMinVols)
!hpf$ distribute water(block)


!  Function to calculate diagonal entry index.

       ml(i) = i * (i+1) / 2

       tolerance = 1.e-3
       ns = 1
       print *, 'region =',region

       do i = 1,im
            alam    = dlon * (i-1) * pi180
            clam(i) = cos ( alam )
            slam(i) = sin ( alam )
       end do  

 
c  Set some counters

      matrices = 0
      isizes   = 0
      inreps   = 0
 
!  Loop over grid volumes ( horizontal and vertical ) forming and
!  solving matrices.

      istart=lbound(nhgp,dim=1)
      ifinish=ubound(nhgp,dim=1)

      do 1 igp = istart, ifinish    ! <-- BEGIN looping over volumes.

      print *, igp,mv_begin+igp-1
         ireg = region(igp) + 0.1
      print *, igp,mv_begin+igp-1,region(igp),ireg
         force2(1:270000) = 0.
         isize = 0
 
!  Calculate covariance matrix.

      print *, 'nreports= ',nreports(igp)
         if( nreports(igp).ne.0 ) then
           print *, 'ready to fill local arrays'

           testing_loc1(igp)=nreports(igp)
           wdrg(igp,:) =0.0
           plonl=plon(igp,1:700)
           platl=plat(igp,1:700)
           pstdl=pstd(igp,1:700)
           wlatl=wlat(igp,1:700)
           wlonl=wlon(igp,1:700)
           wstdl=pstd(igp,1:700)
           wdrgl=wdrg(igp,1:700)
           udell=udel(igp,1:700)
           vdell=vdel(igp,1:700)
           pdell=pdel(igp,1:700)
           pkxl=pkx(igp,1:700)
           wkxl=wkx(igp,1:700)
           npl=np(igp)
           nwl=nw(igp)
           nhgpl=nhgp(igp)
           waterl=water(igp)
           print *, 'ready to call puvcov, np,nw= ', npl,nwl
           call PUVCOV(npl,nwl,nhgpl,
     &              nv,IREG,LEV,waterl,crscov,OBSERR,
     &              hcorfe,radei,pi180,rhobar,omega2,
     &              plonl,platl,
     &              pstdl,wlatl,
     &              wlonl,wstdl,
     &              wdrgl,udell,
     &              vdell,pdell,
     &              pkxl, wkxl,
     &              isize, bc, bf, scale,
     &              symat, delta, force)


           matrices = matrices + 1
           print *, 'matrix # ',matrices,' size= ',isize

         endif

         solved = .false.
         itry = 0

         icount = 1
         do i = 1,nq
           do j = 1,nhgp(igp)*nv
             do k = 1,isize
                  force2(icount+k-1) = force(k,j,i)

             end do
              icount = icount + nr
           end do
         end do

         symat_temp(1:12000) = symat(1:12000)

!  Solve the matrices.

        do while ((.not. solved).and. (isize .gt. 0) ) 

         npd_mark = 0.
         weights(1:150,1:600,1:3) = 0.

         print *, 'ready to call vv2'
         call vv2( symat, force2, npd_mark, nhgp(igp), 1, nq, nr,
     &               ml(nr), ns, nv, weights, tolerance )
         print *, 'finished  vv2'
 
         npd_flag(igp,lev)=0
         if ( npd_mark .ne. float(isize) ) then
               print *, 'try augmentation'
               itry = itry + 1
               symat(1:12000) = symat_temp(1:12000)
               do 2 augment = 1,isize
                 idiag = augment*(augment+1)*.5
                 symat(idiag) = symat(idiag)*aug(itry)   
 2             continue
               if ( itry .ge. 3 ) then
                 npd_flag(igp,lev)=1
               print *, ' augmentation failed for ',igp,isize,npd_mark,
     &                    npd_flag(igp,lev)
               solved=.true.
               endif

         else
 
               solved = .true.

         end if
         print *, 'npd_flag ', npd_flag(igp,lev),' igp= ',igp

        end do

  1   continue   ! <--  END looping over grid volumes.
 
 
      return
      end
