      PROGRAM RESPON
C     
C..   A PROGRAM TO FORMULATE COMPLEX TRANSFER FUNCTION MATRIX H IN
C     EQUATION:  X(OMEGA) = H(OMEGA)*F(OMEGA)
C     AND TO CALCULATE THE MAXIMUM X-MAGNITUDE PER GIVEN EXCITATION
C     FREQUENCY UNDER UNIT LOADING
C
C     Original version by Sami Saarinen, CSC/1988
C
C     Version 2 for OPMVAX Replacement Project 1991 by Sami Saarinen
C
C     - I/O to 20-unit removed
C     - source code reformatted using emacs 'indent-region'
C     - modified to SINGLE PRECISION (32-bit)
C     
CDEP  
      IMPLICIT REAL (A-H,O-Z)
      parameter (nproc=4)
      PARAMETER (NDIM=882,NMODES=10,MAXSAM=200)
      DIMENSION FREQR(NMODES),XSIIR(NMODES),TABLE(MAXSAM,2)
c-old      DIMENSION FREQS(MAXSAM),XRESP(NDIM)
      DIMENSION FREQS(MAXSAM)
      dimension jfreq(nproc+1),idof(maxsam)
C-old:DEP  --- NOTE: ALL COMPLEX*16 MEANS A 64-BIT REAL AND IMAGINARY PART
C-old      COMPLEX*16 H,XRESP
      COMPLEX H,XRESP
c-old      COMMON /FSCC1/ H(NDIM,NDIM)
      COMMON /FSCC1/ H(NDIM,NDIM,nproc)
      COMMON /FSCC2/ FII(NDIM,NMODES)
      common /fscc3/ xresp(ndim,maxsam)
c-old      CHARACTER*40 FILNAM
C     
      call mesage('** JOB STARTED **   ',1,1)
C
      write(6,*)'Enter number of processors to be used:'
      read(5,*) kproc
      kproc = max(1,min(kproc,nproc))
      write(6,*) 'Actual number of processors equals to ',kproc
      WRITE(6,*)'ENTER NUMBER OF EXCITATION FREQUENCY POINTS:'
      READ(5,*) NFREQ
      NFREQ=MIN(NFREQ,MAXSAM)
      WRITE(6,*) 'NFREQ=',NFREQ
      WRITE(6,*)'ENTER LIST OF FREQUENCY VALUES:'
      READ(5,*) (FREQS(I),I=1,NFREQ)                      
      WRITE(6,*)'** LIST OF FREQUENCIES TO BE USED:'
      WRITE(6,FMT='(1P,(1X,6G13.5))') (FREQS(I),I=1,NFREQ)
      WRITE(6,*)'ENTER SAMPLE DEGREE OF FREEDOM:'
      READ(5,*) KDOF
      WRITE(6,*)'KDOF=',KDOF
C     WRITE(6,*)'ENTER FILENAME FOR MODAL DATA INPUT:'
C     READ(5,'(A)') FILNAM
C     WRITE(6,*)'FILENAME='//FILNAM
CDEP  
      OPEN(UNIT=1,file='fort.1',STATUS='OLD')
      READ(1,*) KDIM,KMODES
      KDIM=MIN(KDIM,NDIM)
      KMODES=MIN(KMODES,NMODES)
      WRITE(6,*)'KDIM=',KDIM,' KMODES=',KMODES
      call mesage('**READ MODES BEGIN  ',2,1)
      CALL READI(FII,FREQR,XSIIR,KDIM,KMODES)
      call mesage('**READ MODES DONE   ',2,0)
      CLOSE(1)
C     
      call mesage('**SCALING MODES     ',4,1)
      CALL SCALE(FII,KDIM,KMODES)
      call mesage('**SCALING DONE      ',4,0)
C     
C..   FORMULATE H-MATRIX IN EVERY FREQUENCY POINTS AND CALCULATE RESPONSE
C     
CDEP  
C-old      OPEN(20,STATUS='UNKNOWN')
      call mesage('CALCULATING RESPONSE',3,1)          
      PII180=57.29577951

      jinc = (nfreq+1) / kproc
      jfreq(1) = 1
      do 120 ii=2,nproc
         jfreq(ii) = jfreq(ii-1) + jinc
 120  continue
      jfreq(kproc+1) = nfreq

c---  candidate for 'do parallel' by using nproc's

      do 110 ii=1,kproc

c-old         DO 100 IFREQ=1,NFREQ
         DO 100 IFREQ=jfreq(ii),jfreq(ii+1)

            CALL HMATRX(
     +           FREQS(IFREQ),FREQR,XSIIR,FII,H(1,1,ii),KDIM,KMODES)
            CALL FINDMX(H(1,1,ii),KDIM,XRESP(1,ifreq),FREQS(IFREQ),
     +           IDOF(ifreq))

 100     CONTINUE
 110  continue

      DO 200 IFREQ=1,NFREQ
         FREQ=FREQS(IFREQ)
         RREAL=REAL(XRESP(IDOF(ifreq),ifreq))
C-old RIMAG=DIMAG(XRESP(IDOF))
C-old RMAGN=DSQRT(RREAL**2+RIMAG**2)
C-old ANGLE=DATAN2(RREAL,RIMAG)*PII180               
         RIMAG=AIMAG(XRESP(IDOF(ifreq),ifreq))
         RMAGN=SQRT(RREAL**2+RIMAG**2)
         ANGLE=ATAN2(RREAL,RIMAG)*PII180               
         WRITE(6,1000) IFREQ,FREQ,RMAGN,ANGLE,IDOF(ifreq)
C-old WRITE(20,1000) IFREQ,FREQ,RMAGN,ANGLE,IDOF(ifreq)
 1000    FORMAT(5X,I5,': ',1P,3(G13.5,2X),' AT DOF ',I5)
         RREAL=REAL(XRESP(KDOF,ifreq))
C-old RIMAG=DIMAG(XRESP(KDOF))
         RIMAG=AIMAG(XRESP(KDOF,ifreq))
C-old TABLE(IFREQ,1)=DSQRT(RREAL**2+RIMAG**2)
C-old TABLE(IFREQ,2)=DATAN2(RREAL,RIMAG)*PII180               
         TABLE(IFREQ,1)=SQRT(RREAL**2+RIMAG**2)
         TABLE(IFREQ,2)=ATAN2(RREAL,RIMAG)*PII180               
C-old call mesage('**INTERMEDIATE TIME ',3,0)
 200  continue

      call mesage('CALCULATION DONE    ',3,0)
C     
      RMAX=0.0
      DO 300 I=1,NFREQ
         RMAX=MAX(RMAX,TABLE(I,1))
 300  CONTINUE
      WRITE(6,1101) 
     +     KDOF,'#','FREQ [Hz]',
     +     'Magnitude','Phase Angle','   Plot Bar'

      DO 310 I=1,NFREQ
         NUMSTR=INT(30.0*TABLE(I,1)/RMAX+0.5)
         WRITE(6,1100) I, FREQS(I),TABLE(I,1),TABLE(I,2),
     +        ('*',J=1,NUMSTR)
 310  CONTINUE         
 1100 FORMAT(1X,I5,F12.2,2X,E13.5,2X,F12.2,2X,40A1)
 1101 FORMAT(//1X,'*** RESULTS AT SAMPLE DEGREE OF FREEDOM ',I5,' :'//
     +       1X,A5,A12,2X,A13,2X,A12,2X,A/)
C     
      call mesage('** END OF JOB **    ',1,0)
      STOP 
      END
