C                                                                               
C                                                                               
C*******************************************************************************
C******* G E O S T R O P H I C - W I N D - COMPUTATIONS ************************
C*                       WRITTEN BY: M.M. NASSAR                              **
C*                          SUMMER, 1976.                                     **
C*******************************************************************************
C                                                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      REAL*8 NGRAD                                                              
      DIMENSION F1(700),F2(700),F3(700)                                         
      DIMENSION NGRAD(700),EGRAD(700),IVEC(12)                                  
      DOUBLE PRECISION DSIN,DCOS,DSQRT                                          
      DATA INORTH/'NG'/,IEAST/'EG'/                                             
C                                                                               
C READ-IN THE TOTAL NUMBER OF YEARS FOR WHICH DATA RECORDS ARE AVAILABLE        
C AND ALSO READ-IN THE STARTING YEAR OF AVAILABLE RECORDS.                      
C                                                                               
      READ(5,2) NY,IYSTAR                                                       
  2   FORMAT(I4,I4)                                                             
C                                                                               
C COMPUTE THE SIZE OF THE INVOLVED VECTORS = TOTAL NUMBER OF MONTHLY MEANS      
C                                                                               
      N=NY*12                                                                   
C                                                                               
C SET-UP SOME CONVERSION CONSTANTS.                                             
C                                                                               
      PI=3.141592653589793D 0                                                   
      RHODEG=180.D 0/PI                                                         
      RHOSEC=RHODEG*3600.0D 0                                                   
C                                                                               
C INITIALIZE ZERO VALUES FOR THE ELEMENTS OF THE INVOLVED VECTORS.              
C                                                                               
      DO 75 I=1,N                                                               
      F1(I)=0.0D 0                                                              
      F2(I)=0.0D 0                                                              
      F3(I)=0.0D 0                                                              
      NGRAD(I)=0.0D 0                                                           
      EGRAD(I)=0.0D 0                                                           
  75  CONTINUE                                                                  
      DO 1 I=1,12                                                               
  1   IVEC(I)=0                                                                 
C                                                                               
C READ-IN THE COORDINATES (LAT. & LONG.) OF THE THREE POINTS CLOSELY            
C SURROUNDING THE TIDE GUAGE, AND CONTAINING THE PRESSURE                       
C DATA, TO BE USED FOR INTERPOLATION AT THE TIDE-GUAGE ITSELF.                  
C NOTE THAT THE LAT. IS +VE NORTH & THE LONG. IS +VE WEST.                      
C                                                                               
      READ(5,3) PHI1,ALONG1,PHI2,ALONG2,PHI3,ALONG3                             
  3   FORMAT(6F6.2)                                                             
C                                                                               
C COMPUTE THE MEAN LAT. AND LONG. OF THE THREE POINTS FORMING THE TRIANGLE      
C                                                                               
      PHI0=(PHI1+PHI2+PHI3)/3.D 0                                               
      ALONG0=(ALONG1+ALONG2+ALONG3)/3.D 0                                       
C                                                                               
C COMPUTE THE CONVERSION FACTORS FROM CURVILINEAR (LAT. & LONG.)                
C TO LOCAL RECTANGULAR COORDINATES (X, Y).                                      
C                                                                               
      R=6378.D 0                                                                
      FACT1=R/RHODEG                                                            
      PHIRAD=PHI0/RHODEG                                                        
      FACT2=R*DCOS(PHIRAD)/RHODEG                                               
C                                                                               
C COMPUTE THE X (+VE NORTH) AND Y (+VE EAST) LOCAL RECTANGULAR COORDINATES      
C RELATIVE TO THE CENTRE OF THE TRIANGLE AS AN ORIGIN.                          
C IN KILOMETERS UNITS- SINCE THE MEAN RADIUS OF THE EARTH IS IN KM.             
C                                                                               
      X1=(PHI1-PHI0)*FACT1                                                      
      X2=(PHI2-PHI0)*FACT1                                                      
      X3=(PHI3-PHI0)*FACT1                                                      
      Y1=(ALONG0-ALONG1)*FACT2                                                  
      Y2=(ALONG0-ALONG2)*FACT2                                                  
      Y3=(ALONG0-ALONG3)*FACT2                                                  
C                                                                               
C COMPUTE THE DETERMINANT OF THE LEAST-SQUARES APPROXIMATION NORMAL             
C EQUATIONS MATRIX.                                                             
C                                                                               
      XISQ=X1**2+X2**2+X3**2                                                    
      YISQ=Y1**2+Y2**2+Y3**2                                                    
      XIYI=X1*Y1+X2*Y2+X3*Y3                                                    
      DET=XISQ*YISQ-XIYI**2                                                     
C                                                                               
C COMPUTE THE THREE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL FOR            
C THE PRESSURE SURFACE (PLANE) WITHIN THE TRIANGLE UNDER CONSIDERATION.         
C NOTE THAT THESE COEFFICIENTS ARE NOT DEPENDENT ON TIME ,AND ARE               
C CONSTANTS FOR EACH TRIANGLE UNDER CONSIDERATION.                              
C                                                                               
      A1=(YISQ*X1-XIYI*Y1)/DET                                                  
      A2=(YISQ*X2-XIYI*Y2)/DET                                                  
      A3=(YISQ*X3-XIYI*Y3)/DET                                                  
C                                                                               
      B1=(XISQ*Y1-XIYI*X1)/DET                                                  
      B2=(XISQ*Y2-XIYI*X2)/DET                                                  
      B3=(XISQ*Y3-XIYI*X3)/DET                                                  
C                                                                               
C READ-IN AND STORE THE VALUES OF THE OBSERVED PRESSURE AT THE THREE            
C STATIONS FORMING A TRIANGLE AROUND THE TIDE-GUAGE , TO BE USED FOR            
C INTERPOLATION.                                                                
C                                                                               
      K=1                                                                       
      DO 100  I=1,NY                                                            
      JOE=K+11                                                                  
      READ(5,6)ICODE,PHI,ALONG,IYEAR,(F1(J),J=K,JOE)                            
  6   FORMAT(A1,2F3.2,I4,12F5.2)                                                
      K=K+12                                                                    
 100  CONTINUE                                                                  
      K=1                                                                       
      DO 200  I=1,NY                                                            
      JOE=K+11                                                                  
      READ(5,6)ICODE,PHI,ALONG,IYEAR,(F2(J),J=K,JOE)                            
      K=K+12                                                                    
 200  CONTINUE                                                                  
      K=1                                                                       
      DO 300  I=1,NY                                                            
      JOE=K+11                                                                  
      READ(5,6)ICODE,PHI,ALONG,IYEAR,(F3(J),J=K,JOE)                            
      K=K+12                                                                    
 300  CONTINUE                                                                  
C                                                                               
      WRITE(6,51)                                                               
  51  FORMAT('1')                                                               
C                                                                               
C COMPUTE THE MEAN (AVERAGE) VALUE OF THE PRESSURE AT EACH INDIVIDUAL STATION   
C FOR THE WHOLE PERIOD OF RECORDS.                                              
C                                                                               
      NALL=N                                                                    
      NZEROS=0                                                                  
      SUMALL=0.0                                                                
      DO 45  J=1,N                                                              
      SUMALL=SUMALL+F1(J)                                                       
      IF(F1(J).EQ.0.0)  NZEROS=NZEROS+1                                         
  45  CONTINUE                                                                  
      NVALUE=NALL-NZEROS                                                        
      FUNBAR=SUMALL/NVALUE                                                      
      FBAR1=FUNBAR                                                              
      NZEROS=0                                                                  
      SUMALL=0.0                                                                
      DO 245  J=1,N                                                             
      SUMALL=SUMALL+F2(J)                                                       
      IF(F2(J).EQ.0.0)  NZEROS=NZEROS+1                                         
 245  CONTINUE                                                                  
      NVALUE=NALL-NZEROS                                                        
      FUNBAR=SUMALL/NVALUE                                                      
      FBAR2=FUNBAR                                                              
      NZEROS=0                                                                  
      SUMALL=0.0                                                                
      DO 345  J=1,N                                                             
      SUMALL=SUMALL+F3(J)                                                       
      IF(F3(J).EQ.0.0)  NZEROS=NZEROS+1                                         
 345  CONTINUE                                                                  
      NVALUE=NALL-NZEROS                                                        
      FUNBAR=SUMALL/NVALUE                                                      
      FBAR3=FUNBAR                                                              
C                                                                               
C COMPUTE THE VARIATIONS OF PRESSURE AT EACH INDIVIDUAL STATION FROM ITS        
C OWN COMPUTED MEAN VALUE.                                                      
C                                                                               
      DO 145  J=1,N                                                             
      IF(F1(J).EQ.0.0)  GO TO 41                                                
      FUNBAR=FBAR1                                                              
      F1(J)=F1(J)-FUNBAR                                                        
  41  IF(F2(J).EQ.0.0)  GO TO 42                                                
      FUNBAR=FBAR2                                                              
      F2(J)=F2(J)-FUNBAR                                                        
  42  IF(F3(J).EQ.0.0)  GO TO 43                                                
      FUNBAR=FBAR3                                                              
      F3(J)=F3(J)-FUNBAR                                                        
  43  CONTINUE                                                                  
 145  CONTINUE                                                                  
C                                                                               
      DO 500  J=1,N                                                             
C                                                                               
C COMPUTATIONS OF       THE COEFFICIENTS OF THE ASSUMED PRESSURE PLANE          
C IN THE AREA OF THE TRIANGLE JOINING THE THREE GIVEN DATA POINTS.              
C                                                                               
      IF(F1(J).EQ.0.0)  GO TO 500                                               
      IF(F2(J).EQ.0.0)  GO TO 500                                               
      IF(F3(J).EQ.0.0)  GO TO 500                                               
      C0=(F1(J)+F2(J)+F3(J))/3.0D 0                                             
      C1=A1*F1(J)+A2*F2(J)+A3*F3(J)                                             
      C2=B1*F1(J)+B2*F2(J)+B3*F3(J)                                             
      NGRAD(J)=C1*10**7                                                         
      EGRAD(J)=C2*10**7                                                         
 500  CONTINUE                                                                  
C                                                                               
C PRINT-OUT AND PUNCH (ON CARDS) THE NORTH AND EAST GRADIENTS OF THE            
C ATMOSPHERIC PRESSURE PLANE IN THE REGON OF THE TIDE-GAUGE IN QUESTION.        
C                                                                               
      WRITE(6,65)                                                               
  65  FORMAT(///,5X,'VECTOR OF MONTHLY MEANS OF OF NORTH GRADIENT OF THE        
     1 PRESSURE',/,5X,'UNITS OF (INCHES-HG/KM.)*10**7 :')                       
      K=0                                                                       
      ICOUNT=IYSTAR                                                             
      DO 61  I=1,NY                                                             
      DO 62  J=1,12                                                             
      K=K+1                                                                     
  62  IVEC(J)=NGRAD(K)                                                          
      WRITE(6,63) INORTH,ICOUNT,(IVEC(JOE),JOE=1,12)                            
  63  FORMAT(2X,A2,1X,I4,1X,12I8)                                               
      WRITE(7,64) INORTH,ICOUNT,(IVEC(JOE),JOE=1,12)                            
  64  FORMAT(A2,1X,I4,1X,12I6)                                                  
      ICOUNT=ICOUNT+1                                                           
  61  CONTINUE                                                                  
      WRITE(6,66)                                                               
  66  FORMAT(////,5X,'VECTOR OF MONTHLY MEANS OF THE EAST-GRADIENT OF PR        
     1ESSURE',/,5X,'IN UNITS OF (INCHES-HG/KM.)*10**7 :')                       
      K=0                                                                       
      ICOUNT=IYSTAR                                                             
      DO  67  I=1,NY                                                            
      DO  68  J=1,12                                                            
      K=K+1                                                                     
  68  IVEC(J)=EGRAD(K)                                                          
      WRITE(6,63) IEAST,ICOUNT,(IVEC(JOE),JOE=1,12)                             
      WRITE(7,64) IEAST,ICOUNT,(IVEC(JOE),JOE=1,12)                             
      ICOUNT=ICOUNT+1                                                           
  67  CONTINUE                                                                  
      STOP                                                                      
      END                                                                       
