
c               program tilt
c  surface generation with specified edge gradients
c  f(x,y) = a*dx + b*dy + c*dx*dy + dc
	dimension s(4),g(4),p(4)
	dimension r(1500),r2(1500)
	character*56 f1,f2,f3,id
	character pgm*8,ext1*4,ext2*4
        character*3 ans
	data dval,fltbig/0.1701412e39, 1.e29/, ext1,ext2/'.tlt', '.srf'/
	write(*,10)
10	format(' enter input grid filename :'$)
	read(*,37) f1
37	format(a56)
	open (10,file=f1,form='unformatted',status='old',mode='read')
	write(*,11)
11	format(' enter output prefix :'$)
82	read(*,37) f1
        do 80 ii=1,56
        if (f1(ii:ii) .ne. ' ') goto 81
80      continue
        write(*,*) ' Null prefix, re-enter...'
        goto 82
c81	n = leftj(f1)
81      continue
        kk=ii
        do 84 ii=kk,56
        if (f1(ii:ii).eq.' ') goto 86
84      continue
86        jj=ii-1
        f2=f1(kk:jj)
        f1=f2
        n=jj-kk+1
	np = index(f1(1:n),'.') 
	if (np.gt.0) n = np-1
	f2 = f1(1:n)//ext1
	open(11,file=f2,form='unformatted',status='unknown')
	f3 = f1(1:n)//ext2
	open(12,file=f3,form='unformatted',status='unknown')
	write(*,*)'  tilted data grid filename is ',f2(1:n+4)
	write(*,*)' generated surface filename is ',f3(1:n+4)
c
	read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
	write(*,*)'title = ',id
	write(*,*)'ncol,nrow = ',nc,nr
	write(*,*)'xo,yo, dx,dy = ',xo,yo,dx,dy
	wx = dx*float(nc-1)
	wy = dy*float(nr-1)
c
19	write(*,20)
20	format(' help specifying gradients(n/y) ? ',$)
        read(*,'(a)')ans
        if (index(ans,'y').gt.0.or.index(ans,'Y').gt.0) then 
c	if(noyes(ldum).eq.1) then
	write(*,*)' Gradient locations are counterclockwise from the'
	write(*,*)
     1' lower-left: g(1) is bottom, g(2) is right side, etc.'
	write(*,*)' in z units per x units (ie. gammas per kilometer).'
	write(*,*)' The dependent gradient is 999.'
	write(*,*)' DC is the value of the lower-left corner.'
	endif
	write(*,*)' enter dc, and gradients 1 to 4 :'
	read(5,*,end=19,err=19) p(1),g
c
c  derive dependent gradient
	i2=0
	do 90 i=1,4
	  if(abs(g(i)-999.0).lt.0.01) i2=i
90      continue
	if(i2.eq.0) then
	  write(*,*)' one gradient must equal 999'
	  go to 19
	endif
	g(i2) = 0.
	s(1) = g(1)*wx
	s(2) = g(2)*wy
	s(3) = g(3)*wx
	s(4) = g(4)*wy
	s(i2) = (s(1)+s(2)) - (s(3)+s(4))
	if (i2.le.2) s(i2) = -s(i2)
	g(i2) = s(i2)/wx
	if (i2.eq.2 .or. i2.eq.4) g(i2) = s(i2)/wy
60	write(*,70)g
70	format(' edge gradients: ',1p4g11.4)
c
c  calculate corner values
	p(2) = g(1)*wx + p(1)
	p(3) = g(4)*wy + p(1)
	p(4) = g(3)*wx + p(3)
	chk  = g(2)*wy + p(2)
	if (abs(p(4)-chk).gt.1.e-2)
     1 write(*,*)'closure error', p,chk
c  surface coefficients
	a = (p(2)-p(1))/wx
	b = (p(3)-p(1))/wy
	c = (p(4)+p(1)-p(2)-p(3)) / (wx*wy)
	write(*,*)
     1' surface coordinate system origin (0,0) at lower left'
	write(*,*) 
     1' grid corner, x & y in horizontal data units'
	write(*,*)
     1 'surface coefficients f(x,y) = a*x + b*y + c*x*y + dc'
	write(*,*) a,b,c,p(1)
c
	pgm='tilt'
	write(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
	write(id,150) a,b,c,p(1)
150	format('ax+by+cxy+d ',1p4g11.3)
	write(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
c
	do 110 j = 1, nr
	call rowio(nc,r,-1,10,10,ie)
	y = dy*float(j-1)
	byp = b*y + p(1)
	cy  = c*y
	do 100 i = 1, nc
	  x = dx*float(i-1)
	  r2(i) =  x*(a+cy) + byp  	
	  if (r(i).lt.fltbig) then
	    r(i) = r(i) + r2(i)
	    else
	    r(i) = dval
	  endif
100     continue
	call rowio(nc,r,0,11,11,ie)
	call rowio(nc,r2,0,12,12,ie)
110     continue
c
	close(11)
	close(12)
	close(13)
	stop
	end
	subroutine rowio(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
c******************************************************************************
	function leftj(a)
c  left justifies a string and returns the position
c  of the last nonblank character
	character a*(*)
	n = len(a)
	ich = ichar(a(1:1))
	if (ich.ne.0 .and. ich.ne.32) go to 15
	do 1 m = 2, n
	ich = ichar(a(m:m))
1	if (ich.ne.0 .and. ich.ne.32) go to 5
	leftj = 0
	return
c
5	i2 = 0
	do 10 i = m, n
	i2 = i2 + 1
	a(i2:i2) = a(i:i)
10	continue
	do 11 i3 = i2+1, n
11	a(i3:i3) = ' '
	n = n-m+1
c
15	do 20 leftj = n, 1, -1
	ich = ichar(a(leftj:leftj))
20	if (ich.ne.0 .and. ich.ne.32) go to 25
25	return
	end
