ADded subroutine to interpolate seismograms from located source.

apps/print_cmt.F90

@@ -19,5 +19,4 @@ program print_cmt
     call print_source(sources(i), 1)
   enddo
 
-
 end program print_cmt

src/constants.F90

@@ -125,4 +125,12 @@ module constants
   integer, parameter :: NLINES_PER_CMTSOLUTION_SOURCE = 13
   integer, parameter :: NLINES_PER_FORCESOLUTION_SOURCE = 11
 
+  ! Perturbation values
+  double precision,parameter :: dmom  = 1.0d23
+  double precision,parameter :: dlat  = 0.0001d0
+  double precision,parameter :: dlon  = 0.0001d0
+  double precision,parameter :: ddep  = 0.01d0
+  double precision,parameter :: dcmt  = -1.d0
+  double precision,parameter :: dhdur = 0.001d0
+
 end module constants

src/gf3d/CMakeLists.txt

@@ -2,6 +2,7 @@ add_library(gf3d SHARED
   ${CMAKE_CURRENT_LIST_DIR}/gf3d.F90
   ${CMAKE_CURRENT_LIST_DIR}/seismograms.F90
   ${CMAKE_CURRENT_LIST_DIR}/seismograms_windows.F90
+  ${CMAKE_CURRENT_LIST_DIR}/interpolate_source.F90
   # ${CMAKE_CURRENT_LIST_DIR}/gf3d_write_seismograms.F90
 )
 

src/gf3d/gf3d.F90

@@ -27,6 +27,17 @@ module gf3d
     init_log, finalize_log, &
     fftpack
 
+  ! interface interpolate_source
+  !   module subroutine interpolate_source(GF, source, seismograms)
+  !     use gf, only: t_GF
+  !     use sources, only: t_source
+  !     use interpolation, only: interpolateMT
+  !     type(t_GF), intent(in) :: GF
+  !     type(t_source), intent(in) :: source
+  !     double precision, dimension(:,:,:) :: seismograms
+  !   end subroutine interpolate_source
+  ! end interface interpolate_source
+
   interface get_seismograms
     module subroutine get_seismograms(GF, sources, superseismograms)
       use gf, only: t_GF

src/gf3d/interpolate_source.F90

@@ -0,0 +1,52 @@
+submodule (gf3d:seismograms) source_interpolation
+
+  implicit none
+
+contains
+
+  module subroutine interpolate_source(GF, source, seismograms)
+
+    use gf, only: t_GF
+    use sources, only: t_source
+    use interpolation, only: interpolateMT
+
+    ! In
+    type(t_GF), intent(in) :: GF
+    type(t_source), intent(in) :: source
+    ! Out
+    double precision, dimension(:,:,:) :: seismograms
+
+    ! Local
+    integer :: i,j,k,iglob
+    double precision, dimension(:,:,:,:,:,:,:), allocatable :: displacement
+
+    ! Allocate temporary displacement array
+    allocate(displacement(size(GF%displacement,1), 3, 3, GF%ngllx,GF%nglly,GF%ngllz,GF%nsteps))
+
+    ! Get displacement array
+    do k=1,GF%ngllz
+      do j=1,GF%nglly
+        do i=1,GF%ngllx
+          iglob = GF%ibool(i,j,k, source%ispec)
+          displacement(:,:,:,i,j,k, :) = GF%displacement(:,:,:,iglob,:)
+        enddo
+      enddo
+    enddo
+
+    if (source%force) then
+      call throwerror(-1, "Force source interpolation not implemented.")
+    else
+      call interpolateMT(&
+          displacement, size(GF%displacement,1), &
+          GF%nsteps, GF%ngllx, GF%nglly, GF%ngllz, GF%xigll, GF%yigll, GF%zigll, &
+          source%Mxx, source%Myy, source%Mzz, &
+          source%Mxy, source%Mxz, source%Myz, &
+          source%xi, source%eta,source%gamma, &
+          source%xix, source%xiy, source%xiz, &
+          source%etax, source%etay, source%etaz, &
+          source%gammax, source%gammay, source%gammaz, &
+          seismograms)
+    endif
+  end subroutine interpolate_source
+
+end submodule