C                                                                               
C                                                                               
C      PRED                                                                     
C      TO COMPUTE MEAN GRAVITY ANOMALIES FOR 1/3 DEGREE BLOCKS                  
C                                                                               
C                                                                               
C      VERSION 1.0 OCTOBER, 1973  ..  C. MERRY                                  
C                                                                               
C      FOR A DESCRIPTION OF THE MODELS USED, REFER TO:  MERRY AND               
C      VANICEK (1974)  A METHOD FOR ASTROGRAVIMETRIC GEOID DETERMINATION        
C      TECH. REP. 27, DEPT. OF SURV. ENGG., UNB, FREDERICTON, CANADA.           
C                                                                               
C                                                                               
C      INPUT CONSISTS OF:                                                       
C                                                                               
C      UNIT 5 (CARD READER) :                                                   
C                                                                               
C       1) HEADER CARD - 20A4                                                   
C       2)  LIMITS TO AREA OF INTEREST - LAT1, LAT2, LONG1, LONG2 - 4F8.3       
C                                                                               
C      UNIT 3 (DISC) :                                                          
C                                                                               
C        OBSERVED DATA - LAT, LONG, ANOMALY, NO. OF POINT VALUES                
C          USED, STD. DEV.  (UNFORMATTED)   NOTE :  THIS DATA SET IS            
C          FORMED USING THE PROGRAM MEAN (Q. V.)                                
C                                                                               
C                                                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      DIMENSION A(5800),B(5800),C(5800),D(5800)                                 
      DIMENSION NAME(20)                                                        
 100  FORMAT(4F8.3)                                                             
 101  FORMAT(20A4)                                                              
 200  FORMAT(8X,F6.2,F10.2,25X,F13.2,F14.2)                                     
 201  FORMAT(//,10X,'NO. OF OBSERVED ANOMALIES=',I8,4X,'NO. OF PREDICTED        
     * ANOMALIES=',I8,///)                                                      
 202  FORMAT(8X,F6.2,F10.2,2F13.2)                                              
 203  FORMAT('1',20X,20A4,///)                                                  
 204  FORMAT(10X,'LAT.',5X,'LONG.',5X,'OBSERVED',5X,'STANDARD',5X,'PREDIC       
     CTED',5X,'STANDARD',/,29X,'ANOMALY',5X,'DEVIATION',6X,'ANOMALY',5X,'DE     
     *'DEVIATION',//)                                                           
C                                                                               
C      SET CONSTANTS                                                            
C                                                                               
      IPX=0                                                                     
      IPY=0                                                                     
      RD=1.5D0                                                                  
      RHO=206264.80625D0/3600.D0                                                
      SIZE=1.D0                                                                 
      I=1                                                                       
C                                                                               
C      READ AND WRITE HEADER INFORMATION                                        
C                                                                               
      READ(5,101) NAME                                                          
      WRITE(6,203) NAME                                                         
      WRITE(6,204)                                                              
C                                                                               
C      READ LIMITS                                                              
C                                                                               
      READ(5,100) AX,AY,BX,BY                                                   
C                                                                               
C      READ OBSERVED ANOMALIES                                                  
C                                                                               
      REWIND3                                                                   
  4   READ(3,END=5) A(I),B(I),C(I),II,D(I)                                      
      I=I+1                                                                     
      GOTO4                                                                     
  5   CONTINUE                                                                  
      I=I-1                                                                     
C                                                                               
C      LOOP TO COMPUTE ANOMALIES FOR ALL POINTS WITHIN LIMITS                   
C                                                                               
      AA=DFLOAT(IDINT(AX))-0.5D0*SIZE                                           
  6   CONTINUE                                                                  
      BB=DFLOAT(IDINT(BX))-0.5D0*SIZE                                           
  7   CONTINUE                                                                  
      DO 10 J=1,I                                                               
      DX=DABS(A( J)-AA)                                                         
      DY=DABS( B(J)-BB)                                                         
      IF(DX.LT.0.6.AND.DY.LT.0.6) GO TO8                                        
  10  CONTINUE                                                                  
      DIS= 0.D0                                                                 
C                                                                               
C      OBSERVED DATA WITHIN 0.5, 1.0, 1.5 DEG. OF 'EMPTY' BLOCK                 
C        IS SELECTED                                                            
C                                                                               
  12  DIS=DIS+1.5D0*SIZE                                                        
C                                                                               
C      IF NO DATA WITHIN 1.5 DEG. OF EMPTY BLOCK, NO VALUE IS COMPUTED          
C                                                                               
      IF(DIS.GT.(5.D0*SIZE)) GOTO9                                              
      PNM=0.D0                                                                  
      SNM=0.D0                                                                  
      ANM=0.D0                                                                  
      PQA=0.D0                                                                  
      IK=1                                                                      
      DO 11 K=1,I                                                               
      DX=DABS(A(K)-AA)                                                          
      DY=DABS(B(K)-BB)                                                          
      DS=DSQRT(DX**2+DY**2)                                                     
C                                                                               
C      DISCARD NON-ADJACENT DATA                                                
C                                                                               
      IF(DS.GT.DIS) GOTO13                                                      
C                                                                               
C      COMPUTE WEIGHT FUNCTION                                                  
C                                                                               
      DSS=DSQRT(DX**2+(DY*DCOS(A(K)/RHO))**2)                                   
      PEX=DEXP(DSS/RD)                                                          
C                                                                               
C      FORM SUMMATIONS                                                          
C                                                                               
      SNM=SNM+C(K)/D(K)**2/PEX                                                  
      ANM=ANM+1.D0/D(K)**2/PEX                                                  
      PNM=PNM+C(K)**2/D(K)**2/PEX                                               
      PQA=PQA+D(K)                                                              
      IK=IK+1                                                                   
  13  CONTINUE                                                                  
  11  CONTINUE                                                                  
C                                                                               
C      AT LEAST TWO ADJACENT OBSERVED ANOMALIES USED                            
C                                                                               
      IK=IK-1                                                                   
      IF(IK.LT.2) GOTO12                                                        
C                                                                               
C      COMPUTE PREDICTED MEAN ANOMALY AND STANDARD DEVIATION                    
C                                                                               
      QNM=SNM/ANM                                                               
      VAR=PNM-SNM**2/ANM                                                        
      VAR=VAR/((IK-1)*ANM)                                                      
      STD=DSQRT(VAR)                                                            
      PQA=PQA/IK                                                                
      PQB=DSQRT(PQA**2+VAR)                                                     
C                                                                               
C      WRITE PREDICTED ANOMALIES                                                
C                                                                               
      WRITE(6,200) AA,BB,QNM,PQB                                                
      WRITE(4)  AA,BB,QNM,PQB                                                   
      IPY=IPY+1                                                                 
      GOTO1                                                                     
   8  CONTINUE                                                                  
C                                                                               
C      WRITE OBSERVED ANOMALIES                                                 
C                                                                               
      WRITE(6,202) A(J),B(J),C(J),D(J)                                          
      WRITE(4) A(J),B(J),C(J),D(J)                                              
      IPX=IPX+1                                                                 
      GOTO1                                                                     
  9   CONTINUE                                                                  
      QNM=0.D0                                                                  
      PQB=11.56D0                                                               
      WRITE(6,200) AA,BB,QNM,PQB                                                
  1   CONTINUE                                                                  
      BB=BB+SIZE                                                                
      IF(BB.LE.BY) GOTO7                                                        
      AA=AA+SIZE                                                                
      IF(AA.LE.AY) GOTO6                                                        
C                                                                               
C      WRITE NUMBER OF OBSERVED AND PREDICTED ANOMALIES                         
C                                                                               
      WRITE(6,201) IPX,IPY                                                      
      STOP                                                                      
      END