      SUBROUTINE WMCOV(icor,crscov,radei,pi180,hcorfe,IREG,NP,
     &                 NWF,HLON,HLAT,HSTD,WLON,WLAT,
     &                 WSTD,BC,BF,SCALE,UX,VX,GI,GJ,DX,DY,Y1,Y2)
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC                                                         CC
CC    Optimum Interpolation Analysis Version 1.2           CC
CC                                                         CC
CC    WMCOV GENERATES THE MASS-MASS COVARIANCES            CC
CC    USING CHORD-LENGTH DISTANCE FORMULA:                 CC
CC                                                         CC
CC .5*S**2 = 1-COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)) CC
CC                                                         CC
CC    HORIZONTAL FCST ERROR CORRELATION OF (1) SEA LEVEL   CC
CC    PRESSURE OR (2) GEOPOTENTIAL HEIGHT ARE INTERPOLATED CC
CC    LINEARLY FROM A TABLE (HCORFE) INDEXED BY HORZ. DIST.CC
CC                                                         CC
CC                                                         CC
CC    INPUT ARGUMENTS:                                     CC
CC                                                         CC
CC            GI    -  ORDERING ARRAY, WITH GETJ, PUTS XX  CC
CC            GETJ     INTO SYMMETRIC STORAGE FORM         CC
CC            IREG  -  GEOGRAPHICAL REGEION (IN TABLE)     CC
CC            NP    -  NUMBER OF MASS OBS                  CC
CC            NWF   -  (NW + NH)                           CC
CC            HLAT  -  LATITUDES OF MASS POINTS            CC
CC                     (INCLUDES GRID POINTS)              CC
CC            HLON  -  LONGITUDES OF MASS POINTS           CC
CC            HSTD  -  STD DEV'S OF ERROR AT MASS POINTS   CC
CC           SCALE  -  LENGTH SCALE FROM CORTABLE          CC
CC     DX,DY,Y1,Y2  -  WORKING SPACE ARRAYS                CC
CC            icor  -  (PRESSURE)=1 OR  (HEIGHT)=2         CC
CC                                                         CC
CC    OUTPUT ARGUMENT:                                     CC
CC                                                         CC
CC            UX    -  U WIND-MASS COVARIANCE              CC
CC            VX    -  V WIND-MASS COVARIANCE              CC
CC                                                         CC
CC                                                         CC
CC    ******  FORTRAN 77 VERSION *******                   CC
CC                                                         CC
CC    BASED ON J.WOOLLEN'S PUVCOV AND WMCORR ROUTINES      CC
CC                                                         CC
CC    DECYBERIZED AND REORGANIZED BY S.BLOOM  4/90         CC
CC                                                         CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
 
      integer icor
      logical first
      data first/.true./
      real crscov,radei,pi180
      real hcorfe(1000,3)
 
c      DIMENSION HLON(*),HLAT(*),WLON(*),WLAT(*),UX(*),VX(*)
c      DIMENSION GI(*),GJ(*),DX(*),DY(*),Y1(*),Y2(*),BC(*),BF(*)
c      DIMENSION HSTD(*),WSTD(*)
      real,dimension(700):: HLON,HLAT,WLON,WLAT
      real,dimension(50000):: GI,GJ,DX,DY,Y1,Y2,UX,VX
      real,dimension(700):: HSTD,WSTD,BC,BF
!hpf$ distribute HLON(block)
!hpf$ distribute HLAT(block)
!hpf$ distribute WLON(block)
!hpf$ distribute WLAT(block)
!hpf$ distribute UX(block)
!hpf$ distribute VX(block) 
!hpf$ distribute  GI(block)
!hpf$ distribute GJ(block)
!hpf$ distribute DX(block)
!hpf$ distribute DY(block)
!hpf$ distribute  Y1(block)
!hpf$ distribute Y2(block)
!hpf$ distribute BC(block)
!hpf$ distribute BF(block)
!hpf$ distribute  HSTD(block)
!hpf$ distribute WSTD(block)

      DATA   ZERO/0.0/
      DATA    ONE/1.0/
      DATA    TWO/2.0/
      DATA  TENTH/0.1/
      DATA   PT25/0.25/
      DATA    PT5/0.5/
      DATA  PT001/0.001/
 
C------------------------------------------------------
C------------------------------------------------------
 
C
C     SET UP COORDINATE PAIRS
C          HERE WE IDENTIFY GI WITH 'IGAT' IN OLD CODE
C                           GJ WITH 'JGAT' IN OLD CODE
C------------------------------------------------------
 
      NW2 = NP*NWF
 
      DO 10 K=1,NWF
        IOS = (K-1)*NP
        DO 8 L=1,NP
              GI(L+IOS) = float(k)
              GJ(L+IOS) = float(l)
 8      CONTINUE
 10   CONTINUE
 
 
C
C....   GENERATE WIND MASS CORRELATIONS UX,VX
C
 
      CALL WMCORR(icor,crscov,radei,pi180,hcorfe,
     &            IREG,NW2,HLON,HLAT,WLON,WLAT,UX,VX,
     1            SCALE,GI,GJ,DX,Y1,Y2)
 
 
C.........................................................
C
C     BRANCH HERE FOR SURFACE (1) OR UPPER-AIR (2)
C
C.........................................................
 
      IF(icor.EQ.1) THEN
 
C
C.... FORM COVARIANCES BY MULTIPLYING CORRELATIONS BY VARIANCES
C     (NOTE - WSTD HERE CONTAINS THE TROPICALLY DECOUPLED WIND
C     ERROR VARAIANCES)
C
 
              CALL AHGATHER(700,50000,BC,GI,DX,NW2)
              CALL AHGATHER(700,50000,BF,GI,DY,NW2)
              CALL AHGATHER(700,50000,WSTD,GI,Y1,NW2)
              CALL AHGATHER(700,50000,HSTD,GJ,Y2,NW2)
 
              DO 20 K=1,NW2
 
                UPCOV = -Y1(K)*Y2(K)*DX(K)*UX(K) -
     1                         Y1(K)*Y2(K)*DY(K)*VX(K)
 
 
                VPCOV =  Y1(K)*Y2(K)*DX(K)*VX(K) -
     1                         Y1(K)*Y2(K)*DY(K)*UX(K)
 
                UX(K) = CRSCOV * UPCOV
                VX(K) = CRSCOV * VPCOV
 
 20           CONTINUE
 
 
      ELSE
 
C
C     BC CONTAINS LATITUDE SIGN AND TROPICAL DECOUPLING FACTORS
C
              CALL AHGATHER(700,50000,BC,GI,DX,NW2)
              CALL AHGATHER(700,50000,WSTD,GI,Y1,NW2)
              CALL AHGATHER(700,50000,HSTD,GJ,Y2,NW2)
 
              DO 60 K=1,NW2
 
                VX(K) =  DX(K) * Y1(K) * Y2(K) * VX(K)
                UX(K) = -DX(K) * Y1(K) * Y2(K) * UX(K)
 
 60           CONTINUE
 
      END IF
 
      RETURN
      END
