program geogrid
      implicit double precision(a-h,o-z)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                        g e o g r i d
c
c  program for gridding of irregular distributed data into
c  a regular rectangular grid, for interpolationg data in profiles,
c  or for interpolating individual points using enhanced, fast
c  weighted means interpolation or collocation/kriging.
c  the program may detrend data prior to gridding, or just fit
c  a trend surface to data
c
c  a quadrant search method is used internally in the program
c  to speed up computations. this means that at each prediction
c  point only the 'nqmax' nearest points is used in each of the
c  four quadrants around the point. to speed up the quadrant
c  search an internal data organization is used, where an
c  internal data grid ensuring 'rdat' (p.t. 3) points per compartment
c  in average is used. this grid is not related to the possible
c  wanted prediction grid.
c
c  program input:
c
c  <ifile>
c  <ofile>
c  <efile>
c  nd, idno, nqmax, itrend, ipred,
c  <predpar>,
c  mode, rkm,
c  <pred points>
c
c  where
c
c  nd     ...  data values in line pr. point in <ifile>
c              not used when <ifile> contain a grid (mode 9, see below)
c
c  idno   ...  data value used out of the 'nd' possibilities.
c              if 'idno' is zero the height data are gridded.
c              if 'idno' is negative then the data number abs(idno)
c              is assumed to be followed by a standard deviation.
c              the standard deviation is used in the collocation
c              prediction only if it is greater than the specified
c              sigma, see predpar below.
c              Standard deviations on heights is specified as idno = -99     
c              Special bouguer/free-air option: if idno = 99,
c              free-air is used at sea, bouguer on land, with free-air
c              anomaly position first on line. sea points have negative h.
c
c  nqmax  ...  number of closest points per quadrant used in prediction,
c              if nqmax<=0 all points are used.
c
c  itrend ...  trend surface removal: 0  none (just grid data) 
c                                     1  remove mean
c                                     2  remove linear function in x and y
c                                     3  remove 2nd order polynomium
c                                     4  remove 3rd order polynomial
c                                     5  remove 4-par (geographic only), 
c                                        corresponds to 7-parameter datum shift
c              if 'itrend' is negative only the residuals from the trending
c              are output (i.e. data trend is not restored after prediction)
c
c  ipred  ...  prediction type:  0  none (detrend only)
c                                1  collocation (kriging)
c                                2  weighted means
c
c  <predpar>   depends on prediction type:
c              ipred = 1: xhalf(km), rms(noise)
c                -     2: prediction power (e.g. 2)
c              note: for collocation rms(noise) is the smallest standard
c              deviation used when data is given with sigmas. 
c
c  mode   ...  specifies data to be predicted:
c              mode = 1: grid in geographical coordinates
c                     2, 3: grid in utm projection
c                     4: profile interpolation
c                     5: point interpolation, data points given in
c                        <efile> in standard gi format, i.e. as
c                           id, lat, lon, height, ...
c                     6: do, with input data assumed to be bouguer
c                        anomalies (GRS80/2.67). output data will be
c                        gravity values, derived from interpolated
c                        anomalies and prediction point heights.
c                     7: do, compare two files by making difference of 
c                        values between values in <efile> minus
c                        predicted values from <ifile>. Prediction is
c                        only done when distance < dmax (km)
c                        <efile> data position no 'idno' is used
c                     8: like 7, but both <ifile> and <efile> of format
c                           t, lat, lon, height, ...
c                        where t is a real number. values only interpolated
c                        from points with abs difference t > tmin
c                     9: fill-in of unknown (9999) grid values. In this case
c                        <ifile> must contain a grid rather than a point list.
c                        Only the 9999-values are actually predicted.
c              negative: Same as positive value, except only positive
c                        prediction values allowed (all negative values  
c                        set to zero). Use for DTM interpolation only.
c
c  rkm    ...  margin (km) for data selection area surrounding the
c              wanted grid (or rectangular envelope of the profile
c              or wanted points for mode > 3).
c
c  <pred points>  specifies prediction points depending on 'mode':
c
c  mode 1: fi1, fi2, la1, la2, dfi, dla        (geographic grid)
c  -    2: n1, n2, e1, e2, dn, de / iell, zone (utm grid)
c  -    3: fi1, la1, dn, de, nn, ne / iell, zone (utm grid with
c                                               sw-corner lat/lon,
c                                               spacing, and number of
c                                               points)
c  -    4: fi1, la1, fi2, la2, n               (profile with 'n' points)
c  -    5, 6: lthr
c  _    7: dmax, lthr
c  _    8: dmax, tmin
c  _    9: (nothing)
c
c  the utm projections are defined by their zone no ('uzone') and
c  ellipsoid no ('ie'): ie = 1: wgs84/nad83, ie = 2: hayf/ed50,
c  ie = 3: clarke/nad27.
c
c  the <ifile> contains data points in standard line format, i.e. given
c  as lines 'id, lat, lon, height, data(1), .. , data(nd)', where
c  'idno' specifies which data to be used.
c  data values <= -9999 or >= 9999 signals missing data and not used.
c
c  output on grid format is written in 'ofile'. for collocation errors
c  are written on 'ofile'. for weighted mean interpolation 'efile'
c  contains distances in km to the closest data point.
c  for profile and point interpolation all output is on 'ofile'.
c
c  for collocation a second order markov model is used, with c0
c  assigned from the data variance.
c
c  current program limits:
c  maximal number of points: 15790000
c  max e-w prediction grid dimension: 5000
c  max selected pr quadrant: 25
c
c  programmer: rene forsberg, danish geodetic institute, nov 85
c  updated  oct 29, 1991, rf/se
c  modified for vax, rf unsw, dec 88
c  inclusion of individual variances and detrending, rf jan 1990
c  grid fill-in included aug 1990, rf
c  negative modes dec 16, 1991, rf
c  sorting arrays expanded and lambert/polar projections jan2002,rf
c  xgeogrid comparison features mar 2004, rf
c  thr options oban july 2004
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      dimension dl(36), sa(22)
      dimension pred(1000), sdev(1000), isel(100), row(5000)
      dimension c(5150), cp(100), ov(11), tcoef(11)
      real*4 z,zsig2,t
      common /dat/x(15790000),y(15790000),z(15790000),zsig2(15790000),
     .t(15790000)
      common /dsort/ij(15790000),ijn(15790000),ifc(8100),ih(8100)
      common /spar/xx,yy,dyy,dxx,yy1,xx1,r2min,dmax,
     .thr,tmin,nqmax,nyy,nxx,ld,lt
      common /tpar/cosfi,radeg,xx0,yy0
      common /g2lpar/
     .re,rfi0,rla0,cosfi0,sinfi0,resin,sinfi2,sinfi12,itype
      logical lutm, lsig, lfaba, lsel, lrst, ligrid, lneg, lt, ld, lthr
      character*72 ifile, ofile, efile
c
c  constants
c
      npmax = 15790000
      namax = 8100
      irwdim = 5000
      radeg = 180/3.1415926535d0
      degkm = 1.852*60
      rdat = 3.0
c
      write(*,1)
1     format(/,
     .' *********************************************************',
     .'************'/
     .' *    GEOGRID - GRAVSOFT data gridding - vers. JUL 04 - (c',
     .') RF/KMS   *'/
     .' *********************************************************',
     .'************')
      read(*,2) ifile
      read(*,2) ofile
      read(*,2) efile
2     format(a72)
c
      open(10,file=ifile,status='old')
      open(20,file=ofile,status='unknown')
      open(21,file=efile,status='unknown')
c
c  input remaining control data
c  ----------------------------
c
      read(*,*) nd,idno,nqmax,itrend,ipred
c
      lrst = (itrend.gt.0)
      itrend = abs(itrend)
      lsig = (idno.lt.0)
      if (idno.eq.-99) idno = 0
      if (lsig) idno = abs(idno)
      lfaba = (idno.eq.99)
      if (lfaba) idno = 1
c
      if (ipred.eq.1) read(*,*) xhalf,rmsn
      if (ipred.eq.2) read(*,*) ipwr
      read(*,*) mode,rkm
      lneg = (mode.lt.0)
      if (lneg) mode = iabs(mode)
      if (mode.lt.1.or.mode.gt.9) then
        write(*, 11)
11      format(' *** mode parameter undefined ')
        stop
      endif
      if (itrend.gt.5) stop 'itrend parameter too big'
      if (itrend.eq.5.and.(mode.eq.2.or.mode.eq.3.or.mode.eq.9)) 
     .stop 'itrend=3 only allowed for usual geographic grid'
      lutm = .false.
      if (mode.ge.2.and.mode.le.3) lutm = .true.
      ld = (mode.eq.7.or.mode.eq.8)
      lt = (mode.eq.8)
      ligrid = (mode.eq.9)
      lthr = .false.
c
      goto (12,12,13,14,15,15,15,15,18),mode
c
c  grid limits
c  mode 1 and 2: grid limit specs
c
12    read(*,*) rfi1,rfi2,rla1,rla2,dfi,dla
      if (lutm) read(*,*) iell,izone
      nn = (rfi2-rfi1)/dfi+1.499999
      ne = (rla2-rla1)/dla+1.499999
      if (lutm) call utmcon(iell, izone, sa)
      goto 20
c
c  mode 3: utm grid spec by sw corner lat/lon
c
13    read(*,*) rfi,rla,dfi,dla,nn,ne
      read(*,*) iell,izone
      call utmcon(iell, izone, sa)
      call utg(rfi/radeg, rla/radeg, rn, re, sa, .false., .true.)
      rfi1 = rn
      rla1 = re
      rfi2 = rfi1 + (nn-1)*dfi
      rla2 = rla1 + (ne-1)*dla
      goto 20
c
c  mode 4: profile data
c
14    read(*,*) rfi1,rla1,rfi2,rla2,nn
      dfi = .00001
      dla = .00001
      if (nn.gt.1) dfi = (rfi2-rfi1)/(nn-1)
      if (nn.gt.1) dla = (rla2-rla1)/(nn-1)
      goto 20
c
c  mode 5,6,7,8: point data. find min and max lat and lon.
c
15    rfi1 = 999999.9
      rla1 = 999999.9
      rfi2 = -999999.9
      rla2 = -999999.9
      if (mode.eq.5.or.mode.eq.6) read(*,*) lthr
      if (ld) then
        if (mode.eq.7) read(*,*) dmax,lthr
        if (mode.eq.8) read(*,*) dmax,tmin
        dmax = dmax*1000
        dmax2 = dmax**2
      endif
c
      nn = 0
16    if (mode.le.6) then
        read(21,*,end=17) thr,rlat,rlon,rh
        id = thr
      elseif (mode.eq.7) then
        read(21,*,end=17) thr,rlat,rlon,rh,(dl(j),j=1,nd)
        id = thr
      elseif (lt) then
        read(21,*,end=17) thr,rlat,rlon,rh,(dl(j),j=1,nd)
      endif
c
      nn = nn + 1
      if (rlat.lt.rfi1) rfi1 = rlat
      if (rlat.gt.rfi2) rfi2 = rlat
      if (rlon.lt.rla1) rla1 = rlon
      if (rlon.gt.rla2) rla2 = rlon
      goto 16
17    rewind(21)
      goto 20
c
c  mode 9: grid input
c
18    read(10,*) rfi1,rfi2,rla1,rla2,dfi,dla
      lutm = (abs(rfi1).gt.100.or.abs(rfi2).gt.100)
      if (lutm) then 
        read(10,*) iell,izone
        call utmcon(iell, izone, sa)
        write(*,1802) rfi1,rfi2,rla1,rla2,dfi,dla,iell,izone
1802    format(/' input grid: ',6f9.0,2i4)
      else
        write(*,1801) rfi1,rfi2,rla1,rla2,dfi,dla
1801    format(/' input grid: ',6f9.4)
      endif
      nn = (rfi2-rfi1)/dfi + 1.5
      ne = (rla2-rla1)/dla + 1.5
      if (ne.gt.irwdim) stop 'too long rows in grid'
c
c  parameter checks
c
20    if (lutm) then
        if (iell.lt.1.or.iell.gt.4.or.
     .  izone.lt.1.or.izone.gt.99)
     .  stop 'utm grid label not ok'
      endif
      lsel = .true.
      if (nqmax.le.0) lsel = .false.
      if (nqmax.gt.25) stop '*** nqmax too large'
      if (mode.le.3.and.(nn.lt.1.or.ne.lt.1.or.ne.gt.irwdim)) then
        write(*, 21)
21      format(' *** grid specification illegal')
        stop
      endif
c
c  define internal map projection for geographic 
c  make special case for polar regions
c
      if (lutm) then
        cosfi = 1.0
      else 
        if (rfi2.eq.90.and.abs(rla2-rla1).gt.180) then
          rlat0 = 90
          cosfi = cos(rfi1/radeg)
        elseif (rfi1.eq.-90.and.abs(rla2-rla1).gt.180) then
          rlat0 = -90
          cosfi = cos(rfi2/radeg)
        else
          rlat0 = (rfi1+rfi2)/2
          cosfi = cos(rlat0/radeg)
        endif
        rlon0 = (rla1+rla2)/2
        call setlamb(rlat0,rlon0)
        write(*,22) rlat0,rlon0
22      format(' internal lambert map projection center: ',2f9.3)
      endif
      if (mode.eq.4) dpr = sqrt(dfi**2+(dla*cosfi)**2)*degkm
c
c   input list of data from ifile
c   -----------------------------
c
      ndp = 0
      np = 0
      nneg = 0
      if (lutm) then
        r = rkm*1000 + 0.1 
        rfimin = min(rfi1,rfi2) - r
        rfimax = max(rfi1,rfi2) + r
        rlamin = min(rla1,rla2) - r
        rlamax = max(rla1,rla2) + r
      else
        r = rkm/degkm + 0.000001
        rfimin = min(rfi1,rfi2) - r
        rfimax = max(rfi1,rfi2) + r
        rlamin = min(rla1,rla2) - r/cosfi
        rlamax = max(rla1,rla2) + r/cosfi
      endif
      zmin = 9999999.9
      zmax = -9999999.9
      zsmin = zmin
      zsmax = zmax
      zpsmin = zmin
      zpsmax = zmax
      zsum = 0
      zsum2 = 0
      xx1 = 9999999.9
      xx2 = -9999999.9
      yy1 = 9999999.9
      yy2 = -9999999.9
      ii = nn+1
      jj = ne+1
c
c  loop entry point for data input
c  -------------------------------
c
25    if (ligrid) then
24      if (jj.eq.ne+1) then 
          if (ii.eq.1) goto 30
          read(10,*) (row(kk),kk=1,ne)
          jj = 1
          ii = ii-1
        endif 
        if (lutm) then
          yy = rfi1 + (ii-1)*dfi
          xx = (rla1 + (jj-1)*dla)
        else
          call geo2lamb(.true.,rfi1+(ii-1)*dfi,rla1+(jj-1)*dla,yy,xx)
        endif
        zz = row(jj)
        jj = jj+1
        ndp = ndp + 1
        if (zz.ge.9999) goto 24
      else
c
c  input point reading
c
        read(10,*,end=30) thr,rlat,rlon,rh,(dl(kk),kk=1,nd)
        id = thr
c
        ndp = ndp + 1
        if (.not.lutm) then
          if (rlat.lt.rfimin.or.rlat.gt.rfimax.or.
     .    rlon.lt.rlamin.or.rlon.gt.rlamax) goto 25
          call geo2lamb(.true.,rlat,rlon,yy,xx)
        else
          call utg(rlat/radeg,rlon/radeg,yy,xx,sa,.false.,.true.)
          if (xx.lt.rlamin.or.xx.gt.rlamax.or.
     .    yy.lt.rfimin.or.yy.gt.rfimax) goto 25
        endif
        zz = dl(idno)
        if (zz.le.-9999.or.zz.ge.9999) goto 25
        if (lfaba.and.rh.gt.0.and.rh.lt.9999) zz = dl(idno+1)
        if (idno.eq.0) zz = rh
        zzsig = dl(idno+1)
      endif
c
      np = np+1
      if (xx.lt.xx1) xx1 = xx
      if (xx.gt.xx2) xx2 = xx
      if (yy.lt.yy1) yy1 = yy
      if (yy.gt.yy2) yy2 = yy
      if (np.gt.npmax) then
        write(*, 28) npmax
28      format(' *** too many data points - max: ',i7)
        stop
      endif
      x(np) = xx
      y(np) = yy
      z(np) = zz
      if (lsig) zsig2(np) = zzsig**2
      if (lt) t(np) = thr
      if (zz.lt.zmin) zmin = zz
      if (zz.gt.zmax) zmax = zz
      if (lsig.and.zzsig.lt.zsmin) zsmin = zzsig
      if (lsig.and.zzsig.gt.zsmax) zsmax = zzsig
      zsum = zz+zsum
      zsum2 = zz**2 +zsum2
      goto 25
c
c  detrend data if required
c  ------------------------
c
30    if (itrend.eq.0) goto 309
      xx0 = xx1
      yy0 = yy1
      do 301 i = 1, 66
301   c(i) = 0.0
      do 302 k = 1, np
        zz = z(k)
        call setov(x(k),y(k),zz,itrend,nt,ov)
        ii = 0
        do 302 i = 1, nt+1
        do 302 j = 1, i
          ii = ii+1
          c(ii) = c(ii) + ov(i)*ov(j)
302   continue          
      call chol(c,nt,nsing)
      if (nsing.ne.0) write(*,303) nsing
303   format(' *** warning: detrending ',i2,' singularities')
      k = nt*(nt+1)/2
      do 304 i = 1, nt
304   tcoef(i) = c(k+i)
c
      tsum = 0
      tsum2 = 0
      trmin = 999999.9
      trmax = -trmin
      do 305 k = 1, np 
        rr = z(k) - trend(x(k),y(k),itrend,tcoef)
        tsum = rr + tsum
        tsum2 = rr**2 + tsum2
        if (rr.lt.trmin) trmin = rr
        if (rr.gt.trmax) trmax = rr
        z(k) = rr
305   continue
c
c  output info on input, data and selected internal grid
c  -----------------------------------------------------
c
309   if (np.lt.1) then
        write(*, 32)
32      format(' *** too few points in selection area')
        write(*,'(4f14.4)') rfimin,rfimax,rlamin,rlamax
        stop
      endif
      if (np.gt.1) then
        rs = sqrt((zsum2 - zsum**2/np)/(np-1))
      else 
        rs = 0
      endif
      rm = zsum/np
      if (itrend.gt.0) then
        if (np.gt.1) then
          ts = sqrt((tsum2 - tsum**2/np)/(np-1))
        else
          ts = 0
        endif
        tm = tsum/np
      endif
c
c  covariance parameters - second order markov model
c
      if (ipred.ne.1) goto 79
      sqrc0 = rs
      if (itrend.gt.0) sqrc0 = ts
      c0 = sqrc0**2
      cnn = rmsn**2
      alfa = 0.595*xhalf*1000
c
79    continue
c
      if (.not.lutm) then
        write(*,85) rkm,rfimin,rfimax,rlamin,rlamax
85      format(' rkm = ',f6.1,', data selection area: ',4f9.4)
      else
        write(*,86) rkm,
     .  rfimin,rfimax,rlamin,rlamax
86      format(' rkm = ',f6.1,', data selection area: ',4f10.0)
        call utg(rfi2,rla1,yy,xx,sa,.true.,.true.)
        call utg(rfi2,rla2,yi,xi,sa,.true.,.true.)
        call utg(rfi1,rla1,yj,xj,sa,.true.,.true.)
        call utg(rfi1,rla2,yk,xk,sa,.true.,.true.)
        write(*,88) yy*radeg,xx*radeg,yi*radeg,xi*radeg,
     *  yj*radeg,xj*radeg,yk*radeg,xk*radeg
88      format(' geographical coordinates of grid corners:',
     .  2(/,'   ',2f9.3,'   ',2f9.3))
      endif
c
      if (ligrid) then
        ndp = nn*ne
        nd = 1
        idno = 1
      endif
      write(*,87) nd,idno,ndp,np,zmin,zmax,rm,rs
87    format(/' data values per point: ',i2,', used no.: ',i2,/,
     *' total points in file:',i7,', selected:',i7,/,
     *' min max mean stddev: ',4f10.3)
      if (lsig) write(*,871) zsmin, zsmax
871   format(' minimal and maximal standard deviations of data: ',2f9.2)
c
      if (itrend.gt.0) write(*,801) itrend,nt,(tcoef(j),j=1,nt)
      if (itrend.gt.0) write(*,802) trmin,trmax,tm,ts
801   format(/' detrending done on data, itrend = ',i1,/
     .' no of trend parameters estimated: ',i2,/
     .' solution: ',5e12.4,/,' ',5e12.4)
802   format(' detrended data (min,max,mean,stddev): ',4f9.3)
c
      if (ipred.eq.1) write(*,81) sqrc0, xhalf, rmsn
81    format(/' collocation prediction - sqrc0,xhalf(km),rmsn = ',3f7.2)
      if (ipred.eq.2) write(*,82) ipwr
82    format(/' weighted means prediction, ipwr =',i2)
      if (.not.lsel) write(*,83)
83    format(' all points used in each prediction')
      if (lsel) write(*,84) nqmax
84    format(' selection:',i3,' closest points per quadrant')
c
      if (lt) write(*,841) tmin
841   format(' - predictions only using data with time dif ',
     .'greater than: ',f8.3)
c
c  skip data organization if all points used
c
      if (.not.lsel) goto 100
c
33    r = sqrt((yy2-yy1)*(xx2-xx1)/np*rdat)
      if (r.eq.0) r = 1.0
      nyy = (yy2-yy1)/r+.4999
      nxx = (xx2-xx1)/r+.4999
      if (nyy.lt.1) nyy = 1
      if (nxx.lt.1) nxx = 1
      if (nxx*nyy.le.namax) goto 35
        write(*, 34)
34      format(' - large organization, rdat increased -')
        rdat = 2.0*rdat
        goto 33
35    continue
      dyy = (yy2-yy1)/nyy
      dxx = (xx2-xx1)/nxx
      n = nxx*nyy
c
      if (.not.lutm) write(*,37) yy1,yy2,xx1,xx2
37    format(/' data organization limits in lambert proj:',4f9.0)
      if (lutm) write(*,38) yy1,yy2,xx1,xx2
38    format(/' data organization limits: ',4f10.0)
      write(*,39) nyy,nxx,n,dyy/1000,dxx/1000,np*1.d0/n
39    format(' subrectangles (n,e,total): ',3i5,
     */' size (km): ',2f9.2,', average pts per rect (rdat):',f9.3)
c
c  organization of irregular data
c  ------------------------------
c
c  the following piece of code corresponds closely to subroutine
c  'oaf' of 'gspp' by hans sunkel, see osu internal report: gspp
c  manual, ohio state university, 1980
c
c  ij(1...np)   stores element index
c  ijn(1...np)  stores data index in increasing sort element order
c  ifc(1...n1)  stores number of data per sort element
c  ih(1...n1)   pointer vector, data in sort elem k has index ijn(ih(k))
c
      do 41 i = 1,n
41    ifc(i) = 0
      do 45 i = 1,np
        if (dyy.eq.0) then
          iy = 1
        else
          iy = (y(i)-yy1)/dyy
        endif
        if (dxx.eq.0) then
          ix = 1
        else
          ix = (x(i)-xx1)/dxx
        endif
        if (iy.ge.nyy) iy = nyy-1
        if (ix.ge.nxx) ix = nxx-1
        ixy = iy*nxx + ix+1
        ij(i) = ixy
        ifc(ixy) = ifc(ixy) + 1
45    continue
c
      ih(1) = 1
      do 46 i = 2,n
        i1 = i-1
        ih(i) = ifc(i1)+ih(i1)
46    continue
c
      do 47 i = 1,np
        k = ij(i)
        kk = ih(k)
        ijn(kk) = i
        ih(k) = ih(k) + 1
47    continue
c
      j = 0
      jj = 0
      do 48 i = 1,n
        k = ifc(i)
        ih(i) = ih(i)-k
        if (k.eq.0) j=j+1
        if (k.gt.jj) jj=k
48    continue
c
      write(*, 49) jj, 100.0*j/n
49    format(' max points in subrects:',i5,
     *', percentage with no points: ',f5.1)
c
c  prediction loop
c  ---------------
c
c  the point numbers to be used in prediction is stored in
c  array 'isel', selected by subroutine 'select'.
c
100   zmin = 9999999.9
      zmax = -9999999.9
      zsum = 0
      zsum2 = 0
      zvalmin = 9999999.9
      zvalmax = -9999999.9
      zvalsum = 0
      zvalsum2 = 0
      zdifmin = 9999999.9
      zdifmax = -9999999.9
      zdifsum = 0
      zdifsum2 = 0
      no = 0
c
      if (mode.ge.4.and.mode.le.8) ne = 1
      if (lsel) goto 102
        nsel = np
        if (nsel.gt.100) stop '** max 100 pts in collocation'
        do 101 i = 1, nsel
101     isel(i) = i
c
102   if (ligrid) then
        rewind(10)
        read(10,*) r,r,r,r,r,r
        if (lutm) read(10,*) iell,izone
      endif
c
c  main loop entry
c
      do 240 ii = nn, 1, -1
      if (ligrid) read(10,*) (row(kk),kk=1,ne)
      do 200 jj = 1, ne
        goto (105,105,105,106,107,107,107,107,105),mode
c
c  grid points
c
105     if (lutm) then
          xx = (jj-1)*dla+rla1
          yy = (ii-1)*dfi+rfi1
        else
          call geo2lamb(.true.,(ii-1)*dfi+rfi1,(jj-1)*dla+rla1,yy,xx)
        endif
        if (ligrid) then
          if (row(jj).lt.9999) then
            zpred = row(jj)
            zstd = 0
            goto 190
          endif  
        endif
        goto 108
c
c  profile points
c
106     if (lutm) then
          xx = ((nn-ii)*dla+rla1)
          yy = (nn-ii)*dfi+rfi1
        else
          call geo2lamb(.true.,(nn-ii)*dfi+rfi1,(nn-ii)*dla+rla1,yy,xx)
        endif
        goto 108
c
c  individual points
c
107     if (mode.le.6) then
          read(21,*) thr,rlat,rlon,rh
          id = thr
        elseif (mode.ge.7) then
          read(21,*) thr,rlat,rlon,rh,(dl(j),j=1,nd)
          id = thr
          if (nd.eq.0) then
            zval = rh
          else
            zval = dl(idno)
          endif 
        endif
c
        if (lutm) then
          yy = rlat
          xx = rlon
        else
          call geo2lamb(.true.,rlat,rlon,yy,xx)
        endif
c
c  skip predictions for ipred = 0, otherwise select closest data
c
108     if (ipred.eq.0) then
          zpred = 0.0
          zstd = 0.0
          goto 190
        endif
        if (lsel) call select(isel,nsel)
        if (ld) then
          if (r2min.gt.dmax2) goto 240
        endif
c
        if (ipred.eq.1) goto 150
c
c  prediction by weighted means
c  ----------------------------
c
        wsum = 0
        wfsum = 0
        do 120 i = 1, nsel
          k = isel(i)
          dx = x(k)-xx
          dy = y(k)-yy
          r2 = dx**2 + dy**2
c
c  non-zero value if dist lt 10 cm to avoid singularity
c
          if (r2.lt.0.01) r2 = 0.01     
c
c  distance weights
c
          if (ipwr.eq.2) then
            w = 1/r2
          else          
            w = 1/sqrt(r2)**ipwr
          endif
          wsum = w + wsum
          wfsum = w*z(k) + wfsum
120     continue
        zpred = wfsum/wsum
        zstd = sqrt(r2min)/1000
        goto 190
c
c  collocation prediction
c  ----------------------
c
150     k = 1
        do 156 i = 1,nsel+1
          if (i.le.nsel) goto 153
          xi = xx
          yi = yy
          i1 = i-1
          goto 154
153       xi = x(isel(i))
          yi = y(isel(i))
          i1 = i
154       do 155 j = 1,i1
            kk = isel(j)
            r = sqrt((x(kk)-xi)**2 + (y(kk)-yi)**2)
c
c  second order markov covariance function
c
            cov = c0*(1 + r/alfa)*exp(-r/alfa)
            if (i.eq.j) cov = cov + max(cnn,dble(zsig2(kk)))
c
            c(k) = cov
            k = k+1
          if (i.gt.nsel) cp(j) = cov
155       continue
156     continue
c
        call chol(c,nsel,nsing)
        if (nsing.gt.0) write(*,157) ii,jj,nsing
157     format(' *** cholesky: ',i3,' singularities point: ',2i5)
c
        zpred = 0
        zstd = 0
        i1 = nsel*(nsel+1)/2
        do 158 i = 1,nsel
          r = c(i+i1)
          zpred = zpred + r*z(isel(i))
          zstd = zstd + r*cp(i)
158     continue
        r = c0 - zstd
        if (r.lt.0) zstd = -sqrt(-r)
        if (r.ge.0) zstd = sqrt(r)
c
c  statistics for prediction
c  -------------------------
c
190     if (lrst) zpred = zpred + trend(xx,yy,itrend,tcoef)
        if (lneg) then
          if (zpred.lt.0) zpred = 0
          nneg = nneg+1
        endif
        no = no+1
        if (zpred.lt.zmin) zmin = zpred
        if (zpred.gt.zmax) zmax = zpred
        zsum = zsum + zpred
        zsum2 = zsum2 + zpred**2
        if (zstd.lt.zpsmin) zpsmin = zstd
        if (zstd.gt.zpsmax) zpsmax = zstd
        pred(jj) = zpred
        sdev(jj) = zstd
200   continue
c
c  output result and possible heading
c  ----------------------------------
c
      mode10 = mode
      if (mode.eq.9.and.lutm) mode10 = 10
      goto (201,211,211,221,231,231,231,231,201,211),mode10
c
c  output of geographical grid data
c
201   if (ii.eq.nn) write(20,204) rfi1,rfi1+(nn-1)*dfi,
     *rla1,rla1+(ne-1)*dla,dfi,dla
      write(20,205) (pred(jj),jj=1,ne)
204   format(' ',4f12.6,2f12.8)
205   format(/,50(' ',7f10.3,/))
c
      if (ii.eq.nn.and.ipred.eq.1) write(21,207)
     *rfi1,rfi1+(nn-1)*dfi,rla1,rla1+(ne-1)*dla,dfi,dla
      if (ii.eq.nn.and.ipred.eq.2) write(21,208)
     *rfi1,rfi1+(nn-1)*dfi,rla1,rla1+(ne-1)*dla,dfi,dla
207   format(' ',4f12.6,2f12.7)
208   format(' ',4f12.6,2f12.7)
      if (ipred.eq.1) write(21,205) (sdev(jj),jj=1,ne)
      if (ipred.eq.2) write(21,206) (sdev(jj),jj=1,ne)
206   format(/,40(' ',8f9.3,/))
      goto 240
c
c  utm grid data
c
211   yi = rfi1
      yj = (rfi1+(nn-1)*dfi)
      xi = rla1
      xj = (rla1+(ne-1)*dla)
      yy = dfi
      xx = dla
      if (ii.lt.nn) goto 217
      write(20,212) yi,yj,xi,xj,yy,xx,iell,izone
212   format(' ',6f11.0,/,' ',2i4)
217   write(20,205) (pred(jj),jj=1,ne)
c
      if (ii.lt.nn) goto 220
      write(21,212) yi,yj,xi,xj,yy,xx,iell,izone
220   if (ipred.eq.1) write(21,205) (sdev(jj),jj=1,ne)
      if (ipred.eq.2) write(21,206) (sdev(jj),jj=1,ne)
      goto 240
c
c  profile data
c
221   if (ii.eq.nn) write(*,224) rfi1,rla1,rfi2,rla2,dpr
224   format(/' profile ',2f8.3,'  to ',2f8.3,', spacing ',f9.2,
     *' km')
      write(20,223) (nn-ii)*dpr, zpred, zstd
223   format(' ', f9.2,2f12.3)
      goto 240
c
c  point data
c  conversion grs80 bouguer anomaly to gravity value
c
231   if (mode.eq.6) then
        r = sin(rlat/radeg)**2
        gamma80 = 978032.677d0*(1+.00193185135d0*r)/
     .  sqrt(1 - .00669438002d0*r)
        if (rh.ge.0) hfakt = 0.1967
        if (rh.lt.0) hfakt = -0.0687
        zpred = zpred + gamma80 - hfakt*rh
      endif
      if (ld) then
c
        zdif = zval - zpred 
c
        if (zval.lt.zvalmin) zvalmin = zval
        if (zval.gt.zvalmax) zvalmax = zval
        zvalsum = zvalsum + zval
        zvalsum2 = zvalsum2 + zval**2
        if (zdif.lt.zdifmin) zdifmin = zdif
        if (zdif.gt.zdifmax) zdifmax = zdif
        zdifsum = zdifsum + zdif
        zdifsum2 = zdifsum2 + zdif**2
c
        if (lt.or.lthr) then
          write(20,234) thr,rlat,rlon,rh,zpred,zstd,zdif
        else
          write(20,233) id,rlat,rlon,rh,zpred,zstd,zdif
        endif
      else
        if (lthr) then        
          write(20,234) thr,rlat,rlon,rh,zpred,zstd
        else
          write(20,233) id,rlat,rlon,rh,zpred,zstd
        endif
      endif
233   format(' ',i9,2f11.5,f9.2,f11.2,2f9.2)
234   format(' ',f12.5,2f11.5,f9.2,f11.3,2f9.3)
c
240   continue
c ------------------------------------------------
c  prediction loop exit
c
      write(*,*)
      if (ld) then
        if (no.le.1) rs = 0.0
        if (no.gt.1) rs = sqrt((zvalsum2 - zvalsum**2/no)/(no-1))
        rm = 0
        if (no.gt.0) rm = zsum/no
        write(*,249) nn,dmax/1000,rm,rs,zvalmin,zvalmax
249     format(' total number of points in datafile (efile): ',i7,/ 
     .  ' prediction pts selected closer than ',f7.3,' km'/
     .  ' selected pts   mean std.dev. min max: ',4f9.2)
      endif

      if (no.le.1) rs = 0.0
      if (no.gt.1) rs = sqrt((zsum2 - zsum**2/no)/(no-1))
      rm = 0
      if (no.gt.0) rm = zsum/no
      write(*,250) no,rm,rs,zmin,zmax,zpsmin,zpsmax
250   format(' predicted: ',i7,' points'/
     *' prediction pts mean std.dev. min max: ',4f9.2/,
     .' prediction error values min max:   ',21x,2f9.2)
c
      if (ld) then
        if (no.le.1) rs = 0.0
        if (no.gt.1) rs = sqrt((zdifsum2 - zdifsum**2/no)/(no-1))
        rm = 0
        if (no.gt.0) rm = zdifsum/no
        write(*,251) rm,rs,zdifmin,zdifmax
251     format(' differences    mean std.dev. min max: ',4f9.2)
      endif
c
      if (.not.lrst.and.itrend.gt.0) write(*, 255)
255   format(' - trend in data not restored in output -'/)
      if (lneg) write(*,256) nneg
256   format(' - number of negative predictions set to zero:',i7)
c
      close(10)	
      close(20)	
      close(21)	
      end
c
      subroutine select(isel,nsel)
      implicit double precision(a-h,o-z)
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                       s e l e c t
c
c  subroutine for selection of 'nqmax' points in each quadrant
c
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension isel(100)
      real*4 z,zsig2,t
      common /dat/x(15790000),y(15790000),z(15790000),zsig2(15790000),
     .t(15790000)
      common /dsort/ij(15790000),ijn(15790000),ifc(8100),ih(8100)
      common /spar/xx,yy,dyy,dxx,yy1,xx1,r2min,dmax,
     .thr,tmin,nqmax,nyy,nxx,ld,lt
      common /qpar/rq2(4,25),iq(4,25),ns(4),r2max(4),iimax(4)
      logical lq1,lq2,lq3,lq4,ld,lt
c
c  initialization of parameters for rect elem sort
c
      lq1 = .false.
      lq2 = .false.
      lq3 = .false.
      lq4 = .false.
      r2min = 9.999d9
      do 11 j = 1,4
        r2max(j) = 0.0
        ns(j) = 0
11    continue
c
      if (dyy.eq.0) then
        i0 = 1
      else
        i0 = (yy-yy1)/dyy+1
      endif
      if (dxx.eq.0) then
        j0 = 1
      else
        j0 = (xx-xx1)/dxx+1
      endif
c
c  smallest scan area is 3 x 3
c
      i1 = i0-1
      i2 = i0+1
      j1 = j0-1
      j2 = j0+1
      if (i1.lt.1) i1=1
      if (i2.gt.nyy) i2=nyy
      if (j1.lt.1) j1=1
      if (j2.gt.nxx) j2=nxx
      do 8 i = i1,i2
      do 8 j = j1,j2
8     call rsort(i,j)
c
c  check section
c  check quadrants one by one and scan
c
10    if (ns(1).ge.nqmax) lq1 = .true.
      if (ns(2).ge.nqmax) lq2 = .true.
      if (ns(3).ge.nqmax) lq3 = .true.
      if (ns(4).ge.nqmax) lq4 = .true.
      i1 = i1-1
      i2 = i2+1
      j1 = j1-1
      j2 = j2+1
      if (i2.gt.nyy.and.j1.lt.1) lq1 = .true.
      if (i2.gt.nyy.and.j2.gt.nxx) lq2 = .true.
      if (i1.lt.1.and.j1.lt.1) lq3 = .true.
      if (i1.lt.1.and.j2.gt.nxx) lq4 = .true.
c
c  check if points in new cells farther away than dmax
c
      if (ld) then
        if ((i0-i1-1)*dyy.gt.dmax.and.(j0-j1-1)*dxx.gt.dmax) goto 90
      endif
      if (lq1.and.lq2.and.lq3.and.lq4) goto 90
c
c  scan edges
c
      if (.not.(lq1.and.lq2)) call rsort(i2,j0)
      if (.not.(lq2.and.lq4)) call rsort(i0,j2)
      if (.not.(lq3.and.lq4)) call rsort(i1,j0)
      if (.not.(lq1.and.lq3)) call rsort(i0,j1)
      if (lq1) goto 30
        do 25 i = i0+1,i2
25      call rsort(i,j1)
        do 26 j = j1+1,j0-1
26      call rsort(i2,j)
30    if (lq2) goto 40
        do 35 i = i0+1,i2
35      call rsort(i,j2)
        do 36 j = j0+1,j2-1
36      call rsort(i2,j)
40    if (lq3) goto 50
        do 45 i = i1,i0-1
45      call rsort(i,j1)
        do 46 j = j1+1,j0-1
46      call rsort(i1,j)
50    if (lq4) goto 60
        do 55 i = i1,i0-1
55      call rsort(i,j2)
        do 56 j = j0+1,j2-1
56      call rsort(i1,j)
c
60    goto 10
c
90    nsel = 0
      k = 0
      do 91 i = 1,4
        k = ns(i)
        do 92 j = 1,k
92      isel(j+nsel) = iq(i,j)
        nsel = nsel + k
91    continue
      return
      end
c
      subroutine rsort(i,j)
      implicit double precision(a-h,o-z)
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                       r s o r t
c
c  subroutine sorts 'nqmax' closest points within each quadrant
c  in a subrectangle
c
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*4 z,zsig2,t
      common /dat/x(15790000),y(15790000),z(15790000),zsig2(15790000),
     .t(15790000)
      common /dsort/ij(15790000),ijn(15790000),ifc(8100),ih(8100)
      common /spar/xx,yy,dyy,dxx,yy1,xx1,r2min,dmax,
     .thr,tmin,nqmax,nyy,nxx,ld,lt
      common /qpar/rq2(4,25),iq(4,25),ns(4),r2max(4),iimax(4)
      logical ld,lt
c
      if (i.lt.1.or.i.gt.nyy) return
      if (j.lt.1.or.j.gt.nxx) return
c
      ixy = (i-1)*nxx+j
      k = ifc(ixy)
      if (k.eq.0) return
c
      k1 = ih(ixy)
      k2 = k1+k-1
      do 50 k = k1,k2
        kk = ijn(k)
        dx = x(kk)-xx
        dy = y(kk)-yy
        if (lt) then
          if (abs(t(kk)-thr).lt.tmin) goto 50
        endif
        if (dy.lt.0) goto 15
        iiq = 1
        if (dx.ge.0) iiq = 2
        goto 20
15      iiq = 3
        if (dx.ge.0) iiq = 4
20      r2 = dx**2+dy**2
        if (r2.lt.r2min) r2min = r2
c
c  add point if not enough in quadrant
c
        ii = ns(iiq)+1
        if (ii.gt.nqmax) goto 30
        ns(iiq) = ii
        iq(iiq,ii) = kk
        rq2(iiq,ii) = r2
        if (r2.lt.r2max(iiq)) goto 50
        r2max(iiq) = r2
        iimax(iiq) = ii
        goto 50
c
c  check if point is among the closest points
c
30      if (r2.gt.r2max(iiq)) goto 50
        ii = iimax(iiq)
        iq(iiq,ii) = kk
        rq2(iiq,ii) = r2
        do 31 jj = 1,nqmax
          rr = rq2(iiq,jj)
          if (rr.le.r2) goto 31
          r2 = rr
          iimax(iiq) = jj
31      continue
        r2max(iiq) = r2
c
50    continue
      end
c
      subroutine setov(xx,yy,obs,itrend,nt,ov)
      implicit double precision(a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                      s e t o t v
c
c  sets observation equation coefficients for detrending
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension ov(11)
      common /tpar/cosfi,radeg,xx0,yy0
      common /g2lpar/
     .re,rfi0,rla0,cosfi0,sinfi0,resin,sinfi2,sinfi12,itype
c
      if (itrend.eq.5) goto 50
c
c  just mean value or polynomials
c
      dx = xx-xx0
      dy = yy-yy0
      nt = 1  
      ov(1) = 1.0
      if (itrend.eq.1) goto 90
      nt = 3
      ov(2) = dx
      ov(3) = dy
      if (itrend.eq.2) goto 90
      nt = 6
      ov(4) = dx**2
      ov(5) = dy**2
      ov(6) = dx*dy
      if (itrend.eq.3) goto 90
      nt = 10
      ov(7) = dx**3
      ov(8) = dy**3
      ov(9) = dx**2*dy
      ov(10) = dx*dy**2
      goto 90
c
c  4-parameter transformation corresponding to 7-par spatial transformation
c
50    nt = 4
      call geo2lamb(.false.,yy,xx,rfi,rla)
      cosf = cos(rfi/radeg)
      ov(1) = cosf*cos(rla/radeg)
      ov(2) = cosf*sin(rla/radeg)
      ov(3) = sin(rfi/radeg)
      ov(4) = 1.0
c
90    ov(nt+1) = obs
      return
      end
c
      function trend(xx,yy,itrend,tsol)
      implicit double precision(a-h,o-z)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                         t r e n d
c
c  finds value of trend function for given point
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension tsol(11),ov(11)
c
      rr = 0.0
      call setov(xx,yy,rr,itrend,nt,ov)
      do 10 i = 1, nt
10    rr = rr + ov(i)*tsol(i)
      trend = rr
      return
      end
c
      subroutine chol(c,n,nsing)
      implicit double precision(a-h,o-z)
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                      c h o l
c
c  subroutine solves positive definite symmetric linear equations
c  using cholesky decomposition. coefficients stored columnwise
c  in c, followed by righthand side. n is number of unknowns.
c  solution is returned as last column in c.
c  'nsing' is number of singularities. it should be zero for a
c  succesful solution.
c
c  rf, nov 85
c
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension c(5150)
      nsing = 0
c
      do 50 nr = 1,n+1
        i=nr*(nr-1)/2
        ir=i
        do 40 nc = 1,nr
          if (nc.gt.n) goto 50
          sum=0
          ic=nc*(nc-1)/2
          i=i+1
          nc1=nc-1
          do 30 np = 1,nc1
30        sum = sum - c(ir+np)*c(ic+np)
          ci = c(i)+sum
          if (nr.eq.nc) goto 31
          cc = c(ic+nc)
          if (cc.eq.0) then
            write(*,*) '*** Cholesky reduction zero-division'
            cc = 1.0e9
          endif
          c(i) = ci/cc
          goto 40
31        if (nr.gt.n) goto 40
          if (ci.gt.0) goto 32
            nsing = nsing+1
            c(i) = 1.0e9
            goto 40
32        c(i) = sqrt(ci)
40      continue
50    continue
c
c  back substitution
c
      do 80 nc = n,1,-1
        ir=i
        ic=nc*(nc+1)/2
        cc = c(ic)
        if (cc.eq.0) then
          write(*,*) '*** Cholesky subsitution zero-division'
          cc = 1.0e9
        endif
        c(i) = c(i)/cc
        do 70 np = nc-1,1,-1
          ir=ir-1
          ic=ic-1
          c(ir) = c(ir) - c(i)*c(ic)
70      continue
        i = i-1
80    continue
      return
      end
c
      subroutine utmcon(isys, izone, sa)
      implicit double precision (a-h,o-z)
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                         u t m c o n
c
c the procedure produces the constants needed in the transfor-
c mations between transversal mercator and geographical coordina-
c tes and the tolerances needed for the check of the transfor-
c mations. the transformation constants are generated from a
c reg_label defining a transversal mercator system. the formulae
c are taken from konig und weise : mathematische grundlagen der
c hoheren geodasie und kartographie, erster band, berlin 1951.
c
c parameters
c __________
c
c isys, izone         (call)            integers
c specifies ellipsoid and utm zone. the following ellipsoids are
c currently implemented:
c
c     1: wgs84,  2: hayford (ed50),  3: clarke (nad27)
c     4: bessel, if utmzone = 99 then bessel and national swedish
c                projection is assumed (nb: this version is only
c                good to approximatively 10 m for sweden, meridian
c                exact longitude unknown)   
c     5, 6: irish and british systems (utmzone 97 and 98)
c
c sa                  (return)          array
c the constants needed in the transformation.
c sa(1) =           normalized meridian quadrant  (geotype),
c sa(2) =           easting at the central meridian (geotype),
c sa(3) =           longitude of the central meridian (geotype),
c sa(4)  -  sa(7) = const for ell. geo -> sph. geo (real)
c sa(8)  - sa(11) = const for sph. geo -> ell. geo (real),
c sa(12) - sa(15) = const for sph. n, e -> ell. n, e (real),
c sa(16) - sa(19) = const for ell. n, e -> sph. n, e (real),
c sa(20) =          toler. for utm input, 5 mm.
c sa(21) =          toler. for geo input, do.
c sa(22) =          north offset 
c                  
c the user may change sa(20) - sa(21) for special checks.
c
c prog: knud poder, danish geodetic institute, 7 nov 1977,
c updated 18 sep 1983;
c rc fortran version alp/rf oct 86, last updated dec 1990
c updated nov 2001/jan 2002 qatar and british systems
c
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension sa(22)
      double precision n,m,northpr
      dimension aa(6),ff(6)
c
c  ellipsoid parameters
c                  WGS84      Hayford      Clarke       Bessel
c                  Mod.Airy      Airy
      data aa / 6378137.0d0, 6378388.0d0, 6378206.4d0, 6377397.155d0,
     .          6377340.189d0, 6377563.396d0/
      data ff / 298.2572236d0, 297.0d0, 294.9786982d0, 299.153d0,
     .          299.3250d0,  299.3250d0/
      radeg = 180/3.1415926536d0
c
      if (isys.lt.1.or.isys.gt.6) stop 'illegal ellipsoid in utmcon'
      if (izone.lt.1.or.izone.gt.99) stop 'illegal UTM zone in utmcon'
c
      eastpr = 500000.d0
      northpr = 0.d0
      dm = 4.0d-4
      if (izone.eq.95) then
c  qbc
        dm = 1 - 0.999996d0
        eastpr =  400000.d0
        northpr = -2766043.105d0+500000d0
        clong = (50+49/60.d0)/radeg
      elseif (izone.eq.96) then
c  qng
        dm = 1 - 0.99999d0
        eastpr = 200000.d0
        northpr = -2705104.270d0+300000d0
        clong = (51+13/60.d0)/radeg
      elseif (izone.eq.97) then
c  irish national grid
        isys = 5
        dm = -0.000035d0
        eastpr = 200000.d0
        northpr = -5929822.893d0+250000d0
        clong = -8.d0/radeg
      elseif (izone.eq.98) then
c  british osgb36
        isys = 6
        dm = 1 - 0.9996012717d0
        eastpr = 400000.d0 
        northpr = -5427063.818d0-100000d0
        clong = -2.d0/radeg
      elseif (izone.eq.99) then
c  swedish
        isys = 4
        dm = 0.d0
        eastpr = 1500000.d0
        clong = 15.8067/radeg
      else
c  UTM
        clong = ((izone - 30)*6 - 3)/radeg
      endif
c
      a = aa(isys)
      f = 1.d0/ff(isys)
c
c  normalized meridian quadrant
c  see k|nig und weise p.50 (96), p.19 (38b), p.5 (2)
c
      n=f/(2.0-f)
      m=n**2*(1.0/4.0+n**2/64.0)
      w= a*(-n - dm+m*(1.0-dm))/(1.0+n)
      sa(1)=a + w
c
c  central easting and longitude
c
      sa(2)=eastpr
      sa(3)=clong
c
c  check-tol for transformation
c  5.0 mm on earth
c
      sa(20) = 0.0050
      sa(21) = sa(20)/a
c
c  coef of trig series
c
c  ell. geo -> sph. geo., kw p186 - 187 (51) - (52)
c
      sa(4) = n*(-2 + n*(2.0/3.0 + n*(4.0/3.0 + n*(-82.0/45.0))))
      sa(5) = n**2*(5.0/3.0 + n*(-16.0/15.0 + n*(-13.0/9.0)))
      sa(6) = n**3*(-26.0/15.0 + n*34.0/21.0)
      sa(7) = n**4*1237.0/630.0
c
c   sph. geo - ell. geo., kw p190 - 191 (61) - (62)
c
      sa(8)=n*(2.0+n*(-2.0/3.0 +n*(-2.0+n*116.0/45.0)))
      sa(9)=n**2*(7.0/3.0+n*(-8.0/5.0+n*(-227.0/45.0)))
      sa(10)=n**3*(56.0/15.0+n*(-136.0)/35.0)
      sa(11)=n**4* (4279.0/630.0)
c
c  sph. n, e -> ell. n, e,  kw p196 (69)
c
      sa(12)=n*(1.0/2.0+n*(-2.0/3.0+
     .                            n*(5.0/16.0+n*41.0/180.0)))
      sa(13)=n**2*(13.0/48.0+
     .                            n*(-3.0/5.0+n*557.0/1440.0))
      sa(14)=n**3*(61.0/240.0+n*(-103.0/140.0))
      sa(15)=n**4*(49561.0/161280.0)
c
c  ell. n, e -> sph. n, e,  kw p194 (65)
c
      sa(16)=n*(-1.0/2.0+n*(2.0/3.0+
     .                           n*(-37.0/96.0+n*1.0/360.0)))
      sa(17)=n**2*(-1.0/48.0+
     .                           n*(-1.0/15.0+n*437.0/1440.0))
      sa(18)=n**3* (-17.0/480.0+n*37.0/840.0)
      sa(19)=n**4*(-4397.0/161280.0)
c
      sa(22)=northpr
      return
      end
c
      subroutine utg(rn, re, b, l, sa, direct, tcheck)
      implicit double precision(a-h,o-z)
      double precision b, l, sa(22)
      logical direct, tcheck
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                       u t g
c
c dual autochecking transformation procedure for transformation
c between geographical coordinates and transversal mercator co-
c ordinates. the procedure transforms utm->geo when direct is
c true and the reverse when direct is false.
c an alarm is produced when the check by the inverse transforma-
c tion exceeds the tolerance of 5.0 mm or an other value set by
c the user in sa(19) for utm->geo or sa(20) for geo->utm.
c
c n, e              (call)             real
c the utm- or geographical coordinates input for trans-
c formation in meters or radians.
c
c b, l              (return)           real
c the geographical or utm-coordinates output from the procedure
c as transformed and checked coordinates in radians or meters
c
c sa                (call)             array
c transformation constants for direct and inverse transf.
c see fields in sa or set_utm_const for a description
c
c direct            (call)             logical
c direct = true => transformation utm -> geogr.
c direct = false => transformation geogr -> utm
c
c tcheck            (call)             logical
c tcheck = true => check by back transformation
c tcheck = false => no check. this possibility doubles the
c                   speed of the subroutine, with the risk of
c                   obtaining bad results at large (> 60 deg)
c                   distances from the central meridian.
c
c programmer: knud poder, dgi, nov 1977
c rcfortran alp/rf oct 86
c
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision np, lg, l0, ndif, ncheck
      integer h, s
      double precision n, e, n0
c     radeg = 180/3.1415926536d0
c
      n = rn
      e = re
c
      qn = sa(1)
      e0 = sa(2)
      n0 = sa(22)
      l0 = sa(3)
c
c  transformation sequence
c
      if (direct) i=1
      if (.not.direct) i=3
      h = 4-i
      s = 2-i
c
c  check-values
c
      ncheck=n
      echeck=e
c
c  transformation cases
c
      do 100 i=i,h,s
      goto (10,20,30),i
c
c  case 1, utm -> geo
c  ------------------
c
10    np = (n - n0)/qn
      ep = (e - e0)/qn
c
c  ellip. n, e -> sph. n, e
c
      np = np + clcsin(sa, 15, 4, np*2, ep*2, dn, de)
      ep = ep + de
c
c  sph. n, e = compl. sph. lat -> sph lat, lng
c
      cosbn = cos(np)
      if (cosbn.eq.0) cosbn = 1.0d-33
      snh = (exp(ep) - exp(-ep))/2
      lg   = atan(snh/cosbn)
      bbg  = atan(sin(np)*cos(lg)/cosbn)
c
c  sph. lat, lng -> ell. lat, lng
c
      bbg = bbg + clsin(sa, 7, 4, bbg*2)
      lg = lg + l0
      goto 100
c
c  case 2, transf results
c  ----------------------
c
20    b     = bbg
      n     = bbg
      l     = lg
      e     = lg
      if (tcheck) goto 100
      return
c
c  case 3, geo -> utm
c  ------------------
c
30    bbg   = n + clsin(sa, 3, 4, n*2)
      lg    = e - l0
c
c  sph. lat, lng -> compl. sph. lat = sph n, e
c
      cosbn = cos(bbg)
      if (cosbn.eq.0) cosbn = 1.0d-33
      np    = atan(sin(bbg)/(cos(lg)*cosbn))
      rr = sin(lg)*cosbn
      if (abs(rr).ge.0.95) goto 40
      ep = log((1+rr)/(1-rr))/2
      goto 41
40    ep = 1.0d38
41    continue
c
c  sph. normalized n, e -> ell. n, e
c
      np = np + clcsin(sa, 11, 4, np*2, ep*2, dn, de)
      ep = ep + de
      bbg = qn*np + n0
      lg = qn*ep + e0
c
  100 continue
c
c  in/rev-dif for check
c
      ndif  = bbg - ncheck
      edif  = lg - echeck
      edcos = edif
      if(.not.direct) edcos = edcos*cos(ncheck)
c
c  error actions
c
      if (direct) tol = sa(20)
      if (.not.direct) tol = sa(21)
      if (abs(ndif).lt.tol.and.abs(edcos).lt.tol) return
c
      n = rn
      e = re
      if (direct) n = b
      if (direct) e = l
      ep = 6371000.0
      if (direct) ep = 1.0
      write(*, 90) n*57.29578, e*57.29578, ndif*ep, edcos*ep
90    format(' *** utg coor ',2f7.1,' checkdiff too large: ',
     .2f8.3, ' m')
      return
      end
c
      double precision function clsin(a, i0, g, arg)
      implicit double precision(a-h,o-z)
      dimension a(22)
      integer g
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                      c l s i n
c
c computes the sum of a series of a(i+i0)*sin(i*arg) by clenshaw
c summation from g down to 1.
c the sum is the value of the function.
c
c prog.: knud poder 1978
c
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      cosarg = 2*cos(arg)
c
      hr1 = 0.0
      hr  = a(g+i0)
c
      do 10 it = g - 1,1,-1
        hr2 = hr1
        hr1 = hr
        hr = -hr2 + cosarg*hr1 + a(it+i0)
   10 continue
c
      clsin = hr*sin(arg)
c
      return
      end
c
      double precision function clcsin(a, i0, g, argr, argi, r, i)
      implicit double precision(a-h,o-z)
      dimension a(22)
      double precision r, i
      integer g
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                        c l c s i n
c
c computes the sum of a series a(i+i0)*sin(i*arg_r + j*i*arg_i)
c (where j is the imaginary unit) by clenshaw summation
c from g down to 1. the coefficients are here real and
c the argument of sin is complex. the real part of the
c sum is the value of the function.
c
c prog.: knud poder 1978
c
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision ii
c
      sinar = sin(argr)
      cosar = cos(argr)
      ex = exp(argi)
      emx = exp(-argi)
      sinhar = (ex - emx)/2
      coshar = (ex + emx)/2
      rr = 2*cosar*coshar
      ii = -2*sinar*sinhar
c
      hr1 = 0.0
      hi1 = 0.0
      hi  = 0.0
      hr  = a(g+i0)
c
      do 10 it = g-1, 1, -1
        hr2 = hr1
        hr1 = hr
        hi2 = hi1
        hi1 = hi
        hr  = -hr2 + rr*hr1 - ii*hi1 + a(it+i0)
        hi  = -hi2 + ii*hr1 + rr*hi1
   10 continue
c
      rr = sinar*coshar
      ii = cosar*sinhar
c
      r = rr*hr - ii*hi
      clcsin = r
      i = rr*hi + ii*hr
c
      return
      end
c
      subroutine setlamb(rlat0,rlon0)
      implicit double precision(a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 
c  set constants needed for geo2lamb
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      common /g2lpar/ 
     .re,rfi0,rla0,cosfi0,sinfi0,resin,sinfi2,sinfi12,itype
c
      radeg = 180/3.1415926536d0
      re = 6371000.d0
c
      if (abs(rlat0).eq.90) then
        itype = 1
      elseif(rlat0.eq.0) then
        itype = 3
      elseif (abs(rlat0).lt.90) then
        itype = 2
      else
        stop 'too large latitude in setlamb'
      endif
      rfi0 = rlat0/radeg
      rla0 = rlon0/radeg
      if (itype.ne.2) return
      cosfi0 = cos(rfi0)
      sinfi0 = sin(rfi0)
      resin = re/sinfi0
      sinfi2 = 2*sinfi0
      sinfi12 = 1 + sinfi0**2
      return
      end
c
      subroutine geo2lamb(ldir,xin,yin,xout,yout)
      implicit double precision(a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  program implements lambert equivalent projections
c  including special cases on equator and pole.
c  Reference: Ricardus and Adler: Map Projections, 1971.     
c
c  (c) Rene Forsberg, KMS, Dec 2001
c
c  ldir: true: geo->lam, false: lam->geo
c 
c  xin, yin:  lat,lon   or N,E  input
c  xout, yout:  N,E     or lat,lon  output
c
c  (c) Rene Forsberg, KMS, Dec 2001
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      logical ldir
      common /g2lpar/ 
     .re,rfi0,rla0,cosfi0,sinfi0,resin,sinfi2,sinfi12,itype
      data pi / 3.1415926536d0 /
c
      radeg = 180/pi
      goto (100,200,300),itype
c
c  polar equidistant - R&A p167-13
c
100   if (ldir) then
        rfi = xin/radeg
        rla = yin/radeg-rla0
        if (rfi0.gt.0) then
          r = 2*re*sin(pi/4-rfi/2)
        else
          r = 2*re*sin(pi/4+rfi/2)
        endif
        xout = r*cos(rla)
        yout = r*sin(rla)
      else
        r = sqrt(xin**2+yin**2)/(2*re)
        if (rfi0.gt.0) then
          rfi = pi/2 - 2*asin(r)
        else
          rfi = -(pi/2 - 2*asin(r))
        endif
        rla = atan2(yin,xin) + rla0
        xout = rfi*radeg 
        yout = rla*radeg
      endif
      return
c
c  lambert conical equivalent (albers) - R&A p166-11
c 
200   if (ldir) then
        rfi = xin/radeg
        rla = yin/radeg - rla0
        r = sqrt(sinfi12-sinfi2*sin(rfi))
        xout = resin*(cosfi0 - cos(rla*sinfi0)*r)
        yout = resin*sin(rla*sinfi0)*r
      else
        rx = cosfi0-xin/resin
        ry = yin/resin
        rla = atan2(ry,rx)/sinfi0 + rla0
        rfi = asin((sinfi12-(rx**2+ry**2))/sinfi2) 
        xout = rfi*radeg 
        yout = rla*radeg
      endif
      return
c
c  lambert equator cylindrical equivalent - R&A p166-12
c
300   if (ldir) then
        rfi = xin/radeg
        rla = yin/radeg - rla0
        xout = re*sin(rfi)
        yout = re*rla
      else
        rfi = asin(xin/re)
        rla = yin/re + rla0
        xout = rfi*radeg 
        yout = rla*radeg
      endif
      return
      end