c  HAMMER
c  This program calculates terrain corrections according to Hammer's 
c  method, and is designed to accomodate elevations not included on 
c  Hammer's chart.
c  The Hammer zones are as follows (where r1 and r2 are in the inner
c  and outer radii of the zones, n is the no. of compartments in that
c  zone, units all in feet)
c
c	zone      r1       r2       n
c	 A        0       6.56      1
c	 B       6.56     54.6      4
c        C       54.6     175       6
c	 D       175      558       6
c	 E       558      1280      8
c	 F       1280     2936     12
c	 G       2936     5018     12
c	 H       5018     8578     12
c
c  For any other zone user must enter r1, r2  and n
c Tien Grauch
c
	character zone*1,word*16,sta*8
	dimension tc(50),sumtc(8)
	g=6.67e-03
	pi=3.141592654
	fac=.3048006e0
	write(*,'(a,$)')' Want help?  '
	read 800, zone
	if(zone.eq.'y'.or.zone.eq.'Y') then
		print 798
		print 799
	endif
798	format(' This program calculates the tc according to ',
     1 'Hammer''s method, and'/
     2 ' is designed to accomodate heights not included on ',
     3 'Hammer''s chart.'/' Program gives values of specified ',
     4 'zones and keeps a running total.'/
     5' The Hammer zones are as follows (where r1 and r2 are ',
     6'the inner and'/
     7' outer radii of the zones, n is the no. of ',
     8'compartments in that'/
     9' zone, units all in feet)'/)
799	format(10x,'	zone      r1       r2       n'/
     110x,'	 A        0       6.56      1'/
     210x,'	 B       6.56     54.6      4'/
     310x,'	 C       54.6     175       6'/
     410x,'	 D       175      558       6'/
     510x,'	 E       558      1280      8'/
     610x,'	 F       1280     2936     12'/
     710x,'	 G       2936     5018     12'/
     810x,'	 H       5018     8578     12'//
     9'    For any other zone user must enter r1, r2 and n.'/)

	print 901
901	format(' Output TC''s to tc.tbl file? ',$)
	read 800,zone
	iout=1
	if(zone.eq.'y'.or.zone.eq.'Y') then
	print 922
922	format(1x,' (Only Hammer zones A-H are supported by this option)')
		iout=2
		open(12,file='tc.tbl',form='formatted',status='unknown')
	endif
	write(*,'(a,$)')' Enter density for tc''s  '
	read*,den
	if(iout.eq.2) write(12,911) den
911	format('sta id  elev(ft) TC(',f4.2,')  A',6x,'B',6x,'C',
     1 6x,'D',6x,'E',6x,'F',6x,'G',6x,'H')
c reset for new station
11	do 5 i=1,8
	sumtc(i)=0.0
5	continue
	write(*,'(a,$)')' Enter station elev in ft (car ret to exit) '
	read 900,word
900	format(a16)
	if(word.eq.' ') go to 300   
        call shiftr(word) 
	read(word,8000) elev
	elev=elev*fac
	c1=2.0*pi*g*den
	sumall=0.0
	if(iout.eq.2) then
		print 902
902		format(' Enter station id (8 chars max): ',$)
		read 903, sta
903		format(a8)
	endif
1	print 813
813	format(' Enter zone (car ret for new station): '$)
	read 800, zone
800	format(a1)
	if(zone.eq.' ') then
		if(iout.eq.1) go to 11
		write(12,915) sta,elev/fac,sumall,sumtc
915	format(a8,f8.0,9f7.3)
		go to 11
	endif
c  internally work in meters, user enters feet
	if(zone.eq.'a'.or.zone.eq.'A') then
	  r1=0.
	  r2=1.99949
	  n=1
	  izone=1
	  go to 50	
	endif
	if(zone.eq.'b'.or.zone.eq.'B') then
	  r1=1.99949
	  r2=16.64
	  n=4
	  izone=2
	  go to 50
	endif
	if(zone.eq.'c'.or.zone.eq.'C') then
	  r1=16.64
	  r2=53.34
	  n=6
	  izone=3
	  go to 50
	endif
	if(zone.eq.'d'.or.zone.eq.'D') then
	  r1=53.34
	  r2=170.1
	  n=6
	  izone=4
	  go to 50
	endif
	if(zone.eq.'e'.or.zone.eq.'E') then
	  r1=170.1
	  r2=390.1
	  n=8
	  izone=5
	  go to 50
	endif
	if(zone.eq.'f'.or.zone.eq.'F') then
	  r1=390.1
	  r2=895.
	  n=8
	  izone=6
	  go to 50
	endif
	if(zone.eq.'g'.or.zone.eq.'G') then
	  r1=895.
	  r2=1530.
	  n=12
	  izone=7
	  go to 50
	endif
	if(zone.eq.'h'.or.zone.eq.'H') then
	  r1=1530.
	  r2=2615.
	  n=12
	  izone=8
	  go to 50
	endif
	if(iout.eq.2) then
		print*,'Can''t output customized zones in tc.tbl'
		go to 1
	endif
	print*,'Enter inner and outer radius in ft and no. of
     & compartments'
	read*,r1,r2,n
	izone=1
	r1=r1*fac
	r2=r2*fac
50	r12=r1*r1
	r22=r2*r2
	c2=r2-r1
	c3=c1/float(n)
	print 810,n
810	format(' Enter ave. elev for the following ',i2,' compartments.'/
     1  ' Give car ret to correct previous entry'/)
	k=0
52	k=k+1
	if(k.gt.n) go to 80
55	print 801,k
801	format(' Ave. elev of compartmt ',i2,': ',$)
	read 900,word
	if(word.eq.' ') then
		if(k.eq.1) go to 1
		sumtc(izone)=sumtc(izone)-tc(k-1)
		k=k-1
		go to 55
	endif
        call shiftr(word)
	read(word,8000) h
	h=abs(h*fac-elev)
	h2=h*h	
	tcn=c2+sqrt(r12+h2)-sqrt(r22+h2)
	tc(k)=c3*tcn
	print 802,tc(k)
802	format(10x,f7.3,' mgals')
	sumtc(izone)=sumtc(izone)+tc(k)
	go to 52
80	sumall=sumall+sumtc(izone)
	print*
	print 803,sumtc(izone),sumall
803	format(' TC of this zone = ',f7.2,' in mgals'/' Running total
     1 for this station = ',f7.2,' in mgals'/)
	go to 1
300	if(iout.eq.2) close(12)
	stop          
8000    format(g16.0)
	end                 
        subroutine shiftr(a)
c       shift list directed real input to right end of field for internal read
        character*16 a,c
        character*1 b(16)
        equivalence(c,b)
        c=a
        do 10 i=16,1,-1
        if(b(i) .ne. ' ') goto 11
10      contiNUE
        I=1
11      ishft=16-i
        a=' '
        a(ishft+1:16)=c(1:i)
        return
        end

