
C	to compute dispersion curve using dispcurve.f

C	Sven Ivansson 1993




	PROGRAM main   !!!  RMODFGR

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C	nol <= noldim
	INTEGER NFRDIM
	PARAMETER ( NFRDIM = 1000 )
C       nfr <= nfrdim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'DPCFSUB.INC'
	INCLUDE 'SOURCFW.INC'  ! for convenience only

	INTEGER MHQEFF, JF,NFR,NFROK, IUROT0, IGRCALL, IASFAIL
	REAL*8 HSTEPFSUBPRE, FREQMIN,FREQINT, ALADAUER,PREALA,ALA, 
     $	DSFMIN,DSFMAX, DSFFAKTMAX, DSUMAX, DRELGOALFAKT, RSUG2, 
     $	REDUCFAKT,DLOGREDUCFAKT,DDELT,XYTOL, PI, SLASK, 
     $  OFREQ(nfrdim)
	COMPLEX*16 UPSTART,USTART, U0LOKINV, 
     $	CFREQARR(nfrdim), UCURVE(nfrdim), DUDFCURVE(nfrdim), 
     $	UGRCURVE(nfrdim)
	CHARACTER*32 FILM, FILFU

	DATA PI /3.141592653589793D0/

	U0 = (1.0D0,0.0D0)   ! internal(!!!) computations without 'u rotation'
 	UPREF = (0.0D0,0.0D0)

	WRITE(*,*)
	WRITE(*,*) 'Give infile FILM:'
	FILM = ' '
	DO WHILE (FILM.EQ.' ') 
	  READ(*,'(A)') FILM
	END DO
	CALL BLANKTOEND (FILM)
	IAZIMI = 0
	WRITE(*,*) 'Give MHQEFF:'
	READ(*,*) MHQEFF
	WRITE(*,*) 'Give ICUTTYP ( 1="hyperbolic for up" 2="vertical for up" '
	WRITE(*,*) '3="mixture with aimcut" ) :'
	WRITE(*,*) '( "-5 or -15 or -45 or -60" for Muir )'
	WRITE(*,*) '( irrelevant for a rigid or free bottom )'
	WRITE(*,*) '( rk=omeg*u )'
	READ(*,*) ICUTTYP
	IF (ICUTTYP.EQ.3) THEN
	  WRITE(*,*) 'Give AIMCUT:'
	  READ(*,*) AIMCUT
	END IF
	IF ( ICUTTYP.NE.1 .AND. ICUTTYP.NE.2 .AND. ICUTTYP.NE.3 .AND. 
     $	ICUTTYP.NE.-5 .AND. ICUTTYP.NE.-15 .AND. 
     $	ICUTTYP.NE.-45 .AND. ICUTTYP.NE.-60 ) STOP 340
	WRITE(*,*) 'Give ISHEETTYP ( 1="++" 2="-+" 3="+-" 4="--" ) for up" :'
	WRITE(*,*) '( irrelevant for a rigid or free bottom )'
	WRITE(*,*) '( "1 or 3" "2 or 4" irrelevant for a fluid half-space )'
	READ(*,*) ISHEETTYP
	WRITE(*,*) 'Give EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,',
     $	'EPSMAGNIFMAX :'
	WRITE(*,*) '( standard is  epsfakt=0.1, hfaktmax=1.05, hoff=6d-16 )'
	READ(*,*) EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,EPSMAGNIFMAX
	WRITE(*,*) 'Give EPSFSUB,HSTEPFSUBPRE :'
	WRITE(*,*) '( hstepfsub will be "hstepfsubpre/freq" , which '
	WRITE(*,*) '  is a "compromise" regarding rkra & rkrb )'
	READ(*,*) EPSFSUB,HSTEPFSUBPRE

	WRITE(*,*)
	WRITE(*,*) 'Give FREQMIN(>0.0),FREQINT (both in Hz) , NFR :'
	WRITE(*,*) '( " freqmin<0.0 & freqint=0.0 "  for  ',
     $	'"explicit ofreq(1:nfr) " )'
	WRITE(*,*) '( "freq in mainmodu" if nfr<0 )'
	READ(*,*) FREQMIN,FREQINT, NFR
	IF ( FREQMIN.EQ.0.0D0 .OR. NFR.GT.nfrdim ) STOP 6666
	IF (FREQMIN.GE.0.0D0) THEN
	  SLASK = FREQINT/NFR
	  DO JF = 1, NFR
	    OFREQ(JF) = FREQMIN + (JF-1)*SLASK
	  END DO
	ELSE
	  IF (FREQINT.NE.0.0D0) STOP 6660
	  WRITE(*,*) 'give ofreq(1:nfr) increasing'
	  READ(*,*) (OFREQ(JF),JF=1,NFR)
	END IF
	WRITE(*,*) 'Give ALADAUER (s), PREALA :'
	WRITE(*,*) '( 0.0<preala<=1.0, preala=1.0 for real "omeg & rk" )'
	WRITE(*,*) '( "preala" is "temporary reduct factor during aladauer" )'
	WRITE(*,*) '( "omeg" will be " 2*pi*freq + i*ala " '
	WRITE(*,*) '   where "ala" = "-ln(preala)/aladauer" )'
	READ(*,*) ALADAUER,PREALA
	IF (PREALA.EQ.1.0D0) THEN
	  ALA = 0.0D0
	ELSE
	  ALA = -DLOG(PREALA)/ALADAUER
	END IF
	DO JF = 1, NFR
	  CFREQARR(JF) = DCMPLX ( OFREQ(JF), ALA/(2*PI) )
	END DO
	WRITE(*,*) 'Give IUROT0 ( 0 for "u as specified", 1 for ',
     $	'"u rotation with -arg(omeg)" :'
	WRITE(*,*) '( u = "up rotated by -arg(omeg)" )'
	WRITE(*,*) '( irrelevant if "preala"=1.0 )'
	READ(*,*) IUROT0

	WRITE(*,*)
	WRITE(*,*) 'Give UPSTART (mode-up at cfreqarr(1), not on a ',
     $	'branch-cut) :'
	READ(*,*) UPSTART
	IF (IUROT0.EQ.0) THEN
	  USTART = UPSTART
	ELSE
	  USTART = (CDABS(CFREQARR(1))/CFREQARR(1))  * UPSTART
	END IF

	WRITE(*,*)
	WRITE(*,*) 'Give outfile name FILFU :'
	FILFU = ' '
	DO WHILE (FILFU.EQ.' ') 
	  READ(*,'(A)') FILFU
	END DO
	CALL BLANKTOEND (FILFU)

	WRITE(*,*)
	WRITE(*,*) 'Give DSFMIN (non-negative), DSFMAX, DSFFAKTMAX :'
	READ(*,*) DSFMIN,DSFMAX,DSFFAKTMAX
	WRITE(*,*) 'Give DSUMAX, DRELGOALFAKT (the latter less than one) :'
	READ(*,*) DSUMAX,DRELGOALFAKT
	WRITE(*,*) 'Give RSUG2 (max u search-rectangle half-side) :'
	READ(*,*) RSUG2
	WRITE(*,*) 'Give REDUCFAKT,DDELT,XYTOL (for dzeroholo) :'
	READ(*,*) REDUCFAKT,DDELT,XYTOL
	DLOGREDUCFAKT = DLOG(REDUCFAKT)

	OPEN(UNIT=20,FILE=FILM,STATUS='OLD')
	CALL RMODEL ( 20,0, MHQEFF )
	NSOURCFW = 0  ! irrelevant here
	CLOSE(20)
	IF (RHO(NOL).GT.0.0D0) THEN
	  UPBPA = CDSQRT(UA2(NOL))      ! u0=(1,0), freq-indep here
	  UPBPA2 = UPBPA * UPBPA
	  IF (ISOLID(NOL).NE.0) THEN
	    UPBPB = CDSQRT(UB2(NOL))    ! u0=(1,0), freq-indep here
	    UPBPB2 = UPBPB * UPBPB
	  ELSE
	    UPBPB = (0.0D0,0.0D0)
	    UPBPB2 = (0.0D0,0.0D0)
	  END IF
	ELSE
	  UPBPA = (0.0D0,0.0D0)
	  UPBPA2 = (0.0D0,0.0D0)
	  UPBPB = (0.0D0,0.0D0)
	  UPBPB2 = (0.0D0,0.0D0)
	END IF
	IF (RHO(0).GT.0.0D0) THEN
	  UPBPASURF = CDSQRT(UA2(0))    ! u0=(1,0), freq-indep here
	  UPBPA2SURF = UPBPASURF * UPBPASURF
	  IF (ISOLID(0).NE.0) THEN
	    UPBPBSURF = CDSQRT(UB2(0))  ! u0=(1,0), freq-indep here
	    UPBPB2SURF = UPBPBSURF * UPBPBSURF
	  ELSE
	    UPBPBSURF = (0.0D0,0.0D0)
	    UPBPB2SURF = (0.0D0,0.0D0)
	  END IF
	ELSE
	  UPBPASURF = (0.0D0,0.0D0)
	  UPBPA2SURF = (0.0D0,0.0D0)
	  UPBPBSURF = (0.0D0,0.0D0)
	  UPBPB2SURF = (0.0D0,0.0D0)
	END IF

	IGRCALL = 0
	CALL DISPCURVE ( HSTEPFSUBPRE, NFR,CFREQARR, USTART, 
     $	DSFMIN,DSFMAX,  DSFFAKTMAX, DSUMAX, DRELGOALFAKT, RSUG2, 
     $	DLOGREDUCFAKT,DDELT,XYTOL, NFROK,UCURVE,DUDFCURVE,UGRCURVE, 
     $	IGRCALL,IASFAIL, .TRUE.)
	IF (IASFAIL.NE.0) THEN
	  WRITE(*,*) ' %%% iasfail =',IASFAIL,' from dispcurve'
	  WRITE(*,*) '     nfrok =',NFROK
	END IF

	OPEN(UNIT=1,FILE=FILFU,STATUS='NEW')
	WRITE(1,*) 12,'  ! lines to skip, from rmodfgr'
	WRITE(1,*) ' FILM,MHQEFF: ',FILM,' ',MHQEFF
	IF (ICUTTYP.NE.3) THEN
	  WRITE(1,*) ' ICUTTYP: ',ICUTTYP
	ELSE
	  WRITE(1,*) ' ICUTTYP,AIMCUT: ',ICUTTYP,SNGL(AIMCUT)
	END IF
	WRITE(1,*) ' ISHEETTYP: ',ISHEETTYP
	WRITE(1,*) ' EPSGLOBLN,EPSGLOB: ',SNGL(EPSGLOBLN),SNGL(EPSGLOB)
	WRITE(1,*) ' EPSFAKT,HFAKTMAX: ',SNGL(EPSFAKT),SNGL(HFAKTMAX)
	WRITE(1,*) ' HOFF,EPSMAGNIFMAX: ',SNGL(HOFF),SNGL(EPSMAGNIFMAX)
	WRITE(1,*) ' EPSFSUB,HSTEPFSUBPRE: ',SNGL(EPSFSUB),SNGL(HSTEPFSUBPRE)
	WRITE(1,*) ' UPSTART: ',UPSTART
	WRITE(1,*) ' DSFMIN,DSFMAX, DSFFAKTMAX :', SNGL(DSFMIN),
     $	SNGL(DSFMAX),SNGL(DSFFAKTMAX)
	WRITE(1,*) ' DSUMAX,DRELGOALFAKT :', SNGL(DSUMAX),SNGL(DRELGOALFAKT)
	WRITE(1,*) ' RSUG2 :', SNGL(RSUG2)
	WRITE(1,*) ' REDUCFAKT,DDELT,XYTOL :', 
     $	SNGL(REDUCFAKT),SNGL(DDELT),SNGL(XYTOL)

	WRITE(1,*)
	WRITE(1,*) 3, ICUTTYP,ISHEETTYP, ALA,IUROT0

	WRITE(1,*)
	WRITE(1,*) NFROK,' ! nfrok (nfr=',NFR,' ) "freq,up"'
	IF (IUROT0.EQ.0) THEN
	  DO JF = 1, NFROK
	     WRITE(1,*) SNGL(OFREQ(JF)), UCURVE(JF)          ! freq,up
	  END DO
	ELSE
	  DO JF = 1, NFROK
	    U0LOKINV = CFREQARR(JF)/CDABS(CFREQARR(JF))
	    WRITE(1,*) SNGL(OFREQ(JF)), U0LOKINV*UCURVE(JF)  ! freq,up
	  END DO
	END IF
	WRITE(1,*)
	WRITE(1,*) NFROK,' ! nfrok (nfr=',NFR,' ) "freq,du/df"'
	DO JF = 1, NFROK
	  WRITE(1,*) SNGL(OFREQ(JF)), DUDFCURVE(JF)          ! freq,du/df
	END DO
	WRITE(1,*)
	WRITE(1,*) NFROK,' ! nfrok (nfr=',NFR,' ) "freq,group-u"'
	DO JF = 1, NFROK
	  WRITE(1,*) SNGL(OFREQ(JF)), UGRCURVE(JF)           ! freq,gr-u
	END DO
	WRITE(1,*)
	WRITE(1,*) ' IGRCALL was ',IGRCALL
	IF (IASFAIL.NE.0) WRITE(1,*) ' %%% iasfail =',IASFAIL

	STOP
	END

