      SUBROUTINE WWCOV(icor,crscov,radei,pi180,hcorfe,
     &                 IREG,NW,NH,NP,WLON,WLAT,WSTD,HSTD,BC,BF,     
     &                 XU,XV,UU,UV,VU,VV,GI,GJ,DX,DY,Y1,Y2)             
                                                                        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC           
CC                                                         CC           
CC                                                         CC           
CC    Program Name: WWCOV            Date: April, 1990     CC
CC    Programmers: J. Woollen, S. Bloom                    CC
CC                                                         CC           
CC    WWCOV GENERATES THE WIND-WIND 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            NH    -  NUMBER OF GRIDPTS IN MINIVOLUME     CC           
CC            NP    -  NUMBER MASS OBSERVATION LOCATIONS   CC           
CC            NW    -  NUMBER WIND OBSERVATION LOCATIONS   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            ireg  -  (PRESSURE)=1 OR (HEIGHT)=2 or mix=3 CC           
CC                                                         CC           
CC    OUTPUT ARGUMENT:                                     CC           
CC                                                         CC           
CC            XU    -  U MASS-WIND COVARIANCE              CC           
CC            XV    -  V MASS-WIND COVARIANCE              CC           
CC                                                         CC           
CC            UU    -  U WIND - U WIND COVARIANCE          CC           
CC            UV    -  U WIND - V WIND COVARIANCE          CC           
CC            VU    -  V WIND - U WIND COVARIANCE          CC           
CC            VV    -  V WIND - V WIND COVARIANCE          CC           
CC                                                         CC           
CC                                                         CC           
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC           
                                                                        
                                                                        
      integer icor
      real crscov,radei,pi180
      real HCORFE(1000,3)

      logical first                                                     
      data first /.true./                                               
                                                                        
      real,dimension(700):: WLON,WLAT,WSTD,HSTD,BC,BF
!hpf$ distribute  WLON(block)
!hpf$ distribute WLAT(block)
!hpf$ distribute WSTD(block)
!hpf$ distribute  HSTD(block)
!hpf$ distribute XU(block)
!hpf$ distribute XV(block)             
      real,dimension(50000):: GI,GJ,DX,DY,Y1,Y2,XU,XV                  
!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)
      real,dimension(50000):: UU,UV,VU,VV                                 
!hpf$ distribute  UU(block)
!hpf$ distribute UV(block)
!hpf$ distribute VU(block)
!hpf$ distribute VV(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------------------------------------------------------                 
      IF(FIRST) THEN                                                    
c              CALL GETCONH('RADE    ',RADE)                             
c              CALL GETCONH('PI180   ',PI180)                            
              FIRST = .FALSE.                                           
      END IF                                                            
C------------------------------------------------------                 
                                                                        
      NWF = NW + NH                                                     
                                                                        
      NW2 = NWF * NW                                                    
                                                                        
C                                                                       
C     SET UP COORDINATE PAIRS                                           
C          HERE WE IDENTIFY GI WITH 'IGAT' IN OLD CODE                  
C                           GJ WITH 'JGAT' IN OLD CODE                  
C------------------------------------------------------                 
                                                                        
      DO 10 K=1,NWF                                                     
        IOS = (K-1)*NW                                                  
        DO 8 L=1,NW                                                     
              GI(L+IOS) = float(k)                                      
              GJ(L+IOS) = float(l)                                      
 8      CONTINUE                                                        
 10   CONTINUE                                                          
                                                                        
C                                                                       
C....   GENERATE WIND MASS CORRELATIONS UX,VX                           
C                                                                       
                                                                        
      CALL WWCORR(icor,crscov,radei,pi180,hcorfe,
     &            IREG,NW2,WLON,WLAT,XU,XV,UU,UV,VU,VV,             
     1            GI,GJ,DX,Y1,Y2)                                       
                                                                        
C.........................................................              
C                                                                       
C     BRANCH HERE FOR SURFACE (PUV) OR UPPER-AIR (ZUV)                  
C                                                                       
C.........................................................              
                                                                        
c      IF(XUV.EQ.'PUV') THEN                                             
      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,GJ,DX,NW2)                               
              CALL AHGATHER(700,50000,BF,GJ,DY,NW2)                               
                                                                        
              DO 13 K=1,NW                                              
                                                                        
                IOS = (K-1)*NW                                          
                                                                        
                DO 12 L=1,NW                                            
                                                                        
                  GI(L+IOS) = float(k)                                  
                                                                        
 12             CONTINUE                                                
 13           CONTINUE                                                  
                                                                        
              DO 15 K=1,NH                                              
                                                                        
                IOS = (K-1)*NW + NW*NW                                  
                                                                        
                DO 14 L=1,NW                                            
                                                                        
                  GI(L+IOS) = float(NP+K)                               
                                                                        
 14             CONTINUE                                                
 15           CONTINUE                                                  
                                                                        
              CALL AHGATHER(700,50000,HSTD,GI,Y1,NW2)                             
              CALL AHGATHER(700,50000,WSTD,GJ,Y2,NW2)                             
                                                                        
              DO 20 K=1,NW2                                             
                                                                        
C                                                                       
C...         ***HERE Y1 NOW CONTAINS THE VARIANCE                       
                                                                        
                Y1(K) = Y1(K) * Y2(K)                                   
                                                                        
                PUCOV = -Y1(K)*DX(K)*XU(K) -                            
     1                         Y1(K)*DY(K)*XV(K)                        
                                                                        
                PVCOV =  Y1(K)*DX(K)*XV(K) -                            
     1                         Y1(K)*DY(K)*XU(K)                        
                                                                        
                XU(K) = CRSCOV * PUCOV                                  
                XV(K) = CRSCOV * PVCOV                                  
                                                                        
 20           CONTINUE                                                  
                                                                        
              DO 23 K=1,NWF                                             
                                                                        
                IOS = (K-1)*NW                                          
                                                                        
                DO 22 L=1,NW                                            
                                                                        
                  GI(L+IOS) = float(k)                                  
                                                                        
 22             CONTINUE                                                
 23           CONTINUE                                                  
                                                                        
              CALL AHGATHER(700,50000,WSTD,GI,Y1,NW2)                             
                                                                        
              DO 25 KI=1,NW2                                            
                Y1(KI) = Y1(KI) * Y2(KI)                                
 25           CONTINUE                                                  
                                                                        
              CALL AHGATHER(700,50000,BC,GI,Y2,NW2)                               
              CALL AHGATHER(700,50000,BF,GI,GJ,NW2)                               
                                                                        
              DO 30 K=1,NW2                                             
                                                                        
                BCC   = Y2(K) * DX(K) * Y1(K)                           
                BCF   = Y2(K) * DY(K) * Y1(K)                           
                BFC   = GJ(K) * DX(K) * Y1(K)                           
                BFF   = GJ(K) * DY(K) * Y1(K)                           
                                                                        
                UUT   = BCC*UU(K) + BCF*UV(K) + BFC*VU(K) + BFF*VV(K)   
                UVT   = BCF*UU(K) - BCC*UV(K) + BFF*VU(K) - BFC*VV(K)   
                VUT   = BFC*UU(K) + BFF*UV(K) - BCC*VU(K) - BCF*VV(K)   
                VVT   = BFF*UU(K) - BFC*UV(K) - BCF*VU(K) + BCC*VV(K)   
                                                                        
                UU(K) = UUT                                             
                UV(K) = UVT                                             
                VU(K) = VUT                                             
                VV(K) = VVT                                             
                                                                        
 30           CONTINUE                                                  
                                                                        
      ELSE                                                              
                                                                        
C-------------------------------------------------                      
C   MUST REDEFINE GI TO SELECT PROPERLY THE MASS                        
C   GRID POINTS TO CORRELATE WITH THE WIND OBS PTS                      
C-------------------------------------------------                      
                                                                        
              DO 43 K=1,NW                                              
                                                                        
                IOS = (K-1)*NW                                          
                                                                        
                DO 42 L=1,NW                                            
                                                                        
                  GI(L+IOS) = float(k)                                  
                                                                        
 42             CONTINUE                                                
 43           CONTINUE                                                  
                                                                        
C                                                                       
C ... NOTE, HERE 'NH' IS REALLY NH*NV                                   
C     AND 'NP'IS REALLY NZ (UPPER-AIR)                                  
C                                                                       
              DO 45 K=1,NH                                              
                                                                        
                IOS = (K-1)*NW + NW*NW                                  
                                                                        
                DO 44 L=1,NW                                            
                                                                        
                  GI(L+IOS) = float(NP+K)                               
                                                                        
 44             CONTINUE                                                
 45           CONTINUE                                                  
                                                                        
                                                                        
              CALL AHGATHER(700,50000,HSTD,GI,Y1,NW2)                             
              CALL AHGATHER(700,50000,WSTD,GJ,Y2,NW2)                             
              CALL AHGATHER(700,50000,BC,GJ,DX,NW2)                               
                                                                        
              DO 60 K=1,NW2                                             
                                                                        
                XU(K) = -DX(K) * Y1(K) * Y2(K) * XU(K)                  
                XV(K) =  DX(K) * Y1(K) * Y2(K) * XV(K)                  
                                                                        
 60           CONTINUE                                                  
                                                                        
C-------------------------------------------------                      
C   RETURN GI TO ORIGINAL DEFINITION, WIND                              
C   GRID POINTS TO CORRELATE WITH THE WIND OBS PTS                      
C-------------------------------------------------                      
                                                                        
              DO 53 K=1,NWF                                             
                                                                        
                IOS = (K-1)*NW                                          
                                                                        
                DO 52 L=1,NW                                            
                                                                        
                  GI(L+IOS) = float(k)                                  
                                                                        
 52             CONTINUE                                                
 53           CONTINUE                                                  
                                                                        
              CALL AHGATHER(700,50000,WSTD,GI,Y1,NW2)                             
              CALL AHGATHER(700,50000,BF,GI,DX,NW2)                               
              CALL AHGATHER(700,50000,BF,GJ,DY,NW2)                               
                                                                        
              DO 70 K=1,NW2                                             
                                                                        
                SGN   = DX(K) * DY(K)                                   
                if ( sgn .eq. 0. ) sgn = 1.
                                                                        
                UU(K) =  SGN  *  Y1(K) * Y2(K) * UU(K)                  
                UV(K) = -SGN  *  Y1(K) * Y2(K) * UV(K)                  
                VU(K) = -SGN  *  Y1(K) * Y2(K) * VU(K)                  
                VV(K) =  SGN  *  Y1(K) * Y2(K) * VV(K)                  
                                                                        
 70           CONTINUE                                                  
                                                                        
      END IF                                                            
      RETURN                                                            
      END                                                               
