CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of PUVCOV from puvcov.f
CMAX: ----------------------------------------------------
      SUBROUTINE puvcov(NumMinVols,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 PUVCOV forms a covariance matrix for p-u-v sea level analysis using th
C damped cosine formulas.
C
C The order of the matrix is returned in common SYMETRIC.
C
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
C  Record of modifications:
C              J.WOOLLEN 4/84
C              M.NESTLER 7/84
C              J.WOOLLEN 9/85
C              J.WOOLLEN 9/86
C              M.SEABLOM 6/89    * TROPICAL DECOUPLING ADDED
C              S.Bloom   9/90    * Original Cray version.
C              M.Seablom 8/91    * Upgrade to 1.2; variable minivolumes.
C              M.Seablom 8/93    * Parallelization of matrix solves.
c
c     Maximal observation points
c
c
c
      integer nrep(NumMinVols)
      REAL LEVLINDX (MAXNREP,NumMinVols) 
      real OBSERR (NumMinVols,200, 4, 19), FCERR (*), data (MaxObs
     &ervPoints * 10), symat (12000), delta (150), 
     &force (NumMinVols,150, 600, 3), 
     &gp (600)
      COMMON / DATAPTRS / IKX, ILAT, ILON, IPRS, ITYP, IDEL, ISTD
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / SFCDRAG / DRAG (imr,jmr+1), IWATER (13104)
      COMMON / zuvmat / plat (700), plon (700), pstd (700), pkx (700), p
     &del (700), wlat (700), wlon (700), wdrg (700), wstd (700), wkx (70
     &0), siny (700), cosy (700), plst (700), udel (700), vdel (700), zl
     &st (700), wlst (700), type (700), bc (700), bf (700), gdec (200), 
     &ufac (700), vfac (700), twstd (700), glat (700), glon (700), gprs 
     &(700), gstd (700), cv1 (50000), cv2 (50000), cv3 (23569), cv4 (235
     &69), cv5 (23569), cv6 (23569), w1 (50000), w2 (50000), w3 (50000),
     & w4 (50000), w5 (50000), w6 (50000), dtemp (3000)
      PARAMETER (maxdat = 150, zero = 0.0, one = 1.0, two = 2.0, onept1 
     &= 1.10)
      PARAMETER (ninety = 90.0, one80 = 180.0, tenth = 0.10, pt05 = 0.05
     &)
      PARAMETER (tenm4 = 0.0001)
      COMMON / PUV / CRSCOV
      COMMON / sparam / lsolve, ns, nr, nv, nh (13104), nq, kq (3), volc
     &entre_x (13104), ngv, volcentre_y (13104), glatup (13104), glonup 
     &(13104)
      CHARACTER * 1 QI, QJ
      CHARACTER * 60 message
      DIMENSION IPOINT (150), NUMW (150), NUMP (150), NUM (150)
      LOGICAL WATER, FIRST, UNISLP, onestage
!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 drag(block)
!hpf$ distribute dtemp(block)
!hpf$ distribute fcerr(block)
!hpf$ distribute force(block,block,block)
!hpf$ distribute gdec(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 ipoint(block)
!hpf$ distribute iwater(block)
!hpf$ distribute kq(block)
!hpf$ distribute levlindx(block)
!hpf$ distribute nh(block)
!hpf$ distribute num(block)
!hpf$ distribute nump(block)
!hpf$ distribute numw(block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute pdel(block)
!hpf$ distribute pkx(block)
!hpf$ distribute plat(block)
!hpf$ distribute plon(block)
!hpf$ distribute plst(block)
!hpf$ distribute pstd(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 vdel(block)
!hpf$ distribute vfac(block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(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 wdrg(block)
!hpf$ distribute wkx(block)
!hpf$ distribute wlat(block)
!hpf$ distribute wlon(block)
!hpf$ distribute wlst(block)
!hpf$ distribute wstd(block)
!hpf$ distribute zlst(block)
      INTEGER  KX100(700)
!hpf$ distribute kx100(block)
      INTEGER  KX101(700)
!hpf$ distribute kx101(block)
      DATA FIRST / .true. /
c      save first, unislp, maxrep, im, pi180,
c     &     rhobar, omega2, radei2, radei, rade
C  --- Statement functions -----------------------!
C-------------------------------------------------!
      INTEGER Ti1, Ti2, Ti3
      i1(i) = i * (i - 1) / 2
      i2(i) = i * (i + 1) / 2
      jx(xlat,dlat) = (xlat + ninety) / dlat + onept1
      ix(xlon,dlon) = (xlon + one80) / dlon + onept1
      message(1:60) = ' '
      CALL getconh('RADE    ',rade)
      CALL getconh('RADEI   ',radei)
      CALL getconh('RADEI2  ',radei2)
      CALL getconh('OMEGA2  ',omega2)
      CALL getconh('RHOBAR  ',rhobar)
      CALL getconh('PI180   ',pi180)
      CALL getconi('IM      ',im)
      CALL getconi('MAXREP  ',maxrep)
      CALL getconl('SLP     ',unislp)
      IF (unislp) THEN
        crscov = zero
      ELSE
        crscov = one
      ENDIF
      IF (maxrep * 3 .GT. maxdat) THEN
        message = ' PUVCOV: Value of MAXREP exceeds MAXDAT.'
        CALL flash(5,message,1)
      ENDIF
      CALL getconh('DLAT    ',dlat)
      CALL getconh('DLON    ',dlon)
      force=0.0
      symat = -1.
      delta = 0.
      FORALL (i = 1:maxdat) num(:,i) = i
      irank(:) = 0
C  Locate grid points within a volume.
      glon=0.0
      glat=0.0
      nhlon=nh/2
      nhgp_max=maxval(nh)
      nhgp=nh
      nhmask=.false.
      nhindex=0
      forall(igp=1:NumMinVols,j=1,nhlon(igp)) nhmask(igp,j)=.true.
      forall(igp=1:NumMinVols,j=1,nh(igp))    nhmask2(igp,j)=.true.
      where(nhmask) 
          nhindex_s=1
      elsewhere 
          nhindex_s=0
      endwhere
      call cmf_scan_add(nhindex,nhindex_s,NO_SEGMENT,AXIS_2,CMF_UPWARD,
     &                  CMF_INCLUSIVE,CMF_NONE,nhmask)
      forall(igp=1:NumMinVols,j=1:nhlon(igp))
     &   glon(igp,nhindex(igp,j))=glonup(igp)+(j-1)*dlon
      forall(igp=1:NumMinVols,j=1:nhlon(igp))
     &   glon(igp,nhlon(igp)+nhindex(igp,j))=glonup(igp)+(j-1)*dlon
      forall(igp=1:NumMinVols,j=1:nhlon(igp))
     &   glat(igp,nhindex(igp,j))=glatup(igp)
      forall(igp=1:NumMinVols,j=1:nhlon(igp))
     &   glat(igp,nhlon(igp)+nhindex(igp,j))=glatup(igp)+dlat
      where(glat.gt.90.0.and.nhmask2) glat=90.0
      where(nhmask2)
        iigp=((glat+90.0)/dlat + onept1 - 1)*im +
     &        ((glon+180.0)/dlon + onept1)
      elsewhere
        iigp=0
      endwhere
      forall(igp=1:NumMinVols,j=1:nh(igp))
     &    gstd(igp,j) = fcerr(iigp(igp,j))

C  MAKE LISTS OF PRESSURE AND WIND DATA
C  ------------------------------------
      mv_data=0.0
      forall(igp=1:NumMinVols,j=1:nreps(igp,lev))
     &       mv_data(igp,j,:)=data_mm(datalist(j,lev,igp),:)

      where(mv_data(:,:,ityp).ne.two) 
          pw_mask=.TRUE.
      elsewhere
          pw_mask=.FALSE.
      endwhere
      where(pw_mask) 
          pw_indx_s=1
      elsewhere
          pw_indx_s=0
      endwhere
      np(:) = count(MASK=pw_mask,DIM=2)
      np_max=maxval(np)
      call cmf_scan_add(pw_indx,pw_indx_s,NOSEGMENT,AXIS_2,CMF_UPWARD,
     &                  CMF_INCLUSIVE,CMF_NONE,pw_mask)
      forall(igp=1:NumMinVols, j=1:nreps(igp,lev),pw_mask(igp,j))
     &       plst(igp,p_indx(igp,j))=datalist(j,lev,igp)
      forall(igp=1:NumMinVols, j=1:nreps(igp,lev),pw_mask(igp,j))
     &       nump(igp,pw_indx(igp,j))=j

      where(mv_data(:,:,ityp).ne.one) 
          pw_mask=.TRUE.
      elsewhere
          pw_mask=.FALSE.
      endwhere
      where(pw_mask) 
          pw_indx_s=1
      elsewhere
          pw_indx_s=0
      endwhere
      nw(:) = count(MASK=pw_mask,DIM=2)
      nw_max=maxval(nw)
      call cmf_scan_add(pw_indx,pw_indx_s,NOSEGMENT,AXIS_2,CMF_UPWARD,
     &                  CMF_INCLUSIVE,CMF_NONE,pw_mask)
      forall(igp=1:NumMinVols, j=1:nreps(igp,lev),pw_mask(igp,j))
     &       wlst(igp,pw_indx(igp,j))=datalist(j,lev,igp)
      forall(igp=1:NumMinVols, j=1:nreps(igp,lev),pw_mask(igp,j))
     &       numw(igp,pw_indx(igp,j))=j

C
C     IF THE MINI-VOLUME IS NOT A WATER POINT
C     THEN ONLY PERFORM PRESSURE ANALYSIS
C
C.......................................................
      kxx =((glatup+90.0)/dlat - 1)*im +
     &        (glonup+180.0)/dlon
 
      water = iwater(kxx) 
      where (.NOT.water) nw = 0
      where (np .EQ. 0 .AND. nw .EQ. 0) irank = 0

C  GATHER STUFF FOR PRESSURE-PRESSURE COVARIANCE
C  ---------------------------------------------
      plat=0.0
      plon=0.0
      pstd=0.0
      pdel=0.0
      pkx= 0.0
      where (np .NE. 0) 
              plat=mv_data(:,:,1)
              plon=mv_data(:,:,2)
              pstd=mv_data(:,:,3)
              pdel=mv_data(:,:,4)
              pkx =mv_data(:,:,5)
      endwhere
      plat(:,1 + np_max:nhgp_max + np_max) = glat(:,:nhgp_max)
      plon(:,1 + np_max:nhgp_max + np_max) = glon(:,:nhgp_max)
      pstd(:,1 + np_max:nhgp_max + np_max) = gstd(:,:nhgp_max)

C
C..... LOAD IN THE GRIDPOINT INFO FOR THE MINIVOLUME PTS
C
C  OBTAIN   PRES - PRES ERROR COVARIANCES
C  --------------------------------------
C
C     PRESSURE-PRESSURE COVARIANCE WORKING ARRAY
C
C     CV1 = PPCOV
C
        np2 = (np + nhgp) * np
        nps = np * (np + 1) / 2
        npf = np + nhgp
        npf_max=np_max+nphg_max
        CALL mmcovh('PUV',ireg,np_max,nphg_max,npf_max,plon,
     &              glon,plat,glat,pstd,gstd,cv1)
C  LOAD OBS-OBS COVARIANCES INTO SYMMETRIC MATRIX 'SYMAT'
C  ------------------------------------------------------
C  (NB..HERE OBS-OBS MEANS F.E.'S AT OBS LOCATIONS)
        DO kz = 1 , np
          izz = i1(kz)
          ioffset = (kz - 1) * np
          DO kzz = 1 , kz
            symat(izz + kzz) = cv1(ioffset + kzz)
          ENDDO
        ENDDO
C  LOAD GRID-OBS COVARIANCES INTO ARRAY 'FORCE'
C  --------------------------------------------
        ios = 0
        nos = np * np
        DO kg = 1 , nhgp * nv
          ios = nos + (kg - 1) * np
          FORALL (kz = 1:np) force(kz,kg,1) = cv1(ios + kz)
        ENDDO
      ENDIF
C  DECIDE IF A WIND ANALYSIS IS REQUIRED
C  -------------------------------------
      IF (.NOT.water) THEN
        GOTO 100
      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(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 + nw) = glat(:nhgp)
      wlon(1 + nw:nhgp + nw) = glon(:nhgp)
      wstd(1 + nw:nhgp + nw) = gstd(:nhgp)
      nwf = nw + nhgp
      siny(:nwf) = sin(wlat(:nwf) * pi180)
      WHERE (siny(:nwf) .GT. -pt05 .AND. siny(:nwf) .LT. zero)
        siny(:nwf) = -pt05
      ENDWHERE
      WHERE (siny(:nwf) .LT. pt05 .AND. siny(:nwf) .GE. zero)
        siny(:nwf) = pt05
      ENDWHERE
      cosy(:nwf) = cos(wlat(:nwf) * pi180)
      WHERE (cosy(:nwf) .EQ. zero)
        cosy(:nwf) = tenm4
      ENDWHERE
C
C
C.... OBATAIN TROPICALLY DECOUPLED WIND ERROR VARIANCES (WITH
C     EKMAN BALANCE ASSUMPTION IN EXTRA-TROPICS)
C
C     IN ADDITION, COMPUTE THE EKMAN BALANCE COEFFICENTS BC & BF
C
      CALL wvarp_v1(ireg,nwf,siny,1,700,cosy,1,700,wlat,1,700,wlon,1,700
     &    ,wstd,1,700,wdrg,1,700,bc,1,700,bf,1,700,scale)
C  COMPUTE THE WIND-PRESSURE COVARIANCES
C  -------------------------------------
C
C     WORKING ARRAYS
C
C     CV1 = UPCOV
C     CV2 = VPCOV
C
      IF (np .NE. 0) THEN
        nwp = nwf * np
ccmic$ guard
        CALL wmcov('PUV',ireg,np,nwf,plon,plat,pstd,wlon,wlat,wstd,bc,bf
     &      ,scale,cv1,cv2,w1,w2,w3,w4,w5,w6)
ccmic$ endguard
C
C.... STORE THE RESULTS IN APPROPRIATE LOCATIONS IN SYMAT AND FORCE
C
        IF (nw .NE. 0) THEN
          DO kw = 1 , nw
            iup = i1(np + kw)
            ivp = i1(np + nw + kw)
            ios = (kw - 1) * np
            DO kp = 1 , np
              symat(iup + kp) = cv1(ios + kp)
            ENDDO
            DO kp = 1 , np
              symat(ivp + kp) = cv2(ios + kp)
            ENDDO
          ENDDO
        ENDIF
        nos = nw * np
        DO kg = 1 , nhgp
          ios = nos + (kg - 1) * np
          FORALL (kp = 1:np) force(kp,kg,2) = cv1(ios + kp)
          FORALL (kp = 1:np) force(kp,kg,3) = cv2(ios + kp)
        ENDDO
      ENDIF
C  COMPUTE THE WIND-WIND COVARIANCES
C  ---------------------------------
C
C     WORKING ARRAYS
C
C     CV1 = PUCOV
C     CV2 = PVCOV
C
C     CV3 = UUCOV
C     CV4 = UVCOV
C     CV5 = VUCOV
C     CV6 = VVCOV
C
      IF (nw .NE. 0) THEN
ccmic$ guard
        CALL wwcov('PUV',ireg,nw,nhgp,np,wlon,wlat,wstd,pstd,bc,bf,cv1,
     &      cv2,cv3,cv4,cv5,cv6,w1,w2,w3,w4,w5,w6)
ccmic$ endguard
        DO kw = 1 , nw
          iuu = i1(np + kw) + np
          ivu = i1(np + nw + kw) + np
          ivv = i1(np + nw + kw) + np + nw
          ios = (kw - 1) * nw
          DO kww = 1 , kw
            symat(iuu + kww) = cv3(ios + kww)
          ENDDO
          DO kww = 1 , kw
            symat(ivv + kww) = cv6(ios + kww)
          ENDDO
          DO kww = 1 , nw
            symat(ivu + kww) = cv5(ios + kww)
          ENDDO
        ENDDO
        nos = nw * nw
        DO kg = 1 , nhgp
          ios = nos + (kg - 1) * nw
          FORALL (kw = 1:nw) force(np + kw,kg,1) = cv1(ios + kw)
          FORALL (kw = 1:nw) force(np + nw + kw,kg,1) = cv2(ios + kw)
          FORALL (kw = 1:nw) force(np + kw,kg,2) = cv3(ios + kw)
          FORALL (kw = 1:nw) force(np + nw + kw,kg,2) = cv4(ios + kw)
          FORALL (kw = 1:nw) force(np + kw,kg,3) = cv5(ios + kw)
          FORALL (kw = 1:nw) force(np + nw + kw,kg,3) = cv6(ios + kw)
        ENDDO
      ENDIF
C  ADD THE OBSEVATIONAL ERROR TO THE FORECAST ERROR
C  ------------------------------------------------
  100 CONTINUE
      IF (np .NE. 0) THEN
        kx100(:np) = pkx(:np) + tenth
        DO i = 1 , np
          ipp = i2(i)
          oberr = obserr(kx100(i),2,19)
          symat(ipp) = symat(ipp) + oberr * oberr
        ENDDO
      ENDIF
      IF (nw .NE. 0) THEN
        kx101(:nw) = wkx(:nw) + tenth
        DO i = 1 , nw
          iuu = i2(i + np)
          oberr = obserr(kx101(i),3,19)
          symat(iuu) = symat(iuu) + oberr * oberr
        ENDDO
        DO i = 1 , nw
          ivv = i2(i + np + nw)
          oberr = obserr(kx101(i),4,19)
          symat(ivv) = symat(ivv) + oberr * oberr
        ENDDO
      ENDIF
C  COUNT THE RANK AND STORE THE DELTAS
C  -----------------------------------
      irank = np + nw + nw
      FORALL (k = 1:np) delta(k) = pdel(k)
      FORALL (k = 1:nw) delta(np + k) = udel(k)
      FORALL (k = 1:nw) delta(np + nw + k) = vdel(k)
C  FOR PRINTING INDIVIDUAL CASES
C  -----------------------------
      IF (igp .EQ. 0) THEN
        PRINT 912, igp
912   FORMAT (' PUVCOV: Diagnostic printout for gridpoint ',i5,//)
        IF (np .NE. 0) THEN
          PRINT *, ' Values of I, PKX, PLAT, PLON, PDEL, PSTD:'
          DO i = 1 , np
            PRINT 914, i,pkx(i),plat(i),plon(i),pdel(i),pstd(i)
          ENDDO
        ENDIF
        IF (nw .NE. 0) THEN
          WRITE( 6,'(///)' ) 
          PRINT *, ' Values of I,KX,LAT,LON,UDEL,VDEL,STD, AND DRAG:'
          DO i = 1 , nw
            PRINT 914, i,wkx(i),wlat(i),wlon(i),udel(i),vdel(i),wstd(i),
     &          wdrg(i)
          ENDDO
          i = 1 + nw
          PRINT *, ' Values of BF and BC:',bf(i),bc(i)
        ENDIF
914   FORMAT (' ',I5,6F12.2,2X,E12.6)
C        PRINT*, ' VALUES OF OBSERVATIONAL ERROR --------- '
C        DO 1066 IMIKE = 2,4
C           PRINT*, ' VARIABLE TYPE = ', IMIKE
C           DO 1067 JMIKE = 1,200
C             PRINT*, OBSERR(JMIKE,IMIKE,19)
C1067       CONTINUE
C1066    CONTINUE
        ip = np
        iu = np + nw
        DO 7056 i = 1 , irank
          IF (i .LE. ip) THEN
            qi = 'P'
            numi = i
          ENDIF
          IF (i .GT. ip .AND. i .LE. iu) THEN
            qi = 'U'
            numi = i - np
          ENDIF
          IF (i .GT. iu) THEN
            qi = 'V'
            numi = i - iu
          ENDIF
          DO 7054 j = 1 , i
            IF (j .LE. ip) THEN
              qj = 'P'
              numj = j
            ENDIF
            IF (j .GT. ip .AND. j .LE. iu) THEN
              qj = 'U'
              numj = j - np
            ENDIF
            IF (j .GT. iu) THEN
              qj = 'V'
              numj = j - iu
            ENDIF
            ij = i * (i - 1) / 2 + j
  101       CONTINUE
            PRINT 911, qi,numi,qj,numj,ij,symat(ij)
 7054     CONTINUE
 7056   CONTINUE
911   FORMAT (' ',A2,I2,A5,I2,I10,G15.4)
        PRINT *, ' Printout of forcing terms:'
        PRINT 913, (force(i,1,1),i = 1,irank)
        PRINT 913, (force(i,1,2),i = 1,irank)
        PRINT 913, (force(i,1,3),i = 1,irank)
913   FORMAT (10G12.4)
      ENDIF
      RETURN
      END
