!--------( EIMGA : preparation for 3D imaging  with Auto-scaling )--------
  implicit none
  integer,parameter :: KVER = 2018
  integer :: lwkdir, lrtrim;  real calma
  real,allocatable :: g(:,:), ds(:), dt(:)
  real,allocatable :: t(:,:), s(:,:), w(:,:), th(:,:), f(:,:), h(:,:)
  character(400) :: buf
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, fnam5
  character(72) :: out;  character(40) :: strs
  character(25) :: strt;  character(20) :: sdtm
  character(8) :: area, tarea, area1, salt
  integer :: ldr, i, j, ia, ja, l, ixc, iyc, nhx, nhy, kds, nds
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: nct, ncct, iorgt, korgt, ispat, ispbt
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixl, iyl, mszxm, mszym, mxm, mym
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  integer :: ixst, iyst, mszxt, mszyt, mxt, myt
  integer :: mszxh, mszyh, kic, kjc, ki, kj, ir, jr
  real :: vnul, alt, hnul, tnul, ht, ht1, ht2, ha, hb, hc
  real :: xs, ys, xt, yt, fmx, fmy, xst, yst, xtt, ytt, fmxt, fmyt
  real :: xl, yl, xh, yh, fmxm, fmym, dx, dy, x, y, xx, yy
  real :: am, hw, ct, dip, dec, the, phi
  real :: v, vm, sz, za, ss, ggn

  call premsg('EIMGA : preparation for 3D imaging  with Auto-scaling')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1;  flog = wdr
  fnam1 = wdr;  fnam2 = wdr;  fnam3 = wdr;  fnam4 = wdr;  fnam5 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Mag.Anomaly data filename   ==> ', 81-ldr, fnam1(ldr:80))
  open(12,file=fnam1,status='old')
  read(12,'(a)') buf
  do while (buf(1:1) == '#')
    read(12,'(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(12,'(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'
  mszxh = mszx / 2;  xs = float(ixs)/1000.
  mszyh = mszy / 2;  ys = float(iys)/1000.

  allocate(f(mx,my));  allocate(h(mx,my))

  read(12,*) ((f(i,j),i=1,mx),j=1,my)
  if (alt == 0.) then
    read(12,'(a)') buf
    do while (buf(1:1) == '#')
      read(12,'(a)') buf
    enddo
    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
    read(12,*) ixs1, iys1, mszx1, mszy1, mx1, my1, hnul
    if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg) .or.   &
        (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb) .or.   &
        (ixs1 /= ixs) .or. (iys1 /= iys) .or. (mszx1 /= mszx) .or.   &
        (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
       call abendm('Mag.Anom Header inconsistent')
    read(12,*) ((h(i,j),i=1,mx),j=1,my)
    strt = ' given as 2nd set data   '
  else
    if (alt < 0.) call abendm('Undefined Altitude')
    write(strt,'(f7.0,a)') alt, ' m above sea level'
  endif
  close(12)
  call premsg('  obs.surface : ' // strt)

  call gparma('Enter Src.Surf.Alt data filename  ==> ', 81-ldr, fnam2(ldr:80))
  open(10,file=fnam2,status='old')
  read(10,'(a)') buf
  do while (buf(1:1) == '#')
    read(10,'(a)') buf
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    if ((ncct<0) .or. (ncct>=400)) call abendm('invalid coord.no.')
    nct = ncct
    if (ncct < 200) then
      nct = ncct; ncct = nct + 800
    else
      ncct = ncct - 200; nct = ncct
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    nct = ncct
    if ((ncct >= 800) .and. (ncct < 1000)) nct = ncct - 800
  endif
  if ((nct >= 1) .and. (nct <= 62)) then
    iorgt = 0;  ispat = 0;  ispbt = 0
    korgt = (nct-30)*360 - 180
    if (nct > 60) korgt = 0
    if (nct == 61) iorgt = +5400
    if (nct == 62) iorgt = -5400
  endif
  read(10,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul
  xst = float(ixst)/1000.;  yst = float(iyst)/1000.

  allocate(t(mxt,myt))

  read(10,*) ((t(i,j),i=1,mxt),j=1,myt)
  close(10)

  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
     '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  if ((ncct /= ncc) .or. (iorgt /= iorg) .or. (korgt /= korg)   &
      .or. (ispat /= ispa) .or. (ispbt /= ispb)) then
    write(out,'(12x,a8,12x,i4,4x,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    call premsg(out)
    call abendm('Header inconsistent')
  endif
  if (tarea /= area) call premsg('           ' // tarea)
  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)')   &
     '--- mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6)')   &
     '--- source alt. ', xst, yst, mszxt, mszyt, mxt, myt
  call premsg(out)
  xt = xs + float(mszx*(mx-1))/1000.
  yt = ys + float(mszy*(my-1))/1000.
  fmx = float(mszx);  fmy = float(mszy)
  xtt = xst + float(mszxt*(mxt-1))/1000.
  ytt = yst + float(mszyt*(myt-1))/1000.
  fmxt = float(mszxt);  fmyt = float(mszyt)
  if ((xst > xs) .or. (yst > ys) .or. (xtt < xt) .or. (ytt < yt))   &
    call abendm('ranges conflict')

  call premsg('')
  call premsg('Specify Source distribution Parameters')
  call premsg('  Southwest corner Coord. in km ?')
  call gparmf('    Northing : ', xl)
  ixl = nint(xl*1000.)
  xl = float(ixl)/1000.
  if ((xl < xst) .or. (xl > xs)) call abendm('ranges conflict')
  call gparmf('    Easting  : ', yl)
  iyl = nint(yl*1000.)
  yl = float(iyl)/1000.
  if ((yl < yst) .or. (yl > ys)) call abendm('ranges conflict')
  call premsg('  Mesh size and numbers to N ?')
  call gparmi('    Unit Mesh size in m :  ', mszxm)
  if (mszxm <= 0) call abendm('invalid value')
  fmxm = float(mszxm)
  call gparmi('    Number of Mesh to N :  ', mxm)
  xh = xl + float(mszxm*(mxm-1))/1000.
  if ((xh < xt) .or. (xh > xtt)) call abendm('ranges conflict')
  call premsg('  Mesh size and numbers to E ?')
  call gparmi('    Unit Mesh size in m :  ', mszym)
  if (mszym <= 0) call abendm('invalid value')
  fmym = float(mszym)
  call gparmi('    Number of Mesh to E :  ', mym)
  yh = yl + float(mszym*(mym-1))/1000.
  if ((yh < yt) .or. (yh > ytt)) call abendm('ranges conflict')
  call premsg('  Truncation of Source effect ?')
  call gparmf('    Half width (in km) :  ', hw)
  nhx = nint(hw*1000./fmx);  nhy = nint(hw*1000./fmy)
  if ((nhx < 10) .or. (nhy < 10)) call abendm('truncation too short')
  if ((nhx > 999) .or. (nhy > 999)) call abendm('truncation too long')

  allocate(g(-nhx:+nhx,-nhy:+nhy))
  allocate(s(mxm,mym));  allocate(w(mxm,mym));  allocate(th(mxm,mym))

  call premsg('Specify Direction Parameters (in degrees)')
  call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
  if ((dip < -90.) .or. (dip > 90.)) call abendm('invalid value')
  if ((dec < -180.) .or. (dec >= 360.)) call abendm('invalid value')
  call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
  if ((the < -90.) .or. (the > 90.)) call abendm('invalid value')
  if ((phi < -180.) .or. (phi >= 360.)) call abendm('invalid value')
  call premsg('')
  call premsg('Specify Source-Layer configuration')
  call premsg('    Horiz.Flat (0), Const.Thickness below Surf. (1),')
  call gparmi('     or Var.Thickness below Surf. (2) : ', kds)
  if ((kds < 0) .or. (kds > 2)) call abendm('invalid value')
  call gparmi('  Number of layers : ', nds)
  if ((nds<1) .or. (nds>999)) call abendm('invalid value')

  allocate(ds(0:nds));  allocate(dt(nds))
  do i=1,nds; dt(i) = -99999.; enddo

  if (kds == 0) then
    call premsg('  Enter Bottom depths of Layers (in m BSL)')
    call parmin('  (shallow to deep) : ', 400, buf)
    ds(0) = -99999.
    read(buf,*) (dt(i),i=1,nds)
    do i=1,nds
      ds(i) = dt(i)
      if (ds(i) <= ds(i-1)) call abendm('depths conflict')
    enddo
    strs = ' variable depth horizon below sea level '
  else if (kds == 1) then
    call gparmf('  Const. thickness of Layers (in m) : ', ct)
    if (ct <= 0.) call abendm('invalid value')
    ds(0) = 0.
    do i=1,nds
      ds(i) = ds(i-1) + ct
    enddo
    strs = ' constant thickness below terr.surface  '
  else
    call premsg('  Enter Thicknesses of Layers (in m)')
    call parmin('  (shallow to deep) : ', 400, buf)
    ds(0) = 0.
    read(buf,*) (dt(i),i=1,nds)
    do i=1,nds
      if (dt(i) <= 0.) call abendm('invalid thickness')
      ds(i) = ds(i-1) + dt(i)
    enddo
    strs = ' variable thickness below terr.surface  '
  endif
  call gparmf('Enter Initial source magnetization (A/m) ==> ', vm)

  call premsg('')
  call gparma('Enter CFIM matrix output filename ==> ', 81-ldr, fnam3(ldr:80))
  open(20,file=fnam3,form='unformatted',status='new')
  call gparma('Enter AIMG model output filename  ==> ', 81-ldr, fnam4(ldr:80))
  call clspin()
  call strdtm(sdtm)
  open(21,file=fnam4,status='new')
  call gparma('Enter FSCL Par.Scal.Coefs outfile ==> ', 81-ldr, fnam5(ldr:80))
  open(22,file=fnam5,form='unformatted',status='new')
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'EIMGA : preparation for 3D imaging  with Auto-scaling'
  write(6,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Src.Surf.Alt infile : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' CFIM matrix outfile : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' AIMG model outfile  : ', fnam4(1:lrtrim(fnam4))
  write(6,'(a,a)') ' FSCL Par.Scal.Coefs : ', fnam5(1:lrtrim(fnam5))

  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)')   &
    'mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  write(6,'(4x,a12,2f10.2,2i8,2i6)')   &
    'src.surf.alt', xst, yst, mszxt, mszyt, mxt, myt
  write(6,'(4x,a12,2f10.2,2i8,2i6,i5,a3)')   &
    'source model', xl, yl, mszxm, mszym, mxm, mym, nds, 'Lys'
  write(6,'(a,f6.2,a,i4,a,i4,a)')   &
    ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  mag.mesh )'
  write(6,'(a,f6.1,3x,a,f7.1)') '  ambient field :  dip=', dip, 'dec=', dec
  write(6,'(a,f6.1,3x,a,f7.1)') '  magnetization :  dip=', the, 'dec=', phi
  write(6,'(a,a25)') '  obs. surface  : ', strt
  write(6,'(a,a40)') '  layer config. : ', strs
  if (kds == 0) then
    write(6,'(4x,a)') 'bottom depth horizons of layers (m):'
    write(6,'((5x,10f7.0))') (ds(i),i=1,nds)
  else if (kds == 1) then
    write(6,'(4x,i3,a,f7.0,a)')   &
       nds, ' layers with const.thickness of', ct, ' (m)'
  else
    write(6,'(4x,a)') 'thicknesses of layers (m):'
    write(6,'((5x,10f7.0))') (dt(i),i=1,nds)
  endif
  write(6,'(a,f8.2)') ' Initial source magnetization (A/m) :', vm

  call magafd(the, phi, dip, dec)
  write(20) area, KVER, alt, nhx,nhy,nds, ncc,iorg,korg,ispa,ispb,   &
          ixs,iys,mszx,mszy,mx,my, ixl,iyl,mszxm,mszym,mxm,mym
  write(22) area, KVER, ncc, iorg,korg, ispa,ispb, nds,   &
          ixl, iyl, mszxm, mszym, mxm, mym
  am = fmxm*fmym

  do j=1,mym
    iyc = iyl + (j-1)*mszym;  ja = (iyc-iyst) / mszyt
    dy = float(iyc-iyst - ja*mszyt) / fmyt
    do i=1,mxm
      ixc = ixl + (i-1)*mszxm;  ia = (ixc-ixst) / mszxt
      dx = float(ixc-ixst - ia*mszxt) / fmxt
      ht = t(ia+1,ja+1)
      if (ht == tnul) call abendm('source alt. undefined')
      if (dx > 0.) then
        if ((ia+2) > mxt) call abendm('source alt. undefined')
        ht2 = t(ia+2,ja+1)
        if (ht2 == tnul) call abendm('source alt. undefined')
        ht = ht*(1.-dx) + ht2*dx
      endif
      if (dy > 0.) then
        if ((ja+2) > myt) call abendm('source alt. undefined')
        ht1 = t(ia+1,ja+2)
        if (ht1 == tnul) call abendm('source alt. undefined')
        if (dx > 0.) then
          ht2 = t(ia+2,ja+2)
          if (ht2 == tnul) call abendm('source alt. undefined')
          ht1 = ht1*(1.-dx) + ht2*dx
        endif
        ht = ht*(1.-dy) + ht1*dy
      endif
      s(i,j) = ht
    enddo
  enddo

  call dpcini('... gen. CFIM matrix data : ')
  do l=1,nds
    do j=1,mym
      iyc = iyl + (j-1)*mszym;  kjc = (iyc-iys + mszyh) / mszy
      do i=1,mxm
        ixc = ixl + (i-1)*mszxm;  kic = (ixc-ixs + mszxh) / mszx
        if (kds == 0) then
          if ((s(i,j)+ds(l)) <= 0.) then
            write(20) i,j,l, kic,kjc, ((0.,ki=-nhx,nhx),kj=-nhy,nhy)
            w(i,j) = 0.;  th(i,j) = 1.
            call dpcent(((l-1)*mym+(j-1))*mxm+i, mxm*mym*nds);  cycle
          else
            ha = amin1(s(i,j), -ds(l-1));  hb = -ds(l)
          endif
        else
          ha = s(i,j) - ds(l-1);  hb = s(i,j) - ds(l)
        endif
        do kj=-nhy,nhy
          jr = kjc + kj;  y = float(iyc-iys - jr*mszy);  yy = y * y
          do ki=-nhx,nhx
            ir = kic + ki;  x = float(ixc-ixs - ir*mszx);  xx = x * x
            ss = xx + yy;  v = 0.
            if ((ir < 0) .or. (ir >= mx) .or. (jr < 0) .or. (jr >= my)   &
                 .or. (f(ir+1,jr+1) == vnul)) then
              g(ki,kj) = 0.;  cycle
            endif
            hc = alt
            if (alt == 0.) hc = h(ir+1,jr+1)
            za = hc - ha;  sz = ha - hb
            if ((iabs(ki) <= 10) .and. (iabs(kj) <= 10)) then
              call mprism(0.01, x, fmxm, y, fmym, za, sz)
            else
              call mvline(0.01, x, y, za, sz, am)
            endif
            g(ki,kj) = calma(0., 0., 0.)
          enddo
        enddo
        ggn = 0.
        do kj=-nhy,nhy;  do ki=-nhx,nhx
          ggn = ggn + g(ki,kj)*g(ki,kj)
        enddo;  enddo
        if (ggn == 0.) then
          w(i,j) = 0.
        else
          w(i,j) = 1. / sqrt(ggn)
        endif
        th(i,j) = ha - hb
        write(20) i,j,l, kic,kjc, ((g(ki,kj),ki=-nhx,nhx),kj=-nhy,nhy)
        call dpcent(((l-1)*mym+(j-1))*mxm+i, mxm*mym*nds)
      enddo
    enddo
    write(22) ((w(i,j),i=1,mxm),j=1,mym)
    write(22) ((th(i,j),i=1,mxm),j=1,mym)
  enddo
  close(20);  close(22)

  write(21,'(a,5(1x,i5))') '#AIMG:  ', kds, nds, KVER, nhx, nhy
  l = (nds-1) / 8
  do i=1,l
    write(21,'(a,8(1x,f7.1))') '#AIMG:  ', (ds(j),j=i*8-7,i*8)
  enddo
  write(21,'(a,8(1x,f7.1))') '#AIMG:  ', (ds(j),j=l*8+1,nds)
  write(21,'(a)') '#'
  write(21,'(a8,4x,i4,4i8,8x,2i4)') area, ncc, iorg,korg, ispa,ispb, nds,0
  write(21,'(2i12,4i6,1x,f7.1,1x,f7.0,i6)')   &
     ixl, iyl, mszxm, mszym, mxm, mym, 9999.9, -1., 0
  do j=1,mym
    write(21,'((f7.1,9(1x,f7.1)))') (s(i,j),i=1,mxm)
  enddo
  do l=1,nds
    write(21,'(a8,4x,i4,4i8,8x,2i4)') area, ncc, iorg,korg, ispa,ispb, nds,l
    write(21,'(2i12,4i6,1x,f7.1,1x,f7.0,i6)')   &
       ixl, iyl, mszxm, mszym, mxm, mym, 9999.9, -1., 0
    do j=1,mym
      write(21,'((f7.1,9(1x,f7.1)))') ((vm*100.),i=1,mxm)
    enddo
  enddo
  close(21)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'EIMGA : preparation for 3D imaging with Auto-scaling'
    write(99,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Src.Surf.Alt infile : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' CFIM matrix outfile : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' AIMG model outfile  : ', fnam4(1:lrtrim(fnam4))
    write(99,'(a,a)') ' FSCL Par.Scal.Coefs : ', fnam5(1:lrtrim(fnam5))
    write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
    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)')   &
      'mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
    write(99,'(4x,a12,2f10.2,2i8,2i6)')   &
      'src.surf.alt', xst, yst, mszxt, mszyt, mxt, myt
    write(99,'(4x,a12,2f10.2,2i8,2i6,i5,a3)')   &
      'source model', xl, yl, mszxm, mszym, mxm, mym, nds, 'Lys'
    write(99,'(a,f6.2,a,i4,a,i4,a)')   &
      ' Half width of calc.window :', hw, ' km (', nhx,   &
      ' *', nhy, '  mag.mesh )'
    write(99,'(a,f6.1,3x,a,f7.1)') '  ambient field :  dip=', dip, 'dec=', dec
    write(99,'(a,f6.1,3x,a,f7.1)') '  magnetization :  dip=', the, 'dec=', phi
    write(99,'(a,a25)') '  obs. surface  : ', strt
    write(99,'(a,a40)') '  layer config. : ', strs
    if (kds == 1) then
      write(99,'(18x,i3,4h  * ,f8.0)') nds, ct
    else
      write(99,'((8x,8f8.0))') (dt(i),i=1,nds)
    endif
    write(99,'(a,f8.2)') ' Initial source magnetization (A/m) :', vm
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end