      subroutine mmcovh ( icor, crscov, radei, pi180, hcorfe,
     &                    ireg, npf, np2, plon, plat, pstd, xx,
     &                    gi, gj, dx, y1, y2, dy, np )
                                                                        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC           
CC                                                         CC           
CC    MMCOV  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            geti  -  ORDERING ARRAY, WITH GETJ, PUTS XX  CC           
CC            getj     INTO SYMMETRIC STORAGE FORM         CC           
CC            IREG  -  GEOGRAPHICAL REGEION (IN TABLE)     CC           
CC            NP2   -  .5*(NP+NH)(NP+NH+1)                 CC           
CC            PLAT  -  LATITUDES OF MASS POINTS            CC           
CC                     (INCLUDES GRID POINTS)              CC           
CC            PLON  -  LONGITUDES OF MASS POINTS           CC           
CC            PSTD  -  STD DEV'S OF ERROR AT MASS POINTS   CC           
CC  GI,GJ,DX,Y1,Y2  -  WORKING SPACE ARRAYS                CC           
CC            XUV   -  'PUV' (PRESSURE) OR 'ZUV' (HEIGHT)  CC           
CC                                                         CC           
CC    OUTPUT ARGUMENT:                                     CC           
CC                                                         CC           
CC            XX    -  MASS-MASS COVARIANCE                CC           
CC                                                         CC           
CC                                                         CC           
CC    ******  FORTRAN 77 VERSION *******                   CC           
CC                                                         CC           
CC                                                         CC           
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC           
                                                                        
                                                                        
                                                                        
      real HCORFE(1000,3)                              
      real crscov, radei, pi180
      real rade
      integer icor
                                                                        
      logical first                                                     
                                                                        
c      DIMENSION PLON(*),PLAT(*),XX(*)                   
c      DIMENSION Y1(*),Y2(*),GI(*),GJ(*),DX(*),DY(*)                     
c      DIMENSION PSTD(*)                                                 
                                                                        
      real,dimension(50000):: Y1,Y2,GI,GJ,DX,DY,XX
      real,dimension(700):: PLON,PLAT,PSTD
         
      DIMENSION COSA(1200)                                               
                                                                        
      DATA   ZERO/0.0/                                                  
      DATA    ONE/1.0/                                                  
      DATA    TWO/2.0/                                                  
      DATA    TEN/10.0/                                                 
      DATA  TENTH/0.1/                                                  
      DATA   PT25/0.25/                                                 
      DATA    PT5/0.5/                                                  
      DATA SIXHUN/600.0/                                                
      DATA  ONE50/150.0/                                                
      DATA  PT001/0.001/                                                
                                                                        
c      CALL GETCONH('RADE    ',RADE)                             
c      CALL GETCONH('PI180   ',PI180)                            
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------------------------------------------------------                 
C     call zeitbeg('mmvsetup')                                          
c      print *, 'inside mmcovh'
c      print *, 'bounds ', ubound(gi,1), ubound(gj,1)
                                                                        
      rade=1./radei
      do k = 1,npf
         ioffset = ( k-1 )*np
         do kk = 1,np
            gi( kk+ioffset ) = float(k)
            gj( kk+ioffset ) = float(kk)
c            print *, kk+ioffset
         end do
      end do
                                                                        
                                                                        
c      print *, 'ready to call ahgather'
      CALL AHGATHER(700,50000,PLON,GI,Y1,NP2)                                     
      CALL AHGATHER(700,50000,PLON,GJ,Y2,NP2)                                     
                                                                        
      DO 15 N=1,NP2                                                     
              DX(N) = COS((Y2(N)-Y1(N))*PI180)                          
 15   CONTINUE                                                          
                                                                        
      DO 20 N=1,NPF                                                     
              COSA(N) = COS( PLAT(N)*PI180 )                            
 20   CONTINUE                                                          
                                                                        
      CALL AHGATHER(700,50000,PLAT,GI,Y1,NP2)                                     
      CALL AHGATHER(700,50000,PLAT,GJ,Y2,NP2)                                     
                                                                        
      DO 25 N=1,NP2                                                     
              DY(N) = ONE                                               
              IF(Y1(N)*Y2(N) .LT. ZERO) DY(N) = -ONE                    
 25   CONTINUE                                                          
                                                                        
      CALL AHGATHER(1200,50000,COSA,GI,Y1,NP2)                                     
      CALL AHGATHER(1200,50000,COSA,GJ,Y2,NP2)                                     
C..............................................                         
C IN WHAT FOLLOWS, THE REMAINDER OF THE OLD                             
C MMCORR CODE IS PUT INTO ONE 'DO' LOOP                                 
C..............................................                         
C     call zeitend                                                      
                                                                        
c      ICOR = 0                                                          
c      IF(XUV.EQ.'PUV') ICOR=1                                           
c      IF(XUV.EQ.'ZUV') ICOR=2                                           
c      IF(XUV.EQ.'MIX') RETURN                                           
c      IF(ICOR.EQ.0) then                                                
c       print*,' FATAL: MMCOV>> Bad XUV encountered. '                   
c       stop 'Fatal.'                                                    
c      end if                                                            

                                                                        
C     call zeitbeg('mmvloop')                                           
                                                                        
      SQ2 = SQRT(TWO)                                                   
      DO 50 K=1,NP2                                                     
c      print *, ' inside do loop ', k
                                                                        
C             COSDX = COS(DX(K)*PI180)                                  
C             COSY1 = COS(Y1(K)*PI180)                                  
C             COSY2 = COS(Y2(K)*PI180)                                  
C             COSDY = COS((Y1(K)-Y2(K))*PI180)                          
                                                                        
C             S =  SQRT(TWO - TWO*DY(K) + TWO*Y1(K)*Y2(K)*(ONE-DX(K)))  
C             S =  SQRT(TWO - TWO*COSDY + TWO*COSY1*COSY2*(ONE-COSDX))  
                                                                        
              ZZ = SQRT((ONE-Y1(K)*Y1(K))*(ONE-Y2(K)*Y2(K)))            
              S =  SQ2 * SQRT( ABS(ONE - Y1(K)*Y2(K)*DX(K) - DY(K)*ZZ)) 
C             S =  SQ2 * SQRT( ONE - COSY1*COSY2*COSDX -                
C    .                  SGN*SQRT((ONE-COSY1*COSY1)*(ONE-COSY2*COSY2)))  
                                                                        
              D = S * RADE                                              
                                                                        
              IF(D.LE.TEN) D=ZERO                                       
                                                                        
              DX(K) = D                                                 
                                                                        
              IF(D.LT.SIXHUN) S=PT5*D                                   
              IF(D.GE.SIXHUN) S=PT25*D + ONE50                          
                                                                        
              INDX = INT(S + PT001)                                     
                                                                        
              R = S - INDX                                              
                                                                        
              IND1 = INDX + 1                                           
              IND2 = IND1 + 1                                           
c             T1 = HCORFE(IND1,1,IREG,ICOR)                            
c             T2 = HCORFE(IND2,1,IREG,ICOR)                             
              T1 = HCORFE(IND1,1)                             
              T2 = HCORFE(IND2,1)                             
c              print *, 'k= ',k,' ind1,2= ',ind1,ind2
                                                                        
              XX(K) = T1 + (T2-T1)*R                                    
                                                                        
 50   CONTINUE                                                          
                                                                        
C     call zeitend                                                      
                                                                        
C                                                                       
C     FORM COVARIANCES BY MULTIPLYING CORRELATIONS BY VARIANCES         
                                                                        
C     call zeitbeg('mmstor')                                            
      CALL AHGATHER(700,50000,PSTD,GI,Y1,NP2)                                     
      CALL AHGATHER(700,50000,PSTD,GJ,Y2,NP2)                                     
                                                                        
      DO 60 K=1,NP2                                                     
                                                                        
        XX(K) = Y1(K) * Y2(K) * XX(K)                                   
                                                                        
 60   CONTINUE                                                          
C     call zeitend                                                      
      RETURN                                                            
      END                                                               
