Using gdb, a backtrace of the crash can be found:
#0 0x0000000000407b0c in f90sac::f90sac_tshift at f90sac.F90:1126
#1 0x000000000041ae78 in analyser::an_applysplittingoperators at Analyser.F90:138
#2 0x000000000041d23c in evaluatepath at EvaluatePath.f90:84
#3 0x000000000041c1e8 in evaluatemodel at Evaluate.F90:84
#4 0x000000000041d550 in mhsampler at Samplers.f90:95
#5 0x000000000041bcec in mtssampmain at MTSSampMain.f90:47
The program is crashing on line 1126 of f90sac.F90:
1098 !===============================================================================
1099 subroutine f90sac_tshift(trace,dt)
1100 !===============================================================================
1101 !
1102 ! time shift a sac trace by dt (to nearest sample), zero additional
1103 ! samples
1104 !
1105 ! trace : (I/O) SAC trace
1106 ! dt : (I) time to shift by
1107 !
1108 implicit none
1109 integer :: isamp,ishift
1110 real :: dt
1111 type (SACtrace) :: trace
1112
1113 ishift = nint (dt / (trace % delta) )
1114
1115 ! ** check for no shift
1116 if (ishift == 0 .and. f90sac_suppress_warnings /= 1) then
1117 write(0,'(a)') &
1118 'F90SAC_TSHIFT: Warning: no shift applied, dt too small'
1119 endif
1120
1121 ! ** shift array
1122 trace % trace = cshift((trace % trace),-ishift)
1123
1124 ! ** if negative shift, zero last ishift points
1125 if (ishift < 0) &
1126 trace % trace(trace % npts - abs(ishift): trace % npts) = 0.0
1127 ! ** if positive shift, zero first ishift points
1128 if (ishift > 0) trace % trace(1:abs(ishift)) = 0.0
1129
1130 return
1131 end subroutine f90sac_tshiftSo trace % trace(trace % npts - abs(ishift): trace % npts) = 0.0 is the line where the program crashes.
When the program crashes, ishift is equal to -2147483648 and trace%npts is 1001, so this line is trying to access a negative index of the array trace%trace.
ishift is assigned on line 1113: ishift = nint (dt / (trace % delta) ). At the time of the crash, dt is -inf and trace%delta is 0.050000001.
I'm assuming that dt should not be -inf (it seems it can sometimes be a very high value like -5.0674734e+24, but the outcome is still the same) and this is causing ishift to be the minimum value a 4 byte integer can be.
The backtrace shows that the subroutine f90sac_tshift() is called on line 138 of Analyser.F90:
108!=============================================================================
109 subroutineAn_ApplySplittingOperators(Path)
110!=============================================================================
111 !
112 ! Apply splitting operatorsto the triplet of files t1,t2,t3.It is assumed
113 ! that t1 and t2 are theS-wave components
114 !
115 use f90sac
116 use MTSGlobal
117 implicit none
118
119 ! ** input parameters
120 type(MPath) :: Path
121
122 ! ** locals
123 real(r4) :: rotangtsplit
124 integer :: iop
125
126 ! ** loop over the domainson the path, apply inverseoperators to the data
127 ! this is done in reverseorder, to preserve operatorcommutation
128 do iop=Path % nop,1,-1
129
130 ! Path % op(iop) %fast, Path % op(iop) % tlag
131 rotang = Path % o(iop) % fast !+ t1 % baz forgeographical
132 tsplit = -Path % o(iop) % tlag
133
134 ! ** rotate the traces
135 call f90sac_rotate2(t1,t2,rotang, iforce=1)
136
137 ! ** time shift the slowtrace
138 call f90sac_tshif(t2,tsplit)
139
140 ! ** rotate back tooriginal orientation
141 call f90sac_rotate2(t1,t2,-rotang, iforce=1)
142 enddo
143
144 return
145 end subroutineAn_ApplySplittingOperators
146!=============================================================================So tsplit is the value that is passed in to dt. tsplit is assigned on line 132: tsplit = -Path % op(iop) % tlag.
At this point, iop is 1 and Path%op(iop)%tlag is already inf.
Path%op(iop)%tlag is assigned by a call to the subroutine DomainOperator() in EvaluatePath.f90 on line 59:
59 call DomainOperator(Domain, &
60 Path % op(iop) % azi, &
61 Path % op(iop) % inc, &
62 Path % op(iop) % dist, &
63 Path % op(iop) % fast, &
64 Path % op(iop) % tlag)The subroutine DomainOperator() is defined on line 104. Inside this subroutine, on line 140, the tlag variable is assigned a value: tlag=vs1*dist-vs2*dist
vs1, vs2, and dist are assigned on line 133 by a call to the subroutine CIJ_phasevels():
132 ! ** calculate phase vels
133 call CIJ_phasevels(Domain%c,Domain%density/1.e3,azi,inc, &
134 pol=phi,vs1=vs1,vs2=vs2)CIJ_phasevels() is defined on line 121 in EmatrixUtils.f90.
The vs1 and vs2 values are assigned by a call to the subroutine velo() defined on line 1028 in EmatrixUtils.f90. It seems that vs1 an vs2 are being assigned very high values somehow. An example of some values from one crash:
vs1 = 4.1477941969817330E+096
vs2 = 3.8020460548218681E+089
dist = 100.00000000000000
tlag = 4.1477938167771275E+098
dt = -Infinity
ishift = -2147483648
It seems that the velo() subroutine is not the problem, as the C matrix that is passed as an argument to the subroutine seems to be what is causing the vs1 and vs2 values to be so large.
On a successful run, the C matrix contains values in the range of 106 to 10-5. Whereas, when the program crashes, these values are in the range of 1070 to 1090.
The C matrix that is used in CIJ_phasevels() is passed in as an argument and is then converted to "Mbar" units and normalised. The C matrix that is passed in to CIJ_phasevels() is Domain%c.
Here is where I'm a little stuck, why is Domain%c wrong at this point, and where does it come from. On line 57 of EvaluatePath.f90, it can be seen that Domain is taken from an array of domains that the Model contains, but I'm not sure where and how the domain array is populated in the first place.