text:
	subroutine makefile(ch)
c
C	A PROGRAM TO SEARCH FOR LION ROARS IN FFTL
c
C	WRITES OUTPUT TO FFTL.REC,FOR037,FOR080,FOR090,WHISTLER.POWER
C		ONLY WHISTLER.POWER IS ACCESS=APPEND
C
	integer*4 ch,ok,okt,SCETI4(2),NDATA(2048,4),SCETST(2)
	INTEGER*4 W_CHANNEL_CLOSE,W_EVENT,RET_SIZE,W_MESSAGES_OFF
	INTEGER*4 W_ITEM_I4,W_ITEM_R4,W_ITEM_R8
	INTEGER*4 FFT_CHANNEL,ISRC,SUNCLOCK
	INTEGER*4 ISOURCE(10),LSOURCE(4)
	REAL*8 SCET8,XGSE,YGSE,ZGSE,ERAD
	REAL Y(3),ERR(1024),GSEK(3)
	REAL TDATA(1026),TSPEC(1026),ZCROSS(1024),ZINT(1024),FREQZC(4)
	REAL DBSPEC(513,4),VDATA(2050,4),GDATA(2050,4),VMDATA(2050,4)
	real avspec(513),stdspec(513),spcount(513)
	REAL FREQ(513),PWR(4),SPPWR(4),FREQHZ(4),BW(4),RDATA(1024,4)
	REAL VPWR(4),WT(2050)
	REAL THRESH(513),BACKGND(513)
	complex preamp
	character*32 ITEM
	character*4 event
	COMMON /WPOWER/ TPOWER(513),TCOUNT(513)
	COMMON /VARMATX/ EVAL(3),EVECT(3,3)
C	COMMON /WEIGHT/ WT(2050)
	common /rawvolts/ DATA(1026)
	common /telem/ tmspec(1024)
	integer*4 tmspec
	DATA TWOPI /6.2831853/
	DATA FUNDFR /.33356/
	DATA ERAD /6378.D00/
	DATA PWR /4*0./
	DATA PWRLIM /.6E-5/
	DATA IBACKRD/0/
	DATA VPWRLIM /.7/
	data event /'FFTL'/
	DATA NCOUNT /0/
	DATA IPOWER /0/
	DATA TPOWER /513*0./
	DATA TCOUNT /513*1.E-8/
C
C	CALCULATE FREQUENCIES FOR DBSPEC
C
	FREQ(1) = 0.
	DO N = 2,513
	  FREQ(N) = (N-1)*FUNDFR	    ! FREQ(N) GOES WITH DBSPECT(N)
C	  TH = -29. - 19.*ALOG10(FREQ(N))
C	  TL = -2. - 52.*ALOG10(FREQ(N))
C	  THRESH(N) = AMAX1(TL,TH)
C	  BACKGND(N) = THRESH(N) - 5.
	ENDDO
C	THRESH(228) = THRESH(228) + 5.
C	THRESH(229) = THRESH(229) + 5.
	NUMEVTS = 0
C
	IF(IPOWER.EQ.0) THEN
	  IPOWER = 1
C
C	  READ IN ACCUMULATED TOTAL POWER
C
	  OPEN(UNIT=46,FILE='WHISTLER.POWER',STATUS='OLD')
	  DO N = 2,513
	    READ(46,*) FF,TCOUNT(N),TPOWER(N)
	  ENDDO
	  CLOSE(UNIT=46)
	ENDIF
C
	IF(IBACKRD.EQ.0) THEN
	  IBACKRD = 1
	  OPEN(UNIT=45,FILE='FFTL.BCKG',STATUS='OLD',READONLY)
	  DO N = 1,512
	    READ(45,*) NL,CNT,BACKGNDT,STD2T
	    BACKGND(NL) = BACKGNDT
	    STD = SQRT(STD2T)
	    THRESH(NL) = BACKGND(NL) + 1.5*STD
	  ENDDO
	  CLOSE(UNIT=45)
	  PRINT*,NL,CNT,BACKGND(NL),STD2T,THRESH(NL)
	ENDIF
	IFOUND = 0
 100	ok = w_event(ch,event)
C
C	CHECK FOR END OF RECORD
C
	if (ok.ne.1) then
		if(ok.eq.82) then
		   ok = w_channel_close(ch)
		   open(unit=45,file='fftl.rec',status='old',
     1			access='append')
		   write(45,*) sceti4,numevts
		   close(unit=45)
c
C	do n = 1,512
C	  avspec(n) = avspec(n)/spcount(n)
C	  stdspec(n) = stdspec(n)/spcount(n) - avspec(n)**2
C	  np = n
C	  write(72,*) np,spcount(n),avspec(n),stdspec(n),sqrt(stdspec(n))
C	enddo
c
		   return
	        endif
		write(6,*) 'cannot open ',event, ', ok=', ok
	endif
C
C	END OF END OF RECORD CHECK
C
	OKT = W_MESSAGES_OFF(ch)
C
C
	ITEM = 'CHANNEL_NUMBER'
	ok = W_ITEM_I4(ch, item, FFT_CHANNEL, 1, ret_size)
	IF(FFT_CHANNEL.LE.6) GO TO 100      		! NOT FFTL
C
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
C
C	ITEM = 'CAL_PA'
C	ok = W_ITEM_I4(ch, item, ICAL, 1, ret_size)
C	IF(ICAL.NE.0) THEN
C	  PRINT*,'CALIBRATION'
C	  GO TO 100
C	ENDIF
C
	ITEM = 'SOURCE'
	ok = W_ITEM_I4(ch, item, ISOURCE(FFT_CHANNEL), 1, ret_size)
C
	IF(FFT_CHANNEL.EQ.7) THEN		! START OF A NEW SET
	  ITEM = 'EVENT_SCET'
	  ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
	  IF(ISOURCE(FFT_CHANNEL).GT.6) THEN
	    LD = ISOURCE(FFT_CHANNEL) - 6
	  ELSE
	    LD = 4
	  ENDIF
C
	  BFREQ = 3.
	  TFREQ = 140.
	  CALL FFT_PHYS_BPF(CH,IRAW,TDATA,TSPEC,BFREQ,TFREQ)
	  DO N = 1,1024
	    VDATA(N,LD) = TDATA(N)
	    RDATA(N,LD) = DATA(N)
	  ENDDO
	  DO N = 1,513
	    DBSPEC(N,LD) = TSPEC(N)
	  ENDDO
C	print*,'load channel,source',fft_channel,isource(fft_channel),
C     1		'  to index',ld
C
C	PRINT 722,SCETST,FFT_CHANNEL,ISOURCE(FFT_CHANNEL),LD,IRAW
 722	FORMAT(' EVENT AT',I10,I8,'  CH',I3,'  P/A',I3,'   TO INDEX',
     1		I3,'  IRAW=',I2)
C
	  IF(IRAW.NE.1) THEN
C	    ADVANCE BY ONE HOUR AND RETURN
	    ITEM = 'EVENT_SCET_R8'
	    OK = W_ITEM_R8(ch, item, SCET8, 1, ret_size)
C	    PRINT*,'NOT RAW DATA',SCETI4
	    SCET8 = SCET8 + 1.D00/24.D00
	    CALL W_CHANNEL_POSITION(CH,SCET8)
	    GO TO 100
	  ENDIF
C
	  DO N = 8,10
		ISOURCE(N) = 0
	  ENDDO
	  DO N = 1,4
	 	LSOURCE(N) = 0
	  ENDDO
	  DO N = 1,2
		SCETST(N) = SCETI4(N)
	  ENDDO
	  LSOURCE(LD) = ISOURCE(FFT_CHANNEL)
	  IFOUND = 1
	  GO TO 100
	ENDIF			! end of setup in first channel
C
	IF(ISOURCE(FFT_CHANNEL).GT.6) THEN
	  LD = ISOURCE(FFT_CHANNEL) - 6
	ELSE
	  LD = 4
	ENDIF
	LSOURCE(LD) = ISOURCE(FFT_CHANNEL)
C
	IF(FFT_CHANNEL.GT.7.AND.IFOUND.EQ.1) THEN
C
	  CALL FFT_PHYS_BPF(CH,IRAW,TDATA,TSPEC,BFREQ,TFREQ)
	  DO N = 1,1024
	    VDATA(N,LD) = TDATA(N)
	    RDATA(N,LD) = DATA(N)
	  ENDDO
	  DO N = 1,513
	    DBSPEC(N,LD) = TSPEC(N)
	  ENDDO
C
C	print*,'load channel,source',fft_channel,isource(fft_channel),
C     1		'  to index',ld
C	  PRINT*,'LD,SAMPLE DATA',LD,VDATA(1,LD),VDATA(2,LD),VDATA(3,LD)
C
	  IF(FFT_CHANNEL.NE.10) GO TO 100
	ENDIF
	NUMEVTS = NUMEVTS+1
C
C	********	LOOK FOR EZ, BZ, TOGETHER
c	IEZBZ = 0
C	DO I = 7,10
c	  IF(ISOURCE(I).EQ.9) IEZBZ = IEZBZ + 1
c	  IF(ISOURCE(I).EQ.6) IEZBZ = IEZBZ + 1
c	ENDDO
c	IF(IEZBZ.GE.2) WRITE(77,*) 'EZBZ AT',SCETST,(isource(i),i=7,10)
C	********
C	  print*,'load check,FFT_CH,SOURCE,INDEX',I,ISOURCE(I),LD
C	ENDDO
C
C	CHECK FOR 3 B CHANNELS
C
	I3E3B = 0
	DO N = 7,10
C	  print*,'3B check, FFT_CHANNEL',n,' FROM SOURCE',ISOURCE(n)
	  IF(ISOURCE(N).GT.6) I3E3B = I3E3B + 1
	ENDDO
	IF(I3E3B.LT.3) THEN
C		PRINT*,'NOT 3 Bs MODE',SCETST
		GO TO 100
	ENDIF
C
C	CHECK THAT ALL CHANNELS ARE FULL
C
	DO N = 7,10
	  IF(ISOURCE(N).EQ.0) GO TO 100
	ENDDO
C
c	DO N = 1,1024
c	  WRITE(35,*) N,VDATA(N,1),VDATA(N,2),VDATA(N,3)
c	ENDDO
c	DO N = 1,1024
c	  WRITE(34,*) N,rDATA(N,1),rDATA(N,2),rDATA(N,3)
c	ENDDO
c	IF(1) STOP
C
C	THE FOLLOWING IS SUPPOSED TO PREVENT STARTING AT SOME
C		CHANNEL OTHER THAN THE FIRST
C
	IFOUND = 0
C
C	CALCULATE POWERS FROM SPECTRUM
C
C	DO N = 1,513
C	  WRITE(37,*) N,dbspec(N,1),dbspec(N,2),dbspec(N,3)
C	ENDDO
C
C	FIND ZERO CROSSINGS, AFTER REMOVING AVERAGE
C
	DO LC = 1,2
	  AVR = 0.
	  COUNT = 0.
	  DO IL = 1,1024
	    AVR = AVR+RDATA(IL,LC)
	    COUNT = COUNT+1.
 	  ENDDO
	  IF(COUNT.NE.0.) AVR = AVR/COUNT
	  DO IL = 1,1024
	    RDATA(IL,LC) = RDATA(IL,LC) - AVR
 	  ENDDO
	  IZCNT = 0
	  IL = 1
	  IZ = IL
C	  IF(RDATA(IL,LC).EQ.0) PRINT*,'ZERO DATA',
C     1		IL,RDATA(IL,LC),RDATA(IL+1,LC)
	  DO IL = 2,1023
	    IZ = IL
CT	    IF(RDATA(IL).EQ.0) PRINT*,'ZERO DATA',IL,RDATA(IL-1),
CT     1    RDATA(IL),RDATA(IL+1)
C
C		COUNT ONLY POS TO NEG
	    IF(RDATA(IL,LC).GT.0.AND.RDATA(IL+1,LC).LE.0) THEN
	        IZCNT = IZCNT+1
		S1 = RDATA(IL,LC)
		S2 = RDATA(IL+1,LC)
	        DELZ = .5
		IF((S1-S2).NE.0.) DELZ = S1/(S1 - S2)
		IF(DELZ.LE.1.AND.DELZ.GE.0.) THEN
			ZCROSS(IZCNT) = IL + DELZ
		ELSE
			ZCROSS(IZCNT) = IL + .5
		ENDIF
	    ENDIF
	  ENDDO
	  DO N = 1,IZCNT-1
	    ZINT(N) = ZCROSS(N+1) - ZCROSS(N)
	    IF(ZINT(N).EQ.0.) PRINT*,'ZINT=0 AT ',N
	    IF(ZINT(N).EQ.0.) ZINT(N) = 1.E-6
	  ENDDO
	  FREQZC(LC) = FUNDFR*IZCNT
C	  print*,'zero crossings found',izcnt
C	  print*,'first 5',(zcross(kk),kk=1,5)
	ENDDO
	NCOUNT = NCOUNT+1
C	WRITE(37,*) NCOUNT,FREQZC(1),FREQZC(2)
C
C	LIMLO = 9
C	LIMUP = 512
C	DO NC = 1,4
C	  SPPWR(NC) = 0.
C	  AVRF = 0.
C	  SQF = 0.
C	  DO NF = LIMLO,LIMUP
C	    VLT2 =  10.**(.1*DBSPEC(NF,NC))
C	    SPPWR(NC) = SPPWR(NC) + VLT2
C	    AVRF = AVRF + FREQ(NF)*VLT2
C	    SQF = SQF + VLT2*FREQ(NF)**2
C	  ENDDO
C	  FREQHZ(NC) = AVRF/SPPWR(NC)
C	  BW(NC) = SQF/SPPWR(NC) - FREQHZ(NC)**2
C	  BW(NC) = SQRT(AMAX1(BW(NC),0.))
C	  SPPWR(NC) = .5*SPPWR(NC)
C	ENDDO
CTEMP	IF(PWRTOT.LT.PWRLIM) GO TO 100
C
C	HANDLE BZ GLITCH
C
C
	FFUND = .333557
	ITEM = 'WIND_SPIN_RATE_R4'
	OK = W_ITEM_R4(CH,ITEM,SPINRATE,1,ISIZE)
	SAMP_PER_SEC = 1024.*FFUND
	SPSS = SAMP_PER_SEC
	DANG = SPINRATE/SAMP_PER_SEC			! IN RADIANS
	DANGD = DANG*360./TWOPI
	DANGR = DANG/TWOPI				! IN CYCLES PER SAMP.
	SAMP_PER_SPIN = TWOPI/DANG
	NSPPS = SAMP_PER_SPIN + .5			! ROUNDOFF
C
	  bzmax = 0.
	  bzmin = 0.
	  do i = 1,1024
	    if(rdata(i,3).gt.bzmax) then
		bzmax = rdata(i,3)
		ibzmax = i
	    endif
	    if(rdata(i,3).lt.bzmin) then
		bzmin = rdata(i,3)
		ibzmin = i
	    endif
	  enddo
c	the following is based on solving:
c	  bzmxang = antang + dangd*(1024-ibzmax)
c	  bzmnang = antang + dangd*(1024-ibzmin)
c	for antang, given bzmxang = 34. deg.
c
c	antangb = 79. - dangd*(1024-ibzmax)	! this is s/c angle
	antangb = 34. - dangd*(1024-ibzmax)
	if(antangb.lt.-180.) antangb = antangb + 360.
C
	  DO N = 1,1024
		WT(N) = 1.
	  ENDDO
C
C	CHECK THAT THERE WAS A BZ GLITCH IN THE DATA.  IF NOT, USE
C		EX DATA
C
	IF(BZMAX.GT..3) THEN
C
C	  ON 19960528, BZMAX WAS .391 ALL DAY EXCEPT FOR NO GLITCH
C	  FOR MAY 28 AND 29, BZMIN WAS SLIGHTLY VARIABLE, USE -.219
C
	  ANTANG = ANTANGB
	  IAANG = 0
C
C	  ELIMINATE BZ GLITCH
C
	  N1 = IBZMIN - 38
	  IF(N1.LT.1) N1 = N1 + 1024
	  N2 = IBZMAX + 44
	  IF(N2.GT.1024) N2 = N2 - 1024
C	  PRINT*,'ELIM GIVES N1,N2=',N1,N2
C
C	  NOW ELIMINATE GLITCH DATA, WHICH LIES BETWEEN N1 AND N2
C		IF N2.GT.N1, AND FROM 1 TO N2 AND N1 TO 1024 OTHERWISE
C
C	  IF(IABS(N2-N1).LT.88) THEN
	  IF(N2.GT.N1) THEN
	    DO N = N1,N2
c	      VDATA(N,3) = 0.
c	      RDATA(N,3) = 0.
	      WT(N) = 0.
	    ENDDO
	  ELSE
	    DO N = 1,N2
c	      VDATA(N,3) = 0.
c	      RDATA(N,3) = 0.
	      WT(N) = 0.
	    ENDDO
	    DO N = N1,1024
c	      VDATA(N,3) = 0.
c	      RDATA(N,3) = 0.
	      WT(N) = 0.
	    ENDDO
	  ENDIF
C
	ELSEIF(BZMIN.LT.-.219) THEN
	  antangb = 34. - dangd*(1024-IBZMIN-4)
	  if(antangb.lt.-180.) antangb = antangb + 360.
C	  print*,'b angle from min',antangb
	  IAANG = 2
	ELSE
cold	  CALL FFTLOANGLE(CH,VDATA(1,4),SUNANGLE,ANTANG)
	  CALL FFTLOANGLE(CH,RDATA(1,4),SUNANGLE,ANTANG)
	  IF(ANTANG.GT.180.) ANTANG = ANTANG - 360.
C	  PRINT*,'RETURN FROM FFTLOANGLE, SUN,ANT',SUNANGLE,ANTANG
	  IAANG = 1
	  N1 = 1
	  N2 = 1
	ENDIF
c
C	  write(65,*) 'bzmin,max',scetst,ibzmin,bzmin,ibzmax,bzmax
C	  write(65,*) 'bzmin,max angles',bzmnang,bzmxang
c        endif
C	CALCULATE POWERS FROM RAW VOLTAGE SAMPLES

c	print*,'vdata',vdata(1,1),vdata(1,2),vdata(1,3)
C
	DO NC = 1,4
	  PWR(NC) = 0.
	  AVR = 0.
C
C	  REMOVE AVERAGE
C	  DO NF = 1,1024
C	    AVR = AVR + RDATA(NF,NC)
C	  ENDDO
C	  AVR = AVR/1024.
C	  DO NF = 1,1024
C	    PWR(NC) = PWR(NC) + (RDATA(NF,NC)-AVR)**2
C	  ENDDO
C	  PWR(NC) = PWR(NC)/1024.
C
C	  REMOVE SINFIT
C	  IF(NC.NE.3) THEN
C	    CALL SINFIT(DANGR,RDATA(1,NC),ERR,Y,SUMSQ)
C	    DO NF = 1,1024
C	      PWR(NC) = PWR(NC) + ERR(NF)**2
C	      RDATA(NF,NC) = ERR(NF)
C	    ENDDO
C	  ELSE
C	    DO NF = 1,1024
C	      PWR(NC) = PWR(NC) + RDATA(NF,NC)**2
C	    ENDDO
C	  ENDIF
C	  PWR(NC) = PWR(NC)/1024.
C	      WRITE(67,1067) SCETST,NC,Y,SUMSQ
 1067	      FORMAT(I10,I8,I4,3E12.3,E11.4)
C
	ENDDO
C
C	IF(IRAW.EQ.1) THEN
		CALL FFTLOANGLE(CH,VDATA(1,4),SUNANGLE,ANTANGX)
		IF(ANTANGX.GT.180.) ANTANGX = ANTANGX - 360.
C		PRINT*,'RETURN FROM FFTLOANGLE, SUN,ANT',SUNANGLE,ANTANGX
C	ENDIF
C
c	WRITE(66,1066) SCETST,(FREQHZ(NC),NC=1,2),(SPPWR(NC),NC=1,2),
c     1	(BW(NC),NC=1,2),(FREQZC(NC),NC=1,2),BZMIN,BZMAX,ANTANGB,
c     2  ANTANGX,IAANG
C	WRITE(67,1068) DBSPEC(2,1),DBSPEC(2,2),DBSPEC(3,1),DBSPEC(3,2),
C     1	DBSPEC(4,1),DBSPEC(4,2),DBSPEC(5,1),DBSPEC(5,2),DBSPEC(6,1),
C     1  DBSPEC(6,2)
 1066   FORMAT(I10,I8,2F6.2,4E11.3,2F6.1,2E12.3,2F7.1,I2)
 1068   FORMAT(10F7.1)
C	PWRTOT = PWR(1) + PWR(2) 		! ELIM BZ GLITCH
C
C	CALCULATE POWERS FROM FREQ CORRECTED VOLTAGE SAMPLES
C
c	DO N = 1,1024
c	  WRITE(38,*) N,VDATA(N,1),VDATA(N,2),VDATA(N,3)
c	ENDDO
C
	DO NC = 1,4
	  VPWR(NC) = 0.
	  AVR = 0.
	  IF(NC.NE.3) THEN
C	    REMOVE SINFIT
CS	    CALL SINFIT(DANGR,VDATA(1,NC),ERR,Y,SUMSQ)
CSS	    DO NF = 1,1024
CS	      VPWR(NC) = VPWR(NC) + ERR(NF)**2
CS	      VDATA(NF,NC) = ERR(NF)
CS	    ENDDO
CS	    VPWR(NC) = VPWR(NC)/1024.
	  ELSE
C
C	    REMOVE AVERAGE
	    DO NF = 1,1024
	      AVR = AVR + VDATA(NF,NC)
	    ENDDO
	    AVR = AVR/1024.
	    Y(1) = AVR
	    DO NF = 1,1024
	      VPWR(NC) = VPWR(NC) + (VDATA(NF,NC)-AVR)**2
	    ENDDO
	    VPWR(NC) = VPWR(NC)/945.
	  ENDIF
C	      WRITE(68,1067) SCETST,NC,Y,SUMSQ
C
	ENDDO
c	DO N = 1,1024
c	  WRITE(39,*) N,VDATA(N,1),VDATA(N,2),VDATA(N,3)
c	ENDDO
c	IF(1) STOP
C
	PWRTOT = PWR(1) + PWR(2) 		! AVOID BZ GLITCH
C	WRITE(81,777) SCETST,PWR(1),PWR(2),PWR(3),PWRTOT
 777	FORMAT(I10,I8,2X,4E12.4)
C	IF(PWRTOT.LT.PWRLIM.AND.VPWRTOT.LT.VPWRLIM) GO TO 100
C	RECALCULATE BZ POWER
C
C	NC = 3
C	  PWR(NC) = 0.
C	  AVRF = 0.
C	  SQF = 0.
C	  DO NF = LIMLO,LIMUP
C	    VLT2 =  10.**(.1*DBSPEC(NF,NC))
C	    PWR(NC) = PWR(NC) + VLT2
C	    AVRF = AVRF + FREQ(NF)*VLT2
C	    SQF = SQF + VLT2*FREQ(NF)**2
C	  ENDDO
C	  FREQHZ(NC) = AVRF/PWR(NC)
C	  BW(NC) = SQF/PWR(NC) - FREQHZ(NC)**2
C	  BW(NC) = SQRT(AMAX1(BW(NC),0.))
C	PWRTOT = PWRTOT + PWR(NC)
C
C	TRANSFORM TO GSE SYSTEM
C
C	PRINT*,'AT XFGSE,ANTANG,ANTANGB,DANGD=',ANTANG,ANTANGB,DANGD
C
	CALL XFGSE(CH,NSYS,1024,SPINRATE,SPSS,ANTANG,VDATA,GDATA)
C	DO N = 1,1024
C	  WRITE(39,*) N,gdata(N,1),gdata(N,2),gdata(N,3)
C	ENDDO
C	IF(1) STOP
C
C	print*,' VDATa1', VDATa(1,1), VDATa(2,1), VDATa(3,1)
C	print*,' VDATa4', VDATa(1,4), VDATa(2,4), VDATa(3,4)
C	print*,'GDATa1',GData(1,1),GData(2,1),GData(3,1)
C	print*,'GData2',GData(1,2),GData(2,2),GData(3,2)
C
C	CALCULATE VARIANCE MATRIX, EIGENVALUES, AND EIGENVECTORS
C
	NPT = 1024
C
C	CALL TO VMATRIX TO FIND THIRD EIGENVECTOR, FOR DIRECTION WRT B
C
	CALL VMATRIX(GDATA,WT,NPT)
C	print*,'evals',eval
	DO N = 1,3
	  GSEK(N) = EVECT(N,3)
	ENDDO
C
C	TRANSFORM TO VARIANCE MATRIX EIGENSYSTEM
C
	CALL XFORMVM(CH,NPT,NSYS,GDATA,VMDATA,WT)
c	DO N = 1,1024
c	  WRITE(37,*) N,VMDATA(N,1),VMDATA(N,2),VMDATA(N,3)
c	ENDDO
c	if(1) stop
C	print*,' gdata1', gdata(1,1), gdata(2,1), gdata(3,1)
C	print*,'VMdata1',VMdata(1,1),VMdata(2,1),VMdata(3,1)
C	print*,'xfdata2',xfdata(1,2),xfdata(2,2),xfdata(3,2)
C
C	IF(1) GO TO 100
C
C	TEST, IS IT A KEEPER?
C
C	TEST FOR SEVERAL CONSEQUETIVE LARGE SIGNALS
C
	NBIGXSIG = 0
	NBIGYSIG = 0
	DO NF = 9,512
	  IF(DBSPEC(NF,1).GT.THRESH(NF)) THEN
	    NBIGXSIG = NBIGXSIG + 1
	    NFXB = NF
	print*,'big sig,nf',nbigxsig,nf,thresh(nf),dbspec(nf,1)
	    IF(NBIGXSIG.GE.3) GO TO 205
	  ELSE
	    NBIGXSIG = 0
	  ENDIF
 205	  CONTINUE
	  IF(DBSPEC(NF,2).GT.THRESH(NF)) THEN
	    NBIGYSIG = NBIGYSIG + 1
	    NFYB = NF
	    IF(NBIGYSIG.GE.3) GO TO 200
	  ELSE
	    NBIGYSIG = 0
	  ENDIF
	ENDDO
 200	ITP = 0
	IF(NBIGXSIG.GE.3.OR.NBIGYSIG.GE.3) then
		ITP = 1
C
C	TAKE THE WIDEST RANGE
C
	NFLOW = MIN0(NFXB,NFYB) - 2
	IF(NBIGXSIG.LT.3) NFLOW = NFYB - 2
	NFLOW = MAX0(NFLOW,2)
C
C	else
C	  do n=1,512
C	    avspec(n) = avspec(n) + dbspec(n,1) + dbspec(nf,2)
C	    stdspec(n) = stdspec(n)+dbspec(n,1)**2+dbspec(nf,2)**2
C	    spcount(n) = spcount(n) + 2.
C	  enddo
	endif
C
c	print*,'itp,x,y',itp,nbigxsig,nbigysig
	WTCOUNT = 0.
	DO N = 1,1024
	  WTCOUNT = WTCOUNT + WT(N)
	ENDDO
C
c	print*,'wtcount',n1,n2,wtcount

	RATIO3RD = EVAL(3)/(EVAL(1)+EVAL(2))
	write(37,*) sceti4,ratio3rd
	ITR = 0
	IF(RATIO3RD.LT..15) ITR = 1
C
	print*,'ratio',ratio3rd
	IF(ITP+ITR.EQ.0) GO TO 100
C
	PRINT*,'******START EVENT',SCETST(1),SCETST(2),' *****'
C
	STOKESRQ = (EVAL(1)-EVAL(2))/(EVAL(1)+EVAL(2))
	STOKESI = (EVAL(1)+EVAL(2))/WTCOUNT

C
	IF(ITP.NE.0) THEN
C	  FIND UPPER EDGE OF PEAK
	  NBIGXSIG = 0
	  NBIGYSIG = 0
	  NFXBU = 2
	  NFYBU = 2
	  DO NF = 512,9,-1
	    IF(DBSPEC(NF,1).GT.THRESH(NF)) THEN
	      NBIGXSIG = NBIGXSIG + 1
	      NFXBU = NF
	      IF(NBIGXSIG.GE.3) GO TO 210
	    ELSE
	      NBIGXSIG = 0
	    ENDIF
 210	    CONTINUE
	    IF(DBSPEC(NF,2).GT.THRESH(NF)) THEN
	      NBIGYSIG = NBIGYSIG + 1
	      NFYBU = NF
	      IF(NBIGYSIG.GE.3) GO TO 220
	    ELSE
	      NBIGYSIG = 0
	    ENDIF
	  ENDDO
 220	  CONTINUE
C
	  NFHI  = MAX0(NFXBU,NFYBU) + 2
	  IF(NBIGYSIG.LT.3) NFHI = NFXBU + 2
	  print*,'nflow,nfhi',nflow,nfhi
	  IF(NFHI.LT.NFLOW) NFHI = NFLOW + 2
C
	  DO NC = 1,3
	    SPPWR(NC) = 0.
	    AVRF = 0.
	    SQF = 0.
	    IF(ITP.NE.0) THEN
	      DO NF = NFLOW,NFHI
	        VLT2 =  10.**(.1*DBSPEC(NF,NC)) - 10.**(.1*BACKGND(NF))
	        SPPWR(NC) = SPPWR(NC) + VLT2
	        AVRF = AVRF + FREQ(NF)*VLT2
	        SQF = SQF + VLT2*FREQ(NF)**2
	  	TPOWER(NF) = TPOWER(NF) + VLT2
		TCOUNT(NF) = TCOUNT(NF) + 1./3.
	      ENDDO
	    ENDIF
	    FREQHZ(NC) = AVRF/SPPWR(NC)
	    BW(NC) = SQF/SPPWR(NC) - FREQHZ(NC)**2
	    BW(NC) = SQRT(AMAX1(BW(NC),0.))
	    SPPWR(NC) = .5*SPPWR(NC)		! CONVERT TO RMS
	  ENDDO
	  PWRTOT=SPPWR(1)+SPPWR(2)+SPPWR(3)
	ENDIF
C
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(ch, item, XGSE, 1, ret_size)
	XGSE = XGSE/ERAD
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(ch, item, YGSE, 1, ret_size)
	YGSE = YGSE/ERAD
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(ch, item, ZGSE, 1, ret_size)
	ZGSE = ZGSE/ERAD
	ITEM = 'WIND_MFI_BPHI(GSE)_R4'
	ok = W_ITEM_R4(ch, item, AZMAG, 1, ret_size)
	ITEM = 'WIND_MFI_BTHETA(GSE)_R4'
	ok = W_ITEM_R4(ch, item, ELMAG, 1, ret_size)
	ITEM = 'WIND_3DP_E_TEMP_R4'
	ok = W_ITEM_R4(ch, item, TEMPE, 1, ret_size)
	ITEM = 'WIND_3DP_ION_TEMP_R4'
	ok = W_ITEM_R4(ch, item, TEMPI, 1, ret_size)
	ITEM = 'WIND_3DP_ION_DENSITY_R4'
	ok = W_ITEM_R4(ch, item, DENS, 1, ret_size)
	ITEM = 'WIND_3DP_ION_VX(GSE)_R4'
	ok = W_ITEM_R4(ch, item, VX, 1, ret_size)
	ITEM = 'WIND_3DP_ION_VY(GSE)_R4'
	ok = W_ITEM_R4(ch, item, VY, 1, ret_size)
	ITEM = 'WIND_MFI_BX(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BX, 1, ret_size)
	ITEM = 'WIND_MFI_BY(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BY, 1, ret_size)
	ITEM = 'WIND_MFI_BZ(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BZ, 1, ret_size)
	ITEM = 'WIND_MFI_BMAG_R4'
	ok = W_ITEM_R4(ch, item, BMAG, 1, ret_size)
	ITEM = 'WIND_SPIN_RATE_R4'
	ok = W_ITEM_R4(ch, item, SPINRATE, 1, ret_size)
C****
C	WRITE(37,*), SCETST(2),XMAX,YMAX,XCOUNT,YCOUNT
C****
C
	FP = 9.*SQRT(DENS)			  ! PLASMA FREQ IN KHZ
	FCE = .028*BMAG
	BETA = .402*DENS*(TEMPE+TEMPI)/BMAG**2
C
	AXMAJOR = EVAL(1)/WTCOUNT
	AXMINOR = EVAL(2)/WTCOUNT
	COSKB = (BX*GSEK(1)+BY*GSEK(2)+BZ*GSEK(3))/BMAG
C
C	CALCULATE NUMBER OF TURNS THE VECTOR MAKES IN THE VM SYSTEM
C
C	N1T = MIN0(N1,N2)
C	N2T = MAX0(N1,N2)
C	N1 = N1T
C	N2 = N2T
C	IF(IABS(N2-N1).LT.90) THEN
C	    CALL COUNT_REV(VMDATA(1,1),VMDATA(1,2),1,N1,REVS1)
C	    CALL COUNT_REV(VMDATA(1,1),VMDATA(1,2),N2,1024,REVS2)
C	    COUNT = N1 + (1024-N2+1)
C	    REVS = REVS1 + REVS2
C	ELSE
C	    CALL COUNT_REV(VMDATA(1,1),VMDATA(1,2),N1,N2,REVS)
C	    COUNT = N2-N1+1
C	ENDIF
C
C	DO N = 1,1024			! CHECK COUNT_REV
C	  WRITE(37,*) N,VMDATA(N,1),VMDATA(N,2),VMDATA(N,3)
C	ENDDO
C	IF(1) STOP
C
	CALL COUNT_REV(VMDATA(1,1),VMDATA(1,2),WT,REVS)
	FREQTRN = ABS(REVS)*FUNDFR*1024./WTCOUNT
C
c	WRITE(69,169) SCETST,VPWR(1),VPWR(2),VPWR(3),VPWR(1)+VPWR(2),
c     1	  RATIO3RD,COSKB
	VPWRTOT = VPWR(1)+VPWR(2)
 169	FORMAT(I10,I8,5E11.4,F6.3)
C	PRINT*,'REL STOKES Q,TURNS',STOKESI,STOKESRQ,REVS,FREQTRN
C	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
C	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
C	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
C	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
C
C	PLASMA PARAMETERS, ETC.

C	   write(s,'(i8.8,i6.6)',iostat=ios) s_scet(1), s_scet(2)
C	   s_scet = s(1:4)//'/'//s(5:6)//'/'//s(7:8)//' '//
C	1	s(9:10)//':'//s(11:12)//':'//s(13:14)
C
c	print date, time, event no., freq, Px, Py, Vy/Vx, angles of
c		antenna,B, phiB, mag B or fce, fp
c
	  WRITE(90,1090) scetst,XGSE,YGSE,ZGSE,
     1	     AZMAG,ELMAG,FP,FCE,BETA
 1090	format(I10,I7,3F8.2,F8.3,2F7.1,F7.3,F7.4,F8.4,2F7.0)
C
	EVAL1 = EVAL(1)/WTCOUNT
	EVAL2 = EVAL(2)/WTCOUNT
	EVAL3 = EVAL(3)/WTCOUNT
	FREQAV = .5*(FREQHZ(1) + FREQHZ(2))
	BWAV = .5*(BW(1) + BW(2))
	IF(COSKB.LT.0.) THEN
	  REVS = -REVS
	  COSKB = -COSKB
	ENDIF
	WRITE(80,1011) scetST,IRAW,ITP,ITR,nflow,nfhi,FREQTRN,FREQAV,
     1    BWAV,EVAL1,EVAL2,EVAL3,RATIO3RD,PWRTOT,SPPWR(3),STOKESI,
     2 	  STOKESRQ,COSKB,REVS
 1011	format(I10,I7,3I2,2I4,F6.1,2F6.1,3E11.3,F6.3,2F7.3,F7.3,
     1		F6.2,2F6.2)
C
	IF(OK.NE.82) GO TO 100
	RETURN
	END
	SUBROUTINE COUNT_REV(FX,FY,WT,REVS)
C
C	COUNTS NUMBER OF REVOLUTIONS OF THE VECTOR FX,FY
C	WHICH HAS NPT SAMPLES
C
	REAL FX(1),FY(1),WT(1)
C
C	PRINT*,'COUNT_REV CALLED,N1,N2',N1,N2
C	PRINT*,'FIRST SAMPLES',FX(1),FY(1)
C	PRINT*,'LAST SAMPLES',FX(1024),FY(1024)
	ANGTOT = 0.
	ANGLESV = 0.
	WTCOUNT = 0.
	IF(FY(1)*FX(1).NE.0.) ANGLESV = ATAN2D(FY(1),FX(1))
	YSAVE = FY(1)
	DO N = 1,1024
	  IF(WT(N).LT..1) GO TO 100
	  IF(FY(1)*FX(1).NE.0..AND.WTCOUNT.EQ.0.)
     1		ANGLESV = ATAN2D(FY(N),FX(N))
	  ANG = ANGLESV
	  IF(FY(N).NE.0..OR.FX(N).NE.0.) ANG = ATAN2D(FY(N),FX(N))
	  ANGTOT = ANGTOT + ANG - ANGLESV
	  WTCOUNT = WTCOUNT + 1.
	  SFY = YSAVE*FY(N)
	  IF(FX(N).LT.0..AND.SFY.LT.0.) THEN
	    IF(YSAVE.GT.0.) THEN
	      ANGTOT = ANGTOT + 360.
	    ELSE
	      ANGTOT = ANGTOT - 360.
	    ENDIF
	  ENDIF
C
C	print*,'cnt',fx(n),fy(n),ang,anglesv,angtot
C
	  ANGLESV = ANG
	  YSAVE = FY(N)
 100	  CONTINUE
	ENDDO
	REVS = ANGTOT/360.
C
	RETURN
	END
	SUBROUTINE COUNT_REVOLD(FX,FY,N1,N2,REVS)
C
C	COUNTS NUMBER OF REVOLUTIONS OF THE VECTOR FX,FY
C	WHICH HAS NPT SAMPLES
C
	REAL FX(1),FY(1)
C
C	PRINT*,'COUNT_REV CALLED,N1,N2',N1,N2
C	PRINT*,'FIRST SAMPLES',FX(1),FY(1)
C	PRINT*,'LAST SAMPLES',FX(1024),FY(1024)
	ANGTOT = 0.
	ANGLESV = 0.
	IF(FY(1)*FX(1).NE.0.) ANGLESV = ATAN2D(FY(1),FX(1))
	YSAVE = FY(1)
	DO N = N1+1,N2
	  ANG = ANGLESV
	  IF(FY(N).NE.0..OR.FX(N).NE.0.) ANG = ATAN2D(FY(N),FX(N))
	  ANGTOT = ANGTOT + ANG - ANGLESV
	  SFY = YSAVE*FY(N)
	  IF(FX(N).LT.0..AND.SFY.LT.0.) THEN
	    IF(YSAVE.GT.0.) THEN
	      ANGTOT = ANGTOT + 360.
	    ELSE
	      ANGTOT = ANGTOT - 360.
	    ENDIF
	  ENDIF
C
C	print*,'cnt',fx(n),fy(n),ang,anglesv,angtot
C
	  ANGLESV = ANG
	  YSAVE = FY(N)
	ENDDO
	REVS = ANGTOT/360.
C
	RETURN
	END
       SUBROUTINE SINFIT(X,FDATA,ERR,Y,SUMSQ)
C
C	THIS ROUTINE FITS A SINE WAVE TO DATA.  X IS SUPPOSED TO
C	BE VARIED BY HUNTMN, WHILE Y(1) TO Y(3) ARE DETERMINED BY
C	LEAST SQUARES FIT IN CLOSED FORM.
C
C	X IS FREQUENCY IN CYCLES PER SAMPLE  (A SMALL NUMBER)
C
C	Y(1) IS AVERAGE (OFFSET)
C	Y(2) IS COEFF OF COSINE
C	Y(3) IS COEFF OF SINE
C
        DOUBLE PRECISION COS1,SIN1,COSN,SINN,SINNT,THT0
        REAL X,C(3,3),Y(3),ERR(1024),WT(1024),FDATA(1024)
        DATA TWOPI /6.2831853/
	DATA WT /1024*1./
C
	DO I = 1,3
	  Y(I)= 0.
	  DO J = 1,3
	    C(I,J) = 0.
	  ENDDO
	ENDDO
C
	THT0 = TWOPI*X
	COS1 = DCOS(THT0)
	SIN1 = DSIN(THT0)
	SINN = SIN1
	COSN = COS1
	DO N = 1,1024
	  Y(1) = Y(1) + WT(N)*FDATA(N)
	  Y(2) = Y(2) + WT(N)*FDATA(N)*COSN
	  Y(3) = Y(3) + WT(N)*FDATA(N)*SINN
	  C(1,1) = C(1,1) + WT(N)
	  C(1,2) = C(1,2) + WT(N)*COSN
	  C(1,3) = C(1,3) + WT(N)*SINN
	  C(2,2) = C(2,2) + WT(N)*COSN**2
	  C(2,3) = C(2,3) + WT(N)*COSN*SINN
	  C(3,3) = C(3,3) + WT(N)*SINN**2
C
	  SINNT = COSN*SIN1 + SINN*COS1
	  COSN = COS1*COSN - SIN1*SINN
	  SINN = SINNT
	ENDDO
	C(3,1)= C(1,3)
	C(3,2) = C(2,3)
	C(2,1) = C(1,2)
	CALL GAUSSJ(C,3,3,Y,1,1)
C
	SUMSQ = 0.
	SUMN = 1.E-8
	THT0 = TWOPI*X
	COS1 = DCOS(THT0)
	SIN1 = DSIN(THT0)
	SINN = SIN1
	COSN = COS1
	DO N = 1,1024
          ERR(N) =  FDATA(N) - Y(2)*COSN - Y(3)*SINN - Y(1)
	  SUMSQ = SUMSQ + WT(N)*ERR(N)**2
	  SUMN = SUMN + WT(N)
	  SINNT = COSN*SIN1 + SINN*COS1
	  COSN = COS1*COSN - SIN1*SINN
	  SINN = SINNT
	ENDDO
	SUMSQ = SUMSQ/SUMN
C	FIXED_V1 = Y(1) + Y(2)*COS(3.*THT0) + Y(3)*SIN(3.*THT0)
	FIXED_V = FDATA(3) - ERR(3)
C	PRINT*,'SINFIT EXIT,X,ERR,FDATA,SUMSQ',X,ERR(3),FDATA(3),SUMSQ
       RETURN
       END

ext: for


text:
!
!  Use adaptive mesh refinement to solve the hyperbolic 2-d equation:
!
!              u  +  f(u)    + g(u) + h(u)  = 0
!               t         x        y      z
!
! or the more general non-conservation law form:
!              u  +  A u     + B u    + C u = 0
!               t         x        y       z
!
!  using the wave propagation method as in CLAWPACK in combination
!  with the locally uniform embedded grids of AMR.
!
!  Estimate error with Richardson extrap. (in errest.f)
!  + gradient checking (in errsp.f).  Initial conditions set
!  in (qinit.f), b.c.'s in (bc3amr.f).
!
!  Specify rectangular domain from
!           (xlower,ylower,zlower) to (xupper,yupper,zupper).
!
!
! =========================================================================
!
!  This software is made available for research and instructional use only.
!  You may copy and use this software without charge for these non-commercial
!  purposes, provided that the copyright notice and associated text is
!  reproduced on all copies.  For all other uses (including distribution of
!  modified versions), please contact the author at the address given below.
!
!  *** This software is made available "as is" without any assurance that it
!  *** will work for your purposes.  The software may in fact have defects, so
!  *** use the software at your own risk.
!
!  --------------------------------------
!    AMRCLAW Version 5.0,  2012
!    Homepage: http://www.clawpack.org
!  --------------------------------------
!
!   Authors:
!
!             Marsha J. Berger
!             Courant Institute of Mathematical Sciences
!             New York University
!             251 Mercer St.
!             New York, NY 10012
!             berger@cims.nyu.edu
!
!             Randall J. LeVeque
!             Applied Mathematics
!             Box 352420
!             University of Washington,
!             Seattle, WA 98195-2420
!             rjl@uw.edu
!
! =========================================================================
program amr3

    use amr_module, only: dbugunit, parmunit, outunit, inunit, matunit
    use amr_module, only: mxnest, rinfinity, iinfinity
    use amr_module, only: xupper, yupper, zupper
    use amr_module, only: xlower, ylower, zlower
    use amr_module, only: hxposs, hyposs, hzposs
    use amr_module, only: intratx, intraty, intratz, kratio
    use amr_module, only: cfl, cflv1, cflmax, evol

    use amr_module, only: checkpt_style, checkpt_interval, tchk, nchkpt
    use amr_module, only: rstfile, check_a

    use amr_module, only: max1d, maxvar, maxlv

    use amr_module, only: method, mthlim, use_fwaves,dimensional_split
    use amr_module, only: nghost, mwaves, mcapa, auxtype
    use amr_module, only: tol, tolsp, flag_richardson, flag_gradient

    use amr_module, only: nghost, mthbc
    use amr_module, only: xperdom, yperdom,zperdom

    use amr_module, only: nstop, nout, iout, tfinal, tout, output_style
    use amr_module, only: output_format, printout, verbosity_regrid
    use amr_module, only: output_q_components, output_aux_components
    use amr_module, only: output_aux_onlyonce, matlabu

    use amr_module, only: lfine, lentot, iregridcount, avenumgrids
    use amr_module, only: tvoll,tvollCPU, rvoll, rvol, mstart, possk, ibuff
    use amr_module, only: timeRegridding,timeRegriddingCPU
    use amr_module, only: timeBound,timeStepgrid
    use amr_module, only: timeBoundCPU,timeStepgridCPU
    use amr_module, only: timeValout,timeValoutCPU
    use amr_module, only: kcheck, iorder, lendim, lenmax

    use amr_module, only: dprint, eprint, edebug, gprint, nprint, pprint
    use amr_module, only: rprint, sprint, tprint, uprint

    use amr_module, only: t0, tstart_thisrun

    !use regions_module, only: set_regions
    use gauges_module, only: set_gauges, num_gauges

    implicit none

    ! Local variables
    integer :: i, iaux, mw, level
    integer :: ndim, nvar, naux, mcapa1, mindim
    integer :: nstart, nsteps, nv1, nx, ny, nz, lentotsave
    !integer :: nstart, nsteps, nv1, nx, ny, lentotsave
    integer :: num_gauge_SAVE
    integer :: omp_get_max_threads, maxthreads
    real(kind=8) :: time, ratmet, cut, dtinit, dt_max
    logical :: vtime, rest, output_t0

    ! Timing variables
    integer :: clock_start, clock_finish, clock_rate, ttotal
    real(kind=8) :: cpu_start, cpu_finish,ttotalcpu
    integer, parameter :: timing_unit = 48

    ! Common block variables
    real(kind=8) :: dxmin, dymin, dzmin

    common /comfine/ dxmin,dymin, dzmin

    character(len=364) :: format_string
    character(len=*), parameter :: clawfile = 'claw.data'
    character(len=*), parameter :: amrfile  = 'amr.data'
    character(len=*), parameter :: outfile  = 'fort.amr'
    character(len=*), parameter :: dbugfile = 'fort.debug'
    character(len=*), parameter :: matfile  = 'fort.nplot'
    character(len=*), parameter :: parmfile = 'fort.parameters'
    character(len=*), parameter :: timing_file = 'timing.txt'

    ! Open parameter and debug files
    open(dbugunit,file=dbugfile,status='unknown',form='formatted')
    open(parmunit,file=parmfile,status='unknown',form='formatted')

    maxthreads = 1    !! default, if no openmp

    ! Open AMRClaw primary parameter file
    call opendatafile(inunit,clawfile)

    ! Number of space dimensions, not really a parameter but we read it in and
    ! check to make sure everyone is on the same page.
    read(inunit,"(i1)") ndim
    if (ndim /= 3) then
        print *,'Error ***   ndim = 3 is required,  ndim = ',ndim
        print *,'*** Are you sure input has been converted'
        print *,'*** to Clawpack 5.x form?'
        stop
    endif

    ! Domain variables
    read(inunit,*) xlower, ylower, zlower
    read(inunit,*) xupper, yupper, zupper
    read(inunit,*) nx, ny, nz
    read(inunit,*) nvar    ! meqn
    read(inunit,*) mwaves
    read(inunit,*) naux
    read(inunit,*) t0

    ! ==========================================================================
    ! Output Options
    ! Output style
    read(inunit,*) output_style
    if (output_style == 1) then
        read(inunit,*) nout
        read(inunit,*) tfinal
        read(inunit,*) output_t0

        iout = 0
    else if (output_style == 2) then
        read(inunit,*) nout
        allocate(tout(nout))
        read(inunit,*) (tout(i), i=1,nout)
        output_t0 = (tout(1) == t0)
        ! Move output times down one index
        if (output_t0) then
            nout = nout - 1
            do i=1,nout
                tout(i) = tout(i+1)
            enddo
        endif
        iout = 0
        tfinal = tout(nout)
    else if (output_style == 3) then
        read(inunit,*) iout
        read(inunit,*) nstop
        read(inunit,*) output_t0
        nout = 0
        tfinal = rinfinity
    else
        stop "Error ***   Invalid output style."
    endif

    ! Error checking
    if ((output_style == 1) .and. (nout > 0)) then
        allocate(tout(nout))
        do i=1,nout
            tout(i) = t0 + i * (tfinal - t0) / real(nout,kind=8)
        enddo
    endif

    ! What and how to output
    read(inunit,*) output_format
    allocate(output_q_components(nvar))
    read(inunit,*) (output_q_components(i),i=1,nvar)
    if (naux > 0) then
        allocate(output_aux_components(naux))
        read(inunit,*) (output_aux_components(i),i=1,naux)
        read(inunit,*) output_aux_onlyonce
    endif
    ! ==========================================================================

    ! ==========================================================================
    !  Algorithm parameters

    read(inunit,*) possk(1)   ! dt_initial
    read(inunit,*) dt_max     ! dt_max
    read(inunit,*) cflv1      ! cfl_max
    read(inunit,*) cfl        ! clf_desired
    read(inunit,*) nv1        ! steps_max

    if (output_style /= 3) then
        !nstop = nv1
        nstop = iinfinity ! basically disables this test
    endif

    read(inunit,*) vtime      ! dt_variable
    if (vtime) then
        method(1) = 2
    else
        method(1) = 1
    endif

    read(inunit,*) method(2)  ! order
    iorder = method(2)
    read(inunit,*) method(3)  ! order_trans

    read(inunit,*) dimensional_split
    if (dimensional_split > 1) then
        print *, '*** ERROR ***  dimensional_split = ', dimensional_split
        print *, ' Strang splitting not supported in amrclaw'
        stop
    endif

    read(inunit,*) method(4)   ! verbosity
    read(inunit,*) method(5)   ! src_split
    read(inunit,*) mcapa1

    read(inunit,*) use_fwaves
    allocate(mthlim(mwaves))
    read(inunit,*) (mthlim(mw), mw=1,mwaves)

    ! Boundary conditions
    read(inunit,*) nghost
    read(inunit,*) mthbc(1),mthbc(3),mthbc(5)
    read(inunit,*) mthbc(2),mthbc(4),mthbc(6)

    ! 1,2,3 = low x,y,z  4,5,6 = hi x,y,z
    xperdom = (mthbc(1) == 2 .and. mthbc(2) == 2)
    yperdom =  (mthbc(3) == 2 .and. mthbc(4) == 2)
    zperdom =  (mthbc(5) == 2 .and. mthbc(6) == 2)
    !! spheredom =  (mthbc(3) == 5 .and. mthbc(4) == 5) NOT IMPLEMENTED

    if ((mthbc(1).eq.2 .and. mthbc(2).ne.2) .or. &
        (mthbc(2).eq.2 .and. mthbc(1).ne.2)) then

        print *, '*** ERROR ***  periodic boundary conditions: '
        print *, '  mthbc(1) and mthbc(2) must BOTH be set to 2'
        stop
    endif

    if ((mthbc(3).eq.2 .and. mthbc(4).ne.2) .or. &
        (mthbc(4).eq.2 .and. mthbc(3).ne.2)) then

        print *, '*** ERROR ***  periodic boundary conditions: '
        print *, '  mthbc(3) and mthbc(4) must BOTH be set to 2'
        stop
    endif

    if ((mthbc(5).eq.2 .and. mthbc(6).ne.2) .or. &
        (mthbc(6).eq.2 .and. mthbc(5).ne.2)) then

        print *, '*** ERROR ***  periodic boundary conditions: '
        print *, '  mthbc(5) and mthbc(6) must BOTH be set to 2'
        stop
    endif

   ! if ((mthbc(3).eq.5 .and. mthbc(4).ne.5) .or. &
   !     (mthbc(4).eq.5 .and. mthbc(3).ne.5)) then
   !
   !     print *, '*** ERROR ***  sphere bcs at top and bottom: '
   !     print *, '  mthbc(3) and mthbc(4) must BOTH be set to 5'
   !     stop
   ! endif

   ! if (spheredom .and. .not. xperdom) then
   !
   !        print *,'*** ERROR ***  sphere bcs at top and bottom: '
   !     print *,'must go with periodic bcs at left and right  '
   !     stop
   ! endif

    ! ==========================================================================
    !  Restart and Checkpointing

    read(inunit,*) rest
    read(inunit,*) rstfile

    read(inunit,*) checkpt_style
    if (checkpt_style == 0) then
        ! Never checkpoint:
        checkpt_interval = iinfinity

    else if (abs(checkpt_style) == 2) then
        read(inunit,*) nchkpt
        allocate(tchk(nchkpt))
        read(inunit,*) (tchk(i), i=1,nchkpt)

    else if (abs(checkpt_style) == 3) then
        ! Checkpoint every checkpt_interval steps on coarse grid
        read(inunit,*) checkpt_interval
    endif

    close(inunit)

    ! ==========================================================================
    !  Refinement Control
    call opendatafile(inunit, amrfile)

    read(inunit,*) mxnest
    if (mxnest <= 0) then
        stop 'Error ***   mxnest (amrlevels_max) <= 0 not allowed'
    endif

    if (mxnest > maxlv) then
        stop 'Error ***   mxnest > max. allowable levels (maxlv) in common'
    endif

    ! Anisotropic refinement always allowed in 5.x:
    read(inunit,*) (intratx(i),i=1,max(1,mxnest-1))
    read(inunit,*) (intraty(i),i=1,max(1,mxnest-1))
    read(inunit,*) (intratz(i),i=1,max(1,mxnest-1))
    read(inunit,*) (kratio(i), i=1,max(1,mxnest-1))
    read(inunit,*)

    do i=1,mxnest-1
        if ((intratx(i) > max1d) .or. (intraty(i) > max1d) .or. &
                                      (intratz(i) > max1d)) then
            print *, ""
            format_string = "(' *** Error: Refinement ratios must be no " // &
                            "larger than max1d = ',i5,/,'     (set max1d" // &
                            " in amr_module.f90)')"
            print format_string, max1d
            stop
        endif
    enddo

    if (naux > 0) then
        allocate(auxtype(naux))
        read(inunit,*) (auxtype(iaux), iaux=1,naux)
    endif
    read(inunit,*)

    read(inunit,*) flag_richardson
    read(inunit,*) tol            ! for richardson
    read(inunit,*) flag_gradient
    read(inunit,*) tolsp          ! for gradient
    read(inunit,*) kcheck
    read(inunit,*) ibuff
    read(inunit,*) cut
    read(inunit,*) verbosity_regrid

    ! read verbose/debugging flags
    read(inunit,*) dprint
    read(inunit,*) eprint
    read(inunit,*) edebug
    read(inunit,*) gprint
    read(inunit,*) nprint
    read(inunit,*) pprint
    read(inunit,*) rprint
    read(inunit,*) sprint
    read(inunit,*) tprint
    read(inunit,*) uprint

    close(inunit)
    ! Finished with reading in parameters
    ! ==========================================================================

    ! Look for capacity function via auxtypes:
    mcapa = 0
    do iaux = 1, naux
        if (auxtype(iaux) == "capacity") then
            if (mcapa /= 0) then
                print *, " only 1 capacity allowed"
                stop
            else
                mcapa = iaux
            endif
        endif

        ! Change to Version 4.1 terminology:
        if (auxtype(iaux) == "leftface") auxtype(iaux) = "xleft"
        if (auxtype(iaux) == "bottomface") auxtype(iaux) = "yleft"
        if (auxtype(iaux) == "frontface") auxtype(iaux) = "zleft"
        if (.not. (auxtype(iaux) .eq."xleft" .or. &
                   auxtype(iaux) .eq. "yleft".or. &
                   auxtype(iaux) .eq. "zleft".or. &
                   auxtype(iaux) .eq. "capacity".or. &
                   auxtype(iaux) .eq. "center"))  then
            print *," unknown type for auxiliary variables"
            print *," use  xleft/yleft/zleft/center/capacity"
            stop
        endif
    enddo

    ! Error checking of input data
    if (mcapa /= mcapa1) then
        stop 'Error ***  mcapa does not agree with auxtype'
    endif
    if (nvar > maxvar) then
        stop 'Error ***   nvar > maxvar in common'
    endif
    if (2*nghost > min(nx,ny) .and. ny /= 1) then
        mindim = 2 * nghost
        print *, 'Error ***   need finer domain >', mindim, ' cells'
        stop
    endif
    if (mcapa > naux) then
        stop 'Error ***   mcapa > naux in input file'
    endif

    if (.not. vtime .and. nout /= 0) then
        print *,        ' cannot specify output times with fixed dt'
        stop
    endif

    ! Write out parameters
    write(parmunit,*) ' '
    write(parmunit,*) 'Running amrclaw with parameter values:'
    write(parmunit,*) ' '

    print *, ' '
    print *, 'Running amrclaw ...  '
    print *, ' '

    hxposs(1) = (xupper - xlower) / nx
    hyposs(1) = (yupper - ylower) / ny
    hzposs(1) = (zupper - zlower) / nz

    ! initialize frame number for output.
    ! Note: might be reset in restrt if this is a restart
    if (output_t0) then
        matlabu   = 0
    else
        matlabu   = 1
    endif

    ! Boolean check_a tells which checkpoint file to use next if alternating
    ! between only two files via check_twofiles.f, unused otherwise.
    ! May be reset in call to restrt, otherwise default to using aaaaa file.
    check_a = .true.

    if (rest) then

        open(outunit, file=outfile, status='unknown', position='append', &
                      form='formatted')

        call restrt(nsteps,time,nvar,naux)
        nstart  = nsteps
        tstart_thisrun = time
        print *, ' '
        print *, 'Restarting from previous run'
        print *, '   at time = ',time
        print *, ' '
        ! Call user routine to set up problem parameters:
        call setprob()

        call set_gauges(rest, nvar, naux)
    else

        open(outunit, file=outfile, status='unknown', form='formatted')

        tstart_thisrun = t0

        ! Call user routine to set up problem parameters:
        call setprob()
        call set_gauges(rest, nvar, naux)

        cflmax = 0.d0   ! otherwise use previously heckpointed val

        lentot = 0
        lenmax = 0
        lendim = 0
        rvol   = 0.0d0
        do i   = 1, mxnest
            rvoll(i) = 0.0d0
        enddo
        evol = 0.0d0
        call stst1()

        ! changed 4/24/09: store dxmin,dymin for setaux before
        ! grids are made, in order to average up from finest grid.
        dxmin = hxposs(mxnest)
        dymin = hyposs(mxnest)
        dzmin = hzposs(mxnest)

        call domain(nvar,vtime,nx,ny,nz,naux,t0)

        ! Hold off on gauges until grids are set.
        ! The fake call to advance at the very first timestep
        ! looks at the gauge array but it is not yet built
        num_gauge_SAVE = num_gauges
        num_gauges = 0
        call setgrd(nvar,cut,naux,dtinit,t0)
        num_gauges = num_gauge_SAVE

! commented out to match 4-x version
!!$        if (possk(1) .gt. dtinit*cflv1/cfl .and. vtime) then
!!$            ! initial time step was too large. reset to dt from setgrd
!!$            print *, "*** Initial time step reset for desired cfl"
!!$            possk(1) = dtinit
!!$            do i = 2, mxnest-1
!!$                possk(i) = possk(i-1)*kratio(i-1)
!!$            end do
!!$        endif

        time = t0
        nstart = 0
    endif

    write(parmunit,*) ' '
    write(parmunit,*) '--------------------------------------------'
    write(parmunit,*) ' '
    write(parmunit,*) '   rest = ', rest, '   (restart?)'
    write(parmunit,*) '   start time = ',time
    write(parmunit,*) ' '

!$   maxthreads = omp_get_max_threads()
     write(outunit,*)" max threads set to ",maxthreads
     print *," max threads set to ",maxthreads

    !
    !  print out program parameters for this run
    !
    format_string = "(/' amrclaw parameters:',//," // &
                      "' error tol            ',e12.5,/," // &
                      "' spatial error tol    ',e12.5,/," // &
                      "' order of integrator     ',i9,/," // &
                      "' error checking interval ',i9,/," // &
                      "' buffer zone size        ',i9,/," // &
                      "' nghost                  ',i9,/," // &
                      "' volume ratio cutoff  ',e12.5,/," // &
                      "' max. refinement level   ',i9,/," // &
                      "' user sub. calling times ',i9,/," // &
                      "' cfl # (if var. delt) ',e12.5,/)"
    write(outunit,format_string) tol,tolsp,iorder,kcheck,ibuff,nghost,cut, &
                                 mxnest,checkpt_interval,cfl
    format_string = "(' xupper(upper corner) ',e12.5,/," // &
                     "' yupper(upper corner) ',e12.5,/," // &
                     "' zupper(upper corner) ',e12.5,/," // &
                     "' xlower(lower corner) ',e12.5,/," // &
                     "' ylower(lower corner) ',e12.5,/," // &
                     "' zlower(lower corner) ',e12.5,/," // &
                     "' nx = no. cells in x dir.',i9,/," // &
                     "' ny = no. cells in y dir.',i9,/," // &
                     "' nz = no. cells in z dir.',i9,/," // &
                     "' refinement ratios       ',6i5,/,/)"
    write(outunit,format_string) xupper,yupper,zupper,      &
                                 xlower,ylower,zlower,nx,ny,nz
    write(outunit,"(' refinement ratios:       ',6i5,/)"  ) &
                                                        (intratx(i),i=1,mxnest)
    write(outunit,"(' refinement ratios:       ',6i5,/)"  ) &
                                                        (intraty(i),i=1,mxnest)
    write(outunit,"(' refinement ratios:       ',6i5,/)"  ) &
                                                        (intratz(i),i=1,mxnest)
    write(outunit,"(' no. auxiliary vars.     ',i9)") naux
    write(outunit,"('       var ',i5,' of type ', a10)") &
                                                (iaux,auxtype(iaux),iaux=1,naux)
    if (mcapa > 0) write(outunit,"(' capacity fn. is aux. var',i9)") mcapa

    print *, ' '
    print *, 'Done reading data, starting computation ...  '
    print *, ' '

    call outtre (mstart,printout,nvar,naux)
    write(outunit,*) "  original total mass ..."
    call conck(1,nvar,naux,time,rest)
    if (output_t0) then
        call valout(1,lfine,time,nvar,naux)
    endif
    close(parmunit)

    ! Timing
    call system_clock(clock_start,clock_rate)
    call cpu_time(cpu_start)

    ! --------------------------------------------------------
    !  Tick is the main routine which drives the computation:
    ! --------------------------------------------------------
    call tick(nvar,cut,nstart,vtime,time,naux,t0,rest,dt_max)
    ! --------------------------------------------------------

    call system_clock(clock_finish,clock_rate)
    call cpu_time(cpu_finish)

    !output timing data
    open(timing_unit, file=timing_file, status='unknown', form='formatted')
    write(*,*)
    write(timing_unit,*)
    format_string="('============================== Timing Data ==============================')"
    write(timing_unit,format_string)
    write(*,format_string)

    write(*,*)
    write(timing_unit,*)

    !Integration time
    format_string="('Integration Time (stepgrid + BC + overhead)')"
    write(timing_unit,format_string)
    write(*,format_string)

    !Advanc time
    format_string="('Level           Wall Time (seconds)    CPU Time (seconds)   Total Cell Updates')"
    write(timing_unit,format_string)
    write(*,format_string)
    ttotalcpu=0.d0
    ttotal=0
    do level=1,mxnest
        format_string="(i3,'           ',1f15.3,'        ',1f15.3,'    ', e17.3)"
        write(timing_unit,format_string) level, &
             real(tvoll(level),kind=8) / real(clock_rate,kind=8), tvollCPU(level), rvoll(level)
        write(*,format_string) level, &
             real(tvoll(level),kind=8) / real(clock_rate,kind=8), tvollCPU(level), rvoll(level)
        ttotalcpu=ttotalcpu+tvollCPU(level)
        ttotal=ttotal+tvoll(level)
    end do

    format_string="('total         ',1f15.3,'        ',1f15.3,'    ', e17.3)"
    write(timing_unit,format_string) &
             real(ttotal,kind=8) / real(clock_rate,kind=8), ttotalCPU, rvol
    write(*,format_string) &
             real(ttotal,kind=8) / real(clock_rate,kind=8), ttotalCPU, rvol

    write(*,*)
    write(timing_unit,*)

    format_string="('All levels:')"
    write(*,format_string)
    write(timing_unit,format_string)

    !stepgrid
    format_string="('stepgrid      ',1f15.3,'        ',1f15.3,'    ',e17.3)"
    write(timing_unit,format_string) &
         real(timeStepgrid,kind=8) / real(clock_rate,kind=8), timeStepgridCPU
    write(*,format_string) &
         real(timeStepgrid,kind=8) / real(clock_rate,kind=8), timeStepgridCPU

    !bound
    format_string="('BC/ghost cells',1f15.3,'        ',1f15.3)"
    write(timing_unit,format_string) &
         real(timeBound,kind=8) / real(clock_rate,kind=8), timeBoundCPU
    write(*,format_string) &
         real(timeBound,kind=8) / real(clock_rate,kind=8), timeBoundCPU

    !regridding time
    format_string="('Regridding    ',1f15.3,'        ',1f15.3,'  ')"
    write(timing_unit,format_string) &
            real(timeRegridding,kind=8) / real(clock_rate,kind=8), timeRegriddingCPU
    write(*,format_string) &
            real(timeRegridding,kind=8) / real(clock_rate,kind=8), timeRegriddingCPU

    !output time
    format_string="('Output (valout)',1f14.3,'        ',1f15.3,'  ')"
    write(timing_unit,format_string) &
            real(timeValout,kind=8) / real(clock_rate,kind=8), timeValoutCPU
    write(*,format_string) &
            real(timeValout,kind=8) / real(clock_rate,kind=8), timeValoutCPU

    write(*,*)
    write(timing_unit,*)

    !Total Time
    format_string="('Total time:   ',1f15.3,'        ',1f15.3,'  ')"
    write(timing_unit,format_string) &
            real(clock_finish - clock_start,kind=8) / real(clock_rate,kind=8), &
            cpu_finish-cpu_start
    write(*,format_string) &
            real(clock_finish - clock_start,kind=8) / real(clock_rate,kind=8), &
            cpu_finish-cpu_start

    format_string="('Using',i3,' thread(s)')"
    write(timing_unit,format_string) maxthreads
    write(*,format_string) maxthreads

    write(*,*)
    write(timing_unit,*)

    write(*,"('Note: The CPU times are summed over all threads.')")
    write(timing_unit,"('Note: The CPU times are summed over all threads.')")
    write(*,"('      Total time includes more than the subroutines listed above')")
    write(timing_unit,"('      Total time includes more than the subroutines listed above')")

    !end of timing data
    write(*,*)
    write(timing_unit,*)
    format_string="('=========================================================================')"
    write(timing_unit,format_string)
    write(*,format_string)
    write(*,*)
    write(timing_unit,*)
    close(timing_unit)

    !format_string = "('Total time to solution = ',1f16.8,' s')"
    !write(outunit,format_string) &
    !        real(clock_finish - clock_start,kind=8) / real(clock_rate,kind=8)
    !write(*,format_string) &
    !        real(clock_finish - clock_start,kind=8) / real(clock_rate,kind=8)

    !do level = 1, mxnest
    !  format_string = "('Total advanc time on level ',i3,' = ',1f16.8,' s')"
    !  write(outunit,format_string) level, &
    !         real(tvoll(level),kind=8) / real(clock_rate,kind=8)
    !  write(*,format_string) level, &
    !         real(tvoll(level),kind=8) / real(clock_rate,kind=8)
    !end do

    ! Done with computation, cleanup:
    lentotsave = lentot
    call cleanup(nvar,naux)
    if (lentot /= 0) then
        write(outunit,*) lentot," words not accounted for in memory cleanup"
        print *,         lentot," words not accounted for in memory cleanup"
    endif

    !
    ! report on statistics
    !
    open(matunit,file=matfile,status='unknown',form='formatted')
    write(matunit,*) matlabu-1
    write(matunit,*) mxnest
    close(matunit)

    write(outunit,*)
    write(outunit,*)
    do i = 1, mxnest
      if (iregridcount(i) > 0) then
        write(outunit,801) i,avenumgrids(i)/iregridcount(i),iregridcount(i)
 801    format("for level ",i3, " average num. grids = ",f10.2," over ",i10,  &
               " regridding steps")
      endif
    end do

    write(outunit,*)
    write(outunit,*)
    write(outunit,"('current  space usage = ',i12)") lentotsave
    write(outunit,"('maximum  space usage = ',i12)") lenmax
    write(outunit,"('need space dimension = ',i12,/)") lendim

    write(outunit,"('number of cells advanced for time integration = ',f20.6)")&
                    rvol
    do level = 1,mxnest
        write(outunit,"(3x,'# cells advanced on level ',i4,' = ',f20.2)") &
                    level, rvoll(level)
    enddo

    write(outunit,"('number of cells advanced for error estimation = ',f20.6,/)") &
                     evol
    if (evol + rvol > 0.d0) then
        ratmet = rvol / (evol + rvol) * 100.0d0
    else
        ratmet = 0.0d0
    endif
    write(outunit,"(' percentage of cells advanced in time  = ', f10.2)") ratmet
    write(outunit,"(' maximum Courant number seen = ', f10.2)") cflmax

    write(outunit,"(//,' ------  end of AMRCLAW integration --------  ')")

    ! Close output and debug files.
    close(outunit)
    close(dbugunit)

end program amr3

max_stars_repo_name: BrisaDavis/amrclaw

max_issues_repo_name: BrisaDavis/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!  *********************************************************************
!  *    SUBROUTINE TO DETERMINE LIMITS ON THE EXTINCTION COEFFICIENT   *
!  *********************************************************************
!
      SUBROUTINE CONSTR(SURF,DMIX,EMIN,ROOT,NUMGR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(2),B(2),ROOT(2)
      INCLUDE 'blmdim.inc'
      INCLUDE 'arran.inc'
!
!  DETERMINES LIMITS ON THE EXTINCTION COEFFICIENT
!
    1 ROOT(1)=-1.0
      ROOT(2)=-1.0
      S0=DEXP(-ZVEC(NZ))
      IF (SURF .LE. S0) RETURN
!
!  FIND INTERVAL CONTAINING ALL ROOTS
!
      X=-DLOG(SURF)
      CALL EBCALC(X,FX,DERX,NUMGR)
      PHI=(FUN(NZ,NUMGR)-FX)/EMIN
      Y=X+PHI
      IF (DERX .GE. EMIN) ROOT(1)=0.0
      IF (Y .GE. ZVEC(NZ)) ROOT(2)=PHI/DMIX
      IF (ROOT(1) .GE. 0.0 .AND. ROOT(2) .GE. 0.0) RETURN
!
!  SPLIT INTERVAL TO ISOLATE EACH ROOT.
!
      B(1)=X
      B(2)=Y
      DO 110 K=1,10
      A(1)=0.5*(B(1)+B(2))
      CALL EBCALC(A(1),FM,FPR,NUMGR)
      EBAR=(FM-FX)/(A(1)-X)
      IF (EBAR .GE. EMIN) GO TO 120
      DERIV=(FPR-EBAR)/(A(1)-X)
      IF (DERIV .LT. 0.0) B(2)=A(1)
  110 IF (DERIV .GE. 0.0) B(1)=A(1)
!
!  NO ROOTS EXIST
!
      RETURN
!
!  WE HAVE SEPARATED THE ROOTS.
!
  120 A(2)=A(1)
      DO 130 I=1,2
      IF (ROOT(I) .GE. 0.0) GO TO 130
      DO 125 K=1,10
      XM=0.5*(A(I)+B(I))
      CALL EBCALC(XM,FM,FP,NUMGR)
      EBAR=(FM-FX)/(XM-X)
      IF (EBAR .GT. EMIN) A(I)=XM
  125 IF (EBAR .LE. EMIN) B(I)=XM
      ROOT(I)=(0.5*(A(I)+B(I))-X)/DMIX
  130 CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE PLTRIM(LINE,L)
      CHARACTER*(*) LINE
      CHARACTER CH

      L1 = LEN(LINE)
 2240 IF (.NOT. (L1.GT.0)) GO TO 2260
      CH = LINE(L1:L1)
      IF (CH.EQ.' ') THEN
         GO TO 2250

      END IF

      GO TO 2260

 2250 L1 = L1 - 1
      GO TO 2240

 2260 CONTINUE
      L = L1
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
! Main reading routine, with the read_inputdata subroutine used to obtain input data
module check_input_mod
    use reading_utilities_mod
    use output_mod
    use types_mod
    use convert_mpar_mod
    use gen_util_mod
    implicit none

    private                         ! All private by default

    ! Assign public parts
    public  :: bad_input            ! Routine to notify that bad input entered, and end the program
    public  :: check_glob           ! Check the global settings input

    public  :: check_sim_ext        ! Check the simulation settings input for external simulation
    public  :: check_sim_atp        ! Check the simulation settings input for atp simulation
    public  :: check_sim_mps        ! Check the simulation settings input for mp simulation (mps)
    public  :: check_sim_atp_mr!Check the simulation settings input for atp material removal simulation
    public  :: check_sim_usr        ! Check the settings for user simulation

    public  :: check_opt_nlopt      ! Check the optimization settings input for nlopt
    public  :: check_opt_pspace     ! Check the optimization settings input for parameter space simulations
    public  :: check_opt_usr

    public  :: check_overall        ! Check stuff that can only be checked after all info exists

    ! Values for comparisons
    double precision, parameter     :: numtol = 1.d-20

    ! Default values for glob arrays
    integer, parameter              ::  DEF_ipar_sctype = 1

    ! Default values for sim:atp arrays
    double precision, parameter     :: DEF_err_scale = 1.d0
    integer, parameter, dimension(1):: DEF_result_steps = (/0/)
    integer, parameter, dimension(1):: DEF_error_steps = (/0/)
    double precision, parameter     :: DEF_statev_init = 0.d0

    ! Default values for opt arrays
    double precision, parameter, dimension(1)   :: DEF_initial_step = (/0.d0/)  !0=> nlopt default
    double precision, parameter, dimension(1)   :: DEF_xtol_rel = (/0.d0/)      !0=> no tolerance
    double precision, parameter, dimension(1)   :: DEF_xtol_abs = (/0.d0/)      !0=> no tolerance

    ! Variables
    logical     :: nlgeom
    integer     :: nstatv, nchannels
    integer     :: nnod = 1
    integer     :: nel = 1

    contains

subroutine bad_input(message)
implicit none
    character(len=*)    :: message

    call close_input()
    call write_output(message, 'error', 'inp')

end subroutine

subroutine check_glob(glob)
implicit none
    type(glob_typ)                  :: glob
    integer                         :: nparam, noptim, nsets   ! Number of specified parameters
    integer                         :: k1
    logical                         :: lower_bound_nok, upper_bound_nok !Logical variable for checking bounds
    double precision, allocatable   :: mpar(:)
    ! Subroutine tasks:
    ! 0) Assign default run_type if no assigned (needed for later stages)
    ! 1) Check that all mandatory variables have been defined correctly
    ! 2) Check that input is reasonable
    ! 3) Assign default values to optional allocatable variables.
    ! 4) Set variables for module to be used later

    if (.not.allocated(glob%run_type)) then
        allocate(glob%run_type(3))
        glob%run_type = [1,2,3]
    endif

    ! 1) Check that all mandatory variables have been defined correctly
    if (allocated(glob%ipar_init)) then
        nparam = size(glob%ipar_init)
    else
        call bad_input('ipar_init must be defined')
    endif

    if (allocated(glob%ipar_optim)) then
        noptim = size(glob%ipar_optim)
        ! Check if the indicies are within the bounds of the ipar_init array:
        if ((maxval(glob%ipar_optim)>nparam).or.(minval(glob%ipar_optim)<1)) then
            call bad_input('ipar_optim must contain indicies of the ipar_init array')
        endif
    elseif (all(glob%run_type/=1)) then  ! If not a single simulation
        call bad_input('ipar_optim must be defined unless for a single run')
    else                            ! Single simulation
        allocate(glob%ipar_optim(1))
        glob%ipar_optim(1) = 1
        noptim = 1
    endif

    if (allocated(glob%ipar_max)) then
        if (noptim/=size(glob%ipar_max)) then
            call bad_input('ipar_max must have the same number of elements as ipar_optim')
        endif
    elseif (all(glob%run_type/=1)) then  ! If not a single simulation
        call bad_input('ipar_max must be defined unless for a single run')
    else                            ! Single simulation
        allocate(glob%ipar_max(1))
        glob%ipar_max(1) = abs(glob%ipar_init(1))*2.d0 + 1.d0 ! Ensure larger (multiply and add to avoid numerical issues when scaling)
    endif

    if (allocated(glob%ipar_min)) then
        if (noptim/=size(glob%ipar_min)) then
            call bad_input('ipar_min must have the same number of elements as ipar_optim')
        endif
    elseif (all(glob%run_type/=1)) then  ! If not a single simulation
        call bad_input('ipar_min must be defined unless for a single run')
    else                            ! Single simulation
        allocate(glob%ipar_min(1))
        glob%ipar_min(1) = glob%ipar_init(1)  ! Don't decrement in case log or inverse scaling chosen and initial parameter is close to zero. max increased so they are not equal.
    endif

    ! 2) Check that input is reasonable
    ! Check that parameters are within bounds
    lower_bound_nok = any(glob%ipar_init(glob%ipar_optim)<(glob%ipar_min-numtol))
    upper_bound_nok = any(glob%ipar_init(glob%ipar_optim)>glob%ipar_max+numtol)
    if (lower_bound_nok.or.upper_bound_nok) then
        call bad_input('initial parameters must be within range')
    endif

    if (any(dbl_comp_array2(glob%ipar_max, glob%ipar_min))) then
        call bad_input('ipar_max cannot be equal to ipar_min, use ipar_optim to not optimize this/these parameters')
    endif

    if (glob%nstatv<1) then
        call bad_input('Number of state variables must be greater than 0')
    endif
    nstatv = glob%nstatv

    if (glob%opt_resnr<0) then
        call bad_input('opt_resnr cannot be negative')
    endif

    ! 3) Assign default values to optional allocatable variables. (Other
    !    optional variables are set to default values in the type definition)

    if (.not.allocated(glob%ipar_sctype)) then
        allocate(glob%ipar_sctype(size(glob%ipar_optim)))
        glob%ipar_sctype = DEF_ipar_sctype
    endif

    if (.not.allocated(glob%xvar_sets)) then
        allocate(mpar, source=glob%ipar_init)
        if (allocated(glob%ipar_sets)) then
            nsets = size(glob%ipar_sets, 2)
            allocate(glob%xvar_sets(noptim, nsets))
            do k1 = 1,size(glob%xvar_sets,2)
                mpar(glob%ipar_optim) = glob%ipar_sets(:,k1)
                call ipar_to_xvar(mpar, glob%xvar_sets(:,k1), glob)
            enddo
        else
            allocate(glob%ipar_sets(nparam, 1), glob%xvar_sets(noptim, 1))
            call ipar_to_xvar(mpar, glob%xvar_sets(:,1), glob)
        endif
    endif

    ! 4) Set variables for the module to be used later
    nlgeom = glob%nlgeom

    ! 5) Check user subroutines
    if (glob%user_scale) then
        call check_user_scaling(glob)
    endif

end subroutine

subroutine check_user_scaling(glob)
implicit none
    type(glob_typ)                  :: glob
    double precision, allocatable   :: ipar(:), xvar(:), mpar(:)

    allocate(ipar, source=glob%ipar_init)
    allocate(mpar(size(ipar)))
    allocate(xvar(size(glob%ipar_optim)))

    call ipar_to_xvar(ipar, xvar, glob)
    call xvar_to_ipar(ipar, xvar, glob)

    if (any(.not.dbl_comp_array2(ipar, glob%ipar_init))) then
        call bad_input('User scaling is not reversible')
    endif

    ipar(glob%ipar_optim) = glob%ipar_max
    call ipar_to_xvar(ipar, xvar, glob)
    if (any(.not.dbl_comp_array(xvar, 1.d0))) then
        call bad_input('User scaling from ipar_max doesn''t result in 1.d0')
    endif
    call xvar_to_mpar(mpar, xvar, glob)
    if (any(abs(mpar)>huge(1.d0)/10)) then
        call bad_input('User scaling at upper bound gives infinite parameters to umat')
    endif

    ipar(glob%ipar_optim) = glob%ipar_min
    call ipar_to_xvar(ipar, xvar, glob)
    if (any(.not.dbl_comp_array(xvar, 0.d0))) then
        call bad_input('User scaling from ipar_min doesn''t result in 0.d0')
    endif
    call xvar_to_mpar(mpar, xvar, glob)
    if (any(abs(mpar)>huge(1.d0)/10)) then
        call bad_input('User scaling at lower bound gives infinite parameters to umat')
    endif

end subroutine

subroutine check_sim_usr(usr_sim, simnr)
implicit none
    type(usr_sim_typ)   :: usr_sim
    integer         :: simnr
    if (usr_sim%lib=='') then
        call bad_input('A user library path name (*lib) must be specified when running user simulation. (sim nr. '//int2str(simnr)//')')
    endif

    if (.not.allocated(usr_sim%user_data)) then
        allocate(usr_sim%user_data(1))
        usr_sim%user_data = 1.d0
    endif

end subroutine

subroutine check_sim_ext(sim, simnr)
implicit none
    type(sim_typ)   :: sim
    integer         :: simnr
    if (sim%ext_cmd%script=='') then
        call bad_input('A script command is required when using external simulation. (sim nr. '//int2str(simnr)//')')
    endif

end subroutine

subroutine check_sim_atp(sim, simnr)
implicit none
    type(sim_typ)   :: sim
    integer         :: simnr
    nchannels = 4
    call check_sim_init(sim%init, simnr)

    call check_sim_exp(sim%exp)

    call check_sim_mesh1d(sim%mesh1d, sim%exp%ctrl, continued_analysis=sim%init%cont_analysis.ne.0)

    call check_sim_iter(sim%iter)

    call check_sim_err(sim%err)

    call check_sim_outp(sim%outp)

end subroutine

subroutine check_sim_atp_mr(sim, simnr)
implicit none
    type(sim_typ)       :: sim
    integer             :: simnr
    double precision    :: dummy_ctrl(5,1)

    nchannels = 4

    call check_sim_iter(sim%iter)

    dummy_ctrl(:,1) = (/1.d0, 1.d0, 1.d0, 1.d0, 1.d0/)
    call check_sim_mesh1d(sim%mesh1d, dummy_ctrl, continued_analysis=.true.)

    ! 2) Check that input is reasonable
    if (sim%init%cont_analysis==0) then
        call bad_input('cont_analysis input is mandatory for material removal analysis type. It can not be zero.')
    endif

    if (.not.allocated(sim%iter%time_incr)) then
        allocate(sim%iter%time_incr(1,5))
        sim%iter%time_incr(:,1) = (/1.d0, sim%atp_mr%time_relx, sim%atp_mr%time_relx, sim%atp_mr%time_relx, sim%atp_mr%time_relx/)
    endif

    call check_sim_err(sim%err)

    call check_sim_outp(sim%outp)

    if (.not.allocated(sim%exp%exp_scale)) then
        allocate(sim%exp%exp_scale(1))
        sim%exp%exp_scale = 1.d0
    endif

end subroutine

subroutine check_sim_mps(sim, simnr)
implicit none
type(sim_typ)   :: sim
integer         :: simnr

    if (nlgeom) then
        nchannels = 9
    else
        nchannels = 6
    endif

    call check_sim_exp(sim%exp)
    call check_sim_iter(sim%iter)

    ! 2) Check that input is reasonable
    ! Check initial state variables
    if (allocated(sim%init%statev_init)) then
        if (size(sim%init%statev_init,1)/=1) then
            call bad_input('statev_init can only have one row for mps simulation')
        elseif (size(sim%init%statev_init,2)/=nstatv) then
            call bad_input('statev_init must have the same number of elements as specified by nstatv')
        endif
    endif

    ! 3) Assign default values to optional allocatable variables.
    call check_sim_err(sim%err)

    call check_sim_outp(sim%outp)

end subroutine

subroutine check_sim_init(init, simnr)
implicit none
    type(init_typ)  :: init
    integer         :: simnr

    if (abs(init%cont_analysis)>=simnr) then
        call bad_input('continued analysis must be from an analysis defined before itself')
    endif

end subroutine check_sim_init

subroutine check_sim_exp(exp)
implicit none
type(exp_typ)   :: exp
integer         :: k1

    ! Check *ctrl
    if (allocated(exp%ctrl)) then
        if (size(exp%ctrl, 1)/=(nchannels+1)) then
            call bad_input('ctrl should have '//int2str(nchannels+1)//' columns')
        elseif (any(exp%ctrl(1,:)<0)) then
            call bad_input('The first column in ctrl represents steps and must be strictly positive')
        elseif (.not.dbl_compare(exp%ctrl(1,1),1.d0)) then
            if (exp%exp_info(1)==0) then
                call bad_input('The first row in ctrl should give the settings for the first load step')
            endif
        elseif (size(exp%ctrl, 2)>1) then
            k1 = size(exp%ctrl, 2)
            if (any( exp%ctrl(1,2:k1)<=exp%ctrl(1,1:(k1-1)) )) then
                call bad_input('The steps described by ctrl must be in increasing order')
            endif
        endif
    else
        call bad_input('ctrl must be defined')
    endif

    ! Check *exp_info
    if (allocated(exp%exp_info)) then
        if (size(exp%exp_info)/=2+2*nchannels+1) then
            call bad_input('exp_info should have '//int2str(3+2*nchannels)//' columns')
        endif
    else
        call bad_input('exp_info must be defined')
    endif

    ! Check *exp_data
    if (exp%exp_file=='') then
        call bad_input('exp_data file must be defined')
    endif

    ! Check *exp_scale
    if (.not.allocated(exp%exp_scale)) then
        allocate(exp%exp_scale(1))
        exp%exp_scale = 1.d0
    endif

end subroutine

subroutine check_sim_mesh1d(mesh1d, ctrl, continued_analysis)
implicit none
type(mesh1d_typ)    :: mesh1d
double precision    :: ctrl(:,:)
logical             :: continued_analysis

    if (allocated(mesh1d%node_pos)) then
        if (any(mesh1d%node_pos<0.d0)) then
            call bad_input('node_pos must be positive')
        elseif (size(mesh1d%node_pos)<2) then
            call bad_input('node_pos: At least two nodes must be specified')
        elseif (any(mesh1d%node_pos(2:size(mesh1d%node_pos))<=mesh1d%node_pos(1:(size(mesh1d%node_pos)-1)))) then
            call bad_input('node_pos must be monotonically strictly increasing')
        endif
        allocate(mesh1d%node_pos_undef, source=mesh1d%node_pos)
    elseif (.not.continued_analysis) then
        call bad_input('node_pos must be defined if analysis is not continued')
    endif

    if (dbl_compare(mesh1d%node_pos(1),0.d0,1.d-20)) then
        if (.not.all(dbl_comp_array(ctrl(4,:), 1.d0))) then
            call write_output('Displacement control to zero should be specified on inner radius if it is zero', 'warning', 'inp')
        endif
    endif

    ! Full gauss integration requires element_order=2*ngp-1
    if (mesh1d%element_order>(2*mesh1d%ngp - 1)) then
        call write_output('Reduced element integration is specified', 'warning', 'inp')
    elseif (mesh1d%element_order<(2*mesh1d%ngp - 1)) then
        call write_output('Over-integration of element is specified', 'status', 'inp')
    endif

    nnod = size(mesh1d%node_pos) + (mesh1d%element_order-1)*(size(mesh1d%node_pos)-1)
    nel = size(mesh1d%node_pos) - 1

end subroutine

subroutine check_sim_iter(iter)
implicit none
type(iter_typ)  :: iter
integer         :: k1

    if (allocated(iter%time_incr)) then
        do k1=1,size(iter%time_incr, 2)
            if (any(iter%time_incr(2:,k1)<1/huge(1.d0))) then
                call bad_input('time_incr: time increment must be greater than zero')
            elseif (any(iter%time_incr(3:4,k1)>iter%time_incr(5,k1))) then
                call bad_input('time_incr: dtmax cannot be smaller than dtmin or dt0')
            elseif (any(iter%time_incr(2:5,k1)<iter%time_incr(3,k1))) then
                call bad_input('time_incr: dtmin cannot be larger than any other increment input')
            endif
        enddo
    else
        call bad_input('time_incr must be defined')
    endif

end subroutine

subroutine check_sim_err(err)
implicit none
type(err_typ)   :: err

    if (allocated(err%err_scale)) then
        if (size(err%err_scale,1)/=(nchannels+1)) then
            call bad_input('*err_scale should have '//int2str(nchannels+1)//' columns')
        endif
    else
        allocate(err%err_scale((nchannels+1),1))
        err%err_scale(1,1) = 0.d0
        err%err_scale(2:,1) = DEF_err_scale
    endif

    if (allocated(err%err_scale_ctrl)) then
        if (size(err%err_scale_ctrl,1)/=(nchannels+1)) then
            call bad_input('*err_scale_ctrl should have '//int2str(nchannels+1)//' columns')
        endif
    else
        allocate(err%err_scale_ctrl((nchannels+1),1))
        err%err_scale_ctrl(1,1) = 0.d0
        err%err_scale_ctrl(2:,1) = 0.d0
    endif

    if (.not.allocated(err%err_steps)) then
        allocate(err%err_steps(size(DEF_error_steps)))
        err%err_steps = DEF_error_steps
    endif

    if (.not.allocated(err%user_settings)) then
        allocate(err%user_settings(1))
        err%user_settings(1) = 0.d0
    endif

    if (.not.(err%error_lib=='')) then
        err%error_type = 0  ! Set to user error if library has been specified
    endif

end subroutine

subroutine check_sim_outp(outp)
implicit none
    type(outp_typ)      :: outp
    character(len=50)   :: tmpstr
    integer             :: outpstatus
    integer             :: k1
    integer             :: default_comp(4), default_comp_sigma(4)

    if (.not.allocated(outp%result_steps)) then
        allocate(outp%result_steps(size(DEF_result_steps)))
        outp%result_steps = DEF_result_steps
    endif

    if (outp%dbl_format=='') then
        outp%dbl_format = '(ES25.15E3)'
    else
        write(tmpstr,outp%dbl_format, iostat=outpstatus) 1.0
        if (outpstatus/=0) then
            call write_output('The dbl_format specified "'//trim(outp%dbl_format)//'" has incorrect syntax.', 'error', 'inp', halt=.false.)
            call write_output('It should follow fortran output formats, e.g. "(ES15.5)"', 'error', 'inp', loc=.false.)
        endif
        call write_output('The dbl_format specified gives the number 1.0 written as:', 'status', 'inp')
        call write_output('"'//trim(tmpstr)//'"', 'status', 'inp', loc=.false.)
    endif

    if (allocated(outp%output_energies)) then
        if (any(outp%output_energies>4).or.any(outp%output_energies<=0)) then
            call write_output('output_energies cannot be smaller than 1 or larger than 4 (only sse, spd, scd and rpl possible)', 'error', 'inp')
        endif
    endif

    if (nlgeom) then
        default_comp = [1,2,3,5]    ! 23 shear component at position 5
    else
        default_comp = [1,2,3,6]    ! 23 shear component at position 6
    endif
    default_comp_sigma = [1,2,3,6]  ! Sigma is always in abaqus voigt ordering, i.e. linear geometry

    call add_outp_allocation(outp%output_nodes, outp%ur, [-nnod])
    call add_outp_allocation(outp%stress_comp, outp%stress, default_comp_sigma)
    call add_outp_allocation(outp%strain_comp, outp%strain, default_comp)
    call add_outp_allocation(outp%dfgrd_comp, outp%dfgrd, default_comp)
    call add_outp_allocation(outp%statev_comp, outp%statev, [-nstatv])

    if (.not.allocated(outp%output_elems)) then
        if (any([outp%stress, outp%strain, outp%dfgrd, outp%statev])) then
            allocate(outp%output_elems(nel))
            do k1=1,nel
                outp%output_elems(k1)=k1
            enddo
        endif
    endif

end subroutine

subroutine add_outp_allocation(outp_comp, outp_logical, default_outp_comp)
implicit none
    integer, allocatable    :: outp_comp(:)
    logical                 :: outp_logical
    integer                 :: default_outp_comp(:)

    integer                 :: k1

    if (allocated(outp_comp)) then
        outp_logical = .true.
    elseif (outp_logical) then
        if (default_outp_comp(1)<0) then
            allocate(outp_comp(abs(default_outp_comp(1))))
            do k1 = 1,abs(default_outp_comp(1))
                outp_comp(k1) = k1
            enddo
        else
            allocate(outp_comp(size(default_outp_comp)))
            outp_comp = default_outp_comp
        endif
    endif

end subroutine

subroutine check_opt_nlopt(opt)
implicit none
    type(opt_typ) :: opt

        ! Subroutine tasks:
    ! 1) Check that all mandatory variables have been defined correctly: no mandatory variables
    ! 2) Check that input is reasonable: Nothing to check
    ! 3) Assign default values to optional allocatable variables:

    if (.not.allocated(opt%start%initial_step)) then
        allocate(opt%start%initial_step(1))
        opt%start%initial_step = DEF_initial_step !Set to default 0 value to use nlopt default
    endif

    if (.not.allocated(opt%end_cond%xtol_rel)) then
        allocate(opt%end_cond%xtol_rel(1))
        opt%end_cond%xtol_rel = DEF_xtol_rel
    endif

    if (.not.allocated(opt%end_cond%xtol_abs)) then
        allocate(opt%end_cond%xtol_abs(1))
        opt%end_cond%xtol_abs = DEF_xtol_abs
    endif
end subroutine

subroutine check_opt_pspace(opt)
implicit none
type(opt_typ)       :: opt

end subroutine check_opt_pspace

subroutine check_opt_usr(usr_opt)
implicit none
    type(usr_opt_typ)   :: usr_opt

    if (usr_opt%lib=='') then
        call bad_input('A user library path name (*lib) must be specified when running user optimization')
    endif

    if (.not.allocated(usr_opt%user_data)) then
        allocate(usr_opt%user_data(1))
        usr_opt%user_data = 1.d0
    endif

end subroutine

subroutine check_overall(f_data, opt)
implicit none
    type(fdata_typ) :: f_data
    type(opt_typ)   :: opt(:)
    ! Internal variables
    integer         :: k1

    ! Check if required categories has been read
    if (any(f_data%glob%run_type==2)) then
        if (size(opt)==0) then
            call bad_input('Optimization settings must be defined when a run_type=2 is specified')
        endif
    endif

    ! Check if umat is read, if required (i.e. any stype >0)
    if (f_data%glob%umat_lib=='') then
        do k1=1,size(f_data%sim)
            if (f_data%sim(k1)%stype>0) then
                call bad_input('umat must be defined for simulation with stype = ' // int2str(f_data%sim(k1)%stype))
            endif
        enddo
    endif

end subroutine

end module

max_stars_repo_name: KnutAM/matmodfit

max_issues_repo_name: KnutAM/matmodfit

max_forks_repo_name: KnutAM/matmodfit


text:
      subroutine intrsc(x1,y1,x2,y2,x3,y3,x4,y4,x0,y0)
      a1=y2-y1
      b1=x1-x2
      a2=y4-y3
      b2=x3-x4
C
      d=a1*b2-a2*b1
C.....test: are the lines collinear to each other  ?
      if(abs(d).le.1.e- 8) goto 1
      c1=x2*y1-x1*y2
      c2=x4*y3-x3*y4
C.....compute the intersection point "x0,y0".
      x0=(c2*b1-b2*c1)/d
      y0=(c1*a2-c2*a1)/d
C.....test: does the point "x0,y0" belong to both line segments ?
      if(x0.lt.amax1(amin1(x1,x2),amin1(x3,x4))) goto 1
      if(x0.gt.amin1(amax1(x1,x2),amax1(x3,x4))) goto 1
      if(y0.lt.amax1(amin1(y1,y2),amin1(y3,y4))) goto 1
      if(y0.gt.amin1(amax1(y1,y2),amax1(y3,y4))) goto 1
      return
1     continue
C....in unnormal situation: a half of the element of the "screen".
      x0=0.5*(x3+x4)
      y0=0.5*(y3+y4)
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


! { dg-do run }
! { dg-options "-fdump-tree-original -fcoarray=lib -lcaf_single" }
! { dg-additional-options "-latomic" { target libatomic_available } }
!
! Some dependency-analysis check for coarray communication
!
integer, target, save :: A(10)[*]
integer, pointer :: P(:)
integer, save :: B(10)[*]

A = [1,2,3,4,5,6,7,8,9,10]
B = [1,2,3,4,5,6,7,8,9,10]
A(10:2:-1) = A(9:1:-1)[1] ! 0
B(10:2:-1) = B(9:1:-1)
if (any (A-B /= 0)) call abort

A = [1,2,3,4,5,6,7,8,9,10]
B = [1,2,3,4,5,6,7,8,9,10]
A(9:1:-1) = A(10:2:-1)[1] ! 1
B(9:1:-1) = B(10:2:-1)
if (any (A-B /= 0)) call abort

A = [1,2,3,4,5,6,7,8,9,10]
B = [1,2,3,4,5,6,7,8,9,10]
allocate(P(10))
P(:) = A(:)[1] ! 1
if (any (A-B /= 0)) call abort

A = [1,2,3,4,5,6,7,8,9,10]
B = [1,2,3,4,5,6,7,8,9,10]
allocate(P(10))
P(:) = B(:)[1] ! 0

A = [1,2,3,4,5,6,7,8,9,10]
B = [1,2,3,4,5,6,7,8,9,10]
A(1:5)[1] = A(3:7)[1] ! 1
B(1:5) = B(3:7)
if (any (A-B /= 0)) call abort
end

! { dg-final { scan-tree-dump-times "_gfortran_caf_get \\\(caf_token.0, \\\(integer\\\(kind=\[48\]\\\)\\\) parm.\[0-9\]+.data - \\\(integer\\\(kind=\[48\]\\\)\\\) a, 1, &parm.\[0-9\]+, 0B, &parm.\[0-9\]+, 4, 4, 0, 0B\\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_get \\\(caf_token.0, \\\(integer\\\(kind=\[48\]\\\)\\\) parm.\[0-9\]+.data - \\\(integer\\\(kind=\[48\]\\\)\\\) a, 1, &parm.\[0-9\]+, 0B, &parm.\[0-9\]+, 4, 4, 1, 0B\\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_get \\\(caf_token.0, \\\(integer\\\(kind=\[48\]\\\)\\\) parm.\[0-9\]+.data - \\\(integer\\\(kind=\[48\]\\\)\\\) a, 1, &parm.\[0-9\]+, 0B, &p, 4, 4, 1, 0B\\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_get \\\(caf_token.1, \\\(integer\\\(kind=\[48\]\\\)\\\) parm.\[0-9\]+.data - \\\(integer\\\(kind=\[48\]\\\)\\\) b, 1, &parm.\[0-9\]+, 0B, &p, 4, 4, 0, 0B\\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_sendget \\\(caf_token.0, \\\(integer\\\(kind=\[48\]\\\)\\\) parm.\[0-9\]+.data - \\\(integer\\\(kind=\[48\]\\\)\\\) a, 1, &parm.\[0-9\]+, 0B, caf_token.0, \\\(integer\\\(kind=\[48\]\\\)\\\) parm.\[0-9\]+.data - \\\(integer\\\(kind=\[48\]\\\)\\\) a, 1, &parm.\[0-9\]+, 0B, 4, 4, 0, 0B\\\);" 1 "original" } }


text:
Address(Cypress Lane) is a residential street that runs from Simmons Way at its southern end to Claremont Drive at its northern end.

Along the way it passes Tulane Drive, Alice Street, Union Drive, and Drexel Drive. It also has a bike path connecting it to Spruce Lane.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
 !
 ! This file is part of the ExaHyPE project.
 ! Copyright (c) 2016  http://exahype.eu
 ! All rights reserved.
 !
 ! The project has received funding from the European Union's Horizon
 ! 2020 research and innovation programme under grant agreement
 ! No 671698. For copyrights and licensing, please consult the webpage.
 !
 ! Released under the BSD 3 Open Source License.
 ! For the full license text, see LICENSE.txt
 !

SUBROUTINE ADERVolumeIntegralNonlinear(lduh,lFhi_x,lFhi_y,lFhi_z,lShi,dx)
    USE typesDef
    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    ! Argument list
    REAL, INTENT(IN)              :: lFhi_x(nVar,nDOF(1),nDOF(2),nDOF(3))    ! nonlinear flux tensor in each space-time DOF in x direction
    REAL, INTENT(IN)              :: lFhi_y(nVar,nDOF(2),nDOF(1),nDOF(3))    ! nonlinear flux tensor in each space-time DOF in y direction
    REAL, INTENT(IN)              :: lFhi_z(nVar,nDOF(3),nDOF(1),nDOF(2))    ! nonlinear flux tensor in each space-time DOF in z direction
    REAL, INTENT(IN)              :: lShi(nVar,nDOF(1),nDOF(2),nDOF(3))      ! nonlinear source term
    REAL, INTENT(OUT)             :: lduh(nVar,nDOF(1),nDOF(2),nDOF(3))      ! spatial degrees of freedom
    DOUBLE PRECISION, INTENT(IN)  :: dx(d)                                   ! mesh spacing
    ! Local variables
    INTEGER           :: i,j,k,l
    REAL              :: aux(d)
    !
    ! Initialize the update DOF
    lduh = 0.
    ! x - direction
    DO k = 1, nDOF(3)
        DO j = 1, nDOF(2)
            aux = (/ 1.d0, wGPN(j), wGPN(k) /)
            lduh(:,:,j,k) = lduh(:,:,j,k) + MATMUL( lFhi_x(:,:,j,k), TRANSPOSE(Kxi) )*PRODUCT(aux(1:nDim))/dx(1)
        ENDDO
    ENDDO
    !
    ! y - direction (not needed in 1D)
    IF(nDim>=2) THEN
        DO k = 1, nDOF(3)
            DO i = 1, nDOF(1)
                aux = (/ 1.d0, wGPN(i), wGPN(k) /)
                lduh(:,i,:,k) = lduh(:,i,:,k) + MATMUL( lFhi_y(:,:,i,k), TRANSPOSE(Kxi) )*PRODUCT(aux(1:nDim))/dx(2)
            ENDDO
        ENDDO
    ENDIF
    !
    ! z - direction (node needed in 1D and 2D)
    IF(nDim>=3) THEN
        DO j = 1, nDOF(2)
            DO i = 1, nDOF(1)
                aux = (/ 1.d0, wGPN(i), wGPN(j) /)
                lduh(:,i,j,:) = lduh(:,i,j,:) + MATMUL( lFhi_z(:,:,i,j), TRANSPOSE(Kxi) )*PRODUCT(aux(1:nDim))/dx(3)
            ENDDO
        ENDDO
    ENDIF
    !
    DO k = 1, nDOF(3)
        DO j = 1, nDOF(2)
            DO i = 1, nDOF(1)
                aux = (/ wGPN(i), wGPN(j), wGPN(k) /)
                lduh(:,i,j,k) = lduh(:,i,j,k) + PRODUCT(aux(1:nDim))*lShi(:,i,j,k)
            ENDDO
        ENDDO
    ENDDO
    !
    !PRINT *, 'lFhi_x=', lfhi_x
    !PRINT *, 'lFhi_y=', lfhi_y
    !PRINT *, 'lFhi_z=', lfhi_z
    !PRINT *, 'lFhi_t=', lshi
    CONTINUE
    !
    !OPEN(UNIT=12, FILE="aoutput_lduh.txt", ACTION="write", STATUS="replace")
    !WRITE(12, '(ES24.16,1x)') , lduh
    !CALL EXIT

END SUBROUTINE ADERVolumeIntegralNonlinear

max_stars_repo_name: linusseelinger/ExaHyPE-Tsunami

max_issues_repo_name: linusseelinger/ExaHyPE-Tsunami

max_forks_repo_name: linusseelinger/ExaHyPE-Tsunami


! { dg-do compile }
! { dg-options "-std=f2008" }
!
! PR fortran/45170
!
! Character deferred type parameter
!

subroutine one(x, y) ! { dg-error "Entity .y. at .1. has a deferred type parameter" }
  implicit none
  character(len=:), pointer :: x
  character(len=:) :: y
  character(len=:), allocatable, target :: str2
  character(len=:), target :: str ! { dg-error "deferred type parameter" }
end subroutine one

subroutine two()
  implicit none
  character(len=:), allocatable, target :: str1(:)
  character(len=5), save, target :: str2
  character(len=:), pointer :: pstr => str2
  character(len=:), pointer :: pstr2(:)
end subroutine two

subroutine three()
!  implicit none  ! Disabled because of PR 46152
  character(len=:), allocatable, target :: str1(:)
  character(len=5), save, target :: str2
  character(len=:), pointer :: pstr
  character(len=:), pointer :: pstr2(:)

  pstr => str2
  pstr2 => str1
  str1 = ["abc"]
  pstr2 => str1

  allocate (character(len=77) :: str1(1))
  allocate (pstr, source=str2)
  allocate (pstr, mold=str2)
  allocate (pstr) ! { dg-error "requires either a type-spec or SOURCE tag" }
  allocate (character(len=:) :: str1(1)) ! { dg-error "cannot contain a deferred type parameter" }

  str1 = [ character(len=2) :: "abc" ]
  str1 = [ character(len=:) :: "abc" ] ! { dg-error "cannot contain a deferred type parameter" }
end subroutine three

subroutine four()
  implicit none
  character(len=:), allocatable, target :: str
  character(len=:), pointer :: pstr
  pstr => str
  str = "abc"
  if(len(pstr) /= len(str) .or. len(str)/= 3) STOP 1
  str = "abcd"
  if(len(pstr) /= len(str) .or. len(str)/= 4) STOP 2
end subroutine four

subroutine five()
character(len=4) :: str*(:)
allocatable :: str
end subroutine five