Hi Everyone, <div><br></div><div>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.  </div>
<div><br></div><div>As I said, code is below.  Thanks! </div><div>Sumedh Joshi</div><div>Center for Applied Math</div><div>Cornell University</div><div><br></div><div><br></div><div>-------</div><div><br></div><div><div>subroutine define_external_model(x,y,iflag_element,rho,vp,vs,QKappa_attenuation, &</div>
<div>       Qmu_attenuation,c11,c13,c15,c33,c35,c55)</div><div><br></div><div>  implicit none</div><div><br></div><div>  include "constants.h"</div><div><br></div><div>! user can modify this routine to assign any different external Earth model (rho, vp, vs)</div>
<div>! based on the x and y coordinates of that grid point and the flag of the region it belongs to</div><div><br></div><div>  integer, intent(in) :: iflag_element</div><div><br></div><div>  double precision, intent(in) :: x,y</div>
<div><br></div><div>  double precision, intent(out) :: rho,vp,vs</div><div>  double precision, intent(out) :: QKappa_attenuation,Qmu_attenuation</div><div>  double precision, intent(out) :: c11,c15,c13,c33,c35,c55</div><div>
<br></div><div>  double precision :: zc, zb, e, z1, h0, h1<span class="Apple-tab-span" style="white-space:pre">     </span></div><div><span class="Apple-tab-span" style="white-space:pre">     </span></div><div><span class="Apple-tab-span" style="white-space:pre">     </span>!</div>
<div><span class="Apple-tab-span" style="white-space:pre">      </span>! Encode the Munk profile; no seafloor. </div><div><span class="Apple-tab-span" style="white-space:pre">             </span>!</div><div><span class="Apple-tab-span" style="white-space:pre">            </span>! Make this a fluid layer with a Munk profile. </div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>rho = 1000.d0</div><div><span class="Apple-tab-span" style="white-space:pre">                </span>vs = 0.0d0 </div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre">         </span>!</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>! Set the SOFAR channel depth, going down from the surface. </div><div><span class="Apple-tab-span" style="white-space:pre">         </span>zc = 200.d0</div><div>
<br></div><div><span class="Apple-tab-span" style="white-space:pre">                </span>!</div><div><span class="Apple-tab-span" style="white-space:pre">            </span>! Map to the appropriate z-axis. </div><div><span class="Apple-tab-span" style="white-space:pre">            </span>h0 = 0.d0<span class="Apple-tab-span" style="white-space:pre">                   </span>! Bottom layer thickness.</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>h1 = 600.d0<span class="Apple-tab-span" style="white-space:pre">         </span>! Water layer thickness.</div><div><span class="Apple-tab-span" style="white-space:pre">             </span>z1 = h0 + h1 - y</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span></div><div><span class="Apple-tab-span" style="white-space:pre">             </span>!</div><div><span class="Apple-tab-span" style="white-space:pre">            </span>! Set the sound speed value based on the Munk profile.</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>e = 0.00737d0</div><div><span class="Apple-tab-span" style="white-space:pre">                </span>zb = 2.d0*(z1 - zc)/zc</div><div><span class="Apple-tab-span" style="white-space:pre">               </span>vp = 1500.d0*(1.d0 + e*(zb - 1.d0 + exp(-zb))) </div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span></div><div><span class="Apple-tab-span" style="white-space:pre">             </span>!</div><div><span class="Apple-tab-span" style="white-space:pre">            </span>! Set attenuation. </div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>Qkappa_attenuation = 9999.</div><div><span class="Apple-tab-span" style="white-space:pre">           </span>Qmu_attenuation = 9999.</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre">             </span>!</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>! Make anisotropic. </div><div><span class="Apple-tab-span" style="white-space:pre">         </span>c11 = 0.d0</div><div><span class="Apple-tab-span" style="white-space:pre">           </span>c13 = 0.d0</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>c15 = 0.d0</div><div><span class="Apple-tab-span" style="white-space:pre">           </span>c33 = 0.d0</div><div><span class="Apple-tab-span" style="white-space:pre">           </span>c35 = 0.d0</div>
<div><span class="Apple-tab-span" style="white-space:pre">              </span>c55 = 0.d0<span class="Apple-tab-span" style="white-space:pre">          </span></div><div><span class="Apple-tab-span" style="white-space:pre">             </span></div><div>
  end subroutine define_external_model</div></div><div><br></div>