C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *    
C     SUBROUTINE FOR FREQUENCY RESPONSE FUNCTION OF LAYERED SOIL        
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *    
C                                                                       
C                                              CODED BY Y.OHSAKI        
C                                                                       
C     PURPOSE                                                           
C       TO COMPUTE ACCELERATION FREQUENCY-RESPONSE-FUNCTION AT THE TOP  
C       OF SPECIFIED OBJECTIVE LAYER OR STRAIN FREQUENCY-RESPONSE-      
C       FUNCTION AT THE MIDDLE POINT OF SPECIFIED OBJECTIVE LAYER IN    
C       TERMS OF INPUT ACCELERATION AT THE TOP OF SPECIFIED REFERENCE   
C       LAYER IN A HORIZONTALLY LAYERED SOIL DEPOSIT                    
C                                                                       
C     USAGE                                                             
C       CALL FESP(L,TH,UW,G,ALPHA,BETA,ND1,LOBJ,LREF,DF,N,H,ND2,        
C                IND1,IND2)                                             
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       L          - TOTAL NUMBER OF LAYERS, INCLUDING THE BASE LAYER   
C                    L.LE.60                                            
C       TH(ND1)    - THICKNESS IN METERS                                
C       UW(ND1)    - UNIT WEIGHT IN TONS/CU.METER                       
C       G(ND1)     - SHEAR MODULUS IN TONS/SQ.METER                     
C       ALPHA(ND1) - SCATTER DAMPING COEFFICIENT IN 1/SEC               
C       BETA(ND1)  - MATERIAL DAMPING COEFFICIENT IN DECIMAL FRACTION   
C       ND1        - DIMENSION OF TH,UW,G,ALPHA,BETA IN CALLING PROGRAM 
C       LOBJ       - NUMBER OF OBJECTIVE LAYER                          
C       LREF       - NUMBER OF REFERENCE LAYER                          
C       DF         - FREQUENCY INCREMENT IN FREQUENCY RESPONSE FUNCTION 
C                    IN HZ                                              
C       N          - TOTAL NUMBER OF DATA IN FREQUENCY RESPONSE         
C                    FUNCTION                                           
C       H(ND2)     - FREQUENCY RESPONSE FUNCTION IN COMPLEX NUMBER      
C       ND2        - DIMENSION OF H IN CALLING PROGRAM                  
C       IND1       - 0 FOR ACCELERATION-TO-ACCELERATION RESPONSE        
C                    1 FOR STRAIN-TO-ACCELERATION RESPONSE              
C       IND2       - 0 IF REFERENCE LAYER IS INTERNAL                   
C                    1 IF REFERENCE LAYER IS EXPOSED                    
C                                                                       
C     REMARKS                                                           
C       THE THICKNESS OF BASE LAYER TH(L) IS NOT REQUIRED               
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       NONE                                                            
C                                                                       
      SUBROUTINE FESP(L,TH,UW,G,ALPHA,BETA,ND1,LOBJ,LREF,DF,N,H,ND2,    
     *               IND1,IND2)                                         
C                                                                       
      DIMENSION TH(ND1),UW(ND1),G(ND1),ALPHA(ND1),BETA(ND1)             
      COMPLEX   H(ND2),GB(60),P(60),R(59),A(60),B(60),EA,EB,EC,         
     *          HOBJ,HREF                                               
      PARAMETER (P2=6.283185)                                           
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      W=0.                                                              
      DW=P2*DF                                                          
      H(1)=CMPLX(REAL(1-IND1),0.)                                       
      A(1)=(1.,0.)                                                      
      B(1)=(1.,0.)                                                      
      DO 170 K=2,N                                                      
      W=W+DW                                                            
C                                                                       
C     IMPEDANCE RATIOS                                                  
C                                                                       
      DO 110 I=1,L                                                      
      CALL SUB(W,ALPHA(I),BETA(I),DAMP)                                 
      GB(I)=G(I)*CMPLX(1.,2.*DAMP)                                      
      P(I)=CSQRT(UW(I)/9.8/GB(I))*(0.,1.)                               
  110 CONTINUE                                                          
      DO 120 I=1,L-1                                                    
      R(I)=GB(I)/GB(I+1)*P(I)/P(I+1)                                    
  120 CONTINUE                                                          
C                                                                       
C     FREQUENCY RESPONSE FUNCTIONS                                      
C                                                                       
      DO 160 I=1,L                                                      
      IF(I.NE.LOBJ) GO TO 140                                           
      IF(IND1.EQ.1) GO TO 130                                           
      HOBJ=A(I)+B(I)                                                    
      GO TO 140                                                         
  130 EA=CEXP(P(I)*W*TH(I)/(2.,0.))                                     
      EB=(1.,0.)/EA                                                     
      HOBJ=-P(I)/W*(A(I)*EA-B(I)*EB)*0.01                               
  140 IF(I.NE.LREF) GO TO 150                                           
      HREF=A(I)+B(I)                                                    
      IF(IND2.NE.0) HREF=(2.,0.)*A(I)                                   
  150 IF(I.EQ.L) GO TO 160                                              
      EC=CEXP(P(I)*W*TH(I))                                             
      EA=A(I)*EC                                                        
      EB=B(I)/EC                                                        
      A(I+1)=((1.+R(I))*EA+(1.-R(I))*EB)/(2.,0.)                        
      B(I+1)=((1.-R(I))*EA+(1.+R(I))*EB)/(2.,0.)                        
  160 CONTINUE                                                          
      H(K)=HOBJ/HREF                                                    
  170 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C                                                                       
C     SUBROUTINE FOR DAMPING FACTOR                                     
C                                                                       
      SUBROUTINE SUB(W,ALPHA,BETA,DAMP)                                 
C                                                                       
      DAMP=ALPHA/W+BETA                                                 
      RETURN                                                            
      END                                                               
                                                                        
