C  * * * * * * * * * * * * * * * * * * * *                              
C     SUBROUTINE FOR OHSAKI'S SPECTRUM                                  
C  * * * * * * * * * * * * * * * * * * * *                              
C                                                                       
C                    CODED BY Y.OHSAKI                                  
C                                                                       
C     PUBPOSE                                                           
C       TO COMPUTE, FOR A GIVEN PERIOD, THE SPECTRL VELOCITY IN         
C       OHSAKI'S SPECTRUM WHEN THE MAGNITUDE AND EPICENTRAL DISTANCE    
C       OF AN EARTHQUAKE AND THE DAMPING FACTOR OF THE SPECTRUM ARE     
C       SPECIFIED                                                       
C                                                                       
C     USAGE                                                             
C       CALL OHSP(EM,R,H,T,SV,ICALL)                                    
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       EM    - MAGNITUDE                                               
C       R     - EPICENTRAL DISTANCE IN KILOMETERS                       
C       H     - DAMPING FACTOR IN DECIMAL FRACTION                      
C       T     - PERIOD IN SEC                                           
C       SV    - SPECTRAL VELOCITY IN KINES                              
C       ICALL - ICALL.EQ.0 FOR FIRST CALL                               
C               ICALL.NE.0 FOR SUBSEQUENT CALLS WITH THE SAME EM,R,H    
C                          AND DIFFERENT T                              
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       VELK                                                            
C                                                                       
      SUBROUTINE OHSP(EM,R,H,T,SV,ICALL)                                
C                                                                       
      DIMENSION   D(3,3),W(3,3,8),U(3,3,8),UK(8)                        
      EQUIVALENCE (UK(1),TB),(UK(2),TC),(UK(3),TD),(UK(4),SVA),         
     *            (UK(5),SVB),(UK(6),SVC),(UK(7),SVD),(UK(8),SVE)       
      SAVE        UK                                                    
      DATA        D/5.,15.,60.,10.,45.,150.,25.,120.,350./,             
     *            TA/.02/,TE/2./                                        
      DATA        W/0.10,.12,.14,.20,.28,.28,.35,.45,.50,               
     *              0.14,.25,.32,.40,.45,.50,.60,.70,.70,               
     *              0.80,.65,.60,1.2,.90,.80,1.5,1.2,1.0,               
     *              1.22,.76,.62,.52,.42,.46,.28,.30,.36,               
     *               15., 9., 8.,12.,10.,11., 4.,10.,18.,               
     *               21.,25.,27.,25.,27.,32.,30.,33.,38.,               
     *               21.,25.,27.,25.,27.,32.,30.,33.,38.,               
     *                9.,12.,15.,20.,24.,27.,28.,32.,37./               
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      IF(ICALL.NE.0) GO TO 210                                          
      DO 130 K=1,8                                                      
      DO 120 J=1,3                                                      
      DO 110 I=1,3                                                      
      U(I,J,K)=W(I,J,K)                                                 
  110 CONTINUE                                                          
  120 CONTINUE                                                          
  130 CONTINUE                                                          
C                                                                       
C     EFFECT OF DAMPING FACTOR                                          
C                                                                       
      T0=10.**(0.31*EM-1.2)                                             
      DO 160 K=5,8                                                      
      DO 150 J=1,3                                                      
      DO 140 I=1,3                                                      
      U(I,J,K)=U(I,J,K)/SQRT(1.+17.*(H-0.05)*EXP(-2.5*U(I,J,K-4)/T0))   
  140 CONTINUE                                                          
  150 CONTINUE                                                          
  160 CONTINUE                                                          
C                                                                       
C     INTERPOLATION                                                     
C                                                                       
      J=INT((EM+0.0001)/7.)+1                                           
      EMR=EM-REAL(J+5)                                                  
      R1=ALOG(D(2,J+1))*EMR+ALOG(D(2,J))*(1.-EMR)                       
      RLOG=ALOG(R)                                                      
      IF(RLOG-R1) 170,180,180                                           
  170 I=1                                                               
      R2=R1                                                             
      R1=ALOG(D(1,J+1))*EMR+ALOG(D(1,J))*(1.-EMR)                       
      GO TO 190                                                         
  180 I=2                                                               
      R2=ALOG(D(3,J+1))*EMR+ALOG(D(3,J))*(1.-EMR)                       
  190 RR=(RLOG-R1)/(R2-R1)                                              
      IF(RR.LT.0.) RR=0.                                                
      VMULT=VELK(EM,R)/10.                                              
      DO 200 K=1,8                                                      
      U1=ALOG(U(I,J+1,K))*EMR+ALOG(U(I,J,K))*(1.-EMR)                   
      U2=ALOG(U(I+1,J+1,K))*EMR+ALOG(U(I+1,J,K))*(1.-EMR)               
      UK(K)=EXP(U2*RR+U1*(1.-RR))                                       
      IF(K.LE.3) GO TO 200                                              
      UK(K)=UK(K)*VMULT                                                 
  200 CONTINUE                                                          
C                                                                       
C     SPECTRAL VELOCITY FOR GIVEN PERIOD                                
C                                                                       
  210 IF(T.GT.TA) GO TO 220                                             
      SV=T/TA*SVA                                                       
      RETURN                                                            
  220 IF(T.GE.TB) GO TO 230                                             
      SV=T*(SVA/TA+(T-TA)/(TB-TA)*(SVB/TB-SVA/TA))                      
      RETURN                                                            
  230 IF(T.GE.TC) GO TO 240                                             
      SV=T*(SVB/TB+(T-TB)/(TC-TB)*(SVC/TC-SVB/TB))                      
      RETURN                                                            
  240 IF(T.GE.TD) GO TO 250                                             
      SV=EXP(ALOG(SVC)+ALOG(T/TC)/ALOG(TD/TC)*ALOG(SVD/SVC))            
      RETURN                                                            
  250 SV=EXP(ALOG(SVD)+ALOG(T/TD)/ALOG(TE/TD)*ALOG(SVE/SVD))            
      RETURN                                                            
      END                                                               
