prepare_fine_par_y.f90

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

SUBROUTINE PREPARE_FINE_PAR_Y (L, H, Z0, NRFREQ)

  USE NRTYPE              , ONLY: WP
  USE MOD_ATTEN           , ONLY: FRK, FRJ, PFP21, PF21, CF,           &
                                  I2NRFREQ1
  USE MOD_INTEGRATION     , ONLY: PTS, PTS2,                           &
                                  M1_FINE, YM1_FINE
  USE MOD_PARAMETERIZATION, ONLY: FREF
  USE MOD_INTERFACES      , ONLY: FUNC_YM

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

  IMPLICIT NONE

  INTEGER , INTENT(IN) :: L, NRFREQ
  REAL(WP), INTENT(IN) :: H, Z0

  INTEGER              :: LI, J
  
  REAL(WP)             :: HPTS_H, H2_H, ZSTART,                        &
                          POFRJ, POFRJ2P1,                             &
                          THE1, THE2, R, M1, Z

  REAL(WP)   , DIMENSION(  NRFREQ ) :: YY

  COMPLEX(WP), DIMENSION(I2NRFREQ1) :: SUM1

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


  HPTS_H  =  H/REAL(PTS,WP)
  H2_H    = (H - HPTS_H)/2._WP

  ZSTART = Z0 + REAL(L,WP)*H - H2_H

  DO LI = PTS+PTS2, PTS2+1, -1
    Z = ZSTART + REAL(LI-1,WP)*HPTS_H

    YY  = FUNC_YM (Z, NRFREQ)

    SUM1 = CMPLX(0.,0.)
    THE1 = 0._WP
    THE2 = 0._WP
   
    DO J = 1, NRFREQ

      POFRJ    = FREF / FRJ(J)
      POFRJ2P1 = POFRJ * POFRJ + 1._WP

      SUM1(:)  = SUM1(:) + YY(J)*CF(J,:)

      THE1    = THE1 + YY(J) / PF21 (J)
      THE2    = THE2 + YY(J) * PFP21(J)
  
    END DO

    THE1 = 1._WP - THE1
    R    = SQRT( THE1*THE1 + THE2*THE2 )

    M1=  M1_FINE(LI) / ((R+THE1)/(2._WP*R*R))

    YM1_FINE (LI,:) = M1 / (1._WP - SUM1(:))

     M1_FINE (LI  ) = M1

  END DO

  IF ( L == 0 ) THEN

    DO LI = 1, PTS2
      YM1_FINE (LI,:) = YM1_FINE (PTS+1-LI,:)
       M1_FINE (LI  ) =  M1_FINE (PTS+1-LI  )
    END DO

  ELSE

    DO LI = 1, PTS2
      Z = ZSTART + REAL(LI-1,WP)*HPTS_H

      YY  = FUNC_YM (Z, NRFREQ)

      SUM1 = CMPLX(0.,0.)
      THE1 = 0._WP
      THE2 = 0._WP
      
      DO J = 1, NRFREQ

        POFRJ    = FREF / FRJ(J)
        POFRJ2P1 = POFRJ * POFRJ + 1._WP

        SUM1(:)  = SUM1(:) + YY(J)*CF(J,:)

        THE1    = THE1 + YY(J) / PF21 (J)
        THE2    = THE2 + YY(J) * PFP21(J)
      
      END DO

      THE1 = 1._WP - THE1
      R    = SQRT( THE1*THE1 + THE2*THE2 )

       M1=  M1_FINE(LI) / ((R+THE1)/(2._WP*R*R))

      YM1_FINE (LI,:) = M1 / (1._WP - SUM1(:))

       M1_FINE (LI  ) = M1
     
    END DO

  END IF

END SUBROUTINE PREPARE_FINE_PAR_Y