[CIG-SEISMO] SPECFEM2D: Using define_external_model.f90

[Sumedh Joshi]

Hi Everyone,

I'm trying to use SPECFEM2D with an external model (for those of you that
do underwater acoustics, I'm trying to implement a Munk sound speed profile
in the water column).  To do this, I've modified define_external_model.f90
to include the sound speed variation that I want (the code is pasted
below), but I get no meaningful results.  After a few timesteps, all of the
seismograms are NaNs, although the code continues running for as long as
I'd like.  I'm trying to figure out if I'm doing something obviously wrong,
as I have never gotten define_external_model.f90 to operate correctly, even
when I just set constant material parameters with no spatial variation.

As I said, code is below.  Thanks!
Sumedh Joshi
Center for Applied Math
Cornell University


-------

subroutine
define_external_model(x,y,iflag_element,rho,vp,vs,QKappa_attenuation, &
       Qmu_attenuation,c11,c13,c15,c33,c35,c55)

  implicit none

  include "constants.h"

! user can modify this routine to assign any different external Earth model
(rho, vp, vs)
! based on the x and y coordinates of that grid point and the flag of the
region it belongs to

  integer, intent(in) :: iflag_element

  double precision, intent(in) :: x,y

  double precision, intent(out) :: rho,vp,vs
  double precision, intent(out) :: QKappa_attenuation,Qmu_attenuation
  double precision, intent(out) :: c11,c15,c13,c33,c35,c55

  double precision :: zc, zb, e, z1, h0, h1
 !
! Encode the Munk profile; no seafloor.
!
! Make this a fluid layer with a Munk profile.
rho = 1000.d0
vs = 0.0d0

!
! Set the SOFAR channel depth, going down from the surface.
zc = 200.d0

!
! Map to the appropriate z-axis.
h0 = 0.d0 ! Bottom layer thickness.
h1 = 600.d0 ! Water layer thickness.
z1 = h0 + h1 - y
 !
! Set the sound speed value based on the Munk profile.
e = 0.00737d0
zb = 2.d0*(z1 - zc)/zc
vp = 1500.d0*(1.d0 + e*(zb - 1.d0 + exp(-zb)))
 !
! Set attenuation.
Qkappa_attenuation = 9999.
Qmu_attenuation = 9999.

!
! Make anisotropic.
c11 = 0.d0
c13 = 0.d0
c15 = 0.d0
c33 = 0.d0
c35 = 0.d0
c55 = 0.d0
  end subroutine define_external_model
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://geodynamics.org/pipermail/cig-seismo/attachments/20130610/4546ae70/attachment.html>