program gcomb
      implicit double precision(a-h,o-z)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                          g c o m b
c
c  program for combining two grids into one by various manipulations.
c  the grids must have the same grid spacings and relative positions,
c  but not neccessarily the same coverage area.
c  the output grid will cover the same area as the first grid
c  9999-values (unknown) are not added/subtracted
c
c  input:
c
c  gridfile1
c  gridfile2
c  outfile
c  mode, lint
c
c  where mode:
c
c  mode = 0:  one-grid modification. Additional input: factor, bias
c             (new  = old*factor + bias)
c
c  mode = 1:  subtract 'grid1' minus 'grid2'
c
c  mode = 2:  add 'grid1' plus 'grid2'
c
c  mode = 3:  convert bouguer to free-air
c             'grid1' is a bouguer anomaly grid, 'grid2' a height grid.
c
c  mode = 4:  treshold values. set values in 'grid1' to 9999 when
c             values in 'grid2' are greater than 'treshold' (input after lint).
c
c  mode = 5:  grid overwrite
c             values in 'grid2' is written on top of 'grid1', except when
c             9999 is encountered, then the grid1-values are kept.
c
c  mode = 6:  grid select
c             values in 'grid1' are written only when there is no data in
c             'grid2'.
c
c  mode = 7:  Treshold select. Keep values in grid1 when grid1 > treshold,
c             keep values in grid2 when both grid1 and grid2 < trsh
c
c  mode = 8:  Variable multiplication. Multiply grid 1 by "fak1" when
c             grid 2 .lt. 'treshold', else multiply by "fak2"
c  mode = 9:  Grid division. Divide grid 1 by grid 2 in each cell
c
c  special option:
c  if 'grid1' = '0' then a label must be input, and grid1 will be
c  all zeros. this is useful for e.g. extending a grid with a border
c  of zero values.
c  if 'grid1' = '9999' as above, except unknown-9999 values are used.
c
c  rene forsberg, thessaloniki/unsw nov 1988, last updated sep 1989
c  updated and changed dec 89, feb 92, rf
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension glab1(6),glab2(6)
      dimension grow1(20000),grow2(20000) 
      character*80 ifile1,ifile2,ofile
      logical lzero, l9999, llab, lint, lutm
c
      write(*,3)
      lutm = .false.
3     format(/
     .' ************************************************************',
     .'********'/
     .' *     GCOMB - GRAVSOFT grid combination - vers. feb 92 (c) R',
     .'F/KMS  *'/
     .' ************************************************************',
     .'********')
      write(*,1)
1     format(' input file names - file1/file2/ofile: ')
      read(*,2) ifile1
      read(*,2) ifile2
      read(*,2) ofile
2     format(a80)
c
      write(*,4)
4     format(
     .' input: mode (1 dif, 2 sum, 3 fa, 4 trsh, 5 overwr..), lint:')
      read(*,*) mode, lint
c
      if (mode.eq.4) write(*,5)
5     format(' input grid2 treshold values for setting grid1 to 9999:')
      if (mode.eq.7) write(*,502)
502   format(' input treshold select value: ')
      if (mode.eq.4.or.mode.eq.7) read(*,*) trsh
      if (mode.eq.8) write(*,*) 'input treshold, fak1, fak2:'
      if (mode.eq.8) read(*,*) trsh,fak1,fak2
      if (mode.eq.0) write(*,501)
501   format(' input factor and bias to be added: ')
      if (mode.eq.0) read(*,*) fact,bias
c
      lzero = (ifile1.eq.'0')
      l9999 = (ifile1.eq.'9999')
      llab = (lzero.or.l9999)
      if (llab) write(*,6)
6     format(' input wanted label:')
      if (llab) read(*,*) (glab1(j),j=1,6)
c     if (llab) 
c    .lutm = (abs(glab1(4)).gt.100.or.abs(glab1(5)).gt.100)
      if (llab.and.lutm) write(*,7)
7     format(' input ellipsoid and zone:')
      if (llab.and.lutm) read(*,*) iell,izone
c
      if (.not.llab) open(10,file=ifile1,status='old')
      if (mode.gt.0) open(11,file=ifile2,status='old')
      open(20,file=ofile,status='unknown')
c
      if (.not.llab) then
        read(10,*) glab1
c       write(*,*) glab1
        lutm = (abs(glab1(4)).gt.500.or.abs(glab1(5)).gt.500)
        if (lutm) read(10,*) iell,izone
      endif
      if (mode.gt.0) then
        read(11,*) glab2
        if (lutm) then
          read(11,*) iell2,izone2
          if (iell.ne.iell2.or.izone.ne.izone2)
     .    stop 'utm spec difference'
        endif
        if (abs(glab1(5)-glab2(5)).gt.0.00001.or.
     .  abs(glab1(6)-glab2(6)).gt.0.00001)
     .  stop 'grid spacings different'
      else
        do 9 j = 1, 6
9       glab2(j) = glab1(j)
      endif
c
      if (.not.lutm) write(20,11) (glab1(j),j=1,6)
11    format(' ',6f12.6)
      if (lutm) write(20,12) (glab1(j),j=1,6),iell,izone
12    format(' ',6f12.1,/,' ',2i4)
c
      rfi0 = glab1(1)
      rla0 = glab1(3)
      dfi = glab1(5)
      dla = glab1(6)
      nn = nint((glab1(2)-rfi0)/dfi + 1)
      ne = nint((glab1(4)-rla0)/dla + 1)
      nn2 = nint((glab2(2)-glab2(1))/dfi + 1)
      ne2 = nint((glab2(4)-glab2(3))/dla + 1)
      in1 = nint((glab2(1)-rfi0)/dfi + 1)
      in2 = nint((glab2(2)-rfi0)/dfi + 1)
      ie1 = nint((glab2(3)-rla0)/dla + 1)
      ie2 = nint((glab2(4)-rla0)/dla + 1)
c
      r1min = 9999.99 
      r1max = -r1min
      nr1 = 0
      r1sum = 0.0
      r2sum2 = 0.0
      r2min = 9999.99 
      r2max = -r2min
      nr2 = 0
      r2sum = 0.0
      r2sum2 = 0.0
      rrmin = 9999.99 
      rrmax = -rrmin
      nrr = 0
      rrsum = 0.0
      rrsum2 = 0.0
c
      write(*,20) 1,nn,1,ne,in1,in2,ie1,ie2
 20   format(/' ---  G C O M B  ---',/,
     .' number of points in grid1,  north:',2i7,', east:',2i7,/
     .' relative position of grid2, north:',2i7,', east:',2i7) 
c
c  skip records in grid2 north of grid1
c
      if (in2.gt.nn) then 
        do 21 i = 1, in2-nn
21      read(11,*) (grow2(j),j=1,ne2)
      endif
c
      do 50 i = nn, 1, -1
      if (.not.llab) read(10,*) (grow1(j),j=1,ne)
      if (llab) then
        if (lzero) rr = 0.0
        if (l9999) rr = 9999.0
        do 202 j = 1,ne
202     grow1(j) = rr
      endif
      if (in1.le.i.and.i.le.in2.and.mode.gt.0)
     .read(11,*) (grow2(j),j=1,ne2)
      do 41 j = 1, ne
        r1 = grow1(j)
        if (i.lt.in1.or.i.gt.in2.or.j.lt.ie1.or.j.gt.ie2) then
c      OA - Entered that r2 is 0 when not present - hereby r1 is kept.
          if (mode.eq.5) r2 = r1
          if (mode.lt.5) r2 = 0.0
        elseif (mode.gt.0) then
          r2 = grow2(j-ie1+1)
        endif
c
        goto (30,31,32,33,34,35,36,37,38,29),mode+1
c
c  one-grid modification
c
 29     if (r1.gt.9998.9999.or.abs(r2).lt.0.0001) goto 39
	write(*,*) r1,r2
        rr = r1 / r2
        goto 40

c  one-grid modification
c
 30     if (r1.gt.9998.9999) goto 39
        rr = r1*fact + bias
        goto 40
c
c  difference
c
 31     if (r1.gt.9998.9999.or.r2.ge.9998.9999) goto 39
        rr = r1 - r2
        goto 40
c
c  sum
c
 32     if (r1.gt.9998.9999.or.r2.gt.9998.9999) goto 39 
        rr = r1 + r2
        goto 40
c
c  convert bouguer to free-air
c
 33     if (r1.gt.9998.9999.or.r2.gt.9998.9999) goto 39 
        rr = r1 + 0.1119*r2
        goto 40
c
c  treshold
c
 34     rr = r1
        if (r2.gt.trsh) rr = 9999.0
        goto 40
c
c  overwrite
c
 35     rr = r1
        if (r2.le.9998.9999) rr = r2
        goto 40
c
c  select where no data
c
 36     rr = r1
        if (r2.le.9998.9999) goto 39
        goto 40
c
c  treshold select
c
 37     rr = 0
        if (r1.gt.trsh) rr = r1
        if (r1.le.trsh.and.r2.lt.trsh) rr = r2
        goto 40
c
c  variable multiplication
c
 38     if (r1.gt.9998.9999) goto 39
        if (r2.lt.trsh) rr = r1*fak1
        if (r2.ge.trsh) rr = r1*fak2
        goto 40
c
 39     rr = 9999.0
c
 40     grow1(j) = rr
        if (r1.lt.r1min) r1min = r1
        if (r1.gt.r1max.and.r1.lt.9998.9999) r1max = r1
        if (r1.lt.9998.9999) nr1 = nr1 + 1
        if (r1.lt.9998.9999) r1sum = r1sum + r1
        if (r1.lt.9998.9999) r1sum2 = r1sum2 + r1**2
        if (r2.lt.r2min) r2min = r2
        if (r2.gt.r2max.and.r2.lt.9998.9999) r2max = r2
        if (r2.lt.9998.9999) nr2 = nr2 + 1
        if (r2.lt.9998.9999) r2sum = r2sum + r2
        if (r2.lt.9998.9999) r2sum2 = r2sum2 + r2**2
        if (rr.lt.rrmin) rrmin = rr
        if (rr.gt.rrmax.and.rr.lt.9998.9999) rrmax = rr
        if (rr.lt.9998.9999) nrr = nrr + 1
        if (rr.lt.9998.9999) rrsum = rrsum + rr
        if (rr.lt.9998.9999) rrsum2 = rrsum2 + rr**2
 41   continue
c
      if (lint) write(20,51) (nint(grow1(j)),j=1,ne)
 51   format(30(/,12i6))
      if (.not.lint) write(20,52) (grow1(j),j=1,ne)
 52   format(30(/,8f10.2))
 50   continue
c
      n = nn*ne
      if (nr1.gt.1) r1sum2 = sqrt((r1sum2-r1sum**2/nr1)/(nr1-1))
      if (nr1.le.1) r1sum2 = 0.0
      if (nr1.gt.0) r1sum = r1sum/nr1
      if (nr2.gt.1) r2sum2 = sqrt((r2sum2-r2sum**2/nr2)/(nr2-1))
      if (nr2.le.1) r2sum2 = 0.0
      if (nr2.gt.0) r2sum = r2sum/nr2
      if (nrr.gt.1) rrsum2 = sqrt((rrsum2-rrsum**2/nrr)/(nrr-1))
      if (nrr.le.1) rrsum2 = 0.0
      if (nrr.gt.0) rrsum = rrsum/nrr
      write(*,80) r1sum,r1sum2,r1min,r1max,n-nr1,
     .r2sum,r2sum2,r2min,r2max,n-nr2,
     .rrsum,rrsum2,rrmin,rrmax,n-nrr
80    format(' statistics of grids (in grid1 area, '
     .,'grid 2 poss. extended with 9999-values)',/,
     .'            mean   std.dev.    min     max  unknown',/,
     .' file1:',4f9.2,i7/' file2:',4f9.2,i7/' ofile:',4f9.2,i7)
      close(20)
      end