CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of MIXCOV from mixcov.f
CMAX: ----------------------------------------------------
C......................................................................
C                                                                     .
C  Program:    MIXCOV                                                 .
C  Language:   Cray cft77 fortran version 4.0                         .
C  Computer:   Cray YMP 864                                           .
C                                                                     .
C  Programmer: M. Nestler, D. Lamich, M. Seablom                      .
C  Date:       August, 1991                                           .
C                                                                     .
C  MIXCOV computes moisture covariances.                              .
C                                                                     .
C  * MODIFICATIONS:                                                   .
C                                                                     .
C    FEBRUARY 1991: (M. SEABLOM)  Capability to process > 100k obs.   .
C    AUGUST   1991: (M. SEABLOM)  Variable minivolumes incorporated.  .
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
      SUBROUTINE mixcov(nrep,igp,ireg,lev,rindx,obserr,fcerr,isize,data,
     &    nobs,symat,delta,force,gp)
      include '/usr/include/cm/CMF_defs.h'
      
      INTEGER MaxObservPoints
      PARAMETER (MaxObservPoints = 150000)
C  COMPUTES COVARIANCE MATRIX FOR MIXING RATIO ANALYSIS
c
c     Maximal observation points
c
c
c
      INTEGER MXGRID
      PARAMETER (MXGRID = 700)
      INTEGER MXREPS
      PARAMETER (MXREPS = 150)
      REAL COSDX
      REAL COSDY
      REAL COSY1
      REAL COSY2
      REAL D
      REAL data (MaxObservPoints * 10)
      REAL DLAT
c      SAVE              DLAT
      REAL DLON
c      SAVE              DLON
      DIMENSION FCERR (IM * JNP, MLEV)
      LOGICAL FIRST
!hpf$ distribute data(block)
!hpf$ distribute delta(block)
!hpf$ distribute fcerr(block,block)
!hpf$ distribute force(block,block,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 kq(block)
!hpf$ distribute nh(block)
!hpf$ distribute ober(block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute pmand(block)
!hpf$ distribute prsi(block)
!hpf$ distribute prsj(block)
!hpf$ distribute rindx(block)
!hpf$ distribute symat(block)
!hpf$ distribute vcor(block)
!hpf$ distribute vcorfe(block,block,block,block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
!hpf$ distribute xdel(block)
!hpf$ distribute xkx(block)
!hpf$ distribute xlat(block)
!hpf$ distribute xlon(block)
!hpf$ distribute xprs(block)
!hpf$ distribute xstd(block)
!hpf$ distribute xxcov(block)
      INTEGER  INDX100(mxreps+mxgrid)
!hpf$ distribute indx100(block)
      REAL  COSY1100(mxreps+mxgrid)
!hpf$ distribute cosy1100(block)
      REAL  Y1100(mxreps+mxgrid)
!hpf$ distribute y1100(block)
      REAL  X1100(mxreps+mxgrid)
!hpf$ distribute x1100(block)
      DATA FIRST / .true. /
c      SAVE              FIRST
      REAL GLAT (MXGRID)
      REAL GLATX
      REAL GLON (MXGRID)
      REAL GLONX
      REAL GPRS (MXGRID)
      REAL GSTD (MXGRID)
      REAL HCOR
      INTEGER I
      INTEGER IGP
      INTEGER IIGP
      INTEGER INDX
      INTEGER ISIZE
      INTEGER ix
      INTEGER I4X5
      INTEGER J
      INTEGER jx
      INTEGER J4X5
      INTEGER K
      INTEGER LEV
      REAL OBER (MXREPS)
      REAL PI180
c      SAVE              PI180
      DIMENSION RINDX (NREP)
      INTEGER N
      INTEGER NHROOT
c      SAVE              NHROOT
      INTEGER NL
      INTEGER NNH
      INTEGER NNI
      INTEGER NNJ
      INTEGER NNV
      INTEGER NREG
      INTEGER NX
      DIMENSION OBSERR (200, 4, 19)
      REAL PRSI ((MXREPS + MXGRID) * (MXREPS + MXGRID + 1) / 2)
      REAL PRSJ ((MXREPS + MXGRID) * (MXREPS + MXGRID + 1) / 2)
      REAL RADE
c      SAVE              RADE
      REAL S
      REAL VAR
      REAL VCOR ((MXREPS + MXGRID) * (MXREPS + MXGRID + 1) / 2)
      REAL XDEL (MXREPS)
      REAL XKX (MXREPS)
      REAL XLAT (MXREPS + MXGRID)
      REAL XLON (MXREPS + MXGRID)
      REAL XPRS (MXREPS + MXGRID)
      REAL XSTD (MXREPS + MXGRID)
      REAL XXCOV ((MXREPS + MXGRID) * (MXREPS + MXGRID + 1) / 2)
      REAL X1
      REAL X2
      REAL Y1
      REAL Y2
      COMMON / DATAPTRS / IKX, ILAT, ILON, IPRS, ITYP, IDEL, ISTD
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / GRIDXX / IM, JNP, MLEV
      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)
      INTEGER Ti1, Ti2, Ti3
      ix(xxlon,dlon) = (xxlon + 180.) / dlon + 1.1
      jx(xxlat,dlat) = (xxlat + 90.) / dlat + 1.1
      CALL getconh('DLAT    ',dlat)
      CALL getconh('DLON    ',dlon)
      CALL getconh('RADE    ',rade)
      CALL getconh('PI180   ',pi180)
      DO ti1 = 1 , 150
        delta(ti1) = 0.
      ENDDO
      DO ti3 = 1 , 3
        DO ti2 = 1 , 600
          DO ti1 = 1 , 150
            force(ti1,ti2,ti3) = 0.
          ENDDO
        ENDDO
      ENDDO
      DO ti1 = 1 , 12000
        symat(ti1) = -1.
      ENDDO
C  Set "nhgp" to be the number of horizontal grid points for current vol
      nhgp = nh(igp)
      n = 0
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  GATHER STUFF FOR MAKING THE COVARIANCES
      nx = nrep
      isize = nrep
      FORALL (k = 1:nx) indx100(k) = rindx(k) + 0.1
      FORALL (k = 1:nx) xlat(k) = data(ilat + indx100(k) - 1)
      FORALL (k = 1:nx) xlon(k) = data(ilon + indx100(k) - 1)
      FORALL (k = 1:nx) xprs(k) = data(iprs + indx100(k) - 1)
      FORALL (k = 1:nx) xstd(k) = data(istd + indx100(k) - 1)
      FORALL (k = 1:nx) xdel(k) = data(idel + indx100(k) - 1)
      FORALL (k = 1:nx) xkx(k) = data(ikx + indx100(k) - 1)
      xlat(1 + nx:nh(igp) * nv + nx) = glat(:nh(igp) * nv)
      xlon(1 + nx:nh(igp) * nv + nx) = glon(:nh(igp) * nv)
      xprs(1 + nx:nh(igp) * nv + nx) = gprs(:nh(igp) * nv)
      xstd(1 + nx:nh(igp) * nv + nx) = gstd(:nh(igp) * nv)
C  COMPUTE THE COVARIANCES
      n = 0
      x1100(:nx + nh(igp) * nv) = xlon(:nx + nh(igp) * nv)
      y1100(:nx + nh(igp) * nv) = xlat(:nx + nh(igp) * nv)
      cosy1100(:nx + nh(igp) * nv) = cos(y1100(:nx + nh(igp) * nv) * 
     &    pi180)
      DO i = 1 , nx + nh(igp) * nv
        DO j = 1 , min0(i,nx)
          n = n + 1
          prsi(n) = xprs(i)
          prsj(n) = xprs(j)
          x2 = xlon(j)
          y2 = xlat(j)
          cosy2 = cos(y2 * pi180)
          cosdy = cos((y1100(i) - y2) * pi180)
          cosdx = cos((x1100(i) - x2) * pi180)
          s = 1.0 - cosdy + cosy1100(i) * cosy2 * (1.0 - cosdx)
          s = sqrt(2.0 * s)
          IF (s .LE. 0.002) THEN
            s = 0.0
          ENDIF
          d = s * rade
          hcor = exp(-0.002 * d)
          var = xstd(i) * xstd(j)
          xxcov(n) = var * hcor
        ENDDO
      ENDDO
      CALL vcorrfe_v1('MIX',prsi,1,(mxreps + mxgrid) * (mxreps + mxgrid 
     &    + 1) / 2,prsj,1,(mxreps + mxgrid) * (mxreps + mxgrid + 1) / 2,
     &    vcor,1,(mxreps + mxgrid) * (mxreps + mxgrid + 1) / 2,n)
      CALL obertrp0_v1(xkx,1,mxreps,xprs,1,mxreps + mxgrid,ober,1,mxreps
     &    ,obserr,kq(1),nx)
      n = 0
      DO i = 1 , nx + nh(igp) * nv
        DO j = 1 , min0(i,nx)
          IF (i .LE. nx .AND. i .EQ. j) THEN
            delta(i) = xdel(i)
          ENDIF
        ENDDO
      ENDDO
      DO i = 1 , nx + nh(igp) * nv
        DO j = 1 , min0(i,nx)
          n = n + 1
          xxcov(n) = xxcov(n) * vcor(n)
          IF (i .LE. nx) THEN
            symat(n) = xxcov(n)
            IF (i .EQ. j) THEN
              symat(n) = symat(n) + ober(i) ** 2
            ENDIF
          ENDIF
          IF (.NOT.i .LE. nx) THEN
            force(j,i - nx,1) = xxcov(n)
          ENDIF
        ENDDO
      ENDDO
C  PRINTING INTERESTING CASES
      IF (igp .EQ. 0) THEN
        WRITE( 6,900 ) (i,xkx(i),xlat(i),xlon(i),xprs(i),xdel(i),xstd(i)
     &      ,i = 1,nx)
        WRITE( 6,910 ) igp
        DO 7018 i = 1 , nx
          WRITE( 6,920 ) i,(symat(k),k = i * (i + 1) / 2 - i + 1,i * (i
     &        + 1) / 2)
  400     CONTINUE
 7018   CONTINUE
        WRITE( 6,930 ) igp
        DO 7022 i = 1 , nh(igp)
          DO 7020 j = 1 , nv
            WRITE( 6,940 ) i,j,(force(n,j + (i - 1) * nv,1),n = 1,nx)
  405       CONTINUE
 7020     CONTINUE
          CONTINUE
 7022   CONTINUE
      ENDIF
      RETURN
900   FORMAT (' MIXCOV: ',I5,6F12.2)
910   FORMAT (//' MIXCOV: COV. MATRIX FOR MOISTURE ANALYSIS IGP =',I6)
920   FORMAT ('0MIXCOV:  I = ',I4/(' MIXCOV: ',10F7.3))
930   FORMAT (//' MIXCOV: FORCING TERMS FOR MOISTURE ANALYSIS IGP =',I6)
940   FORMAT ('0MIXCOV: IH =',I2,' IV =',I2/(' MIXCOV: ',10F7.3))
      END
