!--------( TXPROJ : translate grid-data into new projection )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:,:), g(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, buf
  character(72) :: out;  character(20) :: sdtm;  character(9) :: s9
  character(8) :: area, area1, salt, saltn
  integer :: ldr, i, j, i1, i2, n, ip, jp, il, jl, ih, jh, in, jn
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ncn, nccn, iorgn, korgn, ispan, ispbn
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixt, iyt, mszx2, mszy2, ixl, iyl, ixh, iyh, mxm, mym
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1, ixln, iyln, mxn, myn
  real :: vnul, alt, vnul1, alt1, smx, smy, sx, sy, dx, dy, x, y
  real :: xs, ys, xt, yt, xl, yl, xh, yh, xp, yp, xln, yln
  real :: fi0, fi1, fi2, fi3, fk0, fk1, fk2, fk3, fi, fk
  real :: wlat, wlon, walt, tlat, tlon, talt, x0, x1, x2, x3, y0, y1, y2, y3

  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')
  read(1,'(a)') buf
  do while (buf(1:1) == '#')
    read(1,'(a)') buf
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
    if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
    nc = ncc
    if (ncc < 200) then
      nc = ncc; ncc = nc + 800
    else
      ncc = ncc - 200; nc = ncc
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
    nc = ncc
    if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  endif
  if ((nc >= 1) .and. (nc <= 62)) then
    iorg = 0;  ispa = 0;  ispb = 0
    korg = (nc-30)*360 - 180
    if (nc > 60) korg = 0
    if (nc == 61) iorg = +5400
    if (nc == 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 == 0.) salt = '    var.'
  if (alt < 0.) salt = '   undef'
  rewind(1)

  allocate(f(mx,my))

  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) /= '0') .or. (s9(1:1) /= ' '))   &
    call abendm('unable to handle vnul value')
  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  write(out,'(a14,2a10,2x,2a8,2a6,a8)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
    '--- source data ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  saltn = salt
  if (alt < 0.) then
    alt = -1.;  saltn = '   undef'
  endif
  call premsg('')
  call gparmi('Select new proj. Coord. No. (+800 for Old-TOKYO) ==> ', nccn)
  ncn = nccn
  if (nccn >= 800) ncn = nccn - 800
  ispan = 0;  ispbn = 0
  if ((ncn >= 1) .and. (ncn <= 62)) then
    iorgn = 0
    korgn = (ncn-30)*360 - 180
    if (ncn > 60) korgn = 0
    if (ncn == 61) iorgn = +5400
    if (ncn == 62) iorgn = -5400
  else if (ncn == 65) then
    iorgn = 0
    call gparmi2('  Central Longitude ?  (Deg. Min.) : ', i1,i2)
    korgn = i1*60 + i2
  else
    if (ncn == 72) call abendm('nc=72 not implemented')
    call gparmi2('  Origin Latitude ?   (Deg. Min.) : ', i1,i2)
    iorgn = i1*60 + i2
    call gparmi2('  Origin Longitude ?  (Deg. Min.) : ', i1,i2)
    korgn = i1*60 + i2
  endif
  call gparmi2('  Mesh sizes in NS and EW (m) ?  [mszx mszy] : ', mszx2, mszy2)
  if ((mszx2 <= 0) .or. (mszy2 <= 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 >= 800) .and. (nccn < 800)) 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 < 800) .and. (nccn >= 800)) 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 < 0.) ixl = ixl - 1
  ixh = ifix(xh)
  if (xh > 0.) ixh = ixh + 1
  iyl = ifix(yl)
  if (yl < 0.) iyl = iyl - 1
  iyh = ifix(yh)
  if (yh > 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 

  allocate(g(mxm,mym))

  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,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    'Areaname :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'source grid ', xs, ys, mszx, mszy, mx, my, salt
  write(6,'(a20,6x,a6,i4,4x,4i8)')   &
    ' New projection :   ', 'Coord.', nccn, iorgn,korgn, ispan,ispbn
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'gen.new grid', xl, yl, mszx2, mszy2, mxm, mym, saltn

  n = 0
  do while (n < 2)
    read(1,'(a)') buf
    do while (buf(1:1) == '#')
      write(2,'(a)') buf(1:lrtrim(buf))
      read(1,'(a)') buf
    enddo
    n = n + 1
    if (buf(13:16) == '    ') then      !-- v2005
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
      nc1 = ncc1
      if (ncc1 < 200) then
        nc1 = ncc1; ncc1 = nc1 + 800
      else
        ncc1 = ncc1 - 200; nc1 = ncc1
      endif
    else                                !-- v2018
      read(buf,'(a8,4x,i4,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
    endif
    if ((nc1 >= 1) .and. (nc1 <= 62)) then
      iorg1 = 0;  ispa1 = 0;  ispb1 = 0
      korg1 = (nc1-30)*360 - 180
      if (nc1 > 60) korg1 = 0
      if (nc1 == 61) iorg1 = +5400
      if (nc1 == 62) iorg1 = -5400
    endif
    if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg) .or.   &
        (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
      call abendm('Header1 inconsistent')
    read(1,'(a)') buf
    read(buf,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1, alt1
    if ((ixs1 /= ixs) .or. (iys1 /= iys) .or. (mszx1 /= mszx) .or.   &
        (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
      call abendm('Header2 inconsistent')
    read(1,*) ((f(i,j),i=1,mx),j=1,my)

    call dpcini('... calc. new grid data : ')
    do j=1,mym
      y = yl + float((j-1)*mszy2)/1000.
      do 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 >= 800) .and. (nccn < 800)) then
          call xw84t(fi, fk, alt, tlat, tlon, talt)
          fi = tlat;  fk = tlon
        else if ((ncc < 800) .and. (nccn >= 800)) 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 > 0) .and. (ip < mx) .and. (jp > 0) .and. (jp < my)   &
            .and. (f(ip,jp) /= vnul1) .and. (f(ip,jp+1) /= vnul1)   &
            .and. (f(ip+1,jp) /= vnul1) .and. (f(ip+1,jp+1) /= vnul1))   &
          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
        call dpcent((j-1)*mxm+i, mxm*mym)
      enddo
    enddo

    if (n == 1) then
      il = mxm;  ih = 1;  jl = mym;  jh = 1
      do j=1,mym;  do i=1,mxm
        if (g(i,j) == vnul) cycle
        if (i < il) il = i;  if (i > ih) ih = i
        if (j < jl) jl = j;  if (j > jh) jh = j
      enddo;  enddo
      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;  xln = xl;  mxn = mxm
      iyln = iyl;  yln = yl;  myn = mym
      if ((il /= 1) .or. (ih /= mxn) .or. (jl /= 1) .or. (jh /= 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,'(4x,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,4x,i4,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 j=jl,jh
      write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=il,ih)
    enddo
    if (alt /= 0.) n = 2
  enddo
  close(1)
  close(2)

  if (flog(ldr:ldr) /= ' ') 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,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'source grid ', xs, ys, mszx, mszy, mx, my, salt
    write(99,'(a20,6x,a6,i4,4x,4i8)')   &
      ' New projection :   ', 'Coord.', nccn, iorgn,korgn, ispan,ispbn
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'gen.new grid', xln, yln, mszx2, mszy2, mxn, myn, saltn
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end