At this stage the ideal form of the calculation would be to run the simulation
saving the state vector at each step, then analysing the data with a
post-processing program. Unfortunately, the run lengths involved in this
project prohibited this method of operation, as gigabytes of data are produced.
The only alternative was to analyse the data on the fly. The subroutine
SPM_quotOUTPUT" generates a sum of transient responses, and stores the results
in a buffer, which is defined in the following include file:
DOUBLE PRECISION SLTL0J(NGPTS), SUMJ(NGPTS)
DOUBLE PRECISION SUMLT(NGPTS), SLTL0(NGPTS)
DOUBLE PRECISION LTL0J2(NGPTS)
COMMON /BUFFER/SUMLT,SUMJ,SLTL0,SLTL0J,LTL0J2
Output maintains a shift register, in which 
, 
and
are cycled through. The transient response function is then
generated as a function of the time difference between the current value of
and the previous remembered values.
SUBROUTINE OUTPUT
INCLUDE (PARAM)
INTEGER STEP
DOUBLE PRECISION ALPHA, LAMBDA, CALJ
INCLUDE (BUFFER)
DOUBLE PRECISION XZ(D*NP), PZ(D*NP), P
COMMON /GAMMAG/ XZ,PZ,ALPHA,LAMBDA
DOUBLE PRECISION SHIFT(NGPTS), SHIFTJ(NGPTS), SHFTJ2(NGPTS)
COMMON /SHIFTR/ SHIFT, SHIFTJ, SHFTJ2
DOUBLE PRECISION PRESS
CALL SHFTIN(LAMBDA, CALJ(PZ))
P = PRESS()
DO 1 STEP=1,NGPTS
SLTL0(STEP) = SLTL0(STEP) + LAMBDA * SHIFT(STEP)
SLTL0J(STEP) = SLTL0J(STEP) + P * SHIFTJ(STEP)
LTL0J2(STEP) = LTL0J2(STEP) + LAMBDA * SHFTJ2(STEP)
1 CONTINUE
RETURN
END
SUBROUTINE SHFTIN(L,J)
INCLUDE (PARAM)
INTEGER I
DOUBLE PRECISION L,J,SHIFT(NGPTS),SHIFTJ(NGPTS), SHFTJ2(NGPTS)
COMMON /SHIFTR/ SHIFT, SHIFTJ, SHFTJ2
DO 1 I=NGPTS-1,1,-1
SHIFT(I+1)=SHIFT(I)
SHIFTJ(I+1)=SHIFTJ(I)
SHFTJ2(I+1)=SHFTJ2(I)
1 CONTINUE
SHIFT(1)=L
SHIFTJ(1)=L*J
SHFTJ2(1)=L*J**2
RETURN
END
BLOCK DATA BUFNIT
INCLUDE(PARAM)
DOUBLE PRECISION SUMLT(NGPTS), SLTL0(NGPTS)
DOUBLE PRECISION SLTL0J(NGPTS), SUMJ(NGPTS)
DOUBLE PRECISION LTL0J2(NGPTS)
COMMON /BUFFER/SUMLT,SUMJ,SLTL0,SLTL0J,LTL0J2
DATA SLTL0/NGPTS*0.0/
DATA SUMLT/NGPTS*0.0/
DATA SUMJ/NGPTS*0.0/
DATA SLTL0J/NGPTS*0.0/
DATA LTL0J2/NGPTS*0.0/
END