C                                                                               
C***********************************************************************        
C*                                                                     *        
C*  PROGRAM 'T A P I N T':                                             *        
C*  ---------------------                                              *        
C*                                                                     *        
C*                            WRITTEN BY:  M.M. NASSAR,  JULY 1975.    *        
C*                                                                     *        
C*  FOR PREDICTING THE TEMPERATURE AND PRESSURE VARIATIONS AT A        *        
C*  TIDE-GUAGE, FROM THE AVAILABLE INFORMATION AT  ANY THREE POINTS    *        
C*  CLOSELY SURROUNDING THE TIDE-GUAGE: BY APPLYING THE LEAST-SQUARES  *        
C*  APPROXIMATION (INTERPOLATION) TECHNIQUES.                          *        
C*                                                                     *        
C***********************************************************************        
C                                                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      DIMENSION SLG(700),DTG(700),DPG(700),RD(4,700)                            
      DIMENSION GNAME(5)                                                        
      DIMENSION F1(700),F2(700),F3(700)                                         
      DIMENSION DUMVEC(700),IRD(12)                                             
      DIMENSION RDBAR(3)                                                        
      DOUBLE PRECISION DABS,DSIN,DCOS,DSQRT                                     
C                                                                               
C READ-IN THE TOTAL NUMBER OF YEARS IN WHICH DATA IS AVAILABLE,                 
C THE TOTAL NUMBER OF INPUTTED RIVER DISCHARGE DATA SETS, AND THE DATA REFERENCE
C NUMBER THAT CONTAINS THE RELEVANT INFORMATION AT THE TIDE GUAGE.              
C WHICH IS STORED ON DISC.                                                      
C                                                                               
      READ(5,1) NYEARS,NRIVTL,IDG                                               
  1   FORMAT(I2,2X,I2,2X,I2)                                                    
C                                                                               
C COMPUTE THE SIZE OF THE VECTORS OF DATA (MSL, TEMP., ..ETC.).                 
C                                                                               
      NSIZE=NYEARS*12                                                           
C                                                                               
C READ-IN THE HEADER CARD OF THE TIDE GUAGE STATION: NUMBER; NAME;              
C LAT.; LONG.; STARTING YEAR OF OBS.; AND STARTING MONTH IN THIS YEAR.          
C                                                                               
      READ(5,2) NUMG,(GNAME(I),I=1,5),PHIG,ALONGG,IYST,IMST                     
  2   FORMAT(I6,4X,5A4,2F5.2,4X,I4,2X,I2)                                       
C                                                                               
C READ-IN THE COORDINATES (LAT. & LONG.) OF THE THREE POINTS CLOSELY            
C SURROUNDING THE TIDE GUAGE, AND CONTAINING TEMPERATURE AND PRESSURE           
C DATA, TO BE USED FOR INTERPOLATION AT THE TIDE-GUAGE ITSELF.                  
C                                                                               
C   NOTE THAT IF ONE DATA SET (TEMP., PRESSURE) FROM ONLY ONE STATION           
C   WILL BE CONSIDERED, A NUMBER GREATER THAN '10' SHOULD BE INPUT ON           
C   THIS CARD ,ANYWHERE BETWEEN COLS. 1-5, AND THE REST OF THE CARD BLANK.      
C                                                                               
      READ(5,3) PHI1,ALONG1,PHI2,ALONG2,PHI3,ALONG3                             
  3   FORMAT(6F6.2)                                                             
      IF(PHI1.GT.10.0)  GO TO 999                                               
C                                                                               
C SET-UP SOME CONVERSION CONSTANTS.                                             
C                                                                               
      PI=3.141592653589793D 0                                                   
      RHODEG=180.D 0/PI                                                         
      RHOSEC=RHODEG*3600.0D 0                                                   
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 WEST) RECTANGULAR LOCAL COORDINATES      
C RELATIVE TO THE CENTRE OF THE TRIANGLE AS AN ORIGIN.                          
C                                                                               
      X1=(PHI1-PHI0)*FACT1                                                      
      X2=(PHI2-PHI0)*FACT1                                                      
      X3=(PHI3-PHI0)*FACT1                                                      
      XG=(PHIG-PHI0)*FACT1                                                      
      Y1=(ALONG1-ALONG0)*FACT2                                                  
      Y2=(ALONG2-ALONG0)*FACT2                                                  
      Y3=(ALONG3-ALONG0)*FACT2                                                  
      YG=(ALONGG-ALONG0)*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 (TEMP., PRESSURE) AT THE TIDE-GUAGE LOCATION.                                 
C NOTE THAT THESE COEFFICIENTS ARE NOT DEPENDENT ON TIME ,AND ARE               
C CONSTANTS FOR EACH TRIANGLE UNDER CONSIDERATION.                              
C                                                                               
      D1=1.D 0/3.D 0+(YG*Y1*XISQ+XG*X1*YISQ-(XG*Y1+YG*X1)*XIYI)/DET             
      D2=1.D 0/3.D 0+(YG*Y2*XISQ+XG*X2*YISQ-(XG*Y2+YG*X2)*XIYI)/DET             
      D3=1.D 0/3.D 0+(YG*Y3*XISQ+XG*X3*YISQ-(XG*Y3+YG*X3)*XIYI)/DET             
C                                                                               
 999  CONTINUE                                                                  
C                                                                               
C INITIALIZE ZERO VALUES FOR THE ELEMENTS OF DATA VECTORS AT THE TIDE-GUAGE.    
C                                                                               
      DO 10  I=1,NSIZE                                                          
      SLG(I)=0.0D 0                                                             
      DTG(I)=0.0D 0                                                             
      DPG(I)=0.0D 0                                                             
      F1(I)=0.0D 0                                                              
      F2(I)=0.0D 0                                                              
      F3(I)=0.0D 0                                                              
      DUMVEC(I)=0.0D 0                                                          
  10  CONTINUE                                                                  
      DO 20  I=1,4                                                              
      DO 20  J=1,NSIZE                                                          
      RD(I,J)=0.0D 0                                                            
  20  CONTINUE                                                                  
C                                                                               
C READ-IN AND STORE THE MONTHLY MEANS OF M.S.L. DATA AT THE TIDE-               
C GUAGE, FOR THE WHOLE PERIOD OF INVESTIGATION.                                 
C                                                                               
      K=1                                                                       
      DO 30  I=1,NYEARS                                                         
      JACK=K+11                                                                 
      READ(5,4)ICODE,PHIG,ALONGG,IYEAR,(SLG(J),J=K,JACK)                        
  4   FORMAT(A1,2F3.2,I4,12F5.2)                                                
      K=K+12                                                                    
  30  CONTINUE                                                                  
C                                                                               
C READ-IN AND STORE THE MONTHLY MEANS OF THE RIVERS DISCHARGE DATA              
C SETS (UP TO THREE), NEAR THE TIDE-GUAGE.                                      
C                                                                               
C IT SHOULD BE NOTED THAT HERE THERE IS A POSSIBILITY TO ACCEPT MORE THAN       
C THREE SETS OF RIVER DISCHARGES- IN WHICH CASE TWO OR MORE SETS OF THEM        
C HAVE TO BE ADDED (LUMPED) TOGETHER.  THE CHOSEN CRITERION HERE IS :-          
C EACH SET OF RIVER DISCHARGE DATA DECK HAS TO BE ENDED-UP WITH A BLANK         
C CARD, UNLESS THE CURRENT SET IN THE INPUT SEQUENCE IS REQUIRED TO BE          
C ADDED TO THE PREVIOUS ONE, WHERE THIS ENDING BLANK CARD OF THE CURRENT        
C SET MUST INCLUDE A NON-ZERO INTEGER IN COLUMN NUBER ONE.                      
C                                                                               
      IF(NRIVTL.EQ.0)  GO TO 33                                                 
      ILUMP=0                                                                   
      DO 41  II=1,NRIVTL                                                        
      KK=1                                                                      
      IROW=II-ILUMP                                                             
      DO 40  I=1,NYEARS                                                         
      JIM=KK+11                                                                 
      READ(5,5)ICODE,PHI,ALONG,IYEAR,(RD(IROW,J),J=KK,JIM)                      
  5   FORMAT(A1,2F3.2,I4,12G5.0)                                                
      KK=KK+12                                                                  
  40  CONTINUE                                                                  
      READ(5,300) IADD                                                          
 300  FORMAT(I1)                                                                
      IF(IADD.EQ.0)  GO TO 41                                                   
      ISTORE=IROW-1                                                             
      DO 401  J=1,NSIZE                                                         
 401  RD(ISTORE,J)=RD(ISTORE,J)+RD(IROW,J)                                      
      ILUMP=ILUMP+1                                                             
  41  CONTINUE                                                                  
      NRIVER=NRIVTL-ILUMP                                                       
C                                                                               
C CLEANING THE DUMMY VALUES STORED INTO THE RIVER DISCHARGES ARRAY ,            
C WHICH EXCEED THE ACTUAL NUMBER OF DESIRED RIVER DISCHARGES, AND THEN          
C RE-INITIALIZE THEM BACK TO ZEROS.                                             
C                                                                               
      IF(NRIVER.LT.3)  GO TO 625                                                
      GO TO 650                                                                 
 625  ISAMY=NRIVER+1                                                            
      DO 675  I=ISAMY,4                                                         
      DO 675  J=1,NSIZE                                                         
 675  RD(I,J)=0.0D 0                                                            
C                                                                               
 650  CONTINUE                                                                  
C                                                                               
  33  CONTINUE                                                                  
C                                                                               
C PRINT-OUT THE REQUIRED INFORMATION AT THE TIDE-GUAGE STATION.                 
C                                                                               
      WRITE(6,15) IDG                                                           
  15  FORMAT('1',//,5X,'GENERAL INFORMATION ABOUT THE TIDE-GUAGE :-',/,5        
     1X,'----------------------------------------',//,5X,'DATA SET REFER        
     2ENCE NUMBER ON DISC IS =',I4/)                                            
      WRITE(6,16) NUMG,(GNAME(I),I=1,5)                                         
  16  FORMAT(5X,'TIDE-GUAGE NUMBER IS =',I8//,5X,'TIDE-GUAGE NAME IS   =        
     1',5A4//)                                                                  
      WRITE(6,17) IYST,IMST,NSIZE                                               
  17  FORMAT(5X,'DATA STARTING YEAR IS =',I5//,5X,'DATA STARTING MONTH O        
     1F THE ABOVE YEAR IS =',I3//,5X,'SIZE OF ARRAYS OF DIFFERENT TYPES         
     2OF DATA IS =',I4//)                                                       
      WRITE(6,501) PHI1,ALONG1,PHI2,ALONG2,PHI3,ALONG3                          
 501  FORMAT(//,5X,'LAT. & LONG. OF STATION NO. 1 ARE =',2F6.2,/,5X,'LAT        
     1. & LONG. OF STATION NO. 2 ARE =',2F6.2,/,5X,'LAT. & LONG. OF STAT        
     2ION NO. 3 ARE =',2F6.2/)                                                  
      WRITE(6,18)                                                               
  18  FORMAT('1',//,5X,'VECTOR OF MONTHLY MEAN SEA LEVEL IS:')                  
      K=1                                                                       
      DO 90  I=1,NYEARS                                                         
      KL=K+11                                                                   
      WRITE(6,19) (SLG(J),J=K,KL)                                               
  19  FORMAT(10X,12F8.2)                                                        
      K=K+12                                                                    
  90  CONTINUE                                                                  
      WRITE(6,21)                                                               
  21  FORMAT(/////,5X,'FIRST SET OF MONTHLY MEAN RIVER DISCHARGE IS:')          
      K=0                                                                       
      DO 100  I=1,NYEARS                                                        
      DO 101  J=1,12                                                            
      K=K+1                                                                     
      IRD(J)=RD(1,K)                                                            
 101  CONTINUE                                                                  
      WRITE(6,22) (IRD(JOE),JOE=1,12)                                           
  22  FORMAT(10X,12I8)                                                          
 100  CONTINUE                                                                  
      WRITE(6,23)                                                               
  23  FORMAT(//,5X,'SECOND SET OF MONTHLY MEAN RIVER DISCHARGE IS:')            
      K=0                                                                       
      DO 110  I=1,NYEARS                                                        
      DO 111  J=1,12                                                            
      K=K+1                                                                     
      IRD(J)=RD(2,K)                                                            
 111  CONTINUE                                                                  
      WRITE(6,22) (IRD(JOE),JOE=1,12)                                           
 110  CONTINUE                                                                  
      WRITE(6,24)                                                               
  24  FORMAT(//,5X,'THIRD SET OF MONTHLY MEAN RIVER DISCHARGE IS:')             
      K=0                                                                       
      DO 120  I=1,NYEARS                                                        
      DO 121  J=1,12                                                            
      K=K+1                                                                     
      IRD(J)=RD(3,K)                                                            
 121  CONTINUE                                                                  
      WRITE(6,22) (IRD(JOE),JOE=1,12)                                           
 120  CONTINUE                                                                  
C                                                                               
      IF(PHI1.GT.10.0)  GO TO 888                                               
C                                                                               
C CALL THE SUBROUTINE 'F U N I N T', FOR THE                                    
C PREDICTION OF THE MONTHLY MEANS OF TEMPERATURE AT THE                         
C TIDE-GUAGE, FROM THE INFORMATION AVAILABLE AT THE SURROUNDING THREE           
C POINTS, BY THE LEAST-SQUARES APPROXIMATION TECHNIQUES.                        
C                                                                               
      WRITE(6,70)                                                               
  70  FORMAT('1'//,'        THE FOLLOWING PLOT AND ITS ASSOCIATED LEGEND        
     1 PERTAIN TO THE ***TEMPERATURE***')                                       
C                                                                               
      CALL FUNINT(D1,D2,D3,NYEARS,NSIZE,F1,F2,F3,DTG)                           
C                                                                               
C CALL THE SUBROUTINE 'F U N I N T', FOR THE                                    
C PREDICTION OF THE MONTHLY MEANS OF THE PRESSURE AT THE TIDE-                  
C GUAGE STATION FROM THE AVAILABLE INFORMATION AT THE SURROUNDING THREE         
C STATIONS, BY THE LEAST-SQUARES APPROXIMATION TECHNIQUES.                      
C                                                                               
      WRITE(6,80)                                                               
  80  FORMAT('1'//,'        THE FOLLOWING PLOT AND ITS ASSOCIATED LEGEND        
     1PERTAINS TO THE ***PRESSURE***')                                          
C                                                                               
      CALL FUNINT(D1,D2,D3,NYEARS,NSIZE,F1,F2,F3,DPG)                           
C                                                                               
      GO TO 777                                                                 
 888  CONTINUE                                                                  
C                                                                               
C READ-IN THE ONLY ONE AVAILABLE SET OF OBSERVED TEMPERATURE, AT ONE            
C STATION, AND HENCE STORE IT AS IF IT WAS TAKEN AT THE TIDE-GUAGE.             
C                                                                               
      K=1                                                                       
      DO 210  I=1,NYEARS                                                        
      JOE=K+11                                                                  
      READ(5,211)ICODE,PHI,ALONG,IYEAR,(DTG(J),J=K,JOE)                         
 211  FORMAT(A1,2F3.2,I4,12F5.1)                                                
      K=K+12                                                                    
 210  CONTINUE                                                                  
C                                                                               
C READ-IN THE ONLY ONE AVAILABLE SET OF OBSERVED PRESSURE AT ONE                
C STATION, AND HENCE STORE IT AS IF IT WAS TAKEN AT THE TIDE-GUAGE.             
C                                                                               
      K=1                                                                       
      DO 212  I=1,NYEARS                                                        
      JOE=K+11                                                                  
      READ(5,211)ICODE,PHI,ALONG,IYEAR,(DPG(J),J=K,JOE)                         
      K=K+12                                                                    
 212  CONTINUE                                                                  
C                                                                               
 777  CONTINUE                                                                  
C                                                                               
C CALL THE SUBROUTINE 'T S P L O T' ,TO COMPUTE AND PRINT-OUT THE MEAN          
C MAXIMUM AND MINIMUM VALUES OF THE PREDICTED  TEMPERATURE AT THE TIDE-         
C GUAGE, AS WELL AS PLOTTING THE CORRESPONDING TIME-SERIES.                     
C                                                                               
      WRITE(6,460)                                                              
 460  FORMAT('1',//,5X,'***THE FOLLOWING PLOT IS FOR THE PREDICTED TEMPE        
     1RATURE SERIES AT THE TIDE-GUAGE***'/)                                     
C                                                                               
      CALL TSPLOT(DTG,NSIZE,TGBAR)                                              
C                                                                               
C                                                                               
C COMPUTE THE VARIATIONS IN TEMPERATURE VALUES (FROM THEIR MEAN VALUE)          
C FOR THE PREDICTED VALUES AT THE TIDE-GUAGE.                                   
C                                                                               
      DO 216  I=1,NSIZE                                                         
      IF(DTG(I).EQ.0.0)  GO TO 216                                              
      DTG(I)=DTG(I)-TGBAR                                                       
 216  CONTINUE                                                                  
C                                                                               
C PRINT-OUT THE 'PREDICTED' VARIATIONS IN TEMPERATURE AT THE TIDE-GUAGE.        
C                                                                               
      WRITE(6,25)                                                               
  25  FORMAT('1'//,5X,'PREDICTED MONTHLY MEAN TEMPERATURE VARIATIONS AT         
     1THE TIDE-GUAGE ARE:')                                                     
      K=1                                                                       
      DO 130  I=1,NYEARS                                                        
      KL=K+11                                                                   
      WRITE(6,26) (DTG(J),J=K,KL)                                               
  26  FORMAT(10X,12F8.2)                                                        
      K=K+12                                                                    
 130  CONTINUE                                                                  
C                                                                               
C CALL THE SUBROUTINE 'T S P L O T' ,TO COMPUTE AND PRINT-OUT THE MEAN          
C MAXIMUM AND MINIMUM VALUES OF THE PREDICTED  PRESSURE    AT THE TIDE-         
C GUAGE, AS WELL AS PLOTTING THE CORRESPONDING TIME-SERIES.                     
C                                                                               
      WRITE(6,465)                                                              
 465  FORMAT('1',//,5X,'***THE FOLLOWING PLOT IS FOR THE PREDICTED PRESS        
     1URE SERIES AT THE TIDE-GUAGE***'/)                                        
C                                                                               
      CALL TSPLOT(DPG,NSIZE,PGBAR)                                              
C                                                                               
C                                                                               
C COMPUTE THE VARIATIONS IN THE PRESSURE VALUES (FROM THEIR MEAN VALUE)         
C FOR THE PREDICTED VALUES AT THE TIDE-GUAGE.                                   
C                                                                               
      DO 221  I=1,NSIZE                                                         
      IF(DPG(I).EQ.0.0)  GO TO 221                                              
      DPG(I)=DPG(I)-PGBAR                                                       
 221  CONTINUE                                                                  
C                                                                               
C                                                                               
C PRINT-OUT THE 'PREDICTED' VARIATIONS IN PRESSURE AT THE TIDE-GUAGE.           
C                                                                               
      WRITE(6,27)                                                               
  27  FORMAT('1'//,5X,'PREDICTED MONTHLY MEAN PRESSURE VARIATIONS AT THE        
     1 TIDE-GUAGE ARE:')                                                        
      K=1                                                                       
      DO 140  I=1,NYEARS                                                        
      KL=K+11                                                                   
      WRITE(6,26) (DPG(J),J=K,KL)                                               
      K=K+12                                                                    
 140  CONTINUE                                                                  
      IF(NRIVTL.EQ.0)  GO TO 32                                                 
C                                                                               
C CALL THE SUBROUTINE 'T S P L O T' ,TO COMPUTE AND PRINT-OUT THE MEAN          
C MAXIMUM AND MINIMUM VALUES OF THE SETS OF ACCEPTED RIVER DISCHARGE            
C SERIES TO BE UNDER CONSIDERATION (UP TO 3 SETS), AS WELL AS PLOTTING          
C THE CORRESPONDING TIME SERIES FOR EACH SET.                                   
C                                                                               
      DO 201  I=1,NRIVER                                                        
      WRITE(6,470) I                                                            
 470  FORMAT('1',//,5X,'***THE FOLLOWING PLOT IS FOR SET NO.',I2,' OF RI        
     1VER DISCHARGE SERIES***'/)                                                
      DO 202  J=1,NSIZE                                                         
 202  DUMVEC(J)=RD(I,J)                                                         
C                                                                               
      CALL TSPLOT(DUMVEC,NSIZE,DISBAR)                                          
      RDBAR(I)=DISBAR                                                           
C                                                                               
 201  CONTINUE                                                                  
C                                                                               
C COMPUTE THE VARIATIONS IN THE MONTHLY RIVER DISCHARGE (FROM THE CORRESPONDING 
C MEAN VALUE)  AT EVERY AVAILABLE RIVER LOCATION.                               
C                                                                               
      DO 203  I=1,NRIVER                                                        
      DO 203  J=1,NSIZE                                                         
      IF(RD(I,J).EQ.0.0)  GO TO 203                                             
      RD(I,J)=RD(I,J)-RDBAR(I)                                                  
 203  CONTINUE                                                                  
C                                                                               
C PRINT-OUT THE COMPUTED VARIATIONS IN RIVER DISCHARGE VALUES AT EACH           
C GIVEN RIVER LOCATION.                                                         
C                                                                               
      WRITE(6,207)                                                              
 207  FORMAT('1')                                                               
      DO 204  I=1,NRIVER                                                        
      WRITE(6,205) I                                                            
 205  FORMAT(/,/,/,5X,'COMPUTED VARIATIONS IN RIVER DISCHARGE VALUES FOR        
     1 **SET NO.',I2,'** ARE:'/)                                                
      K=0                                                                       
      DO 204  JJ=1,NYEARS                                                       
      DO 206  J=1,12                                                            
      K=K+1                                                                     
      IRD(J)=RD(I,K)                                                            
 206  CONTINUE                                                                  
      WRITE(6,22) (IRD(JOE),JOE=1,12)                                           
 204  CONTINUE                                                                  
  32  CONTINUE                                                                  
C                                                                               
C*******************************************************************************
C                                                                               
C WRITE THE NECESSARY INFORMATION AT THE TIDE GUAGE ON THE DISC--               
C TO BE STORED FOR USE BY ANOTHER PROGRAM.                                      
C                                                                               
      REWIND IDG                                                                
      WRITE(IDG,301) NUMG                                                       
 301  FORMAT(I6)                                                                
      DO 302  I=1,5                                                             
      WRITE(IDG,303) GNAME(I)                                                   
 303  FORMAT(A4)                                                                
 302  CONTINUE                                                                  
      WRITE(IDG,304) IYST                                                       
 304  FORMAT(I4)                                                                
      WRITE(IDG,305) IMST                                                       
 305  FORMAT(I2)                                                                
      WRITE(IDG,306) NSIZE                                                      
 306  FORMAT(I3)                                                                
C                                                                               
      DO 307  I=1,NSIZE                                                         
      WRITE(IDG,308) SLG(I)                                                     
 308  FORMAT(F6.2)                                                              
 307  CONTINUE                                                                  
      DO 309  I=1,NSIZE                                                         
 309  WRITE(IDG,308) DTG(I)                                                     
      DO 310  I=1,NSIZE                                                         
 310  WRITE(IDG,308) DPG(I)                                                     
C                                                                               
      DO 311  I=1,3                                                             
      DO 311  J=1,NSIZE                                                         
      INT=RD(I,J)                                                               
      WRITE(IDG,312) INT                                                        
 312  FORMAT(I7)                                                                
 311  CONTINUE                                                                  
C                                                                               
      STOP                                                                      
      END                                                                       
C                                                                               
C***************************************************************************    
      SUBROUTINE FUNINT(D1,D2,D3,NY,NS,F1,F2,F3,FG)                             
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      DIMENSION F1(NS),F2(NS),F3(NS),FG(NS)                                     
C                                                                               
C READ-IN AND STORE THE VALUES OF THE OBSERVED FUNCTION 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.1)                                                
      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                                                                               
C INITIALIZE ZERO VALUES FOR THE FUNCTIONAL VALUES AT THE TIDE-GUAGE.           
C                                                                               
      DO 75  I=1,NS                                                             
      FG(I)=0.0D 0                                                              
  75  CONTINUE                                                                  
C                                                                               
C COMPUTE THE VALUES OF THE INTERPOLATED FUNCTION AT THE TIDE-GUAGE.            
C                                                                               
      DO 400  I=1,NS                                                            
      FG(I)=D1*F1(I)+D2*F2(I)+D3*F3(I)                                          
 400  CONTINUE                                                                  
C                                                                               
C CALLING THE SUBROUTINE 'F P L O T', TO PLOT THE THREE TIME SERIES             
C OF OBSERVED FUNCTIONAL VALUES AT THE THREE STATIONS ,AND ALSO THE             
C PREDICTED VALUES AT THE TIDE-GUAGE- TO VISUALIZE ITS VARIATIONS.              
C                                                                               
      CALL FPLOT(F1,F2,F3,FG,NS)                                                
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
C***************************************************************************    
      SUBROUTINE FPLOT(F1,F2,F3,FG,NS)                                          
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      DIMENSION F1(NS),F2(NS),F3(NS),FG(NS)                                     
      DIMENSION IPLOT(30),IFRAME(30)                                            
      DOUBLE PRECISION DABS,DMAX1,DMIN1,DFLOAT                                  
      DATA IBLANK/' '/,IDASH/'|'/,IPLUS/'M'/                                    
      DATA IF1/'.'/,IF2/'+'/,IF3/'X'/,IFG/'*'/                                  
C                                                                               
C SPECIFY THE DESIRED LATERAL SPACE FOR THE PLOT ,IN TERMS OF THE               
C EQUIVALENT NUMBER OF CHARACTERS ON THE LINE PRINTER 'NCOL'.                   
C                                                                               
      NCOL=30                                                                   
C                                                                               
C SET-UP VALUES OF 'M' FOR THE ELEMENTS OF THE ARRAY 'IFRAME',                  
C TO BE USED IN FRAMING THE PLOT.                                               
C                                                                               
      DO 9  I=1,NCOL                                                            
      IFRAME(I)=IPLUS                                                           
  9   CONTINUE                                                                  
C                                                                               
C INITIALIZE SOME REASONABLE VALUES - TO HELP EXTRACTING THE MAXIMUM            
C AND MINIMUM OBSERVED (OR PREDICTED) VALUES OF THE FUNCTION UNDER              
C CONSIDERATION- THIS IS REQUIRED MAINLY FOR PLOTTING ROUTINES.                 
C                                                                               
      FMAX=0.0D 0                                                               
      FMIN=1.0D 06                                                              
C                                                                               
C EXTRACT THE MAX. & MIN. VALUES OF THE FUNCTION THAT OCCURRED AT ANY           
C OF THE THREE STATIONS OR THE TIDE-GUAGE.                                      
C                                                                               
      DO 500  I=1,NS                                                            
      FMAX=DMAX1(FMAX,F1(I))                                                    
      FMAX=DMAX1(FMAX,F2(I))                                                    
      FMAX=DMAX1(FMAX,F3(I))                                                    
      FMAX=DMAX1(FMAX,FG(I))                                                    
      IF(F1(I).EQ.0.0)  GO TO 10                                                
      FMIN=DMIN1(FMIN,F1(I))                                                    
  10  IF(F2(I).EQ.0.0)  GO TO 20                                                
      FMIN=DMIN1(FMIN,F2(I))                                                    
  20  IF(F3(I).EQ.0.0)  GO TO 30                                                
      FMIN=DMIN1(FMIN,F3(I))                                                    
  30  IF(FG(I).EQ.0.0)  GO TO 500                                               
      FMIN=DMIN1(FMIN,FG(I))                                                    
 500  CONTINUE                                                                  
C                                                                               
C  COMPUTE THE MEAN VALUE OF THE INTERPOLATED (PREDICTED) FUNCTION              
C AT THE TIDE-GUAGE.                                                            
C                                                                               
C NOTE:- THE ZERO VALUES WHICH REPRESENT MISSING DATA ARE EXCLUDED              
C        FROM THE COMPUTATIONS.                                                 
C                                                                               
      NZEROG=0                                                                  
      SUMFG=0.0D 0                                                              
      DO 600  I=1,NS                                                            
      IF(FG(I).EQ.0.0)  NZEROG=NZEROG+1                                         
      SUMFG=SUMFG+FG(I)                                                         
 600  CONTINUE                                                                  
      NN=NS-NZEROG                                                              
      FGBAR=SUMFG/DFLOAT(NN)                                                    
C                                                                               
C INITIALIZE THE PLOTTING ARRAY.                                                
C                                                                               
      DO 60  I=1,NCOL                                                           
      IPLOT(I)=IBLANK                                                           
  60  CONTINUE                                                                  
      IPLOT(1)=IPLUS                                                            
      IPLOT(NCOL)=IPLUS                                                         
C                                                                               
      X=DFLOAT(NCOL-1)/(FMAX-FMIN)                                              
      Y=FMIN*X-1.0                                                              
      KBAR=FGBAR*X-Y                                                            
      IPLOT(KBAR)=IDASH                                                         
C                                                                               
C PRINT-OUT A LEGEND FOR THE DIFFERENT FUNCTIONS IN THE PLOT                    
C                                                                               
      WRITE(6,5)                                                                
  5   FORMAT(/,5X,'LEGEND :',2X,'.-REPRESENTS THE FUNCTIONAL VALUES AS O        
     1BSERVED AT STATION NO. 1 ',/,15X,'+-REPRESENTS THE FUNCTIONAL VALU        
     2ES AS OBSERVED AT STATION NO. 2',/,15X,'X-REPRESENTS THE FUNCTIONA        
     3L VALUES AS OBSERVED AT STATION NO. 3',/,15X,'*-REPRESENTS THE FUN        
     4CTIONAL VALUES AS PREDICTED AT THE TIDE-GUAGE',/,15X,'M-REPRESENTS        
     5 THE FRAME AROUND THE PLOT')                                              
C                                                                               
C PRINT-OUT THE MEAN, MAX. & MIN. VALUES OF THE TIME SERIES.                    
C                                                                               
      WRITE(6,15) NS,FMAX,FMIN,FGBAR                                            
  15  FORMAT(/,/,10X,'TOTAL NUMBER OF OBSERVED FUNCTIONAL VALUES IS =',I        
     15,/,10X,'MAXIMUM OBSERVED FUNCTIONAL VALUE IS =',F8.3,/,10X,'MINIM        
     1UM OBSERVED FUNCTIONAL VALUE IS =',F8.3,/,10X,'MEAN PREDICTED FUNC        
     2TIONAL VALUE AT THE TIDE-GUAGE IS =',F8.3//)                              
C                                                                               
C COMPUTE AND PRINT-OUT THE SCALE OF THE PLOT, BY CONSIDERING THE LINE          
C PRINTER # 1403 (AT U.N.B.), WHICH PRINTS 10 CHARACTERS PER INCH.              
C                                                                               
      SINCH=10.0                                                                
      SCALE=DABS(FMAX-FMIN)/(DFLOAT(NCOL)/SINCH)                                
      WRITE(6,70) SCALE                                                         
  70  FORMAT(7X,'NO.',5X,'SCALE OF THE PLOT IS: 1 INCH REPRESENTS',E10.3        
     1,' UNITS')                                                                
C                                                                               
C PLOT THE GIVEN TIME SERIES, WITH ITS FUNCTIONAL VALUES.                       
C                                                                               
      WRITE(6,40) (IFRAME(J),J=1,NCOL)                                          
C                                                                               
      DO 72  I=1,NS                                                             
      KF1=F1(I)*X-Y                                                             
      KF2=F2(I)*X-Y                                                             
      KF3=F3(I)*X-Y                                                             
      KFG=FG(I)*X-Y                                                             
      IF(KF1.LT.1)  KF1=1                                                       
      IF(KF2.LT.1)  KF2=1                                                       
      IF(KF3.LT.1)  KF3=1                                                       
      IF(KFG.LT.1)  KFG=1                                                       
      IF(KF1.GT.NCOL)  KF1=NCOL                                                 
      IF(KF2.GT.NCOL)  KF2=NCOL                                                 
      IF(KF3.GT.NCOL)  KF3=NCOL                                                 
      IF(KFG.GT.NCOL)  KFG=NCOL                                                 
C                                                                               
      IPLOT(KF1)=IF1                                                            
      IPLOT(KF2)=IF2                                                            
      IPLOT(KF3)=IF3                                                            
      IPLOT(KFG)=IFG                                                            
C                                                                               
      WRITE(6,80) I,(IPLOT(J),J=1,NCOL)                                         
C                                                                               
C RE-INITIALIZE THE PLOTTING ARRAY.                                             
C                                                                               
      IPLOT(KF1)=IBLANK                                                         
      IPLOT(KF2)=IBLANK                                                         
      IPLOT(KF3)=IBLANK                                                         
      IPLOT(KFG)=IBLANK                                                         
      IPLOT(1)=IPLUS                                                            
      IPLOT(NCOL)=IPLUS                                                         
      IPLOT(KBAR)=IDASH                                                         
  72  CONTINUE                                                                  
      WRITE(6,40) (IFRAME(J),J=1,NCOL)                                          
  40  FORMAT(15X,30A1)                                                          
  80  FORMAT(6X,I4,5X,30A1)                                                     
      RETURN                                                                    
      END                                                                       
C                                                                               
C*******************************************************************************
C                                                                               
C  SUBROUTINE 'T S P L O T':                                                    
C                            FOR COMPUTING AND PRINTING THE MAXIMUM,            
C  MINIMUM AND MEAN VALUES OF A GIVEN TIME SERIES. ALSO IT PLOTS THE            
C  CORRESPONDING TIME SERIES ALONG WITH THE FUNCTIONAL VALUES AND               
C  THE SCALE OF THE PLOT.                                                       
C                                                                               
C*******************************************************************************
      SUBROUTINE TSPLOT(F,N,FBAR)                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      DIMENSION F(N)                                                            
      DIMENSION IPLOT(30),IFRAME(30)                                            
      DOUBLE PRECISION DABS,DMAX1,DMIN1,DFLOAT                                  
      DATA IBLANK/' '/,ISTAR/'*'/,IPLUS/'+'/,IDASH/'|'/                         
C                                                                               
C SPECIFY THE DESIRED LATERAL SPACE FOR THE PLOT ,IN TERMS OF THE               
C EQUIVALENT NUMBER OF CHARACTERS ON THE LINE PRINTER 'NCOL'.                   
C                                                                               
      NCOL=30                                                                   
C                                                                               
C SET-UP VALUES OF '+' FOR THE ELEMENTS OF THE ARRAY 'IFRAME',                  
C TO BE USED IN FRAMING THE PLOT.                                               
C                                                                               
      DO 9  I=1,NCOL                                                            
      IFRAME(I)=IPLUS                                                           
  9   CONTINUE                                                                  
C                                                                               
C COMPUTE THE MEAN, MAXIMUM AND MINIMUM FUNCTIONAL VALUES FOR THE               
C GIVEN TIME-SERIES (REQUIRED TO BE PLOTTED).                                   
C                                                                               
C NOTE:- THE ZERO VALUES WHICH REPRESENT MISSING DATA ARE EXCLUDED              
C        FROM THE COMPUTATIONS.                                                 
C                                                                               
      NZEROF=0                                                                  
      SUMF=0.0                                                                  
      DO 1  I=1,N                                                               
      IF(F(I).EQ.0.0)  NZEROF=NZEROF+1                                          
      SUMF=SUMF+F(I)                                                            
  1   CONTINUE                                                                  
      NSUMF=N-NZEROF                                                            
      FBAR=SUMF/DFLOAT(NSUMF)                                                   
C                                                                               
      FMAX=0.0D 0                                                               
      FMIN=1.0D 10                                                              
      DO 2  I=1,N                                                               
      FMAX=DMAX1(FMAX,F(I))                                                     
      IF(F(I).EQ.0.0)  GO TO 2                                                  
      FMIN=DMIN1(FMIN,F(I))                                                     
  2   CONTINUE                                                                  
C                                                                               
C INITIALIZE THE PLOTTING ARRAY.                                                
C                                                                               
      DO 10  I=1,NCOL                                                           
      IPLOT(I)=IBLANK                                                           
  10  CONTINUE                                                                  
      IPLOT(1)=IPLUS                                                            
      IPLOT(NCOL)=IPLUS                                                         
C                                                                               
      X=DFLOAT(NCOL-1)/(FMAX-FMIN)                                              
      Y=1.0-FMIN*X                                                              
      KBAR=FBAR*X+Y                                                             
      IPLOT(KBAR)=IDASH                                                         
C                                                                               
C PRINT-OUT THE MEAN, MAX. & MIN. VALUES OF THE TIME SERIES.                    
C                                                                               
      WRITE(6,60) FMAX,FMIN,FBAR                                                
  60  FORMAT(5X,'MAX. FUNCTIONAL VALUE IN THE SERIES = ',E11.4,/,5X,'MIN        
     1. FUNCTIONAL VALUE IN THE SERIES = ',E11.4,/,5X,'MEAN FUNCTIONAL V        
     2ALUE IN THE SERIES = ',E11.4/)                                            
C                                                                               
C COMPUTE AND PRINT-OUT THE SCALE OF THE PLOT, BY CONSIDERING THE LINE          
C PRINTER # 1403 (AT U.N.B.), WHICH PRINTS 10 CHARACTERS PER INCH.              
C                                                                               
      SINCH=10.0                                                                
      SCALE=DABS(FMAX-FMIN)/(DFLOAT(NCOL)/SINCH)                                
      WRITE(6,70) SCALE                                                         
  70  FORMAT(6X,'NO.',5X,'VALUE',5X,'SCALE OF THE PLOT IS: 1 INCH REPRES        
     1ENTS',E10.3,' UNITS')                                                     
C                                                                               
C PLOT THE GIVEN TIME SERIES, WITH ITS FUNCTIONAL VALUES.                       
C                                                                               
      WRITE(6,40) (IFRAME(J),J=1,NCOL)                                          
      DO 30  I=1,N                                                              
      KF=F(I)*X+Y                                                               
      IF(KF.LT.1)  KF=1                                                         
      IF(KF.GT.NCOL)  KF=NCOL                                                   
      IPLOT(KF)=ISTAR                                                           
      WRITE(6,20) I,F(I),(IPLOT(J),J=1,NCOL)                                    
C                                                                               
C RE-INITIALIZE THE PLOTTING ARRAY.                                             
C                                                                               
      IPLOT(KF)=IBLANK                                                          
      IPLOT(1)=IPLUS                                                            
      IPLOT(NCOL)=IPLUS                                                         
      IPLOT(KBAR)=IDASH                                                         
  30  CONTINUE                                                                  
      WRITE(6,40) (IFRAME(J),J=1,NCOL)                                          
  20  FORMAT(5X,I4,2X,E11.4,3X,30A1)                                            
  40  FORMAT(25X,30A1)                                                          
      RETURN                                                                    
      END                                                                       
