CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of ZUVCOV from zuvcov.f
CMAX: ----------------------------------------------------
CFPP$ NOCONCUR R
      SUBROUTINE zuvcov(nrep,igp,ireg,lev,levlindx,obserr,fcerr,irank,
     &    data,nobs,symat,delta,force,gp)
      include '/usr/include/cm/CMF_defs.h'

      
      INTEGER MaxObservPoints
      PARAMETER (MaxObservPoints = 150000)
C
c..................................................................c
C                                                                  !
C NASA Goddard Space Flight Center       Data Assimilation Office  !
C                                                                  !
C Goddard Earth Observing System                                   !
C Data Analysis System - Version 1.3.2                             !
C September, 1993                                                  !
C                                                                  !
C------------------------------------------------------------------c
C  Programmers: S. Bloom, J. Woollen, M. Seablom
C  Language   : Cray cft77 fortran version 4.0
C  Computer   : Cray C90
C
C ZUVCOV forms a covariance matrix for low and middle latitude volumes
C in the upper air z-u-v analysis.
C
C  Passed input arguments:
C
C     NREP     - NUMBER OF OBSERVATIONS TO USE IN GRIDPOINT ANALYSIS
C     IGP      - GRID INDEX
C     IREG     - REGION OF CORRELATIONS TO USE
C     LEV      - LEVEL OF ANALYSIS
C     LEVLINDX - ARRAY WITH DATA INDEXES FOR THIS MATRIX
C     OBSERR   - OBSERVATIONAL ERROR ARRAY
C     FCERR    - FORECAST ERROR ARRAY
C
C  OUTPUT ARGUMENTS:
C
C     IRANK    - Order of matrix ( matrix is returned in
C                common block 'symmetric').
C
C
C  Program modification history:
C
C              J.WOOLLEN  4/84
C              M.NESTLER  7/84
C              J.WOOLLEN  9/85
C              J.WOOLLEN  9/86
C              M.SEABLOM  5/89 TROPICAL DECOUPLING ADDED
C              D.LAMICH  11/89 DIAGNOSTIC PRINTS AT MINI-VOLUMES
C              D.LAMICH  11/89 CORRECTED TROPICAL WIND ANALYSIS
C              D.LAMICH   4/90 USING ZSTD*WSTD FOR ZUCOV & ZVCOV INSTEAD
C                              OF WSTD*WSTD. FIXES GP-OB COV IN TROPICS.
C              S.BLOOM    7/90 DECYBERIZATION; MODULARIZATION OF
C                              COVARIANCE COMPUTATIONS
C              M.Seablom  8/91 Incorporation of variable-sized
C                              minivolumes.
C              M.Seablom  5/92 Implementation of hybrid matrix setup.
C              M.Seablom  8/93 Parallelization of matrix solution.
C
C
C--------------------------------------------------------------------
c
c     Maximal observation points
c
c
c
      INTEGER ipoint (150)
      REAL bc (700), bf (700), zoberr (200), uoberr (200), voberr (200),
     & ustdp (200), vstdp (200), levlindx (nrep), fcerr (im * jnp, mlev)
     &, obserr (200, 4, 19), data (MaxObservPoints * 10)
      COMMON / mvdumper / mv_dumper, num_dump, igpmvp (20), levmvp (20)
      COMMON / zuvmat / zlat (700), zlon (700), zprs (700), zstd (700), 
     &zkx (700), wlat (700), wlon (700), wprs (700), wstd (700), wkx (70
     &0), siny (700), cosy (700), zdel (700), udel (700), vdel (700), zl
     &st (700), wlst (700), type (700), ober (700), fcor (700), gdec (20
     &0), ufac (700), vfac (700), twstd (700), glat (700), glon (700), g
     &prs (700), gstd (700), cv1 (50000), cv2 (50000), cv3 (23569), cv4 
     &(23569), cv5 (23569), cv6 (23569), w1 (50000), w2 (50000), w3 (500
     &00), w4 (50000), w5 (50000), w6 (50000), dtemp (3000)
      COMMON / DATAPTRS / IKX, ILAT, ILON, IPRS, ITYP, IDEL, ISTD
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / GETIJ / GETI (20000), GETJ (20000)
      COMMON / GRIDXX / im, jnp, mlev
      COMMON / OBSCORR / VCORLIST (200), HCORLIST (200), VCOROE (13, 13,
     & 4), HCOROE (400, 5)
      COMMON / PMH / PMAND (18)
      COMMON / sparam / lsolve, ns, nr, nv, nh (13104), nq, kq (3), volc
     &entre_x (13104), ngv, volcentre_y (13104), glatup (13104), glonup 
     &(13104)
      REAL symat (12000), delta (150), force (150, 600, 3), gp (600)
      COMMON / VCORS / VCORFE (13, 13, 20, 2)
      LOGICAL FIRST, onestage, mv_dumper, dump_flag, first_obscorrs
!hpf$ distribute bc(block)
!hpf$ distribute bf(block)
!hpf$ distribute cosy(block)
!hpf$ distribute cv1(block)
!hpf$ distribute cv2(block)
!hpf$ distribute cv3(block)
!hpf$ distribute cv4(block)
!hpf$ distribute cv5(block)
!hpf$ distribute cv6(block)
!hpf$ distribute data(block)
!hpf$ distribute delta(block)
!hpf$ distribute dtemp(block)
!hpf$ distribute fcerr(block,block)
!hpf$ distribute fcor(block)
!hpf$ distribute force(block,block,block)
!hpf$ distribute gdec(block)
!hpf$ distribute geti(block)
!hpf$ distribute getj(block)
!hpf$ distribute glat(block)
!hpf$ distribute glatup(block)
!hpf$ distribute glon(block)
!hpf$ distribute glonup(block)
!hpf$ distribute gp(block)
!hpf$ distribute gprs(block)
!hpf$ distribute gstd(block)
!hpf$ distribute hcorlist(block)
!hpf$ distribute hcoroe(block,block)
!hpf$ distribute igpmvp(block)
!hpf$ distribute ipoint(block)
!hpf$ distribute kq(block)
!hpf$ distribute levlindx(block)
!hpf$ distribute levmvp(block)
!hpf$ distribute nh(block)
!hpf$ distribute ober(block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute pmand(block)
!hpf$ distribute siny(block)
!hpf$ distribute symat(block)
!hpf$ distribute twstd(block)
!hpf$ distribute type(block)
!hpf$ distribute udel(block)
!hpf$ distribute ufac(block)
!hpf$ distribute uoberr(block)
!hpf$ distribute ustdp(block)
!hpf$ distribute vcorfe(block,block,block,block)
!hpf$ distribute vcorlist(block)
!hpf$ distribute vcoroe(block,block,block)
!hpf$ distribute vdel(block)
!hpf$ distribute vfac(block)
!hpf$ distribute voberr(block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
!hpf$ distribute vstdp(block)
!hpf$ distribute w1(block)
!hpf$ distribute w2(block)
!hpf$ distribute w3(block)
!hpf$ distribute w4(block)
!hpf$ distribute w5(block)
!hpf$ distribute w6(block)
!hpf$ distribute wkx(block)
!hpf$ distribute wlat(block)
!hpf$ distribute wlon(block)
!hpf$ distribute wlst(block)
!hpf$ distribute wprs(block)
!hpf$ distribute wstd(block)
!hpf$ distribute zdel(block)
!hpf$ distribute zkx(block)
!hpf$ distribute zlat(block)
!hpf$ distribute zlon(block)
!hpf$ distribute zlst(block)
!hpf$ distribute zoberr(block)
!hpf$ distribute zprs(block)
!hpf$ distribute zstd(block)
      INTEGER  IOBSKX100(700)
!hpf$ distribute iobskx100(block)
      INTEGER  IOBSKX101(700)
!hpf$ distribute iobskx101(block)
      DATA FIRST / .true. /, first_obscorrs / .true. /
c     save first, gradei, omega2, pi180, dlat, dlon,
c    &     maxrep, max_nhgp, first_obscorrs
      CHARACTER * 1 QI, QJ
      CHARACTER * 60 message
c      DATA MAXDAT/150/
c      DATA ZERO/0.0/
c      DATA ONE/1.0/
c      DATA TWO/2.0/
c      DATA ONEPT1/1.10/
c      DATA NINETY/90.0/
c      DATA ONE80/180.0/
c      DATA TENTH/0.10/
c      DATA PT05/.050/
c      DATA TENM4/.0001/
      PARAMETER (maxdat = 150, zero = 0.0, one = 1.0, two = 2.0, onept1 
     &= 1.10)
      PARAMETER (ninety = 90.0, one80 = 180.0, tenth = 0.1, pt05 = 0.05)
      PARAMETER (tenm4 = 0.0001)
      LOGICAL QMVPRINT
      DATA QMVPRINT / .false. /
      CHARACTER * 8 OBSTAG (200)
C
C Function statements.
C Branch to alternate subroutine, "zuvcovh", for high
C latitude volumes.
      INTEGER Ti1, Ti2, Ti3
      i1(i) = i * (i - 1) / 2
      i2(i) = i * (i + 1) / 2
      jx(xlat,dlat) = (xlat + 90.) / dlat + 1.1
      ix(xlon,dlon) = (xlon + 180.) / dlon + 1.1
      IF (nh(igp) .GT. 20) THEN
        CALL zuvcovh(nrep,igp,ireg,lev,levlindx,obserr,fcerr,irank,data,
     &      nobs,symat,delta,force,gp)
        RETURN
      ENDIF
C Variable initialization.
      message(1:60) = ' '
      CALL getconh('GRADEI  ',gradei)
      CALL getconh('OMEGA2  ',omega2)
      CALL getconh('PI180   ',pi180)
      CALL getconh('DLAT    ',dlat)
      CALL getconh('DLON    ',dlon)
      CALL getconi('MAXREP  ',maxrep)
      CALL getconi('MAXNHGP ',max_nhgp)
      IF (maxrep * 3 .GT. maxdat) THEN
        message = ' ZUVCOV: Value of MAXREP exceeds MAXDAT. '
        CALL flash(5,message,1)
      ENDIF
C  Set "nhgp" to be the number of horizontal grid points for current
C  volume.
      nhgp = nh(igp)
      n = 0
C Initialize arrays.
      DO ti3 = 1 , 3
        DO ti2 = 1 , 600
          DO ti1 = 1 , maxdat
            force(ti1,ti2,ti3) = 0.0
          ENDDO
        ENDDO
      ENDDO
      DO ti1 = 1 , maxdat
        delta(ti1) = 0.0
      ENDDO
      DO ti1 = 1 , 12000
        symat(ti1) = -1.
      ENDDO
      irank = 0
      n = 0
      DO ii = 1 , maxrep * 3 + nv * nhgp
        DO jj = 1 , ii
          n = n + 1
          geti(n) = float(ii)
          getj(n) = float(jj)
        ENDDO
      ENDDO
C  Locate grid points within a volume.
      nh_count = 0
      nhnv_count = 0
      DO j = 1 , 2
        DO i = 1 , nhgp / 2
          nh_count = nh_count + 1
          glatx = glatup(igp) + (j - 1) * dlat
          IF (glatx .GT. 90.) THEN
            glatx = 90.
          ENDIF
          glonx = glonup(igp) + (i - 1) * dlon
          iigp = (jx(glatx,dlat) - 1) * im + ix(glonx,dlon)
          gp(nh_count) = iigp
          DO k = 1 , nv
            nhnv_count = nhnv_count + 1
            glat(nhnv_count) = glatx
          ENDDO
          DO k = 1 , nv
            glon(nhnv_count) = glonx
          ENDDO
          DO k = 1 , nv
            nl = lev - 1 + k
            gstd(nhnv_count) = fcerr(iigp,nl)
          ENDDO
        ENDDO
      ENDDO
      DO j = 1 , 2
        DO i = 1 , nhgp / 2
          DO k = 1 , nv
            nhnv_count = nhnv_count + 1
            nl = lev - 1 + k
            gprs(nhnv_count) = pmand(nl)
          ENDDO
        ENDDO
      ENDDO
C Set flag for diagnostic print.
      dump_flag = .FALSE.
      IF (mv_dumper) THEN
        DO imvp = 1 , num_dump
          IF (igp .EQ. igpmvp(imvp) .AND. lev .EQ. levmvp(imvp)) THEN
            dump_flag = .TRUE.
          ENDIF
        ENDDO
        DO imvp = 1 , num_dump
          IF (igp .EQ. igpmvp(imvp) .AND. lev .EQ. levmvp(imvp)) THEN
            WRITE( 6,'(/," ZUVCOV: Grid point information at vmv ",     
     &      i6," at level ",i4,//," ",5x,"Index",7x,"Lat",7x,           
     &"Lon",6x,"Pres",2x,"FG H Err",/)' ) igp,lev
          ENDIF
          DO idxgp = 1 , nv * nhgp
            IF (igp .EQ. igpmvp(imvp) .AND. lev .EQ. levmvp(imvp)) THEN
              WRITE( 6,'(" ",i10,f10.2)' ) idxgp,glat(idxgp),glon(idxgp)
     &            ,gprs(idxgp),gstd(idxgp)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
C  MAKE LISTS OF HEIGHT AND WIND DATA
C  ----------------------------------
      nz = 0
      nw = 0
      CALL ahgather(data(ityp),levlindx,type,nrep)
      icount = 0
      DO ipt = 1 , nrep
        IF (type(ipt) .NE. two) THEN
          icount = icount + 1
          zlst(icount) = levlindx(ipt)
        ENDIF
      ENDDO
      nz = icount
      icount = 0
      DO ipt = 1 , nrep
        IF (type(ipt) .NE. one) THEN
          icount = icount + 1
          wlst(icount) = levlindx(ipt)
        ENDIF
      ENDDO
      nw = icount
      IF (nz .EQ. 0 .AND. nw .EQ. 0) THEN
        irank = 0
      ENDIF
      IF (nz .EQ. 0 .AND. nw .EQ. 0) THEN
        RETURN
      ENDIF
      IF (nz + 2 * nw .GT. maxrep * nq) THEN
        PRINT *, ' ZUVCOV: FATAL>> Too much data. '
        PRINT *, ' NZ:',nz,' NW:',nw,' MAXREP:',maxrep,' NQ:',nq
        PRINT *, ' Printout of TYPE array: '
        PRINT *, type
        message = ' ZUVCOV: Excessive amount of data. See prev. msg.'
        CALL flash(5,message,1)
      ENDIF
C  GATHER STUFF FOR HEIGHT-HEIGHT COVARIANCE
C  Load gridpoint information for minivolumes.
      zlat(1 + nz:nhgp * nv + nz) = glat(:nhgp * nv)
      zlon(1 + nz:nhgp * nv + nz) = glon(:nhgp * nv)
      zprs(1 + nz:nhgp * nv + nz) = gprs(:nhgp * nv)
      zstd(1 + nz:nhgp * nv + nz) = gstd(:nhgp * nv)
      IF (nz .NE. 0) THEN
        CALL ahgather_v1(data(ilat),zlst,zlat,nz)
        CALL ahgather_v1(data(ilon),zlst,zlon,nz)
        CALL ahgather_v1(data(iprs),zlst,zprs,nz)
        CALL ahgather_v1(data(istd),zlst,zstd,nz)
        CALL ahgather_v1(data(idel),zlst,zdel,nz)
        CALL ahgather_v1(data(ikx),zlst,zkx,nz)
C  OBTAIN   HEIGHT-HEIGHT ERROR COVARIANCES
C  ----------------------------------------
C
C     HEIGHT - HEIGHT COVARIANCE WORKING ARRAY
C
C     CV1 = ZZCOV
C
        nz2 = (nz + nhgp * nv) * (nz + nhgp * nv + 1) / 2
        nzs = nz * (nz + 1) / 2
        nzf = nz + nhgp * nv
c        call zeitbeg('mmcov')
        CALL mmcov('ZUV',ireg,nzf,nz2,zlon,zlat,zstd,geti,getj,cv1,w1,w2
     &      ,w3,w4,w5,w6)
c        call zeitend
        CALL ahgather_v2(zprs,geti,w1,nz2)
        CALL ahgather_v2(zprs,getj,w2,nz2)
        CALL vcorrfe_v1('ZUV',w1,1,50000,w2,1,50000,w4,1,50000,nz2)
C  LOAD OBS-OBS COVARIANCES INTO SYMMETRIC MATRIX 'SYMAT'
C  ------------------------------------------------------
C  (NB..HERE OBS-OBS MEANS F.E.'S AT OBS LOCATIONS)
C
C     ALSO - MULTIPLY BY VERTICAL CORRELATION 'W3'
        DO k = 1 , nzs
          symat(k) = cv1(k) * w4(k)
        ENDDO
C  LOAD GRID-OBS COVARIANCES INTO ARRAY 'FORCE'
C  --------------------------------------------
        ios = 0
        DO i = 1 , nhgp * nv
          ios = ios + i - 1
          kos = nzs + (i - 1) * nz + ios
          DO k = 1 , nz
            force(k,i,1) = cv1(k + kos) * w4(k + kos)
          ENDDO
        ENDDO
C
C  GENERATE OBSERVATION HEIGHT-HEIGHT COVARIANCES, AND ADD THEM TO
C  THE COVARIANCE MATRIX
C
c        call zeitbeg('obscov')
        CALL obscov(cv1,zkx,geti,getj,obserr,kq,w1,w2,w3,w4,w5,w6,nz)
c        call zeitend
        DO k = 1 , nzs
          symat(k) = symat(k) + cv1(k)
        ENDDO
      ENDIF
C  GATHER STUFF FOR COVARIANCES WITH WIND TERMS
C  --------------------------------------------
      IF (nw .NE. 0) THEN
        CALL ahgather_v1(data(ilat),wlst,wlat,nw)
        CALL ahgather_v1(data(ilon),wlst,wlon,nw)
        CALL ahgather_v1(data(iprs),wlst,wprs,nw)
        CALL ahgather_v1(data(istd),wlst,wstd,nw)
        CALL ahgather_v1(data(idel + nobs),wlst,udel,nw)
        CALL ahgather_v1(data(idel + 2 * nobs),wlst,vdel,nw)
        CALL ahgather_v1(data(ikx),wlst,wkx,nw)
      ENDIF
      wlat(1 + nw:nhgp * nv + nw) = glat(:nhgp * nv)
      wlon(1 + nw:nhgp * nv + nw) = glon(:nhgp * nv)
      wprs(1 + nw:nhgp * nv + nw) = gprs(:nhgp * nv)
      wstd(1 + nw:nhgp * nv + nw) = gstd(:nhgp * nv)
      nwf = nw + nhgp * nv
      siny(:nwf) = sin(wlat(:nwf) * pi180)
      WHERE (siny(:nwf) .GT. -.05 .AND. siny(:nwf) .LT. 0.0)
        siny(:nwf) = -0.05
      ENDWHERE
      WHERE (siny(:nwf) .LT. .05 .AND. siny(:nwf) .GE. 0.0)
        siny(:nwf) = 0.05
      ENDWHERE
      cosy(:nwf) = cos(wlat(:nwf) * pi180)
      WHERE (cosy(:nwf) .EQ. 0.0)
        cosy(:nwf) = 1.E-4
      ENDWHERE
C
C
C.... OBATAIN TROPICALLY DECOUPLED WIND ERROR VARIANCES (WITH
C     GEOSTROPHIC BALANCE ASSUMED IN EXTRA-TROPICS)
C
c     call zeitbeg('wvarz')
      CALL wvarz(ireg,nwf,siny,cosy,wlat,wstd,wprs,bc,bf,scale)
c     call zeitend
C  COMPUTE THE WIND-HEIGHT COVARIANCES
C  -----------------------------------
C
C     WORKING ARRAYS
C
C     CV1 = UPCOV
C     CV2 = VPCOV
C
      IF (nz .NE. 0) THEN
        nwz = nwf * nz
c        call zeitbeg('wmcov')
        CALL wmcov('ZUV',ireg,nz,nwf,zlon,zlat,zstd,wlon,wlat,wstd,bc,bf
     &      ,scale,cv1,cv2,w1,w2,w3,w4,w5,w6)
c        call zeitend
C
C.... STORE THE RESULTS IN APPROPRIATE LOCATIONS IN SYMAT AND FORCE
C
c        IF(NW.NE.0) THEN
        DO k = 1 , nwf
          ios = (k - 1) * nz
          DO l = 1 , nz
            w1(l + ios) = float(k)
          ENDDO
          DO l = 1 , nz
            w3(l + ios) = float(l)
          ENDDO
        ENDDO
        CALL ahgather_v2(wprs,w1,w2,nwz)
        CALL ahgather_v2(zprs,w3,w4,nwz)
        CALL vcorrfe_v1('ZUV',w2,1,50000,w4,1,50000,w5,1,50000,nwz)
        DO kw = 1 , nw
          iup = i1(nz + kw)
          ivp = i1(nz + nw + kw)
          ios = (kw - 1) * nz
          DO kp = 1 , nz
            symat(iup + kp) = cv1(ios + kp) * w5(ios + kp)
          ENDDO
          DO kp = 1 , nz
            symat(ivp + kp) = cv2(ios + kp) * w5(ios + kp)
          ENDDO
        ENDDO
c        END IF
        nos = nw * nz
        DO kg = 1 , nhgp * nv
          ios = nos + (kg - 1) * nz
          DO kp = 1 , nz
            force(kp,kg,2) = cv1(ios + kp) * w5(ios + kp)
          ENDDO
          DO kp = 1 , nz
            force(kp,kg,3) = cv2(ios + kp) * w5(ios + kp)
          ENDDO
        ENDDO
      ENDIF
C  COMPUTE THE WIND-WIND COVARIANCES
C  ---------------------------------
C
C     WORKING ARRAYS
C
C     CV1 = ZUCOV
C     CV2 = ZVCOV
C
C     CV3 = UUCOV
C     CV4 = UVCOV
C     CV5 = VUCOV
C     CV6 = VVCOV
C
      IF (nw .NE. 0) THEN
        nhv = nhgp * nv
        nw2 = nw * nwf
c     call zeitbeg('wwcov')
        CALL wwcov('ZUV',ireg,nw,nhv,nz,wlon,wlat,wstd,zstd,bc,bf,cv1,
     &      cv2,cv3,cv4,cv5,cv6,w1,w2,w3,w4,w5,w6)
c      call zeitend
        DO k = 1 , nwf
          ios = (k - 1) * nw
          DO l = 1 , nw
            w1(l + ios) = float(k)
          ENDDO
          DO l = 1 , nw
            w3(l + ios) = float(l)
          ENDDO
        ENDDO
        CALL ahgather_v2(wprs,w1,w2,nw2)
        CALL ahgather_v2(wprs,w3,w4,nw2)
        CALL vcorrfe_v1('ZUV',w2,1,50000,w4,1,50000,w5,1,50000,nw2)
        DO kw = 1 , nw
          iuu = i1(nz + kw) + nz
          ivu = i1(nz + nw + kw) + nz
          ivv = i1(nz + nw + kw) + nz + nw
          ios = (kw - 1) * nw
          DO kww = 1 , kw
            symat(iuu + kww) = cv3(ios + kww) * w5(ios + kww)
          ENDDO
          DO kww = 1 , kw
            symat(ivv + kww) = cv6(ios + kww) * w5(ios + kww)
          ENDDO
          DO kww = 1 , nw
            symat(ivu + kww) = cv5(ios + kww) * w5(ios + kww)
          ENDDO
        ENDDO
        nos = nw * nw
        DO kg = 1 , nhgp * nv
          ios = nos + (kg - 1) * nw
          DO kw = 1 , nw
            force(nz + kw,kg,1) = cv1(ios + kw) * w5(ios + kw)
          ENDDO
          DO kw = 1 , nw
            force(nz + nw + kw,kg,1) = cv2(ios + kw) * w5(ios + kw)
          ENDDO
          DO kw = 1 , nw
            force(nz + kw,kg,2) = cv3(ios + kw) * w5(ios + kw)
          ENDDO
          DO kw = 1 , nw
            force(nz + nw + kw,kg,2) = cv4(ios + kw) * w5(ios + kw)
          ENDDO
          DO kw = 1 , nw
            force(nz + kw,kg,3) = cv5(ios + kw) * w5(ios + kw)
          ENDDO
          DO kw = 1 , nw
            force(nz + nw + kw,kg,3) = cv6(ios + kw) * w5(ios + kw)
          ENDDO
        ENDDO
      ENDIF
C  ADD THE OBSERVATIONAL ERROR TO THE FORECAST ERROR
C  -------------------------------------------------
  100 CONTINUE
      IF (nz .NE. 0) THEN
        CALL obertrp0_v1(zkx,1,700,zprs,1,700,ober,1,700,obserr,kq(1),nz
     &      )
        zoberr(:nz) = ober(:nz)
        DO i = 1 , nz
          symat(i2(i)) = symat(i2(i)) + ober(i) * ober(i)
        ENDDO
      ENDIF
      IF (nw .NE. 0) THEN
        CALL obertrp0_v1(wkx,1,700,wprs,1,700,ober,1,700,obserr,kq(2),nw
     &      )
        uoberr(:nw) = ober(:nw)
        DO i = 1 , nw
          symat(i2(nz + i)) = symat(i2(nz + i)) + ober(i) * ober(i)
        ENDDO
c        CALL obertrp0(INDX,WPRS,OBER,OBSERR,KQ(3),NW)
        CALL obertrp0_v1(wkx,1,700,wprs,1,700,ober,1,700,obserr,kq(3),nw
     &      )
        voberr(:nw) = ober(:nw)
        DO i = 1 , nw
          symat(i2(nz + nw + i)) = symat(i2(nz + nw + i)) + ober(i) 
     &        * ober(i)
        ENDDO
      ENDIF
C  COUNT THE RANK AND STORE THE DELTAS
C  -----------------------------------
      irank = nz + nw + nw
      FORALL (i = 1:nz) delta(i) = zdel(i)
      FORALL (i = 1:nw) delta(nz + i) = udel(i)
      FORALL (i = 1:nw) delta(nz + nw + i) = vdel(i)
C DIAGNOSTIC PRINTS
C -----------------
      IF (num_dump .GT. 0 .AND. mv_dumper) THEN
c     if ( igp .eq. 461 ) then
C DIAGNOSTIC PRINTS OF HEIGHT DATA
C --------------------------------
        IF (nz .GT. 0) THEN
          PRINT 6701, igp,lev
6701  FORMAT (' ZUVCOV: HEIGHT DATA INFLUENCING MINI-VOLUME',I6,' AT LEV
     &EL',I4,//,' ',5X,'INDEX',7X,'LAT',7X,'LON',6X,'PRES',2X,'FG H ERR'
     &,3X,'OBS ERR',3X,'H DELTA',8X,'KX',6X,'TYPE',/)
          iobskx100(:nz) = zkx(:nz) + 0.1
          DO idxdata = 1 , nz
            PRINT 6703, idxdata,zlat(idxdata),zlon(idxdata),
     &          zprs(idxdata),zstd(idxdata),zoberr(idxdata),
     &          zdel(idxdata),iobskx100(idxdata),
     &          obstag(iobskx100(idxdata))
          ENDDO
6703  FORMAT (' ',I10,6F10.2,I10,2X,A8)
        ENDIF
C DIAGNOSTIC PRINTS OF WIND DATA
C ------------------------------
        IF (nw .GT. 0) THEN
          PRINT 6805, igp,lev
6805  FORMAT (' ZUVCOV: WIND DATA INFLUENCING MINI-VOLUME',I6,' AT LEVEL
     &',I4,//,' ',5X,'INDEX',7X,'LAT',7X,'LON',6X,'PRES',2X,'FG U ERR',2
     &X,'U OB ERR',3X,'U DELTA',2X,'FG V ERR',2X,'V OB ERR',3X,'V DELTA'
     &,8X,'KX',6X,'TYPE',/)
          iobskx101(:nw) = wkx(:nw) + 0.1
          DO idxdata = 1 , nw
            PRINT 6807, idxdata,wlat(idxdata),wlon(idxdata),
     &          wprs(idxdata),ustdp(idxdata),uoberr(idxdata),
     &          udel(idxdata),vstdp(idxdata),voberr(idxdata),
     &          vdel(idxdata),iobskx101(idxdata),
     &          obstag(iobskx101(idxdata))
          ENDDO
6807  FORMAT (' ',I10,9F10.2,I10,2X,A8)
        ENDIF
C  FOR PRINTING INDIVIDUAL CASES
C  -----------------------------
        PRINT 912, igp
912   FORMAT (' GP=',I5,' L = 9  (500 MB)                ',/,' TEST OF D
     &EVECTORIZED UPPER AIR CODE  ',/,'                                 
     &         ',///,'                         ')
        IF (nz .NE. 0) THEN
          PRINT *, ' VALUES OF I, ZKX, ZLAT, ZLON, ZPRS, ZDEL, ZSTD:'
          DO i = 1 , nz
            PRINT 914, i,zkx(i),zlat(i),zlon(i),zprs(i),zdel(i),zstd(i)
          ENDDO
        ENDIF
        IF (nw .NE. 0) THEN
          WRITE( 6,'(///)' ) 
          PRINT *, ' VALUES OF I,KX,LAT,LON,WPRS,UDEL,VDEL,WSTD,HSD:'
          DO i = 1 , nw
            PRINT 914, i,wkx(i),wlat(i),wlon(i),wprs(i),udel(i),vdel(i),
     &          wstd(i)
          ENDDO
        ENDIF
914   FORMAT (' ',I5,6F12.2,2X,2F7.2)
        ip = nz
        iu = nz + nw
        DO 7078 i = 1 , irank
          IF (i .LE. ip) THEN
            qi = 'Z'
            numi = i
          ENDIF
          IF (i .GT. ip .AND. i .LE. iu) THEN
            qi = 'U'
            numi = i - nz
          ENDIF
          IF (i .GT. iu) THEN
            qi = 'V'
            numi = i - iu
          ENDIF
          DO 7076 j = 1 , i
            IF (j .LE. ip) THEN
              qj = 'Z'
              numj = j
            ENDIF
            IF (j .GT. ip .AND. j .LE. iu) THEN
              qj = 'U'
              numj = j - nz
            ENDIF
            IF (j .GT. iu) THEN
              qj = 'V'
              numj = j - iu
            ENDIF
            ij = i * (i - 1) / 2 + j
            PRINT 911, qi,numi,qj,numj,ij,symat(ij)
911   FORMAT (' ',A2,I2,A5,I2,I10,G18.7,G18.7)
  101       CONTINUE
 7076     CONTINUE
 7078   CONTINUE
        PRINT *, ' PRINTOUT OF FORCING TERMS '
        DO 7080 j = 1 , nhgp * nv
          PRINT *, ' Mini volume:',j
          PRINT 913, (force(i,j,1),i = 1,irank)
  112     CONTINUE
 7080   CONTINUE
        DO 7082 j = 1 , nhgp * nv
          PRINT *, ' Mini volume:',j
          PRINT 913, (force(i,j,2),i = 1,irank)
  113     CONTINUE
 7082   CONTINUE
        DO 7084 j = 1 , nhgp * nv
          PRINT *, ' Mini volume:',j
          PRINT 913, (force(i,j,3),i = 1,irank)
  114     CONTINUE
 7084   CONTINUE
913   FORMAT (10G12.4)
      ENDIF
      RETURN
      END
