CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of SOLVEDRV from solvedrv.f
CMAX: ----------------------------------------------------
      SUBROUTINE solvedrv(lu,ld,fcerr,cov,aquan,asdae,data,nobs)
      include '/usr/include/cm/CMF_defs.h'
      
      INTEGER MaxObservPoints
      PARAMETER (MaxObservPoints = 150000)
c.....................................................................c
C                                                                     :
C  Program name:  SOLVEDRV                                            :
C                                                                     :
C  Programmer: M. Seablom                                             :
C  Date:       March, 1993                                            :
C  Computer:   Cray YMP C90                                           :
C  Language:   Cray cft 5.0 fortran                                   :
C                                                                     :
C Program SOLVEDRV is the driver for the matrix solver. Looping is    :
C performed over all grid volumes. The Choleski algorithm is called   :
C for each matrix, while the looping is spread over the number of     :
C available CPUs.                                                     :
C                                                                     :
c.....................................................................c
C     NASA/DAO OPTIMUM INTERPOLATION ANALYSIS  -  VERSION 1.4         :
c.....................................................................c
c
c     Maximal observation points
c
c
c
      INTEGER nreps (50), augment
      LOGICAL first, final, onestage, qobs, offlev, solved
      EXTERNAL cov
      common /obs/      maxdex,  nfgge,   lun,
     &                  vhc(13,13,20),    vcorox(18,18,4),
     &                  obserr(200,4,19), ww(13,13,20),
     &                  qobs(200,5),  offlev(200)
      common / obschar / dataname (5), obsname(200)


      REAL aquan (im, jnp, mlev + 1, nq), asdae (im, jnp, mlev + 1), 
     &clam (144), datalist (150, 50),data (MaxObservPoints * 10),
     & fcerr (im * jnp * mlev), slam (144), weights (150, 600, 3), delta (150), s
     &ymat (12000), force (150, 600, 3), gp (600), force2 (270000), npd_
     &mark, aug (3), symat_temp (12000)
      CHARACTER * 8 dataname, obsname
!hpf$ distribute aquan(block,block,block,block)
!hpf$ distribute asdae(block,block,block)
!hpf$ distribute aug(block)
!hpf$ distribute clam(block)
!hpf$ distribute data(block)
!hpf$ distribute datalist(block,block)
!hpf$ distribute delta(block)
!hpf$ distribute fcerr(block)
!hpf$ distribute force(block,block,block)
!hpf$ distribute force2(block)
!hpf$ distribute glatup(block)
!hpf$ distribute glonup(block)
!hpf$ distribute gp(block)
!hpf$ distribute kq(block)
!hpf$ distribute nh(block)
!hpf$ distribute nreps(block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute offlev(block)
!hpf$ distribute qobs(block,block)
!hpf$ distribute reg(block)
!hpf$ distribute slam(block)
!hpf$ distribute symat(block)
!hpf$ distribute symat_temp(block)
!hpf$ distribute vcorox(block,block,block)
!hpf$ distribute vhc(block,block,block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
!hpf$ distribute weights(block,block,block)
!hpf$ distribute ww(block,block,block)
      REAL  ALAM100(144)
!hpf$ distribute alam100(block)
      DATA aug / 1.01, 1.1, 10. /
      COMMON / gridxx / im, jnp, mlev

      COMMON / region / reg (13104)
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / sparam / lsolve, ns, nr, nv, nh (13104), nq, kq (3), volc
     &entre_x (13104), ngv, volcentre_y (13104), glatup (13104), glonup 
     &(13104)
C  Function to calculate diagonal entry index.
      INTEGER Ti1, Ti2, Ti3
      ml(i) = i * (i + 1) / 2
      tolerance = 1.E-3
      ns = 1
      CALL getconh('DLON    ',dlon)
      CALL getconh('PI180   ',pi180)
      FORALL (i = 1:im) alam100(i) = dlon * (i - 1) * pi180
      clam(:im) = cos(alam100(:im))
      slam(:im) = sin(alam100(:im))
c  Set some counters
      matrices = 0
      isizes = 0
      inreps = 0
c
      WRITE( 6,'(" Now beginning matrix solves.....")' ) 
C  Loop over grid volumes ( horizontal and vertical ) forming and
C  solving matrices.
cmic$ do all shared (aquan, asdae, data, fcerr, im, inreps, isizes)
cmic$*       shared (jnp, ld, lu, matrices, MaxObservPoints, mlev, ngv, 
cmic$*       shared (nh, nq, nv, obserr, reg, nobs, cov, ns, nr )
cmic$*       shared (tolerance)
cmic$*       private (datalist, final, ireg, isize, levels)
cmic$*       private ( npd_mark, force2, solved, weights, icount )
cmic$*       private ( augment, idiag, itry, symat_temp )
cmic$*       private (nreps,symat,delta,force,gp) autoscope
      DO 7014 igp = 1 , ngv
        ireg = reg(igp) + 0.1
        DO ti1 = 1 , 270000
          force2(ti1) = 0.
        ENDDO
        CALL search(igp,lu,ld,datalist,nreps,data,nobs)
C  Loop over the levels at a gridpoint.
        DO 7012 levels = lu , ld , nv
          isize = 0
C  Calculate covariance matrix.
          IF (nreps(levels) .NE. 0) THEN
            CALL cov(nreps(levels),igp,ireg,levels,datalist(1,levels),
     &          obserr,fcerr,isize,data,nobs,symat,delta,force,gp)
            matrices = matrices + 1
          ENDIF
          solved = .FALSE.
          itry = 0
          icount = 1
          DO 7006 i = 1 , nq
            DO 7004 j = 1 , nh(igp) * nv
              DO 7002 k = 1 , isize
                force2(icount + k - 1) = force(k,j,i)
 7002         CONTINUE
              icount = icount + nr
 7004       CONTINUE
 7006     CONTINUE
          CALL cmf_fe_array_to_cm(symat_temp,symat)
C  Solve the matrices.
 7010     IF (((.NOT.solved) .AND. (isize .GT. 0))) THEN
            npd_mark = 0.
            DO ti3 = 1 , 3
              DO ti2 = 1 , 600
                DO ti1 = 1 , 150
                  weights(ti1,ti2,ti3) = 0.
                ENDDO
              ENDDO
            ENDDO
            CALL vv2(symat,force2,npd_mark,nh(igp),1,nq,nr,ml(nr),ns,nv,
     &          weights,tolerance)
            IF (npd_mark .NE. float(isize)) THEN
              WRITE( 6,'(" NPD matrix:",i6,"  Augment #",i2)' ) igp,itry
     &            +1
              itry = itry + 1
              CALL cmf_fe_array_from_cm(symat,symat_temp)
CFPP$ SKIP L
              DO 7008 augment = 1 , isize
                idiag = augment * (augment + 1) * .5
                symat(idiag) = symat(idiag) * aug(itry)
    2           CONTINUE
 7008         CONTINUE
              IF (itry .GE. 3) THEN
                STOP
              ENDIF
            ELSE
              solved = .TRUE.
            ENDIF
            GOTO 7010
          ENDIF
 7011     CONTINUE
          IF (isize .GT. 0) THEN
            CALL store_v1(float(isize),float(igp),float(levels),gp,1,600
     &          ,delta,1,150,force2,weights,aquan,asdae,nh(igp))
          ENDIF
 7012   CONTINUE
        CONTINUE
 7014 CONTINUE
      WRITE( 6,'("...matrix solutions complete.")' ) 
      WRITE( 6,'(//)' ) 
C  Average the wind increments at the poles.
      IF (nq .EQ. 3) THEN
        DO 7016 levels = lu , ld
          CALL avguvinc_v1(im,clam,1,144,slam,1,144,aquan(1,jnp,levels,2
     &        ),aquan(1,jnp,levels,3),aquan(1,1,levels,2),aquan(1,1,
     &        levels,3))
 7016   CONTINUE
      ENDIF
      RETURN
      END
