c--------( TXPROJ : translate grid-data into new projection )--------
      parameter (LXO=4501, LYO=4501, LXN=8001, LYN=8001)
c
      dimension f(LXO,LYO), g(LXN,LYN)
      character*80 wdr, flog, fnam1, fnam2
      character area*8, area1*8, s9*9, salt*8, saltn*8
      character buf*80, out*72, sdtm*20
c
      call premsg('TXPROJ : translate grid-data into new projection')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Input source data filename ==> ',
     *            81-ldr, fnam1(ldr:80))
      open(1,file=fnam1,status='old')
    1 read(1,'(a)') buf
      if(buf(1:1).eq.'#') goto 1
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        ispa = 0
        ispb = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.62) iorg = -5400
      endif
      read(1,'(a)') buf
      read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt
      write(salt,'(1x,f7.0)') alt
      if (alt.eq.0.) salt = '    var.'
      if (alt.lt.0.) salt = '   undef'
      rewind(1)
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      ixt = ixs + (mx-1)*mszx
      iyt = iys + (my-1)*mszy
      smx = float(mszx)
      smy = float(mszy)
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      xt = float(ixt)/1000.
      yt = float(iyt)/1000.
      write(s9,'(f9.2)') vnul
      if(s9(9:9).ne.'0'.or.s9(1:1).ne.' ')
     *     call abendm('unable to handle vnul value')
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      write(out,'(a,2a10,2x,2a8,2a6,a7)')
     *     '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *     'alt'
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6,a8)')
     *     '--- source data ', xs, ys, mszx, mszy, mx, my, salt
      call premsg(out)
      saltn = salt
      if(alt.lt.0.) then
        call premsg('  Altitude [in m] for this data')
        call gparmf('   (Enter -1 to maintain UNDEF) ==> ', alt)
        if (alt.lt.0.5 .and. alt.ne.-1.) call abendm('invalid value')
        write(saltn,'(1x,f7.0)') alt
      endif
      call premsg('')
      call gparmi('Enter new Projection-Number ==> ', nccn)
      ncn = nccn
      if(nccn.ge.200) ncn = nccn - 200
      ispan = 0
      ispbn = 0
      if(ncn.ge.1.and.ncn.le.62) then
        iorgn = 0
        korgn = (ncn-30)*360 - 180
        if(ncn.gt.60) korgn = 0
        if(ncn.eq.61) iorgn = +5400
        if(ncn.eq.62) iorgn = -5400
      else if(ncn.eq.65) then
        iorgn = 0
        call gparmi('  Central Longitude ?  (Deg.) : ', i1)
        korgn = i1*60
      else
        call gparmi2('  Origin Latitude ?   (Deg. Min.) : ', i1,i2)
        iorgn = i1*60 + i2
        call gparmi2('  Origin Longitude ?  (Deg. Min.) : ', i1,i2)
        korgn = i1*60 + i2
        if(ncn.eq.72) then
          call gparmi2('  St.Parallel-1 Lat. ?  (Deg. Min.) : ', i1,i2)
          ispan = i1*60 + i2
          call gparmi2('  St.Parallel-2 Lat. ?  (Deg. Min.) : ', i1,i2)
          ispbn = i1*60 + i2
        endif
      endif
      call gparmi2('  Mesh sizes in NS and EW (m) ?  [mszx mszy] : ',
     *             mszx2, mszy2)
      if (mszx2.le.0 .or. mszy2.le.0) call abendm('illegal size')
      call gparma('Enter new projection grid output filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      call clspin()
      call cvinit(ncc, float(iorg),float(korg),float(ispa),float(ispb))
      call cvenik(ys, xs, fi0, fk0)
      call cvenik(yt, xs, fi1, fk1)
      call cvenik(yt, xt, fi2, fk2)
      call cvenik(ys, xt, fi3, fk3)
      if (ncc.lt.200.and.nccn.ge.200) then
        call xtw84(fi0, fk0, alt, wlat, wlon, walt)
        fi0 = wlat
        fk0 = wlon
        call xtw84(fi1, fk1, alt, wlat, wlon, walt)
        fi1 = wlat
        fk1 = wlon
        call xtw84(fi2, fk2, alt, wlat, wlon, walt)
        fi2 = wlat
        fk2 = wlon
        call xtw84(fi3, fk3, alt, wlat, wlon, walt)
        fi3 = wlat
        fk3 = wlon
      else if (ncc.ge.200.and.nccn.lt.200) then
        call xw84t(fi0, fk0, alt, tlat, tlon, talt)
        fi0 = tlat
        fk0 = tlon
        call xw84t(fi1, fk1, alt, tlat, tlon, talt)
        fi1 = tlat
        fk1 = tlon
        call xw84t(fi2, fk2, alt, tlat, tlon, talt)
        fi2 = tlat
        fk2 = tlon
        call xw84t(fi3, fk3, alt, tlat, tlon, talt)
        fi3 = tlat
        fk3 = tlon
      endif
      call cvinit(nccn,
     *            float(iorgn),float(korgn),float(ispan),float(ispbn))
      call cviken(fi0, fk0, y0, x0)
      call cviken(fi1, fk1, y1, x1)
      call cviken(fi2, fk2, y2, x2)
      call cviken(fi3, fk3, y3, x3)
      xl = amin1(x0, x1, x2, x3)
      xh = amax1(x0, x1, x2, x3)
      yl = amin1(y0, y1, y2, y3)
      yh = amax1(y0, y1, y2, y3)
      ixl = ifix(xl)
      if (xl.lt.0.) ixl = ixl - 1
      ixh = ifix(xh)
      if (xh.gt.0.) ixh = ixh + 1
      iyl = ifix(yl)
      if (yl.lt.0.) iyl = iyl - 1
      iyh = ifix(yh)
      if (yh.gt.0.) iyh = iyh + 1
      xl = float(ixl)
      xh = float(ixh)
      yl = float(iyl)
      yh = float(iyh)
      ixl = ixl * 1000
      ixh = ixh * 1000
      iyl = iyl * 1000
      iyh = iyh * 1000
      mxm = (ixh-ixl)/mszx2 + 1 
      mym = (iyh-iyl)/mszy2 + 1 
      if(mxm.gt.LXN.or.mym.gt.LYN) call abendm('too large array size')
c
      open(2,file=fnam2,status='new')
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'TXPROJ : translate grid-data into new projection'
      write(6,'(a,a)') ' Source infile   : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xs, ys, mszx, mszy, mx,  my, salt,
     *              nccn, iorgn,korgn, ispan,ispbn,
     *              xl, yl, mszx2, mszy2, mxm, mym, saltn
      n = 0
c
    2 read(1,'(a)') buf
      if(buf(1:1).ne.'#') goto 3
      write(2,'(a)') buf(1:lrtrim(buf))
      goto 2
    3 n = n + 1
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if(ncc1.ge.200) nc1 = ncc1 - 200
      if(nc1.ge.1.and.nc1.le.62) then
        iorg1 = 0
        ispa1 = 0
        ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if(nc1.gt.60) korg1 = 0
        if(nc1.eq.61) iorg1 = +5400
        if(nc1.eq.62) iorg1 = -5400
      endif
      if(area1.ne.area.or.ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header1 inconsistent')
      read(1,'(a)') buf
      read(buf,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1, alt1
      if(ixs1.ne.ixs.or.iys1.ne.iys .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mx1.ne.mx.or.my1.ne.my)
     *     call abendm('Header2 inconsistent')
c
      read(1,*) ((f(i,j),i=1,mx),j=1,my)
      call dpcini('... calc. new grid data : ')
      do 100 j=1,mym
        y = yl + float((j-1)*mszy2)/1000.
        do 100 i=1,mxm
          x = xl + float((i-1)*mszx2)/1000.
          call cvinit(nccn,
     *             float(iorgn),float(korgn),float(ispan),float(ispbn))
          call cvenik(y, x, fi, fk)
          if (ncc.lt.200.and.nccn.ge.200) then
            call xw84t(fi, fk, alt, tlat, tlon, talt)
            fi = tlat
            fk = tlon
          else if (ncc.ge.200.and.nccn.lt.200) then
            call xtw84(fi, fk, alt, wlat, wlon, walt)
            fi = wlat
            fk = wlon
          endif
          call cvinit(ncc,
     *             float(iorg),float(korg),float(ispa),float(ispb))
          call cviken(fi, fk, yp, xp)
          sx = (xp-xs)*1000./smx
          ip = ifix(sx + 1.)
          dx = sx - float(ip-1)
          sy = (yp-ys)*1000./smy
          jp = ifix(sy + 1.)
          dy = sy - float(jp-1)
          g(i,j) = vnul1
          if(ip.le.0.or.ip.ge.mx) goto 6
          if(jp.le.0.or.jp.ge.my) goto 6
          if(f(ip,jp).eq.vnul1) goto 6
          if(f(ip,jp+1).eq.vnul1) goto 6
          if(f(ip+1,jp).eq.vnul1) goto 6
          if(f(ip+1,jp+1).eq.vnul1) goto 6
          g(i,j) = f(ip,jp)*(1.-dx)*(1.-dy) + f(ip+1,jp)*dx*(1-dy)
     *              + f(ip,jp+1)*(1.-dx)*dy + f(ip+1,jp+1)*dx*dy
    6     call dpcent((j-1)*mxm+i, mxm*mym)
  100 continue
c
      if (n.eq.1) then
        il = mxm
        ih = 1
        jl = mym
        jh = 1
        do 200 j=1,mym
        do 200 i=1,mxm
          if (g(i,j).eq.vnul) goto 200
          if (i.lt.il) il = i
          if (i.gt.ih) ih = i
          if (j.lt.jl) jl = j
          if (j.gt.jh) jh = j
  200   continue
        in = 1000 / mszx2
        jn = 1000 / mszy2
        il = (il-1)/in*in + 1
        ih = (ih+in-2)/in*in + 1
        jl = (jl-1)/jn*jn + 1
        jh = (jh+jn-2)/jn*jn + 1
        ixln = ixl
        iyln = iyl
        xln = xl
        yln = yl
        mxn = mxm
        myn = mym
        if (il.ne.1.or.ih.ne.mxn.or.jl.ne.1.or.jh.ne.myn) then
          ixln = ixl + (il-1)*mszx2
          iyln = iyl + (jl-1)*mszy2
          xln = float(ixln)/1000.
          yln = float(iyln)/1000.
          mxn = ih - il + 1
          myn = jh - jl + 1
          write(6,'(3x,a,2f10.2,2i8,2i6)')
     *         'reduced into', xln, yln, mszx2, mszy2, mxn, myn
        endif
      endif
      write(2,'(2a)') '# Source : ', fnam1(1:lrtrim(fnam1))
      write(2,'(a)')  '# prog.TXPROJ applied to above file'
      write(2,'(a8,i4,4x,4i8)') area, nccn, iorgn,korgn, ispan,ispbn
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixln,iyln,mszx2,mszy2,mxn,myn,vnul1,alt1
      do 300 j=jl,jh
        write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=il,ih)
  300 continue
      if(n.eq.1.and.alt.eq.0.) goto 2
c
    8 close(1)
      close(2)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a,$)') 'TXPROJ :'
        write(99,'(a)')   ' translate grid-data into new projection'
        write(99,'(a,a)') ' Source infile   : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xs, ys, mszx, mszy, mx,  my, salt,
     *                 nccn, iorgn,korgn, ispan,ispbn,
     *                 xl, yl, mszx2, mszy2, mxn, myn, saltn
        call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        close(99)
      endif
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'source grid ', 2f10.2, 2i8, 2i6, a8/
     *       ' New projection :  ', 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, 'gen.new grid', 2f10.2, 2i8, 2i6, a8)
      end
