C**************************
      SUBROUTINE PRTRB1
C**************************
      include 'para.h'

      xmdot=pxmdot
      xmd=pxmd
      xlin=pxlin
      do 10 j=1,jx-1
         DE   = W0(J,1)
         VX   = W0(J,2)   
         TE   = W0(J,4)
         R    = X1(J)
         RM1  = R-1
         R2   = R*R
         R3   = R2*R
         OK   = (R/RM1)*SQRT(0.5/R3)
         PR   = (AA/3.)*(TE**4)+RMU*DE*TE
         HH   = SQRT(B5*PR/DE)/OK
         SS   = 2*AIN*DE*HH*RG
         vx   = -xmdot/(2.*pi*r*ss)/cc
         W0(J,2)=VX
 10         continue
      return
      end
C**************************
      SUBROUTINE PRTRB2
C**************************
      include 'para.h'
      xmdot=pxmdot
      xmd=pxmd

      do 10 j=jx,jx

         CALL SADMN2(X1(J),ETA,DL,W0(J,1),W0(J,2),W0(J,3),W0(J,4))
 10         continue
      RETURN
      END

C**************************
      SUBROUTINE PRTRB3
C**************************
      include 'para.h'
      dimension ww(4)
c-----------------------
c      rc     : Radius 
c      xlamda : wave length
c      rp     : amplitude of the wave 
c-----------------------
      rc=50.0
      xlamda=1.0
      rp=1.0
      do 10 j=1,jx-1
         do 20 k=1,4
            w0(j,k)=w0(j,k)
     &      +rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda)
 20               continue
 10                  continue
      return
      end

C**************************
      SUBROUTINE PRTRB4
C**************************
      include 'para.h'
      dimension ww(4)
c-----------------------
c      rc     : Radius 
c      xlamda : wave length
c      rp     : amplitude of the wave 
c-----------------------
      rc=30.0
      xkk=30
      rp=2.0

      do 10 j=1,jx-1
         DE   = W0(J,1)
         VX   = W0(J,2)   
         TE   = W0(J,4)
         R    = X1(J)
         RM1  = R-1
         R2   = R*R
         R3   = R2*R
         OK   = (R/RM1)*SQRT(0.5/R3)
         PR   = (AA/3.)*(TE**4)+RMU*DE*TE
         HH   = SQRT(B5*PR/DE)/OK
         SS   = 2*AIN*DE*HH*RG
         vx   = -xmdot/(2.*pi*r*ss)/cc
         W0(J,2)=VX

 20            continue
         dww = dss/xkk/x1(j)*w0(j,1)
         dvr = dss/sqrt(xkk*x1(j))*w0(j,2)
         dvp = dss/sqrt(xkk*x1(j))*w0(j,3)
 10         continue
      write(6,*) "PERTRB 4"
      return
      end
C**************************
