qylsm.f90

1D有限差分波动方程模拟

SUBROUTINE QYLSM(Q,YY,NRFREQ)

  USE NRTYPE   , ONLY: WP
  USE MOD_ATTEN, ONLY: FJJ, FJK, A, B, P, I2NRFREQ1

!----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER                           , INTENT(IN   ) :: NRFREQ
  REAL(WP), DIMENSION (1:  NRFREQ  ), INTENT(OUT) :: YY
  REAL(WP), DIMENSION (1:2*NRFREQ-1), INTENT(INOUT) :: Q

  INTEGER  :: I, J, K
  REAL(WP) :: FJ, FI, SUMN

!----------------------------------------------------------------------

  DO J = 1, NRFREQ
    B(J) = 0.
    A(J,:) = 0.
    DO K = 1, I2NRFREQ1
      FJ = FJK(J,K) + FJJ(J,K)/Q(K)
      B(J) = B(J)+ FJ/Q(K)
      DO I = 1, NRFREQ
        FI = FJK(I,K) + FJJ(I,K)/Q(K)
        A(J,I) = A(J,I) + FI*FJ
      END DO
    END DO
  END DO

  DO I = 1, NRFREQ
    SUMN = A(I,I)-DOT_PRODUCT(A(I,1:I-1),A(I,1:I-1))
    IF (SUMN <= 0.0) THEN
      PRINT *, 'CHOLECKY DECOMPOSITION FAILED'
      STOP
    END IF
    P(I)=SQRT(SUMN)
    A(I+1:NRFREQ,I)= ( A(I,I+1:NRFREQ)                                 &
                     - MATMUL(A(I+1:NRFREQ,1:I-1),A(I,1:I-1)) ) / P(I)
  END DO

  DO I = 1, NRFREQ
    YY(I)=(B(I)-DOT_PRODUCT(A(I,1:I-1),YY(1:I-1)))/P(I)
  END DO
  DO I = NRFREQ, 1, -1
    YY(I)=(YY(I)-DOT_PRODUCT(A(I+1:NRFREQ,I),YY(I+1:NRFREQ)))/P(I)
  END DO
  

END SUBROUTINE QYLSM