c		program addgrd
c  node by node algebraic operations between two grids

	parameter  (maxcol=4000, mnrecl=30 )
	dimension  r(maxcol), r2(maxcol)

	character*56 ifile, i2file, ofile, constant
	character*56 id, id2, id3, b56
        character    opr*1, oopr*1
        data         dval/0.1e+39/, constant/'constant'/,
     1	     inp/12/, inp2/13/, jput/14/
c        cr = char( 13 )
c        lf = char( 10 )
        call askin
        call pfinit('addgrd')

	opr    = ' '
	ifile  = ' '
	i2file = ' '
	ofile  = ' '
	id     = ' '
	id2    = ' '
	id3    = ' '
	b56    = ' '

1       ifile=' '
        call askc('First input grid file',ifile,ie)
        if(ie.eq.-2) go to 999
        if(ie.ne.0) go to 1
        call gopen(inp,ifile,'old','read',ierr)
        if(ierr.ne.0) then
          print*,'Error - try again'
          go to 1
        endif
        call gheader('r',inp,id,nc,nr,xo,dx,yo,dy,ierr)
        if(ierr.ne.0) stop 'Error reading grid header'
        print*, 'title=', id
	print*,     ' nc X nr        =', nc, nr
        print*,     ' xo,yo,dx,dy =', xo, yo, dx, dy
        if ( nc .gt. maxcol ) then
	  print*, ' maximum number of columns=', maxcol
	  stop
        endif
10      opr='+'
        call askc('Operator  h(elp), + - * / m % ',opr,ie)
        if(ie.eq.-2) then
          close(inp)
          go to 1
        endif
        if ( opr .eq. 'h' ) then
          print*, ' arithmetic operators: file1 (opr) file2 = ofile'
	  print*, ' masking:  blank areas in 2 overwrite 1'
	  print*, ' % difference =  100 * ( f1 - f2 ) / f1'
	  go to 10
	endif
        if ( .not. ( opr.eq.'+' .or. opr.eq.'-' .or. opr.eq.'*' .or.
     1	opr.eq.'/' .or. opr.eq.'m' .or. opr.eq.'%' ) )
     1	then
	  print*, ' error, please reenter operator'
	  go to 10
        endif
15      i2file='constant'
        call askc('Second input grid file',i2file,ie)
        if(ie.eq.-2) go to 10
        if ( i2file(1:1) .eq. ' ' ) then
	  print*, ' the 2nd grid will be replaced by a constant value'
	  i2file = 'constant'
        endif
	icon = 0
	if ( i2file .eq. constant ) then
	  icon = -1
          zc=0.0
          if(opr.eq.'/') zc=1.0
          call askf4('Enter constant',zc,ie)
          if(ie.eq.-2) go to 15
          oopr=' '
	  if ( opr .eq. '/' ) then
	    opr = '*'
	    zc  = 1.0 / zc
            oopr='/'
	  endif
          do 50 i = 1, nc
	    r2(i) = zc
50        continue
        else
c  open 2nd file and check
          call gopen(inp2,i2file,'old','read',ierr)
          if(ierr.ne.0) then
            print*,'Error - try again'
            go to 15
          endif
          call gheader('r',inp2,id2,nc2,nr2,xo2,dx2,yo2,dy2,ierr)
          if(ierr.ne.0) stop 'Error reading grid header'
c  different sizes are fatal
          print*, 'title=', id2
          if ( nc .ne. nc2 ) stop ' ncolumns differ'
          if ( nr .ne. nr2 ) stop ' nrows differ'
c  warn of coordinate problems
          if ( abs ( ( dx - dx2 ) / dx ) .gt. .001 )
     1    print*, ' dx differ'
          if ( abs ( ( dy - dy2 ) / dy ) .gt. .001 )
     1    print*, ' dy differ'
          tx1 = xo2 - dx * .1
          tx2 = xo2 + dx * .1
          tmx = amax1( tx1, tx2 )
          tmn = amin1( tx1, tx2 )
          if ( xo .lt. tmn  .or.  xo .gt. tmx ) then
            print*, ' xo1, xo2', xo,xo2
            print*, ' x origin different'
          endif
          ty1 = yo2 - dy * .1
          ty2 = yo2 + dy * .1
          tmx = amax1( ty1, ty2 )
          tmn = amin1( ty1, ty2 )
          if ( yo .lt. tmn  .or.  yo .gt. tmx ) then
            print*, ' yo1, yo2', yo,yo2
            print*, ' y origin different'
          endif
        endif
c  end of check
20      ofile=' '
        if(opr.eq.'m') then
          ofile=ifile(1:index(ifile,'.'))//'msk'
        else if(opr.eq.'+') then
          ofile=ifile(1:index(ifile,'.'))//'sum'
        else if(opr.eq.'-') then
          ofile=ifile(1:index(ifile,'.'))//'dif'
        else if(opr.eq.'*') then
          ofile=ifile(1:index(ifile,'.'))//'prd'
        else if(oopr.eq.'/') then
          ofile=ifile(1:index(ifile,'.'))//'div'
        else if(opr.eq.'%') then
          ofile=ifile(1:index(ifile,'.'))//'pct'
        endif
        call askc('Output grid file',ofile,ie)
        if(ie.eq.-2) then
          close(inp2)
          go to 15
        endif
c        id3=b56
        call askalt()
        call askc('Output title',id,ie)
        call askalt()
        if(ie.eq.-2) go to 20
	length = nc + 1
	if ( length .lt. mnrecl ) length = mnrecl
        call gopen(jput,ofile,'new','write',ierr)
        if(ierr.ne.0) then
          print*,'Error - try again'
          go to 20
        endif
        call gheader('w',jput,id,nc,nr,xo,dx,yo,dy,ierr)
        if(ierr.ne.0) stop 'Error writing output grid header'
	ig = 0
	if ( opr .eq. '+' ) ig = 1
	if ( opr .eq. '-' ) ig = 2
	if ( opr .eq. '*' ) ig = 3
	if ( opr .eq. '/' ) ig = 4
	if ( opr .eq. 'm' ) ig = 5
	if ( opr .eq. '%' ) ig = 6
c
	do 100 j = 1, nr
        call grow('r',inp,dum,r,nc,ierr)
        if(ierr.ne.0) stop 'Error reading grid1 row'
        if ( icon .eq. 0 ) then
          call grow('r',inp2,dum,r2,nc2,ierr)
          if(ierr.ne.0) stop 'Error reading grid2 row'
        endif
	go to( 41, 42, 43, 44, 45, 46 ) ig
c
41      do 51 i = 1, nc
        if ( r(i) .ge. 1.e30  .or.  r2(i) .ge. 1.e30 ) then
          r(i) = dval
        else
          r(i) = r(i) + r2(i)
        endif
51      continue
	go to 90
c
42      do 52 i = 1, nc
        if ( r(i) .ge. 1.e30  .or.  r2(i) .ge. 1.e30 ) then
          r(i) = dval
        else
          r(i) = r(i) - r2(i)
        endif
52      continue
	go to 90
c
43      do 53 i = 1, nc
        if ( r(i) .ge. 1.e30  .or.  r2(i) .ge. 1.e30 ) then
          r(i) = dval
        else
          r(i) = r(i) * r2(i)
        endif
53      continue
	go to 90
c
44      do 54 i = 1, nc
        if ( r(i) .ge. 1.e30.or.  r2(i) .ge. 1.e30 .or.
     1  r2(i) .eq. 0 ) then
          r(i) = dval
        else
          r(i) = r(i) / r2(i)
        endif
54      continue
	go to 90
c
c  masking
45      do 55 i = 1, nc
        if ( r2(i) .ge. 1.e30 ) r(i) = dval
55      continue
	go to 90
c
c  % difference
46      do 56 i = 1, nc
        if ( r(i) .ge. 1.e30      .or.    r2(i) .ge. 1.e30 .or.
     1  r(i) .eq. 0 ) then
          r(i) = dval
        else
          r(i) = 100 * ( r(i) - r2(i) ) / r(i)
        endif
56      continue
c
90      call grow('w',jput,dum,r,nc,ierr)
        if(ierr.ne.0) stop 'Error writing output row'
100	continue
c
        call gclose(inp,'keep')
        call gclose(jput,'keep')
        if ( icon .eq. 0 ) call gclose(inp2,'keep')
999     end
