diff -urN ORG/anam/@mkall anam/@mkall
--- ORG/anam/@mkall	2017-12-23 18:00:00.000000000 +0900
+++ anam/@mkall	2020-11-13 12:00:00.000000000 +0900
@@ -29,6 +29,7 @@
 fort  -Wall adeq.f90     -L/home/SHARE/lib -lgm -o ../bin/adeq
 fort  -Wall adeqc.f90    -L/home/SHARE/lib -lgm -o ../bin/adeqc
 fort  -Wall cdeq.f90     -L/home/SHARE/lib -lgm -o ../bin/cdeq
+fort  -Wall cdeql.f90    -L/home/SHARE/lib -lgm -o ../bin/cdeql
 fort  -Wall edeqs.f90    -L/home/SHARE/lib -lgm -o ../bin/edeqs
 fort  -Wall adeqs.f90    -L/home/SHARE/lib -lgm -o ../bin/adeqs
 fort  -Wall adeqsc.f90   -L/home/SHARE/lib -lgm -o ../bin/adeqsc
diff -urN ORG/anam/adeq.f90 anam/adeq.f90
--- ORG/anam/adeq.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/adeq.f90	2019-12-17 10:44:00.000000000 +0900
@@ -206,6 +206,7 @@
   dv = sqrt(ss / float(ndpt))
   write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, ds, dv
 
@@ -230,6 +231,7 @@
     dv = sqrt(ss / float(ndpt))
     write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, ds, dv
   enddo
   close(20)
diff -urN ORG/anam/adeqc.f90 anam/adeqc.f90
--- ORG/anam/adeqc.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/adeqc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -205,6 +205,7 @@
   lp = 0
   write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
     ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, lp, ds, dv
 
@@ -235,6 +236,7 @@
     endif
     write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
       ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, lp, ds, dv
     if ((lp >= 5) .or. (dv < 0.1) .or. (ilp >= nloop)) exit
   enddo
diff -urN ORG/anam/adeqs.f90 anam/adeqs.f90
--- ORG/anam/adeqs.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/adeqs.f90	2019-12-17 10:44:00.000000000 +0900
@@ -209,6 +209,7 @@
   dv = sqrt(ss / float(ndpt))
   write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, ds, dv
 
@@ -233,6 +234,7 @@
     dv = sqrt(ss / float(ndpt))
     write(6,'(a,i4,a,f8.2,a,f7.2)')   &
       ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, ds, dv
   enddo
   close(20)
diff -urN ORG/anam/adeqsc.f90 anam/adeqsc.f90
--- ORG/anam/adeqsc.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/adeqsc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -209,6 +209,7 @@
   lp = 0
   write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
     ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, lp, ds, dv
 
@@ -239,6 +240,7 @@
     endif
     write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
       ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, lp, ds, dv
     if ((lp >= 5) .or. (dv < 0.1) .or. (ilp >= nloop)) exit
   enddo
diff -urN ORG/anam/amag.f90 anam/amag.f90
--- ORG/anam/amag.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/amag.f90	2019-12-17 10:44:00.000000000 +0900
@@ -245,6 +245,7 @@
   dv = sqrt(ss / float(n))
   write(6,'(a,i4,a,f9.3,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, ds, dv
 
@@ -272,6 +273,7 @@
     dv = sqrt(ss / float(n))
     write(6,'(a,i4,a,f9.3,a,f7.2)')   &
       ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, ds, dv
   enddo
   write(6,'(1x)')
diff -urN ORG/anam/amagc.f90 anam/amagc.f90
--- ORG/anam/amagc.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/amagc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -241,6 +241,7 @@
   dv = sqrt(ss / float(n));  lp = 0
   write(6,'(a,i4,i3,a,f9.3,a,f7.2)')   &
     ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, lp, ds, dv
 
@@ -273,6 +274,7 @@
     endif
     write(6,'(a,i4,i3,a,f9.3,a,f7.2)')   &
       ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, lp, ds, dv
     if ((lp >= 5) .or. (dv < 0.1) .or. (ilp >= nloop)) exit
   enddo
diff -urN ORG/anam/ameq.f90 anam/ameq.f90
--- ORG/anam/ameq.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/ameq.f90	2019-12-17 10:44:00.000000000 +0900
@@ -239,6 +239,7 @@
   dv = sqrt(ss / float(n))
   write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, ds, dv
 
@@ -252,14 +253,14 @@
     do j=1,mym; do i=1,mxm;  p(i,j) = v(i,j) + a*p(i,j);  enddo; enddo
     call conv1(g,nfx,nfy, p,mxm,mym, q,mxf,myf)
     pqq = 0.
-    do j=1,myf; do i=1,myf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
+    do j=1,myf; do i=1,mxf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
     a = pvv / pqq;  ss = 0.
     do j=1,mym; do i=1,mxm
       s(i,j) = s(i,j) + a*p(i,j)
       if (t(i,j) /= tnul) ss = ss + p(i,j)*p(i,j)
     enddo; enddo
     ds = a * sqrt(ss / float(mxm*mym));  n = 0;  ss = 0.
-    do j=1,myf; do i=1,myf
+    do j=1,myf; do i=1,mxf
       if (f(i,j) /= vnul) then
         r(i,j) = r(i,j) - a*q(i,j)
         ss = ss + r(i,j)*r(i,j);  n = n + 1
@@ -268,6 +269,7 @@
     dv = sqrt(ss / float(n))
     write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, ds, dv
   enddo
   close(20)
diff -urN ORG/anam/ameqc.f90 anam/ameqc.f90
--- ORG/anam/ameqc.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/ameqc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -237,6 +237,7 @@
   lp = 0
   write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
     ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, lp, ds, dv
 
@@ -250,14 +251,14 @@
     do j=1,mym; do i=1,mxm;  p(i,j) = v(i,j) + a*p(i,j);  enddo; enddo
     call conv1(g,nfx,nfy, p,mxm,mym, q,mxf,myf)
     pqq = 0.
-    do j=1,myf; do i=1,myf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
+    do j=1,myf; do i=1,mxf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
     a = pvv / pqq;  ss = 0.
     do j=1,mym; do i=1,mxm
       s(i,j) = s(i,j) + a*p(i,j)
       if (t(i,j) /= tnul) ss = ss + p(i,j)*p(i,j)
     enddo; enddo
     ds = a * sqrt(ss / float(mxm*mym));  n = 0;  ss = 0.
-    do j=1,myf; do i=1,myf
+    do j=1,myf; do i=1,mxf
       if (f(i,j) /= vnul) then
         r(i,j) = r(i,j) - a*q(i,j)
         ss = ss + r(i,j)*r(i,j);  n = n + 1
@@ -271,6 +272,7 @@
     endif
     write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
       ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, lp, ds, dv
     if ((lp >= 5) .or. (dv < 0.1) .or. (ilp >= nloop)) exit
   enddo
diff -urN ORG/anam/ameqs.f90 anam/ameqs.f90
--- ORG/anam/ameqs.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/ameqs.f90	2019-12-17 10:44:00.000000000 +0900
@@ -242,6 +242,7 @@
   dv = sqrt(ss / float(n))
   write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, ds, dv
 
@@ -255,14 +256,14 @@
     do j=1,mym; do i=1,mxm;  p(i,j) = v(i,j) + a*p(i,j);  enddo; enddo
     call conv1(g,nfx,nfy, p,mxm,mym, q,mxf,myf)
     pqq = 0.
-    do j=1,myf; do i=1,myf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
+    do j=1,myf; do i=1,mxf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
     a = pvv / pqq;  ss = 0.
     do j=1,mym; do i=1,mxm
       s(i,j) = s(i,j) + a*p(i,j)
       if (t(i,j) /= tnul) ss = ss + p(i,j)*p(i,j)
     enddo; enddo
     ds = a * sqrt(ss / float(mxm*mym));  n = 0;  ss = 0.
-    do j=1,myf; do i=1,myf
+    do j=1,myf; do i=1,mxf
       if (f(i,j) /= vnul) then
         r(i,j) = r(i,j) - a*q(i,j)
         ss = ss + r(i,j)*r(i,j);  n = n + 1
@@ -271,6 +272,7 @@
     dv = sqrt(ss / float(n))
     write(6,'(a,i4,a,f8.2,a,f7.2)')   &
     ' loop', klp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, ds, dv
   enddo
   close(20)
diff -urN ORG/anam/ameqsc.f90 anam/ameqsc.f90
--- ORG/anam/ameqsc.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/ameqsc.f90	2019-12-17 10:44:00.000000000 +0900
@@ -241,6 +241,7 @@
   lp = 0
   write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
     ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+  flush(6)        !--- or  call flush(6)
   open(9,status='scratch')
   write(9,*) klp, lp, ds, dv
 
@@ -254,14 +255,14 @@
     do j=1,mym; do i=1,mxm;  p(i,j) = v(i,j) + a*p(i,j);  enddo; enddo
     call conv1(g,nfx,nfy, p,mxm,mym, q,mxf,myf)
     pqq = 0.
-    do j=1,myf; do i=1,myf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
+    do j=1,myf; do i=1,mxf; pqq = pqq + q(i,j)*q(i,j);  enddo; enddo
     a = pvv / pqq;  ss = 0.
     do j=1,mym; do i=1,mxm
       s(i,j) = s(i,j) + a*p(i,j)
       if (t(i,j) /= tnul) ss = ss + p(i,j)*p(i,j)
     enddo; enddo
     ds = a * sqrt(ss / float(mxm*mym));  n = 0;  ss = 0.
-    do j=1,myf; do i=1,myf
+    do j=1,myf; do i=1,mxf
       if (f(i,j) /= vnul) then
         r(i,j) = r(i,j) - a*q(i,j)
         ss = ss + r(i,j)*r(i,j);  n = n + 1
@@ -276,6 +277,7 @@
     endif
     write(6,'(a,i4,i3,a,f8.2,a,f7.2)')   &
       ' loop', klp, lp, '   rms.Step :', ds, '   rms.Residual :', dv
+    flush(6)        !--- or  call flush(6)
     write(9,*) klp, lp, ds, dv
     if ((lp >= 5) .or. (dv < 0.1) .or. (ilp >= nloop)) exit
   enddo
diff -urN ORG/anam/calmas.f90 anam/calmas.f90
--- ORG/anam/calmas.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/calmas.f90	2020-11-13 12:00:00.000000000 +0900
@@ -6,7 +6,7 @@
   character(5) :: nmdl(4) = (/'Block', 'hRect', 'vLine', 'Point'/)
   character(80) :: wdr, flog, fnam1, fnam2, buf
   character(20) :: sdtm;  character(8) :: area
-  integer :: ldr, n, i, j, kl, ks, IER
+  integer :: ldr, n, i, j, kl, ks
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixl, ixh, iyl, iyh, mszx, mszy, mx, my
   real :: am, cx, cy, dip, dec, hnul, phi, the, vol, rx, ry, rz
@@ -95,10 +95,12 @@
     if ((ks == 1) .or. (ks == 2)) then
       call prompt('   Northing Center (cx), N-S size (rx) ')
       call gparmf2('(in km) ==> ', cx, rx)
+      if (rx <= 0.) call abendm('invalid N-S size')
       if (kl == 1) cx = cx + xl
       pv(1) = cx;  pv(2) = rx
       call prompt('   Easting Center (cy), E-W size (ry)  ')
       call gparmf2('(in km) ==> ', cy, ry)
+      if (ry <= 0.) call abendm('invalid E-W size')
       if (kl == 1) cy = cy + yl
       pv(3) = cy;  pv(4) = ry
     else
@@ -113,6 +115,7 @@
     if ((ks == 1) .or. (ks == 3)) then
       call prompt('   Top depth (za), and Thickness (rz)  ')
       call gparmf2('(in km) ==> ', za, rz)
+      if (rz <= 0.) call abendm('invalid Thickness')
       pv(5) = za;  pv(6) = rz
     else
       call gparmf('    (in km) ==> ', zc)
@@ -120,10 +123,13 @@
     endif
     if (ks == 4) then
       call gparmf('  Equivalent source volume (in cubic-km) ==> ', vol)
+      if (vol <= 0.) call abendm('invalid Eqv.Volume')
     else if (ks == 3) then
       call gparmf('  Equivalent source area (in square-km) ==> ', vol)
+      if (vol <= 0.) call abendm('invalid Eqv.Area')
     else if (ks == 2) then
       call gparmf('  Equivalent source thickness (in km) ==> ', vol)
+      if (vol <= 0.) call abendm('invalid Eqv.Thickness')
     else
       vol = 0.
     endif
@@ -133,8 +139,8 @@
     call gparmf2('    Inclin., Declin. (in degrees) : ', the, phi)
     if ((the < -90.) .or. (the > 90.) .or.   &
         (phi < -180.) .or. (phi >= 360.)) call abendm('invalid value')
-    write(6,'(i4,2x,a5,f8.2,2f8.1)') n, nmdl(ks), am, the, phi
-    write(6,'(6x,3(f9.2,f7.2),f9.4)') (pv(i),i=1,7)
+    write(6,'(i4,1x,a5,2x,f12.3,2x,f10.2,2x,f10.2)') n, nmdl(ks), am, the, phi
+    write(6,'(4x,3(1x,f9.3,1x,f8.3),3x,e11.4)') (pv(i),i=1,7)
     write(9,*) ks, am, the, phi, pv
 
     call magafd(the, phi, dip, dec)
@@ -158,7 +164,22 @@
       endif
     enddo; enddo
   enddo
+  write(9,*) 0, 0., 0., 0., (0.,i=1,7)
   call clspin()
+  endfile(9)
+  rewind(9)
+
+  write(2,'(a)') '### CALMAS model ###'
+  do
+    read(9,*) ks, am, the, phi, pv
+    if (ks == 0) exit
+    write(2,'(a2,a5,1x,i2,2x,f12.3,2x,f10.2,2x,f10.2)')   &
+      '##', nmdl(ks), ks, am, the, phi
+    write(2,'(a1,3x,3(1x,f9.3,1x,f8.3),3x,e11.4)') '#', (pv(i),i=1,7)
+  enddo
+  write(2,'(a7,1x,i2,2x,f12.3,2x,f10.2,2x,f10.2)') '##AmbMF', 0, 0., dip, dec
+  rewind(9)
+
   write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
   write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
     ixl, iyl, mszx, mszy, mx, my, 9999.9, 0.
@@ -172,8 +193,6 @@
     write(2,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
   enddo
   close(2)
-  endfile(9)
-  rewind(9)
 
   if (flog(ldr:ldr) /= ' ') then
     open(99,file=flog,access='append')
@@ -190,13 +209,14 @@
       'model anom. ', xl, yl, mszx, mszy, mx, my, '    var.'
     write(99,'(4x,a,f6.1,a,f6.1)')   &
       'Ambient field direction:    dip=', dip, '    dec=', dec
-    n = 0
+    write(99,'(a)') '### CALMAS model ###'
     do
-      read(9,*,iostat=IER) ks, am, the, phi, pv
-      if (IER < 0) exit
+      read(9,*) ks, am, the, phi, pv
+      if (ks == 0) exit
       n = n + 1
-      write(99,'(i4,2x,a5,f8.2,2f8.1)') n, nmdl(ks), am, the, phi
-      write(99,'(6x,3(f9.2,f7.2),f9.4)') (pv(i),i=1,7)
+      write(99,'(a2,a5,1x,i2,2x,f12.3,2x,f10.2,2x,f10.2)')   &
+        '##', nmdl(ks), ks, am, the, phi
+      write(99,'(a1,3x,3(1x,f9.3,1x,f8.3),3x,e11.4)') '#', (pv(i),i=1,7)
     enddo
     call strdtm(sdtm)
     write(99,'(a,a,a)') '=============== ', sdtm, '==============='
diff -urN ORG/anam/cdeql.f90 anam/cdeql.f90
--- ORG/anam/cdeql.f90	1970-01-01 09:00:00.000000000 +0900
+++ anam/cdeql.f90	2020-11-13 12:00:00.000000000 +0900
@@ -0,0 +1,212 @@
+!--------( CDEQL : calculate the effect of <adeq> source
+!                   on a specified line [StdLIN] )--------
+  implicit none
+  integer :: lwkdir, lrtrim
+  real,parameter :: DPI = 3.141593*2.
+  real,allocatable :: s(:,:), t(:,:)
+  character(80) :: wdr, flog, fnam1, fnam2, fnam3, buf
+  character(72) :: out;  character(20) :: sdtm
+  character(8) :: area, area1, lnam
+  integer :: ldr, klp, i, j, nhx, nhy, ki, kj, ic, jc, is, js
+  integer :: nc, ncc, iorg, korg, ispa, ispb
+  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
+  integer :: ixl, iyl, mszx, mszy, mxm, mym
+  integer :: ixl1, iyl1, mszx1, mszy1, mxm1, mym1
+  integer :: ndpt, idpt, ier
+  real :: smx, smy, amxy, anul, dmy, tnul, hw, hwm
+  real :: ss, x, y, z, a, rr, sr, hp, alt, res, red
+  real (8) :: xl, yl, xp, yp, dfi, dfk, dxn, dye
+
+  call premsg('CDEQL : calculate the effect of ADEQ model source')
+  call premsg('         on a specified line [StdLIN]')
+  call opnpin()
+  ldr = lwkdir(50,wdr) + 1
+  flog = wdr;  fnam1 = wdr;  fnam2 = wdr;  fnam3 = wdr
+  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
+
+  call gparma('Enter ADEQ model input filename  ==> ', 81-ldr, fnam2(ldr:80))
+  open(2,file=fnam2,status='old')
+  do
+    read(2,'(a)') buf
+    if (buf(1:1) /= '#') exit
+  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(2,*) ixl, iyl, mszx, mszy, mxm, mym, anul, dmy, klp
+  if ((mxm < 10) .or. (mym < 10)) call abendm('illegal array size')
+
+  allocate(s(mxm,mym));  allocate(t(mxm,mym))
+
+  read(2,*) ((s(i,j),i=1,mxm),j=1,mym)
+  do
+    read(2,'(a)') buf
+    if (buf(1:1) /= '#') exit
+  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(2,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, tnul
+  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. (mszx1 /= mszx) .or.   &
+      (mszy1 /= mszy) .or. (mxm1 /= mxm) .or. (mym1 /= mym))   &
+    call abendm('Header inconsistent')
+  read(2,*) ((t(i,j),i=1,mxm),j=1,mym)
+  close(2)
+
+  call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
+  xl = dfloat(ixl)/1000.;  yl = dfloat(iyl)/1000.
+  smx = float(mszx);  smy = float(mszy);  amxy = smx * smy
+
+  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
+    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
+  call premsg(out)
+  write(out,'(a,2a10,2x,2a8,2a6,a8)')   &
+    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
+  call premsg(out)
+  call premsg(out)
+  write(out,'(a,2f10.2,2i8,2i6,a8)')   &
+    '--- eq.src.anom.', xl,yl, mszx,mszy, mxm,mym, '    var.'
+  call premsg(out)
+  call premsg('  Truncation of Source effect ?')
+  call gparmf('    Half width of calc.window (in km) : ', hw)
+  hwm = hw * 1000.;  nhx = nint(hwm/smx);  nhy = nint(hwm/smy)
+
+  call premsg('Enter filename of StdLIN on which the effect of')
+  call gparma('    ADEQ source model is calculated ==> ',   &
+              81-ldr, fnam1(ldr:80))
+
+  open(1,file=fnam1,status='old')
+  ndpt = 0;  idpt = 0
+  do
+    read(1,'(a)',iostat=IER) buf
+    if (IER < 0) exit
+    if ((buf(1:1) == '&') .or. (buf(1:1) == '%') .or. (buf(1:1) == '#')) cycle
+    ndpt = ndpt + 1
+  enddo
+  rewind(1)
+
+  call gparma('Enter filename of Output anomaly data ==> ',   &
+              81-ldr, fnam3(ldr:80))
+  call clspin()
+
+  call strdtm(sdtm)
+  open(3,file=fnam3,status='new')
+  write(3,'(a1,a11,1x,a11,4(1x,a9))')   &
+    '#', 'Lat(min)', 'Lon(min)', 'Alt(m)', 'ObsMag', 'Model', 'Diff.'
+
+  write(6,'(a)')   '----------------------------------------'
+  write(6,'(a)')   'CDEQL : calculate the effect of ADEQ source model'
+  write(6,'(a)')   '          on a specified line [StdLIN]'
+  write(6,'(a,a)') ' ADEQ model infile     : ', fnam2(1:lrtrim(fnam2))
+  write(6,'(a)')   ' Filename of StdLIN on which the effect of ADEQ source'
+  write(6,'(a,a)') '   model is calculated : ', fnam1(1:lrtrim(fnam1))
+  write(6,'(a,a)') ' Anomamy data outfile  : ', fnam3(1:lrtrim(fnam3))
+  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,'(a,f6.2,a,i4,a,i4,a)')   &
+    ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  mesh )'
+
+  call dpcini('... calc. Reduced anomaly : ')
+  do
+    read(1,'(a)',iostat=IER) buf
+    if (IER < 0) exit
+    if (buf(1:1) == '#') cycle
+    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) then
+      read(buf,'(1x,a8)') lnam
+      write(3,'(a1,a8)') '&', lnam
+      cycle
+    else
+      idpt = idpt + 1
+      i = index(buf,'N');  if (i /= 0) buf(i:i) = ' '
+      i = index(buf,'E');  if (i /= 0) buf(i:i) = ' '
+      i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
+      i = index(buf,'nT');  if (i /= 0) buf(i:i+1) = '  '
+      read(buf,*) dfi, dfk, alt, res
+      call cvdiken(dfi, dfk, dye, dxn)
+      xp = (dxn-xl)*1000.;  yp = (dye-yl)*1000.;  hp = alt
+      ic = nint(xp/smx) + 1;  jc = nint(yp/smy) + 1
+      ss = 0.
+      do kj=-nhy,+nhy
+        js = jc + kj;  y = sngl(yp - dfloat(js)*smy)
+        do ki=-nhx,+nhx
+          is = ic + ki;  x = sngl(xp - dfloat(is)*smx)
+          if ((is <= 0) .or. (is > mxm)) cycle
+          if ((js <= 0) .or. (js > mym)) cycle
+          if (t(is,js) == tnul) cycle
+          z = hp - t(is,js)
+          rr = x*x + y*y + z*z;  sr = sqrt(rr)
+          a = z*amxy / (rr*sr*DPI);  ss = ss  +  a*s(is,js)
+        enddo
+      enddo
+      red = ss
+      write(3,'(1x,f11.5,1x,f12.5,4(1x,f9.2))')   &
+        dfi, dfk, alt, res, red, res-red
+    endif
+    call dpcent(idpt, ndpt)
+  enddo
+  close(3)
+  close(1)
+
+  if (flog(ldr:ldr) /= ' ') then
+    open(99,file=flog,access='append')
+    write(99,'(//a)') '----------------------------------------'
+    write(99,'(a)')   'CDEQL : calculate the effect of ADEQ source model'
+    write(99,'(a)')   '          on a specified line [StdLIN]'
+    write(99,'(a,a)') ' ADEQ model infile     : ', fnam2(1:lrtrim(fnam2))
+    write(99,'(a)')   ' Filename of StdLIN on which the effect of ADEQ source'
+    write(99,'(a,a)') '   model is calculated : ', fnam1(1:lrtrim(fnam1))
+    write(99,'(a,a)') ' Anomamy data outfile  : ', fnam3(1:lrtrim(fnam3))
+    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,'(a,f6.2,a,i4,a,i4,a)')   &
+      ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  mesh )'
+    call strdtm(sdtm)
+    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
+    close(99)
+  endif
+  stop
+  end
diff -urN ORG/anam/galtf.f90 anam/galtf.f90
--- ORG/anam/galtf.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/galtf.f90	2020-11-13 12:00:00.000000000 +0900
@@ -9,7 +9,7 @@
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixl, iyl, mszx, mszy, mx, my
   real :: sr, srr, alt, alth, altl
-  real :: xl, yl, smx, smy, xn, ye, xc, yc, x, y, rr, wt
+  real :: xl, yl, smx, smy, xp, yp, xc, yc, x, y, rr, wt
   real(8) :: dfi, dfk, dxn, dye
 
   call premsg('GALTF : gen. grid-data of flying altitude')
@@ -100,8 +100,8 @@
     i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
     read(buf,*) dfi, dfk, alt
     call cvdiken(dfi, dfk, dye, dxn)
-    xn = sngl(dxn);  ip = nint((xn-xl)/smx);  xc = xn-xl - float(ip)*smx
-    ye = sngl(dye);  jp = nint((ye-yl)/smy);  yc = ye-yl - float(jp)*smy
+    xp = sngl(dxn-xl);  ip = nint(xp/smx);  xc = xp - float(ip)*smx
+    yp = sngl(dye-yl);  jp = nint(yp/smy);  yc = yp - float(jp)*smy
     do j=-jr,+jr
       if (((jp+j) < 0) .or. ((jp+j) >= my)) cycle
       y = float(j)*smy - yc
@@ -153,8 +153,8 @@
     i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
     read(buf,*) dfi, dfk, alt
     call cvdiken(dfi, dfk, dye, dxn)
-    xn = sngl(dxn);  ip = nint((xn-xl)/smx);  xc = xn-xl - float(ip)*smx
-    ye = sngl(dye);  jp = nint((ye-yl)/smy);  yc = ye-yl - float(jp)*smy
+    xp = sngl(dxn-xl);  ip = nint(xp/smx);  xc = xp - float(ip)*smx
+    yp = sngl(dye-yl);  jp = nint(yp/smy);  yc = yp - float(jp)*smy
     do j=-jr,+jr
       if (((jp+j) < 0) .or. ((jp+j) >= my)) cycle
       y = float(j)*smy - yc
diff -urN ORG/anam/galts.f90 anam/galts.f90
--- ORG/anam/galts.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/galts.f90	2020-11-13 12:00:00.000000000 +0900
@@ -9,7 +9,7 @@
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: ixl, iyl, mszx, mszy, mx, my
   real :: sr, srr, alt, sa, sb, sc
-  real :: xl, yl, smx, smy, xn, ye, xc, yc, x, y, rr, wt
+  real :: xl, yl, smx, smy, xp, yp, xc, yc, x, y, rr, wt
   real(8) :: dfi, dfk, dxn, dye
 
   call premsg('GALTS : generate smooth surface of flight level')
@@ -100,8 +100,8 @@
     i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
     read(buf,*) dfi, dfk, alt
     call cvdiken(dfi, dfk, dye, dxn)
-    xn = sngl(dxn);  ip = nint((xn-xl)/smx);  xc = xn-xl - float(ip)*smx
-    ye = sngl(dye);  jp = nint((ye-yl)/smy);  yc = ye-yl - float(jp)*smy
+    xp = sngl(dxn-xl);  ip = nint(xp/smx);  xc = xp - float(ip)*smx
+    yp = sngl(dye-yl);  jp = nint(yp/smy);  yc = yp - float(jp)*smy
     do j=-jr,+jr
       if (((jp+j) < 0) .or. ((jp+j) >= my)) cycle
       y = float(j)*smy - yc
@@ -140,8 +140,8 @@
     i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
     read(buf,*) dfi, dfk, alt
     call cvdiken(dfi, dfk, dye, dxn)
-    xn = sngl(dxn);  ip = nint((xn-xl)/smx);  xc = xn-xl - float(ip)*smx
-    ye = sngl(dye);  jp = nint((ye-yl)/smy);  yc = ye-yl - float(jp)*smy
+    xp = sngl(dxn-xl);  ip = nint(xp/smx);  xc = xp - float(ip)*smx
+    yp = sngl(dye-yl);  jp = nint(yp/smy);  yc = yp - float(jp)*smy
     do j=-jr,+jr
       if (((jp+j) < 0) .or. ((jp+j) >= my)) cycle
       y = float(j)*smy - yc
diff -urN ORG/anam/plamag.f90 anam/plamag.f90
--- ORG/anam/plamag.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/plamag.f90	2020-02-25 00:00:00.000000000 +0900
@@ -5,14 +5,14 @@
   real,allocatable :: f(:,:), g(:)
   character(80) :: wdr, flog, fnin, fobs, fnam, buf
   character(72) :: out;  character(20) :: sdtm
-  character(8) :: area, area1; character(4) :: s4
+  character(8) :: area, area1; character(6) :: ss
   character(3) :: spl = '   '; character(1) :: yn
   integer :: ldr, i, j, mc, nks, nkn, nkw, nke, npx, npy
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: ixs1, iys1, mszx1, mszy1, nmx1, nmy1
-  integer :: ii, ix, mx, iy, my, n, kpl, lcap, kscl, km, iinc, kinc
+  integer :: ii, ix, mx, iy, my, n, kpl, lcap, kscl, km, ls, iinc, kinc
   integer :: klll, kshl, isf, itf, ksf, ktf, idf, kdf, k
   real :: vnul, xs, ys, vnul1, xs1, ys1, dmy, gsum, gmin, gmax, gdev, dg
   real :: step, astp, wide, high, fscl, ccm, dscl, sclr, sclu, x, y
@@ -29,7 +29,7 @@
     read(10,'(a)') buf
     if (buf(1:1) /= '#') exit
   enddo
-    if (buf(13:16) == '    ') then      !-- v2005
+  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
@@ -127,7 +127,7 @@
     i = (ii-1)*npy
     ix = ixs + mszx*(nks+ii-1)
     mx = nint(float(ix-ixs1)/float(mszx1))
-    if ((mx <= 0) .or. (mx > nmx1)) then
+    if ((mx < 0) .or. (mx >= nmx1)) then
       do j=1,npy
         g(i+j) = vnul
       enddo
@@ -135,7 +135,7 @@
       do j=1,npy
         iy = iys + mszy*(nkw+j-1)
         my = nint(float(iy-iys1)/float(mszy1))
-        if ((my <= 0) .or. (my > nmy1)) then
+        if ((my < 0) .or. (my >= nmy1)) then
           g(i+j) = vnul
         else
           if (f(mx,my) == vnul1) then
@@ -199,17 +199,55 @@
   endif
   fscl = fscl * 100.;  ccm = fscl/100000.;  sclr = 0.;  sclu = 0.
   kscl = 0;  klll = 0;  kshl = 0
-  km = nint(fscl)/25000
-  if (km == 0)  km = 1
-  if (km > 50) km = 50
-  if ((km > 2) .and. (km < 5))   km = 2
-  if ((km > 5) .and. (km < 10))  km = 5
-  if ((km > 10) .and. (km < 20)) km = 10
-  if ((km > 20) .and. (km < 50)) km = 20
-  iinc = km;  kinc = km + km/2
-  if (km == 50) iinc = 60
-  if ((km == 1) .or. (km == 5) .or. (km == 50)) kinc = iinc
-  dscl = float(km)/ccm
+  if (fscl < 25000.) then               !ScaleBar (km [m])
+    km = nint(fscl)/25
+    if (km < 2) then
+      km = 1;  ss = '1 m';  ls = 3
+    else if (km < 5) then
+      km = 2;  ss = '2 m';  ls = 3
+    else if (km < 10) then
+      km = 5;  ss = '5 m';  ls = 3
+    else if (km < 20) then
+      km = 10;  ss = '10 m';  ls = 4
+    else if (km < 50) then
+      km = 20;  ss = '20 m';  ls = 4
+    else if (km < 100) then
+      km = 50;  ss = '50 m';  ls = 4
+    else if (km < 200) then
+      km = 100;  ss = '100 m';  ls = 5
+    else if (km < 500) then
+      km = 200;  ss = '200 m';  ls = 5
+    else
+      km = 500;  ss = '500 m';  ls = 5
+    endif
+      iinc = 1;  kinc = 1               ![LatLonIntv]: 1min x 1min
+    dscl = float(km*100) / fscl
+  else                                  !ScaleBar (km [km])
+    km = nint(fscl)/25000
+    if (km < 2) then
+      km = 1;  ss = '1 km';  ls = 4     ![LatLonIntv]: 1min x 1min
+    else if (km < 5) then
+      km = 2;  ss = '2 km';  ls = 4         !  2min x 3min
+    else if (km < 10) then
+      km = 5;  ss = '5 km';  ls = 4         !  5min x 5min
+    else if (km < 20) then
+      km = 10;  ss = '10 km';  ls = 5       ! 10min x 15min
+    else if (km < 50) then
+      km = 20;  ss = '20 km';  ls = 5       ! 20min x 30min
+    else if (km < 100) then
+      km = 50;  ss = '50 km';  ls = 5       !  1deg x 1deg
+    else if (km < 200) then
+      km = 100;  ss = '100 km';  ls = 6     !  2deg x 3deg
+    else if (km < 500) then
+      km = 200;  ss = '200 km';  ls = 6     !  5deg x 5deg
+    else
+      km = 500;  ss = '500 km';  ls = 6     ! 10deg x 10deg
+    endif
+    iinc = km;  if (km >= 50) iinc = km + km/5;  kinc = iinc
+    if ((km==2) .or. (km==10) .or. (km==20) .or. (km==100))   &
+      kinc = iinc + iinc/2
+    dscl = float(km)/ccm
+  endif
   call prompt('  Draw  Scale Bar  ?  ')
   call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
   if ((yn == 'y') .or. (yn == 'Y')) then
@@ -255,8 +293,8 @@
   if ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
     write(6,'(a,$)') '      with '
     if (kshl /= 0) write(6,'(a,$)') '  Coast'
-    if (mod(kshl,100)/10 /= 0) write(6,'(a,$)') '+River'
-    if (kshl/100 /= 0) write(6,'(a,$)') '+PrefB'
+    if (mod(kshl/10,10) /= 0) write(6,'(a,$)') '+River'
+    if ((kshl/100) /= 0) write(6,'(a,$)') '+PrefB'
     if (klll /= 0) write(6,'(a,$)') '  LatLonLines'
     if (kscl /= 0) then
       write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
@@ -279,9 +317,8 @@
     call plot(sclr+dscl, sclu,     2)
     call plot(sclr,      sclu,     2)
     call plot(sclr,      sclu+0.3, 2)
-    write(s4,'(i2,a2)') km, 'km'
     call lstyle('T', 0.5, 0., 0, -255)
-    call ptext(s4, 4, sclr+dscl/2., sclu+0.5, 2)
+    call ptext(ss, ls, sclr+dscl/2., sclu+0.5, 2)
     call newpen(1)
   endif
   call scisor(0., 0., wide, high)
@@ -352,8 +389,8 @@
     if ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
       write(99,'(a,$)') '      with '
       if (kshl /= 0) write(99,'(a,$)') '  Coast'
-      if (mod(kshl,100)/10 /= 0) write(99,'(a,$)') '+River'
-      if (kshl/100 /= 0) write(99,'(a,$)') '+PrefB'
+      if (mod(kshl/10,10) /= 0) write(99,'(a,$)') '+River'
+      if ((kshl/100) /= 0) write(99,'(a,$)') '+PrefB'
       if (klll /= 0) write(99,'(a,$)') '  LatLonLines'
       if (kscl /= 0) then
         write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
diff -urN ORG/anam/plamagc.f90 anam/plamagc.f90
--- ORG/anam/plamagc.f90	2017-12-23 18:00:00.000000000 +0900
+++ anam/plamagc.f90	2020-02-25 00:00:00.000000000 +0900
@@ -5,14 +5,14 @@
   real,allocatable :: f(:,:), g(:)
   character(80) :: wdr, flog, fnin, fobs, fnam, buf
   character(72) :: out;  character(20) :: sdtm
-  character(8) :: area, area1; character(4) :: s4
+  character(8) :: area, area1; character(6) :: ss
   character(3) :: spl = '   '; character(1) :: yn
   integer :: ldr, i, j, mc, nks, nkn, nkw, nke, npx, npy
   integer :: nc, ncc, iorg, korg, ispa, ispb
   integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
   integer :: ixs, iys, mszx, mszy, nmx, nmy
   integer :: ixs1, iys1, mszx1, mszy1, nmx1, nmy1
-  integer :: ii, ix, mx, iy, my, n, kpl, kscl, km, iinc, kinc
+  integer :: ii, ix, mx, iy, my, n, kpl, kscl, km, ls, iinc, kinc
   integer :: klll, kshl, isf, itf, ksf, ktf, idf, kdf, k
   real :: vnul, xs, ys, vnul1, xs1, ys1, dmy, gsum, gmin, gmax, gdev, dg
   real :: step, vmed, cint, vlo, vhi, wide, high, fscl, ccm, dscl, sclr, sclu
@@ -34,14 +34,16 @@
     if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
     nc = ncc
     if (ncc < 200) then
-      nc = ncc; ncc = nc + 800
+      nc = ncc; ncc = nc + 800;  mc = 0
     else
-      ncc = ncc - 200; nc = ncc
+      ncc = ncc - 200; nc = ncc;  mc = 1
     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
+    nc = ncc;  mc = 1
+    if ((ncc >= 800) .and. (ncc < 1000)) then
+      nc = ncc - 800;  mc = 0
+    endif
   endif
   if ((nc >= 1) .and. (nc <= 62)) then
     iorg = 0;  ispa = 0;  ispb = 0
@@ -125,7 +127,7 @@
     i = (ii-1)*npy
     ix = ixs + mszx*(nks+ii-1)
     mx = nint(float(ix-ixs1)/float(mszx1))
-    if ((mx <= 0) .or. (mx > nmx1)) then
+    if ((mx < 0) .or. (mx >= nmx1)) then
       do j=1,npy
         g(i+j) = vnul
       enddo
@@ -133,7 +135,7 @@
       do j=1,npy
         iy = iys + mszy*(nkw+j-1)
         my = nint(float(iy-iys1)/float(mszy1))
-        if ((my <= 0) .or. (my > nmy1)) then
+        if ((my < 0) .or. (my >= nmy1)) then
           g(i+j) = vnul
         else
           if (f(mx,my) == vnul1) then
@@ -197,17 +199,55 @@
   endif
   fscl = fscl * 100.;  ccm = fscl/100000.;  sclr = 0.;  sclu = 0.
   kscl = 0;  klll = 0;  kshl = 0
-  km = nint(fscl)/25000
-  if (km == 0)  km = 1
-  if (km > 50) km = 50
-  if ((km > 2) .and. (km < 5))   km = 2
-  if ((km > 5) .and. (km < 10))  km = 5
-  if ((km > 10) .and. (km < 20)) km = 10
-  if ((km > 20) .and. (km < 50)) km = 20
-  iinc = km;  kinc = km + km/2
-  if (km == 50) iinc = 60
-  if ((km == 1) .or. (km == 5) .or. (km == 50)) kinc = iinc
-  dscl = float(km)/ccm
+  if (fscl < 25000.) then               !ScaleBar (km [m])
+    km = nint(fscl)/25
+    if (km < 2) then
+      km = 1;  ss = '1 m';  ls = 3
+    else if (km < 5) then
+      km = 2;  ss = '2 m';  ls = 3
+    else if (km < 10) then
+      km = 5;  ss = '5 m';  ls = 3
+    else if (km < 20) then
+      km = 10;  ss = '10 m';  ls = 4
+    else if (km < 50) then
+      km = 20;  ss = '20 m';  ls = 4
+    else if (km < 100) then
+      km = 50;  ss = '50 m';  ls = 4
+    else if (km < 200) then
+      km = 100;  ss = '100 m';  ls = 5
+    else if (km < 500) then
+      km = 200;  ss = '200 m';  ls = 5
+    else
+      km = 500;  ss = '500 m';  ls = 5
+    endif
+      iinc = 1;  kinc = 1               ![LatLonIntv]: 1min x 1min
+    dscl = float(km*100) / fscl
+  else                                  !ScaleBar (km [km])
+    km = nint(fscl)/25000
+    if (km < 2) then
+      km = 1;  ss = '1 km';  ls = 4     ![LatLonIntv]: 1min x 1min
+    else if (km < 5) then
+      km = 2;  ss = '2 km';  ls = 4         !  2min x 3min
+    else if (km < 10) then
+      km = 5;  ss = '5 km';  ls = 4         !  5min x 5min
+    else if (km < 20) then
+      km = 10;  ss = '10 km';  ls = 5       ! 10min x 15min
+    else if (km < 50) then
+      km = 20;  ss = '20 km';  ls = 5       ! 20min x 30min
+    else if (km < 100) then
+      km = 50;  ss = '50 km';  ls = 5       !  1deg x 1deg
+    else if (km < 200) then
+      km = 100;  ss = '100 km';  ls = 6     !  2deg x 3deg
+    else if (km < 500) then
+      km = 200;  ss = '200 km';  ls = 6     !  5deg x 5deg
+    else
+      km = 500;  ss = '500 km';  ls = 6     ! 10deg x 10deg
+    endif
+    iinc = km;  if (km >= 50) iinc = km + km/5;  kinc = iinc
+    if ((km==2) .or. (km==10) .or. (km==20) .or. (km==100))   &
+      kinc = iinc + iinc/2
+    dscl = float(km)/ccm
+  endif
   call prompt('  Draw  Scale Bar  ?  ')
   call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
   if ((yn == 'y') .or. (yn == 'Y')) then
@@ -254,8 +294,8 @@
   if ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
     write(6,'(a,$)') '      with '
     if (kshl /= 0) write(6,'(a,$)') '  Coast'
-    if (mod(kshl,100)/10 /= 0) write(6,'(a,$)') '+River'
-    if (kshl/100 /= 0) write(6,'(a,$)') '+PrefB'
+    if (mod(kshl/10,10) /= 0) write(6,'(a,$)') '+River'
+    if ((kshl/100) /= 0) write(6,'(a,$)') '+PrefB'
     if (klll /= 0) write(6,'(a,$)') '  LatLonLines'
     if (kscl /= 0) then
       write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
@@ -285,9 +325,8 @@
     call plot(sclr+dscl, sclu,     2)
     call plot(sclr,      sclu,     2)
     call plot(sclr,      sclu+0.3, 2)
-    write(s4,'(i2,a2)') km, 'km'
     call lstyle('T', 0.5, 0., 0, -255)
-    call ptext(s4, 4, sclr+dscl/2., sclu+0.5, 2)
+    call ptext(ss, ls, sclr+dscl/2., sclu+0.5, 2)
     call newpen(1)
   endif
   call scisor(0., 0., wide, high)
@@ -359,8 +398,8 @@
     if ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
       write(99,'(a,$)') '      with '
       if (kshl /= 0) write(99,'(a,$)') '  Coast'
-      if (mod(kshl,100)/10 /= 0) write(99,'(a,$)') '+River'
-      if (kshl/100 /= 0) write(99,'(a,$)') '+PrefB'
+      if (mod(kshl/10,10) /= 0) write(99,'(a,$)') '+River'
+      if ((kshl/100) /= 0) write(99,'(a,$)') '+PrefB'
       if (klll /= 0) write(99,'(a,$)') '  LatLonLines'
       if (kscl /= 0) then
         write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
