#52 RES may upset following process modules

[seismick]

If RES (resampling) is followed by FKTR, the job aborts with SIGSEGV.
That is with i686_gfortran_mpich2 build. With x86_64_gfortran_mpich2, there are
4 processing modules that abort after RES.
C4WE: array bound mismatch 'tr_loc'
FKTR: array bound mismatch 'obj'
FXDECON: array ref out of bounds 'cdata1'
STK: array ref out of bounds 'hd'
These are all multi-trace processes. I have not found any single-trace processes
that abort after RES.
If CPSeis is built with Intel compiler, then RES does not seem to cause any
problems.

A workaround would be to follow RES with TROT, then read the resampled traces in
another job.

[seismick]

I looked at RES again after fixing the Bracewell interpolation it uses
(see bug number 92). I was not able to reproduce the problem with
FXDECON after RES. The other 3 processes (C4WE, FKTR and STK) do crash
following RES with mode DEC (decimate).
In each case, the offending line of code is simply copying the trace
array without specifying the limits, implying that the whole array is
used. It seems that the trace array flowing down the pipeline is still
the original length, although the number of samples is halved. With Intel
compiler, the array copy uses the smaller dimension, and carries on, but
gfortran takes the size from the right-hand side, so array bounds mismatch
occurs.

One could patch all these processes by going through the code and adding
"1:obj%ndpt" to all the trace array references. Better to fix the upstream
source of the problem.

I think blind Freddy could see that RES, the resampling process needs
twosets. Then the array length will be consistent with NDPT.
This is not hard. Add a pc_put_control to res_update, add extra
arrays tro and hdo to the execute subroutine, and copy the existing
results into those. And then modify the res_superwrap.f90 to call twosets
instead of onesets.

However, there is one complication. The Bracewell interpolation needs 2
extra samples and linear interpolation 1 extra sample. The original res.f90
zeroed those extra samples, with the trace array being longer than NDPT.

    else if (obj%method_int == 'BRA') then
      tr(obj%ndpt_in+1:obj%ndpt_in+2,IT) = 0.
      call interp_1d_con_bw_real (obj%fctr_res, obj%ndpt_in, 1,   &
                   obj%ft, tr(:,IT), obj%trtmp)
      tr(:obj%ndpt_out,IT) = obj%trtmp(:obj%ndpt_out)
    else
      tr(obj%ndpt_in+1,IT) = 0.
      call interp_1d_con_lin_real (obj%fctr_res, obj%ndpt_in, 1,  &
                   obj%ft, obj%gt, tr(:,IT), obj%trtmp)
      tr(:obj%ndpt_out,IT) = obj%trtmp(:obj%ndpt_out)
    end if

The answer is to use the tro array temporarily:

    else if (obj%method_int == 'BRA') then

! 2019 APR the following lines for twosets
tro(:obj%ndpt_in,IT) = tr(:obj%ndpt_in,IT)
tro(obj%ndpt_in+1:obj%ndpt_in+2,IT) = 0.0
call interp_1d_con_bw_real (obj%fctr_res, obj%ndpt_in, 1, &
obj%ft, tro(:,IT), obj%trtmp)
tro(:obj%ndpt_out,IT) = obj%trtmp(:obj%ndpt_out) ! 2019 APR
else
! 2019 APR the following lines for twosets
tro(:obj%ndpt_in,IT) = tr(:obj%ndpt_in,IT)
tro(obj%ndpt_in+1,IT) = 0.
call interp_1d_con_lin_real (obj%fctr_res, obj%ndpt_in, 1, &
obj%ft, obj%gt, tro(:,IT), obj%trtmp)
tro(:obj%ndpt_out,IT) = obj%trtmp(:obj%ndpt_out) ! 2019 APR
end if

The attached file is updated for twosets. I have tested it with a variety of
processes. Note that C4WE becomes extremely slow after interpolation, so I
doubt that would be used in production.

[seismick]

sorry forgot the attachement

[Bill Menger]

Accepted seismick's input, created a new res_superwrap to incorporate the twosets of arguments, and tested in both interpolate and decimate modes.