      subroutine solvedrvl ( 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
      integer  augment
      real radei, pi180, rhobar, omega, crscov

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

c      real  
c     &      clam(144), 
c     &      slam(144), weights(150,600,3),
c     &      delta(150), symat(12000), force(150,600,3), 
c     &      force2(270000), npd_mark, aug(3), symat_temp(12000)
c     &      force2(2700), npd_mark, aug(3), symat_temp(12000)

      real  
     &      clam(144), 
     &      slam(144), weights(150,60,3),
     &      delta(150), symat(1200), force(150,60,3), 
     &      force2(270), npd_mark, aug(3), symat_temp(1200)

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

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

!hpf$ distribute     udel(block,block)
!hpf$ distribute     vdel(block,block)
!hpf$ distribute     plat(block,block),plon(block,block)
!hpf$ distribute     pstd(block,block)
!hpf$ distribute     pdel(block,block)
!hpf$ distribute     pkx(block,block),wlat(block,block)
!hpf$ distribute     wlon(block,block),wstd(block,block)
!hpf$ distribute     wkx(block,block)

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

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


!  Function to calculate diagonal entry index.

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

       tolerance = 1.e-3
       ns = 1

       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=dlbound(nhgp,dim=1)
      ifinish=dubound(nhgp,dim=1)

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

         ireg = region(igp) + 0.1
c         force2(1:270000) = 0.
         force2(1:2700) = 0.
         isize = 0
 
!  Calculate covariance matrix.

         if( nreports(igp).ne.0 ) then

           testing_loc1(igp)=nreports(igp)
           call PUVCOV(np,nw,nhgp,nv,IREG,LEV,water,crscov,OBSERR,
     &              hcorfe,radei,pi180,rhobar,omega2,plon,plat,
     &              pstd,wlat,wlon,wstd,wdrg,udel,vdel,pdel,pkx,
     &              wkx, isize, bc, bf, scale, symat, delta, force)

           matrices = matrices + 1

         endif

         solved = .false.
         itry = 0

         icount = 1
         do i = 1,nq
           do j = 1,nhgp(igp)*nv
             do k = 1,isize
c                  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.

c         call vv2( symat, force2, npd_mark, nhgp(igp), 1, nq, nr,
c     &               ml(nr), ns, nv, weights, tolerance )
         call vv2( symat, force, npd_mark, nhgp(igp), 1, nq, nr,
     &               ml(nr), ns, nv, weights, tolerance )
 
         npd_flag(igp,lev)=0
         if ( npd_mark .ne. float(isize) ) then

               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
               endif

         else
 
               solved = .true.

         end if

        end do

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