!--------( CIMGS : min.volume imaging )--------
  implicit none
  integer,parameter :: KVER = 2018
  real,parameter :: THRSH = 2.
  integer :: lwkdir, lrtrim
  real,allocatable :: g(:,:), ds(:), b(:,:), r(:,:), q(:,:), t(:,:)
  real,allocatable :: s(:,:,:), p(:,:,:), v(:,:,:), w(:,:,:)
  real,allocatable :: u(:,:,:), z(:,:,:), th(:,:,:), cc(:,:,:), h(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, fnam5, fnam6, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, area1, salt;  character(1) :: yn
  real :: coef(3)
  integer :: ldr, i, j, l, n, nhx, nhy, kds, nds, nver, IER
  integer :: nhx1, nhy1, nds1
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  integer :: ixl, iyl, mszxm, mszym, mxm, mym
  integer :: ixl1, iyl1, mszxm1, mszym1, mxm1, mym1
  integer :: ilp, ip, jp, k, klp, kout, l1, lp, nfx, nfy, nloop, nlp
  real :: anul, alt, alt1, a, af, bf, cf, cmx, cmy, dmy, dv
  real :: ft, hw, p1, p2, rs, ss, vm, vnul, wp, x, y, xl, yl, xs, ys
  real :: del, del2, ee, ratio, rv, ss1, ss2, sza, tr, trp, uvol
  integer :: kc

  call premsg('CIMGS : min.volume imgaging')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr;  fnam3 = wdr
  fnam4 = wdr;  fnam5 = wdr;  fnam6 = 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
  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size')

  allocate(b(mx,my));  allocate(r(mx,my))
  allocate(q(mx,my));  allocate(t(mx,my))

  read(12,*) ((b(i,j),i=1,mx),j=1,my)
  close(12)
  call prompt('  Remove Linear trend ?  ')
  call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
  if ((yn /= 'y') .and. (yn /= 'n')) call abendm('invalid parameter')
  call gparma('Enter CFIM matrix data filename   ==> ', 81-ldr, fnam2(ldr:80))
  open(20,file=fnam2,form='unformatted',status='old')
  read(20) area1, nver, alt1, nhx,nhy,nds, ncc1,iorg1,korg1,ispa1,ispb1,   &
          ixs1,iys1,mszx1,mszy1,mx1,my1, ixl,iyl,mszxm,mszym,mxm,mym
  if (nver /= KVER) call abendm('Version mismatch')
  if ((area1 /= area) .or. (ncc1 /= ncc) .or.   &
      (iorg1 /= iorg) .or. (korg1 /= korg) .or.   &
      (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
    call abendm('Header inconsistent')
  if ((ixs1 /= ixs) .or. (iys1 /= iys) .or.   &
      (mszx1 /= mszx) .or. (mszy1 /= mszy) .or.   &
      (mx1 /= mx) .or. (my1 /= my))   &
    call abendm('Header inconsistent')
  hw = float(mszx*nhx) / 1000.
  nfx = nhx + 1 + nhx;  nfy = nhy + 1 + nhy

  allocate(g(nfx,nfy))
  allocate(s(mxm,mym,nds)); allocate(p(mxm,mym,nds))
  allocate(v(mxm,mym,nds)); allocate(w(mxm,mym,nds))
  allocate(u(mxm,mym,nds)); allocate(z(mxm,mym,nds))
  allocate(th(mxm,mym,nds)); allocate(cc(mxm,mym,nds))
  allocate(h(mxm,mym)); allocate(ds(nds))

  if (alt < 0.) alt = alt1
  write(salt,'(1x,f7.0)') alt
  if (alt == 0.) salt = '    var.'
  if (alt < 0.) salt = '   undef'
  xs = float(ixs)/1000.;  ys = float(iys)/1000.
  xl = float(ixl)/1000.;  yl = float(iyl)/1000.
  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)')   &
    '--- mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  write(out,'(a,2f10.2,2i8,2i6,i5,a)')   &
    '--- source model', xl,yl, mszxm,mszym, mxm,mym, nds, 'Lys'
  call premsg(out)

  call gparma('Enter AIMG model input filename   ==> ', 81-ldr, fnam3(ldr:80))
  open(21,file=fnam3,status='old')
  n = 0;  nver = 0;  l = (nds+7) / 8
  read(21,'(a)') buf
  do while (buf(1:1) == '#')
    if (buf(1:6) == '#AIMG:') then
      if (n == 0) then
        read(buf(7:80),*) kds, nds1, nver, nhx1, nhy1
        if (nver /= KVER) call abendm('AIMG Version mismatch')
        if ((nds1 /= nds) .or. (nhx1 /= nhx) .or. (nhy1 /= nhy))   &
          call abendm('AIMG conflicts with CFIM info')
      else if (n < l) then
        read(buf(7:80),*) (ds(i),i=n*8-7,n*8)
      else if (n == l) then
        read(buf(7:80),*) (ds(i),i=n*8-7,nds)
      endif
      n = n + 1
    endif
    read(21,'(a)') buf
  enddo
  if (n <= l) call abendm('AIMG lacks source config info')
  do l=0,nds
    if (l > 0) read(21,'(a)') buf
    if (buf(13:16) == '    ') then        !-- v2005
      read(buf,'(a8,i4,4x,4i8,8x,2i4)')   &
          area1, ncc1, iorg1,korg1, ispa1,ispb1, nds1,l1
      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,8x,2i4)')   &
          area1, ncc1, iorg1,korg1, ispa1,ispb1, nds1,l1
      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('Header inconsistent')
    read(21,*) ixl1, iyl1, mszxm1, mszym1, mxm1, mym1, anul, dmy, klp
    if ((ixl1 /= ixl) .or. (iyl1 /= iyl) .or.   &
        (mszxm1 /= mszxm) .or. (mszym1 /= mszym) .or.   &
        (mxm1 /= mxm) .or. (mym1 /= mym))   &
      call abendm('Header inconsistent')
    if ((nds1 /= nds) .or. (l1 /= l)) call abendm('AIMG file broken')
    if (l == 0) then
      read(21,*) ((h(i,j),i=1,mxm),j=1,mym)
    else
      read(21,*) ((s(i,j,l),i=1,mxm),j=1,mym)
    endif
  enddo
  call premsg('Enter AIMG model output(new) filename')
  call gparma('   or  Enter ''U'' to update infile ==> ',   &
              81-ldr, fnam6(ldr:80))
  if (fnam6(ldr:lrtrim(fnam6))  ==  'U') then
    rewind(21)
  else
    close(21)
    open(21,file=fnam6,status='new')
  endif
  sza = float(mszxm) * float(mszym)
  uvol = sza * sqrt(sza)

  call gparma('Enter FSCL Par.Scal.Coefs infile  ==> ', 81-ldr, fnam4(ldr:80))
  open(22,file=fnam4,form='unformatted',status='old')
  read(22) area1, nver, ncc1, iorg1,korg1, ispa1,ispb1, nds1,   &
         ixl1, iyl1, mszxm1, mszym1, mxm1, mym1
  if (nver /= KVER) call abendm('Version mismatch')
  nc1 = ncc1
  if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  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('Header inconsistent')
  if ((ixl1 /= ixl) .or. (iyl1 /= iyl) .or.   &
      (mszxm1 /= mszxm) .or. (mszym1 /= mszym) .or.   &
      (mxm1 /= mxm) .or. (mym1 /= mym) .or. (nds1 /= nds))   &
    call abendm('Header inconsistent')
  call gparmf('Specify Par.Scale Weighting power ==> ', wp)
  if (wp < 0.) call abendm('Minus power Invalid')
  do l=1,nds
    read(22) ((w(i,j,l),i=1,mxm),j=1,mym)
    read(22) ((th(i,j,l),i=1,mxm),j=1,mym)
    do j=1,mym; do i=1,mxm
      v(i,j,l) = th(i,j,l) * (sza/uvol)
      if (w(i,j,l) == 0.) then
        v(i,j,l) = 0.
      else
        v(i,j,l) = v(i,j,l) / exp(alog(w(i,j,l)*v(i,j,l))*wp)
      endif
    enddo; enddo
  enddo
  close(22)

  call gparmf('Source Support Threshold (A/m) ==> ', del)
  del2 = (del*100.) * (del*100.)
  call gparmf('Initial Wt.Ratio of SrcVolume / Res. ==> ', ratio)
  if ((ratio < 0.2) .or. (ratio > 5.)) call abendm('abnormal ratio')

  call premsg('Enter Auxiary(2) output filename')
  call gparma('     ( blank to suppress output ) ==> ', 81-ldr, fnam5(ldr:80))
  kout = 0
  if (fnam5(ldr:ldr) /= ' ') then
    open(8,file=fnam5,status='new')
    kout = 1
  endif
  call premsg('How many loops to excute ?')
  call gparmi('   (minus/zero to restart from loop-0) ==> ', nlp)
  nloop = iabs(nlp)
  if (nlp <= 0) then
    call gparmf(' Initial value of magnetiz. (in A/m) ==> ', vm)
  endif
  call clspin()
  call strdtm(sdtm)

  call dpcini('.... read CFIM matrix (0) : ')
  do l=1,nds; do j=1,mym; do i=1,mxm
    read(20) ip, jp, lp, k, k, g
    if ((ip /= i) .or. (jp /= j) .or. (lp /= l))   &
      call abendm('CFIM file broken')
    call dpcent(((l-1)*mym+(j-1))*mxm+i, mxm*mym*nds)
  enddo; enddo; enddo
  rewind(20)

  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'CIMGS : min.volume imaging'
  write(6,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
  if (yn == 'y') then
    write(6,'(22x,a)') 'Linear trend removed'
  else
    write(6,'(22x,a)') 'without removal of Linear trend'
  endif
  write(6,'(a,a)') ' CFIM matrix infile  : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' AIMG model infile   : ', fnam3(1:lrtrim(fnam3))
  if (fnam6(ldr:lrtrim(fnam6))  ==  'U') then
    write(6,'(a)')   ' AIMG model output   : <UPDATE infile>'
  else
    write(6,'(a,a)') ' AIMG model outfile  : ', fnam6(1:lrtrim(fnam6))
  endif
  write(6,'(a,a)') ' FSCL Par.Scal.Coefs : ', fnam4(1:lrtrim(fnam4))
  write(6,'(a,f4.1)') ' Par.Scale Weighting Power : ', wp
  write(6,'(a,f5.2)') ' Source Support Threshold (A/m) : ', del
  write(6,'(a,f5.2)') ' SrcVolume/Res. Initial Wt.Ratio : ', ratio
  if (kout == 1) write(6,'(a,a)')   &
    ' Auxiary(2) outfile  : ', 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,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,i4)') ' Loop Count : ', nloop
  if (nlp <= 0) then
    klp = 0
    do l=1,nds; do j=1,mym; do i=1,mxm
      s(i,j,l) = vm * 100.
    enddo; enddo; enddo
    write(6,'(a,f6.2,a)')   &
      '   Starting from Initial magnetization :', vm, ' A/m'
  endif

  if (yn == 'y') then
    n = 0;  ss1 = 0.
    do j=1,my; do i=1,mx
      if (b(i,j) /= vnul) then
        ss1 = ss1 + b(i,j);  n = n + 1
      endif
    enddo; enddo
    ft = ss1 / float(n)
    cmx = float(mx+1) / 2.;  cmy = float(my+1) / 2.
    call sm2opn(1, 1)
    do j=1,my
      y = float(j) - cmy
      do i=1,mx
        x = float(i) - cmx
        if (b(i,j) /= vnul) call sm2ex(x, y, b(i,j)-ft)
      enddo
    enddo
    call sm2cls(coef, 3, 1)
    cf = coef(1) + ft;  af = coef(2);  bf = coef(3)
    write(6,'(a,f8.2,a,f7.3,a,f7.3,a)')   &
      ' Trend surface in nT :', cf, ' + (', af, ' * x) + (', bf, ' * y)'
    write(6,'(5x,a)')   &
      '( x/y : N/E coord. in mesh unit referred to the center of area )'
    do j=1,my
      y = float(j) - cmy
      do i=1,mx
        x = float(i) - cmx
        if (b(i,j) /= vnul) b(i,j) = b(i,j) - (af*x + bf*y + cf)
      enddo
    enddo
    if (kout == 1) then
      write(8,'(a8,4x,i4,4i8)') 'Obs-Trnd',ncc,iorg,korg,ispa,ispb
      write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
         ixs,iys, mszx,mszy, mx,my, vnul, -1.
      do j=1,my
        write(8,'((f7.1,9(1x,f7.1)))') (b(i,j),i=1,mx)
      enddo
      write(6,'(a)') 'AUX2 <== Obs-Trnd : Trend-removed Obs.'
    endif
  else
    af = 0.;  bf = 0.;  cf = 0.
  endif

  open(9,status='scratch')
  ss1 = 0.
  do l=1,nds; do j=1,mym; do i=1,mxm
    ss = s(i,j,l) * s(i,j,l)
    cc(i,j,l) = (ss + del2) / v(i,j,l)
    ss1 = ss1 + ss / cc(i,j,l)
  enddo; enddo; enddo
  rv = sqrt(ss1 / float(mxm*mym*nds))
  call conv1(g,nfx,nfy, s, mxm,mym,nds, r, mx,my)
  n = 0;  ss2 = 0.
  do j=1,my; do i=1,mx
    if (b(i,j) == vnul) then
      r(i,j) = 0.
    else
      r(i,j) = b(i,j) - r(i,j)
      ss2 = ss2 + r(i,j)*r(i,j)
      n = n + 1
    endif
  enddo; enddo
  dv = sqrt(ss2 / float(n))
  write(6,'(a5,i4,3x,a13,f7.3,a2,f7.4,3x,a6,f7.2,3x,a6,f7.2)')   &
      ' loop', klp, 'rms.Step/Src.', 0., ' /', rv, 'Misfit', dv, 'T.Res.', dv
  write(9,*) klp, 0., 0., rv, dv, dv

  if (nloop /= 0) then
    ilp = 1;  klp = klp + 1
    call conv2(g,nfx,nfy, r, mx,my, u, mxm,mym,nds)
    do l=1,nds; do j=1,mym; do i=1,mxm
      z(i,j,l) = u(i,j,l) * w(i,j,l);  p(i,j,l) = z(i,j,l)
    enddo; enddo; enddo
    call conv1(g,nfx,nfy, p, mxm,mym,nds, q, mx,my)
    p1 = 0.;  p2 = 0.
    do j=1,my; do i=1,mx
      if (b(i,j) /= vnul) then
        p1 = p1 + r(i,j)*q(i,j);  p2 = p2 + q(i,j)*q(i,j)
      endif
    enddo; enddo
    a = p1 / p2
    ss1 = 0.;  ss2 = 0.
    do l=1,nds; do j=1,mym; do i=1,mxm
      s(i,j,l) = s(i,j,l) + a*p(i,j,l)
      ss = s(i,j,l) * s(i,j,l)
      cc(i,j,l) = (ss + del2) / v(i,j,l)
      ss1 = ss1 + ss / cc(i,j,l)
      ss2 = ss2 + p(i,j,l)*p(i,j,l) / cc(i,j,l)
    enddo; enddo; enddo
    rv = sqrt(ss1 / float(mxm*mym*nds))
    rs = a * sqrt(ss2 / float(mxm*mym*nds))
    n = 0;  ss2 = 0.
    do j=1,my; do i=1,mx
      if (b(i,j) /= vnul) then
        r(i,j) = r(i,j) - a*q(i,j)
        ss2 = ss2 + r(i,j)*r(i,j)
        n = n + 1
      endif
    enddo; enddo
    dv = sqrt(ss2 / float(n))
    write(6,'(a5,i4,3x,a13,f7.3,a2,f7.4,3x,a6,f7.2,3x,a6,f7.2)')   &
      ' loop', klp, 'rms.Step/Src.', rs, ' /', rv, 'Misfit', dv, 'T.Res.', dv
    write(9,*) klp, 0., rs, rv, dv, dv
    ee = ratio * ss2/ss1
    tr = sqrt((ss2 + ss1*ee) / float(n + mxm*mym*nds))
    write(6,'(a5,i4,6x,a25,e10.3,10x,a6,f7.2)')   &
      ' loop', klp, 'Updated hyper-parameter =', ee, 'T.Res.', tr
    write(9,*) klp, ee, rs, rv, dv, tr
    kc = 0;  trp = tr
  endif

  do while ((kc < 5) .and. (dv >= 0.1) .and. (ilp < nloop))
    ilp = ilp + 1;  klp = klp + 1
    call conv2(g,nfx,nfy, r, mx,my, u, mxm,mym,nds)
    do l=1,nds; do j=1,mym; do i=1,mxm
      u(i,j,l) = u(i,j,l) - ee*s(i,j,l) / cc(i,j,l)
      z(i,j,l) = u(i,j,l) * w(i,j,l)
    enddo; enddo; enddo
    call conv1(g,nfx,nfy, z, mxm,mym,nds, t, mx,my)
    p1 = 0.;  p2 = 0.
    do j=1,my; do i=1,mx
      if (b(i,j) /= vnul) then
        p1 = p1 + t(i,j)*q(i,j)
        p2 = p2 + q(i,j)*q(i,j)
      endif
    enddo; enddo
    do l=1,nds; do j=1,mym; do i=1,mxm
      p1 = p1 + ee*z(i,j,l)*p(i,j,l) / cc(i,j,l)
      p2 = p2 + ee*p(i,j,l)*p(i,j,l) / cc(i,j,l)
    enddo; enddo; enddo
    a = p1 / p2
    p1 = 0.;  p2 = 0.
    do j=1,my; do i=1,mx
      if (b(i,j) /= vnul) then
        q(i,j) = t(i,j) - a*q(i,j)
        p1 = p1 + r(i,j)*q(i,j)
        p2 = p2 + q(i,j)*q(i,j)
      endif
    enddo; enddo
    do l=1,nds; do j=1,mym; do i=1,mxm
      p(i,j,l) = z(i,j,l) - a*p(i,j,l)
      p1 = p1 - ee*s(i,j,l)*p(i,j,l) / cc(i,j,l)
      p2 = p2 + ee*p(i,j,l)*p(i,j,l) / cc(i,j,l)
    enddo; enddo; enddo
    a = p1 / p2
    ss1 = 0.;  ss2 = 0.
    do l=1,nds; do j=1,mym; do i=1,mxm
      s(i,j,l) = s(i,j,l) + a*p(i,j,l)
      ss = s(i,j,l) * s(i,j,l)
      cc(i,j,l) = (ss + del2) / v(i,j,l)
      ss1 = ss1 + ss / cc(i,j,l)
      ss2 = ss2 + p(i,j,l)*p(i,j,l) / cc(i,j,l)
    enddo; enddo; enddo
    rv = sqrt(ss1 / float(mxm*mym*nds))
    rs = a * sqrt(ss2 / float(mxm*mym*nds))
    n = 0;  ss2 = 0.
    do j=1,my; do i=1,mx
      if (b(i,j) /= vnul) then
        r(i,j) = r(i,j) - a*q(i,j)
        ss2 = ss2 + r(i,j)*r(i,j)
        n = n + 1
      endif
    enddo; enddo
    dv = sqrt(ss2 / float(n))
    tr = sqrt((ss2 + ss1*ee) / float(n + mxm*mym*nds))
    if (((trp-tr)/trp)*100.  >=  THRSH) then
      kc = 0
    else
      kc = kc + 1
    endif
    write(6,'(a5,i4,3x,a13,f7.3,a2,f7.4,3x,a6,f7.2,3x,a6,f7.2)')   &
      ' loop', klp, 'rms.Step/Src.', rs, ' /', rv, 'Misfit', dv, 'T.Res.', tr
    write(9,*) klp, 0., rs, rv, dv, tr

    if (kc /= 0) then
      ee =  ee / 2.
      tr = sqrt((ss2 + ss1*ee) / float(n + mxm*mym*nds))
      write(6,'(a5,i4,6x,a25,e10.3,10x,a6,f7.2)')   &
        ' loop', klp, 'Updated hyper-parameter =', ee, 'T.Res.', tr
      write(9,*) klp, ee, rs, rv, dv, tr
     endif
   trp = tr
  enddo

  write(6,'(1x)')
  endfile(9)
  rewind(9)

  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, anul, -1., klp
  do j=1,mym
    write(21,'((f7.1,9(1x,f7.1)))') (h(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, anul, -1., klp
    do j=1,mym
      do i=1,mxm
        if (s(i,j,l) /= anul) then
          if (s(i,j,l) < -9999.8) s(i,j,l) = -9999.8
          if (s(i,j,l) > +9999.8) s(i,j,l) = +9999.8
        endif
      enddo
      write(21,'((f7.1,9(1x,f7.1)))') (s(i,j,l),i=1,mxm)
    enddo
  enddo
  close(21)
  if (kout == 1) then
    call conv1(g,nfx,nfy, s, mxm,mym,nds, r, mx,my)
    do j=1,my
      y = float(j) - cmy
      do i=1,mx
        x = float(i) - cmx
        if (b(i,j) == vnul) then
          r(i,j) = vnul
          q(i,j) = vnul
        else
          q(i,j) = r(i,j) + (af*x + bf*y + cf)
          b(i,j) = b(i,j) - r(i,j)
        endif
      enddo
    enddo
    write(8,'(a8,4x,i4,4i8)') 'SynthMdl', ncc, iorg,korg, ispa,ispb
    write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
      ixs,iys, mszx,mszy, mx,my, vnul, -1.
    do j=1,my
      write(8,'((f7.1,9(1x,f7.1)))') (r(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'AUX2 <== SynthMdl : Synthetic model anomaly '
    if (yn == 'y') then
      write(8,'(a8,4x,i4,4i8)') 'Mdl+Trnd',ncc,iorg,korg,ispa,ispb
      write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
        ixs,iys, mszx,mszy, mx,my, vnul, -1.
      do j=1,my
        write(8,'((f7.1,9(1x,f7.1)))') (q(i,j),i=1,mx)
      enddo
      write(6,'(a)') 'AUX2 <== Mdl+Trnd : Trend-added Synth. model'
    endif
    write(8,'(a8,4x,i4,4i8)') 'Residual', ncc, iorg,korg, ispa,ispb
    write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
      ixs,iys, mszx,mszy, mx,my, vnul, -1.
    do j=1,my
      write(8,'((f7.1,9(1x,f7.1)))') (b(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'AUX2 <== Residual : Residual <(Obs.)-(Model)>'
    close(8)
  endif
  close(20)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'CIMGS : min.volume imaging'
    write(99,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
    if (yn == 'y') then
      write(99,'(22x,a)') 'Linear trend removed'
      write(99,'(a,f8.2,a,f7.3,a,f7.3,a)')   &
        ' Trend surface in nT :', cf, ' + (', af, ' * x) + (', bf, ' * y)'
      write(6,'(5x,a)')   &
        '( x/y : N/E coord. in mesh unit referred to the center of area )'
    else
      write(99,'(22x,a)') 'without removal of Linear trend'
    endif
    write(99,'(a,a)') ' CFIM matrix infile  : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' AIMG model infile   : ', fnam3(1:lrtrim(fnam3))
    if (fnam6(ldr:lrtrim(fnam6))  ==  'U') then
      write(99,'(a)')   ' AIMG model output   : <UPDATE infile>'
    else
      write(99,'(a,a)') ' AIMG model outfile  : ', fnam6(1:lrtrim(fnam6))
    endif
    write(99,'(a,a)') ' FSCL Par.Scal.Coefs : ', fnam4(1:lrtrim(fnam4))
    write(99,'(a,f4.1)') ' Par.Scale Weighting Power : ', wp
    write(99,'(a,f5.2)') ' Source Support Threshold (A/m) : ', del
    write(99,'(a,f5.2)') ' SrcVolume/Res. Initial Wt.Ratio : ', ratio
    if (kout == 1) write(99,'(a,a)')   &
      ' Auxiary(2) outfile  : ', 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,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,i4)') ' Loop Count : ', nloop
    if (nlp <= 0) then
      write(99,'(a,f6.2,a)')   &
        '   Starting from Initial magnetization :', vm, ' A/m'
    endif
    do
      read(9,*,iostat=IER) klp, ee, rs, rv, dv, tr
      if (IER < 0) exit
      if (ee == 0.) then
        write(99,'(a5,i4,3x,a13,f7.3,a2,f7.4,3x,a6,f7.2,3x,a6,f7.2)')   &
          ' loop', klp, 'rms.Step/Src.', rs, ' /', rv,   &
          'Misfit', dv, 'T.Res.', tr
      else
        write(99,'(a5,i4,6x,a25,e10.3,10x,a6,f7.2)')   &
          ' loop', klp, 'Updated hyper-parameter =', ee, 'T.Res.', tr
      endif
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end
!---
  subroutine conv1(g,nfx,nfy, s, mxm,mym,nds, r, mx,my)
  implicit none
  integer,intent(in) :: nfx, nfy, mxm, mym, nds, mx, my
  real,intent(inout) :: g(nfx,nfy)
  real,intent(in) :: s(mxm,mym,nds)
  real,intent(out) :: r(mx,my)
  character(8) :: area
  real :: ss
  integer :: nhx, nhy, i, j, ip, jp, k, kic, kjc, kis, kjs, is, js, ls

  call dpcini('.... read CFIM matrix (1) : ')
  read(20) area, j, ss, nhx, nhy, (i,k=1,18)
  do jp=1,my; do ip=1,mx;  r(ip,jp) = 0.;  enddo; enddo
  do ls=1,nds; do js=1,mym; do is=1,mxm
    read(20) k, k, k, kic, kjc, g
    kis = kic - nhx;  kjs = kjc - nhy
    do j=1,nfy
      jp = kjs + j;  if ((jp <= 0) .or. (jp > my)) cycle
      do i=1,nfx
        ip = kis + i;  if ((ip <= 0) .or. (ip > mx)) cycle
        r(ip,jp) = r(ip,jp) + g(i,j)*s(is,js,ls)
      enddo
    enddo
    call dpcent(((ls-1)*mym+(js-1))*mxm+is, mxm*mym*nds)
  enddo; enddo; enddo
  rewind(20)
  return
  end subroutine
!---
  subroutine conv2(g,nfx,nfy, r, mx,my, s, mxm,mym,nds)
  implicit none
  integer,intent(in) :: nfx, nfy, mx, my, mxm, mym, nds
  real,intent(inout) :: g(nfx,nfy)
  real,intent(in) :: r(mx,my)
  real,intent(out) :: s(mxm,mym,nds)
  character(8) :: area
  real :: ss
  integer :: nhx, nhy, i, j, ip, jp, k, is, js, ls, kic, kjc, kis, kjs

  call dpcini('.... read CFIM matrix (2) : ')
  read(20) area, j, ss, nhx, nhy, (i,k=1,18)
  do ls=1,nds; do js=1,mym; do is=1,mxm
    read(20) k, k, k, kic, kjc, g
    kis = kic - nhx;  kjs = kjc - nhy;  ss = 0.
    do j=1,nfy
      jp = kjs + j;  if ((jp <= 0) .or. (jp > my)) cycle
      do i=1,nfx
        ip = kis + i;  if ((ip <= 0) .or. (ip > mx)) cycle
        ss = ss + g(i,j)*r(ip,jp)
      enddo
    enddo
    s(is,js,ls) = ss
    call dpcent(((ls-1)*mym+(js-1))*mxm+is, mxm*mym*nds)
  enddo; enddo; enddo
  rewind(20)
  return
  end subroutine