C                                                                               
C                                                                               
C     PROGRAM TO COMPUTE THE GEOID FROM THE GEM10 COEFFICIENTS                  
C                                                                               
C     INPUT: LAT1,LAT2 = MAXIMUM AND MINIMUM LATITUDE OF THE AREA               
C                        OF INTEREST (DEGREES)                                  
C            LON1,LON2 = MAXIMUM AND MINIMUM LONGITUDE OF THE AREA              
C                        OF INTEREST (DEGREES)                                  
C                  INT = GRID INTERVAL (IN DEGREES) ON WHICH GEOID              
C                        HEIGHTS ARE SOUGHT                                     
C                    N = DEGREE AND ORDER OF SPHERICAL HARMONIC                 
C                        EXPANSION                                              
C                                                                               
C     OUTPUT       UND = MATRIX OF GRID GEOID HEIGHTS                           
C                                                                               
C                                                                               
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C   WORKSPACE ALLOCATION                                                        
C                                                                               
C     PN - IS A MATRIX OF DIMENSIONS (NE+1) X (NE+1) WHERE NE IS                
C          THE MAX DEGREE OF SPHERICAL HARMONIC EXPANSION(NE.GE.N)              
C          (CURRENTLY PN IS DIMENSIONED NROW X NROW, NROW=40)                   
C                                                                               
C     VC, VS -ARE VECTORS OF LENGTH .GE. (N+1)*(N+2)/2                          
C                                                                               
C     UND - IS A MATRIX OF SIZE LATR X LONR WHERE                               
C             LATR = (LAT1-LAT2)/INT   AND                                      
C             LONR = (LON1-LON2)/INT                                            
C           CURRENTLY UND IS DIMENSIONED (ILATR X ILONR) =(20 X 40)             
C                                                                               
C                                                                               
      DIMENSION PN(40,40)                                                       
      DIMENSION UND(20,40)                                                      
      DIMENSION VC(500),VS(500)                                                 
C                                                                               
C                                                                               
C   PARAMETERS                                                                  
C                                                                               
C      GM - GRAVITATIONAL CONSTANT                                              
C                                                                               
C       A - SEMI MAJOR AXIS OF THE ELLIPSOID IMPLIED BY THE                     
C          POTENTIAL COEFFICIENTS                                               
C                                                                               
C     FINV - CORRESPONDING FLATTENING (INVERSE)                                 
C                                                                               
C    GE,GP - GRAVITY AT EQUATOR AND POLES                                       
C                                                                               
C      RJ2 - SECOND ZONAL HARMONIC OF THE GEOPOTENTIAL                          
C                                                                               
C                                                                               
      DATA GM/3.986005D14/ , A/6378140.D0/ , FINV/298.257D0/,                   
     #   GE/9.7803/,GP/9.8322/,DR/0.017453293D0/                                
      DATA RJ2/1082.63D0/                                                       
      ILATR=20                                                                  
      NROW=40                                                                   
C                                                                               
      READ(5,8) LAT1,LAT2,LON1,LON2,INT,N                                       
    8 FORMAT(6I4)                                                               
        NP1=N+1                                                                 
      IDP1=NP1                                                                  
       LATR=(LAT1-LAT2)/INT                                                     
      LONR=(LON1-LON2)/INT                                                      
C                                                                               
      F = 1./FINV                                                               
      E2 = 2.*F - F**2                                                          
      B = A*(1.-F)                                                              
C                                                                               
      ID=IDP1-1                                                                 
      NCOL=ID                                                                   
      LMAX=ID                                                                   
C                                                                               
C                                                                               
C     READ POTENTIAL COEFFICIENTS FROM DISK UNIT 8 (ZONAL POTENTIAL             
C     COEFFICIENTS ARE READ FIRST, FOLLOWED BY THE TESSERAL AND                 
C     SECTORIAL ONES).                                                          
C                                                                               
C                                                                               
      LENGTH=IDP1*(IDP1+1)/2                                                    
      READ(8) IDP1,LENGTH,(VC(I),I=1,IDP1)                                      
      DO 10 I=1,IDP1                                                            
   10 VS(I)=0.                                                                  
      IDP11=IDP1+1                                                              
      READ(8) (VC(I),I=IDP11,LENGTH)                                            
      READ(8) (VS(I),I=IDP11,LENGTH)                                            
       WRITE(6,100) (VC(I),VS(I),I=1,LENGTH)                                    
  100 FORMAT(2F20.5)                                                            
C                                                                               
C                                                                               
C   COMPUTE THE FIRST THREE NORMALIZED EVEN ZONAL HARMONICS OF THE              
C   NORMAL GRAVITATIONAL FIELD                                                  
      XJN2 = -RJ2/DSQRT(5.D0)                                                   
      XJN4 = (E2**2/35.D0)*(10.D0*RJ2/E2-1.D0)                                  
      XJN6 = -(E2**3/21.D0)*(15.D0*RJ2/E2-2.D0)/DSQRT(13.D0)                    
C                                                                               
C                                                                               
C     SUBTRACT THE NORMAL POTENTIAL FROM THE GEOPOTENTIAL                       
C     IN ORDER TO OBTAIN THE DISTURBING POTENTIAL.                              
C                                                                               
      VC(3) = VC(3) - XJN2                                                      
      VC(5) = VC(5) - XJN4                                                      
      VC(7) = VC(7) - XJN6                                                      
C                                                                               
      LATST=1                                                                   
      IF(LAT1.EQ.90) LATST=2                                                    
      DO 20 LAT=LATST,LATR                                                      
      LATD=(LAT1-(LAT-1)*INT)                                                   
      XLAT=(LAT1-(LAT-1)*INT)*DR                                                
      XLATG =DATAN((1.D0-E2)*DTAN(XLAT))                                        
      IF(LAT.EQ.1 .OR. LAT.EQ.LATR) XLATG = XLAT                                
      PHI = XLATG/DR                                                            
      CF = DCOS(PHI*DR)                                                         
      SF = DSIN(PHI*DR)                                                         
C                                                                               
C                                                                               
C     COMPUTE NORMAL GRAVITY FROM SOMIGLIANA'S FORMULA                          
C                                                                               
C                                                                               
      GAMMA=(A*GE*CF**2+B*GP*SF**2)/DSQRT(A**2*CF**2+B**2*SF**2)                
      R = A/DSQRT(1.+E2*SF**2/(1.-E2))                                          
      CALL LEGDRE(1,PHI,PN,NROW,NP1)                                            
      DO 21 LON = 1,LONR                                                        
      LOND = LON2 + (LON-1)*INT                                                 
      XLAMDA = DFLOAT(LON2) + (LON-1)*DFLOAT(INT)                               
      SUM = 0.                                                                  
      DO 16 L=2,LMAX                                                            
      SUM1 = 0.                                                                 
      L1 = L + 1                                                                
      DO 15 MM = 1,L1                                                           
      M = MM - 1                                                                
      M1 = MM                                                                   
      LC = (2*(NCOL+1)-M1)*(M1-1)/2 + L1                                        
      C = VC(LC)                                                                
      S = VS(LC)                                                                
      SUM1=SUM1+(C*DCOS(M*XLAMDA*DR)+S*DSIN(M*XLAMDA*DR))*PN(L1,M1)             
   15 CONTINUE                                                                  
      SUM=SUM+(A/R)**L*SUM1*1.D-6                                               
   16 CONTINUE                                                                  
      GEOIDH = (GM/(GAMMA*R))*SUM                                               
      WRITE(6,1000) LATD,LOND,GEOIDH                                            
 1000 FORMAT(5X,'LAT. = ',I6,'  LONG. = ',I6,'  GEOID HT. = ',F8.4,/)           
      UND(LAT,LON) = GEOIDH                                                     
   21 CONTINUE                                                                  
   20 CONTINUE                                                                  
      CALL MOUTD(UND,ILATR,LAT,LON)                                             
      N1=1                                                                      
      N2=10                                                                     
      IF(N2.GT.LON) N2=LON                                                      
   45 CONTINUE                                                                  
      DO 30 I = 1,LAT                                                           
  30  WRITE(6,200) (UND(I,J),J=N1,N2)                                           
  200 FORMAT(//,10F12.5)                                                        
      WRITE(6,201)                                                              
  201 FORMAT(///)                                                               
      IF(N2.EQ.LON) GO TO 50                                                    
      N1=N2+1                                                                   
      N2=N2+10                                                                  
      IF(N2.GT.LON) N2=LON                                                      
      GO TO 45                                                                  
   50 STOP                                                                      
      END