C 
         DOUBLE PRECISION FUNCTION TRAP(X,Y,N,XMAX,NLEAD)               
         IMPLICIT REAL*8(A-H,O-Z)                                       
         DIMENSION X(NLEAD),Y(NLEAD)                                    
         DO 5 I=1,N                                                     
         IF(X(I).GT.XMAX) GO TO 6                                       
         NP=I                                                           
 5       CONTINUE                                                       
 6       SUM=0.0D0                                                      
         TRAP=0.0D0                                                     
         IF(NP.EQ.1) GO TO 20                                           
         DO 10 I=2,NP                                                   
         SUM=SUM+(Y(I)+Y(I-1))*(X(I)-X(I-1))                            
 10      CONTINUE                                                       
         TRAP=0.5D0*SUM                                                 
 20      IF(NP.EQ.N) RETURN                                             
         AUX=(XMAX-X(NP))/(X(NP+1)-X(NP))                               
         YMAX=(1.0D0-AUX)*Y(NP)+AUX*Y(NP+1)                             
         SMALL=(YMAX+Y(NP))*(XMAX-X(NP))*0.5D0                          
         TRAP=TRAP+SMALL                                                
         RETURN                                                         
         END                                                            
C                                                                       
      SUBROUTINE LOCATE(X,R,NMAX,J,NLEAD)                               
      IMPLICIT REAL*8(A-H,O-Z)                                          
      DIMENSION R(NLEAD),IC(13)                                         
      DATA IC/2048,1024,512,256,128,64,32,16,8,4,2,1,1/ 
      J=4096                                                            
      DO 20 IL=1,11                                                     
      IF(J.GT.NMAX) GOTO 12                                             
      IF(X.GT.R(J)) GOTO 15                                             
12    J=J-IC(IL)                                                        
      GOTO 20                                                           
15    J=J+IC(IL)                                                        
20    CONTINUE                                                          
      IF(J.EQ.0) RETURN                                                 
      IF(J.GT.NMAX) J=NMAX                                              
      IF(X.LT.R(J)) J=J-1                                               
      RETURN                                                            
      END                                                               
C                                                                       
         SUBROUTINE CONVOL1(X,F,N,Y,H,FWHM,NLEAD)                       
         IMPLICIT REAL*8(A-H,O-Z)                                       
         PARAMETER(NXXLD=5500)
         PARAMETER(NXYLD=201)
         DIMENSION X(NLEAD),F(NLEAD),Y(NLEAD),H(NLEAD)
         DIMENSION FG(NXYLD),XY(NXYLD)
         DIMENSION XX(NXXLD),FF(NXXLD)                                  
C                                                                       
         DO 5 I=1,N                                                     
         Y(I)=X(I)                                                      
         H(I)=F(I)                                                      
 5       CONTINUE                                                       
         IF(FWHM.LT.0.0D0) THEN 
         WRITE(6,*) 'CONVOL1 WARN: NO CONVOLUTION, FWHM.LT.0.0D0' 
         RETURN 
         ENDIF 
         NN=N+200                                                       
         IF(NN.GT.NXXLD) THEN                                           
         WRITE(6,*) 'CONVOL1 WARN: CHECK XX AND FF DIMENSION'           
         RETURN                                                         
         ENDIF                                                          
         AUXT=2.0D0*FWHM                                                
         DELTA=0.020D0*FWHM                                             
         XMIN=X(1)-AUXT                                                 
         X1=X(1)                                                        
         X2=X(2)                                                        
         F1=DLOG10(F(1))                                                
         F2=DLOG10(F(2))                                                
         DO 10 I=1,100                                                  
         XX(I)=XMIN+DELTA*(I-1)                                         
         AUX=(XX(I)-X1)/(X2-X1)                                         
         FF(I)=10.0D0**((1.0D0-AUX)*F1+AUX*F2)                          
 10      CONTINUE                                                       
         DO 12 I=1,N                                                    
         II=100+I                                                       
         XX(II)=X(I)                                                    
         FF(II)=F(I)                                                    
 12      CONTINUE                                                       
         XMIN=X(N)                                                      
         X1=X(N-1)                                                      
         X2=X(N)                                                        
         F1=DLOG10(F(N-1))                                              
         F2=DLOG10(F(N))                                                
         DO 14 I=1,100                                                  
         II=100+N+I                                                     
         XX(II)=XMIN+DELTA*I                                            
         AUX=(XX(II)-X1)/(X2-X1)                                        
         FF(II)=10.0D0**((1.0D0-AUX)*F1+AUX*F2)                         
 14      CONTINUE                                                       
C                                                                       
         CONST=-4.0D0*DLOG(2.0D0)/FWHM**2                               
         AA=1.0D0/(1.0644644D0*FWHM)                                    
         DO 30 I=1,N                                                    
         XMIN=Y(I)-AUXT                                                 
         XMAX=Y(I)+AUXT                                                 
         IF(XMIN+.05.LT.XX(1)) WRITE(6,*) 'CONVOL1 WARN: XMIN<XX(1)'    
         IF(XMAX-.05.GT.XX(NN)) WRITE(6,*) 'CONVOL1 WARN: XMAX>XX(NN)'  
         DO 20 J=1,201                                                  
         XY(J)=XMIN+DELTA*(J-1)                                         
         CALL LOCATE(XY(J),XX,NN,K,NXXLD)                               
         IF(K.EQ.0) K=1                                                 
         IF(K.EQ.NN) K=NN-1                                             
         CALL CUBIN(XY(J),AUX1,XX,FF,NN,K,NXXLD)                        
         AUX2=AA*DEXP(CONST*(XY(J)-Y(I))**2)                            
         FG(J)=AUX1*AUX2                                                
 20      CONTINUE                                                       
         H(I)=TRAP(XY,FG,NXYLD,XY(NXYLD),NXYLD)                         
 30      CONTINUE                                                       
         RETURN                                                         
         END                                                            
C                                                                       
         SUBROUTINE CUBIN(X,Y,TX,TY,NMAX,I,NLEAD)                       
         IMPLICIT REAL*8(A-H,O-Z)                                       
C                                                                       
C LAGRANGE CUBIC INTERPOLATION ROUTINE                                  
C                                                                       
         DIMENSION TX(NLEAD),TY(NLEAD)                                  
C                                                                       
         KS=1                                                           
         IF(I.EQ.1) KS=2                                                
         NP=NMAX-1                                                      
         KE=4                                                           
         IF(I.EQ.NP) KE=3                                               
C                                                                       
         Y=0.0D0                                                        
         DO 20 K=KS,KE                                                  
         II=I+K-2                                                       
         RLII=1.0D0                                                     
         DO 10 KP=KS,KE                                                 
         JJ=I+KP-2                                                      
         IF(II.EQ.JJ) GO TO 10                                          
         RLII=RLII*(X-TX(JJ))/(TX(II)-TX(JJ))                           
 10      CONTINUE                                                       
         Y=Y+RLII*TY(II)                                                
 20      CONTINUE                                                       
C                                                                       
         RETURN                                                         
         END                                                            
C                                                                       
