prepare_y.f90

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

!-----------------------------------------------------------------------
! this subroutine calculates Y coefficients in every layer in model
!-----------------------------------------------------------------------
SUBROUTINE PREPARE_Y

  USE NRTYPE              , ONLY: WP
  USE MOD_ATTEN           , ONLY: FJK, FJJ, PLO, FRJ, FRK, I2NRFREQ1,  &
                                  CF, PFP21, PF21, Y, A, B, P, Q
  USE MOD_PARAMETERIZATION, ONLY: FRANGE, FRJMAX, FREF, NRFREQ
  
!----------------------------------------------------------------------

  IMPLICIT NONE

  REAL(WP), PARAMETER  :: PI=3.141592654

  INTEGER    :: I, KF, J, K
  REAL(WP)   :: PF
  REAL(WP)   :: RNRFREQ, RNRFREQ1, R2NRFREQ, R2NRFREQ1, R2NRFREQ2

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

  WRITE(*,'(A)')'Preparing Y coefficients...'

  !-------------------------------------- code for calculating Y coef.

  I2NRFREQ1 = (NRFREQ*2)-1

  ALLOCATE (FJJ(NRFREQ,I2NRFREQ1), FJK(NRFREQ,I2NRFREQ1),              &
            CF (NRFREQ,I2NRFREQ1),                                     &
             Y (NRFREQ), B (NRFREQ), P (NRFREQ), FRJ (NRFREQ),         &
             PFP21         (NRFREQ), PF21            (NRFREQ),         &
             A (NRFREQ, NRFREQ),                                       &
             Q (I2NRFREQ1), PLO (I2NRFREQ1), FRK (I2NRFREQ1) )

  RNRFREQ   = REAL(NRFREQ)
  RNRFREQ1  =   RNRFREQ-1.
  R2NRFREQ  = 2*RNRFREQ
  R2NRFREQ1 = 2*RNRFREQ-1.
  R2NRFREQ2 = 2*RNRFREQ-2.

  FREF = FREF*2.*PI

  DO I = 1, NRFREQ
    FRJ(I)=(10**(LOG10(FRJMAX)-(2.*REAL(I)-1.)/R2NRFREQ*FRANGE))*2.*PI
  END DO

  DO I = 1, I2NRFREQ1
    FRK(I)=(10**(LOG10(FRJMAX)-(   REAL(I)   )/R2NRFREQ*FRANGE))*2.*PI
  END DO

  DO KF = 1, I2NRFREQ1
    PLO(KF)= LOG (FRK(KF)/FREF) / PI !FUTTERMAN Q(w) LAW
  END DO

  DO J = 1, NRFREQ
    PF      = FREF / FRJ(J)
    PF21 (J) = PF*PF + 1._WP
    PFP21(J) = PF / PF21(J)
    DO K = 1, I2NRFREQ1
      FJK(J,K) = FRJ(J)*FRK(K) / ( FRJ(J)*FRJ(J) + FRK(K)*FRK(K) )
      FJJ(J,K) = FRJ(J)*FRJ(J) / ( FRJ(J)*FRJ(J) + FRK(K)*FRK(K) )
      CF (J,K) = FRJ(J) / ( CMPLX(0.,FRK(K)) + FRJ(J) )
    END DO
  END DO


END SUBROUTINE PREPARE_Y